Replication for SVAR

Question

I am using Stata and studying VAR, Orthogonalized VAR, Structural VAR. I successfully replicated VAR and OVAR, but I failed to do that for SVAR. Could anyone suggest how to replicate it?

My Replication does not exactly match the R built-in results. How can I fix it?

Below is the built-in R code programmed by Danne, Christian. VARsignR: Estimating VARs Using Sign Restrictions in R, 2015.

###### Install required packages
install.packages("minqa") 
install.packages("HI") 
install.packages("mvnfast")
install.packages("lubridate")  
install.packages("VARsignR")  

###### Import data
rm(list = ls())
set.seed(12345)
library(VARsignR)

SVARdata <- read.csv("https://raw.githubusercontent.com/jayjeo/public/master/LaborShortage/SVARdata_seasonadjusted2_practice.csv")
SVARdata <- ts (SVARdata, frequency = 12, start = c(2012, 1))

###### set sign restrictions
constr <- c(-1,+2)  

###### Uhlig’s (2005) Rejection Method
model <- uhlig.reject(Y=SVARdata, nlags=1, draws=200, subdraws=100, nkeep=1000, KMIN=1, KMAX=6, constrained=constr, constant=TRUE, steps=120)
irfs <- model$IRFS 
vl <- c("y","z")
irfplot(irfdraws=irfs, type="median", labels=vl, save=FALSE, bands=NULL, grid=TRUE, bw=TRUE)

Below is a STATA code that I wrote.

*** SVAR Sign-restrictions. 
clear all

capture program drop subprogram
program subprogram 
    mata mata clear
    putmata ybly yblz zbly zblz
    mata A1=ybly[1],yblz[1]\zbly[1],zblz[1]
    mata sigma7=7600000, 7.76655  \7.76655, .007221
    mata L=cholesky(sigma7)
    mata A=rnormal(1,1,0,1),rnormal(1,1,0,1)\rnormal(1,1,0,1),rnormal(1,1,0,1) 
    mata Q=R=.
    mata qrd(A,Q,R)
    mata pi_svar=L*Q
    mata pi_svar1=pi_svar[1,1]
    mata pi_svar2=pi_svar[2,1]
    forvalues i=2(1)$obs {
        mata: pi_svar=A1*pi_svar    
        mata: pi_svar1=pi_svar1\pi_svar[1,1]
        mata: pi_svar2=pi_svar2\pi_svar[2,1]
    }
    getmata pi_svar1 pi_svar2

    qui tsset t
    qui keep if t<=120
    qui keep t pi_svar1 pi_svar2
    qui rename (pi_svar1 pi_svar2)(p1 p2)
    qui gen d=0
    qui replace d=1 if p1[1]<0&p2[1]>0&p1[2]<0&p2[2]>0&p1[3]<0&p2[3]>0&p1[4]<0&p2[4]>0&p1[5]<0&p2[5]>0&p1[6]<0&p2[6]>0
end 

import delimited "https://raw.githubusercontent.com/jayjeo/public/master/LaborShortage/SVARdata_seasonadjusted2_practice.csv", varnames(1) clear 
global obs=120
gen t=_n
tsset t
sureg (y L.y L.z)(z L.y L.z)
gen ybly=_b[y:L.y]
gen yblz=_b[y:L.z]
gen zbly=_b[z:L.y]
gen zblz=_b[z:L.z]

reg y L.y L.z
predict ey, residual
reg z L.y L.z
predict ez, residual
corr ey ez, covariance 
*mata sigma7=7600000, 7.76655  \7.76655, .007221

forvalues i=1(1)100 { 
    preserve         
        qui subprogram
        qui save pp`i', replace 
        scalar process=`i'
        di process
    restore 
}

use pp1, clear
forvalues i=2(1)100 { 
    append using pp`i'
}
drop if d==0
collapse (median) p1 p2, by(t)

tsset t
tsline p1, name(p1) yline(0)
tsline p2, name(p2) yline(0)
graph combine p1 p2, rows(2)

Below figure is the result using the built-in R program by Danne, Christian. "VARsignR"

Below figure is the result using the STATA code that I wrote (which does not match exactly).