capture program drop svarsirf *! svarirf v1.0.0 GImpavido (gimpavido@imf.org) and Li Tang (ltang@imf.org) 25jul2015 program svarsirf, rclass version 10.0 syntax , [STeps(integer 12) NOSE] ************************* Perform user error checks **************************** if "`e(cmd)'" != "svar" { di as err "svarsirf only works after svar" exit 198 } else { if "`e(lrmodel)'" != "" { di as err "svarsirf only works for short-run model" exit 198 } else { if "`e(cns_a)'`e(cns_b)'" == "" { di as error "cannot find constraints" /* */ "on A and B matrices specified in svar" exit 198 } } } ************************* pass ingredients to locals *************************** local names "`e(depvar_var)'" local namelist "A0est Best bvarest Vvarest Sigma_AB" mat A0est = e(A) mat Best = e(B) mat bvarest = e(b_var) mat Vvarest = e(V_var) mat Sigma_AB = e(V) ***************************** run MATA function ******************************** mata : m_svarsirf("`namelist'", "`names'", "`steps'", "`nose'") ***************************** add matrix headers ******************************* local head = head matrix colnames SIRF = `head' if "`nose'" == "" { matrix colnames SIRFSE = `head' } *********************** return MATA results to STATA *************************** return local head = head return matrix SIRF = SIRF if "`nose'" == "" { return matrix SIRFSE = SIRFSE } end ********************************* start MATA *********************************** capture mata mata drop m_svarsirf() version 10.0 mata: mata set matastrict on // m_svarsirf v1.0.0 GImpavido (gimpavido@imf.org) and Li Tang (ltang@imf.org) 25jul2015 void m_svarsirf( /// string scalar namelist, string scalar names, | /// string scalar steps, string scalar nose) { string rowvector AAA, NAM, head, nse real rowvector b_var real matrix A0, B, V_var, BB, AA, Psr, SIRF, Q1, Q2, /// SSnoc, SSnor, SSr, SS, J, M, MM1, MM2, /// SIGMA_W, SIGMA0, SIGMA_AB, G0, /// Q4, SIGMAH, SIRFSE, STP numeric scalar i, j, s, lags, stp, z, cm, cn, rm, rn pointer rowvector A, PHI, THETA, G, MP, MPP /* start pass the ingredients */ AAA = tokens(namelist) // list of matrices NAM = tokens(names) // list of variables in SVAR nse = tokens(nose) // option for SE stp = strtoreal(steps) // Step for SIRF (default =12) A0 = st_matrix(AAA[1]) B = st_matrix(AAA[2]) b_var = st_matrix(AAA[3]) V_var = st_matrix(AAA[4]) SIGMA_AB = st_matrix(AAA[5]) lags = (cols(b_var)-rows(B))/(rows(B)*cols(B)) /* end pass the ingredients */ /* start the calculation of SIRF */ /* rearrange b_var and eliminate cols for constant */ BB = J(rows(B),cols(B)*lags,0) for (i = 1; i<=rows(B); i++) { cm = 1+(i-1)*(rows(B)*lags+1) cn = cm + rows(B)*lags - 1 BB[i,.] = b_var[|.,cm\.,cn|] } /* rearrange columns of BB (i) in lag order */ AA = J(rows(B),lags*cols(B),.) for (s = 1; s<=lags; s++) { for (i = 1; i<=cols(B); i++) { cm = i+(s-1)*rows(B) cn = s+(i-1)*lags AA[.,cm] = BB[.,cn] } } /* create structure for the lags of the VAR containing n=steps matrices of dimension (rows(B) x cols(B)) by using pointers. These are: a=1 matrix A0; p=lag submatrices of BB (A1,....,Ap); and l = steps - p -1 matrices of zeros */ A = J(1,stp+1,NULL) /* structure of the A0... Ap matrices */ A[1] = &(A0) for (s = 1; s <= stp; s++) { if (s <= lags) A[s+1] = &AA[|.,1+rows(B)*(s-1)\.,rows(B)*s|] else A[s+1] = &J(rows(B),cols(B),0) } Psr = luinv(*A[1])*B PHI = J(1,stp+1,NULL) // structure of the MA representation THETA = J(1,stp+1,NULL) // structure of the SIRF PHI[1] = &(I(rows(B))) // contemporaneous MA THETA[1] = &(I(rows(B))*Psr) // contemporaneous SIRF// s=1,...,step MA and SIRF for (s = 1; s <= stp; s++) { PHI[s+1] = &(J(rows(B),cols(B),0)) for (j = 1; j <= s; j++) { PHI[s+1] = &(*PHI[s+1]+*PHI[s-j+1]**A[j+1]) } THETA[s+1] = &(*PHI[s+1]*Psr) } /* gen SIRF matrix */ SIRF = J(stp+1,rows(B)*cols(B),.) for (s = 1; s <= stp+1; s++) { for (j = 1; j<=cols(B); j++) { for (i = 1; i<=rows(B); i++) { z = (i-1)*rows(B)+j SIRF[s,z] = (*THETA[s])[j,i] } } } /* add first column step */ STP = J(stp+1,1,.) for (s = 1; s<=stp+1; s++) { STP[s,1] = s-1 } SIRF = STP,SIRF /* end the calculation of SIRF */ /* start header matrix */ /* gen matrix header with the names of "impulse_response" */ head = J(1,rows(B)*cols(B)+1,"") head[1,1] = "Step" for (j = 1; j <= cols(B); j++) { for (i = 1; i <= rows(B); i++) { cm = (i-1)*rows(B)+j+1 head[1,cm] = NAM[1,i]+":"+NAM[1,j] } } /* end header matrix */ /* start the calculation of asymptotic SE */ if (nose == "") { /* eliminate cols of const in V_var */ SSnoc = J(lags*rows(B)^2+rows(B),lags*cols(B)^2,.) for (i=1; i <= rows(B); i++) { cm = 1+(i-1)*(rows(B)*lags+1) cn = cm + rows(B)*lags - 1 SSnoc[|.,1+(i-1)*rows(B)*lags\.,i*rows(B)*lags|] = V_var[|.,cm\.,cn|] } /* eliminate rows of const in V_var */ SSnor = J(lags*rows(B)^2,lags*cols(B)^2,.) for (j=1; j <= cols(B); j++) { rm = 1+(j-1)*(rows(B)*lags+1) rn = rm + rows(B)*lags - 1 SSnor[|1+(j-1)*rows(B)*lags,.\j*rows(B)*lags,.|] = SSnoc[|rm,.\rn,.|] } /* rearrange columns (i) in lag order */ z = 1 SSr = J(lags*rows(B)^2,lags*cols(B)^2,.) for (i = 1; i <= rows(B); i++) { for (j = 1; j <= cols(B); j++) { for (s = 1; s <= lags; s++) { SSr[.,i+(j-1)*cols(B)+(s-1)*cols(B)^2] = SSnor[.,z] /* swaps cols */ z++ } } } /* rearrange rows (j) in lag order */ z = 1 SS = J(lags*rows(B)^2,lags*cols(B)^2,.) for (i = 1; i <= rows(B); i++) { for (j = 1; j <= cols(B); j++) { for (s = 1; s <= lags; s++) { SS[i+(j-1)*cols(B)+(s-1)*cols(B)^2,.] = SSr[z,.] /* swaps rows */ z++ } } } /* J matrix */ J = I(rows(B)) for (j = 2; j <= lags; j++) { J = J,J(rows(B),cols(B),0) } /* M matrix */ MM1 = *A[2] for (j = 2; j <= lags; j++) { MM1 = MM1,*A[j+1] } MM2 = I(rows(B)*(lags-1)),J(rows(B)*(lags-1),cols(B),0) M = MM1\MM2 /* SIGMA0 matrix */ Q1 = (I(rows(B))#luinv(B),-luinv(Psr')#luinv(B)) Q2 = Psr'#Psr SIGMA_W = Q1*SIGMA_AB*Q1' SIGMA0 = Q2*SIGMA_W*Q2' /* G_i matrix */ G0 = J(rows(B)^2, rows(B)^2*lags, 0) G = J(1,stp+1,NULL) // structure of the GI G[1] = &(G0) MP = J(1,stp+1,NULL) // power matrix of M MPP = J(1,stp+1,NULL) // power matrix of M' MP[1] = &(I(rows(B)*lags)) MPP[1] = &(I(rows(B)*lags)) for (i = 2; i <= stp+1; i++) { MP[i] = &(*MP[i-1]*M) MPP[i] = &(*MPP[i-1]*M') } for (s = 1; s <= stp; s++) { G[s+1] = &(J(rows(B)^2, rows(B)^2*lags, 0)) for (j = 1; j <= s; j++) { G[s+1] = &(*G[s+1]+(Psr'*J**MPP[s-j+1])#(J**MP[j]*J')) } } /* gen SIRF SE matrix */ for (s = 0; s <= stp; s++) { Q4 = I(rows(B))#(J**MP[s+1]*J') SIGMAH = (*G[s+1])*SS*(*G[s+1])'+Q4*SIGMA0*Q4' if (s == 0) SIRFSE = (diagonal(cholesky(diag(diagonal(SIGMAH)))))' else SIRFSE = SIRFSE \ (diagonal(cholesky(diag(diagonal(SIGMAH)))))' } /* add first column step */ SIRFSE = STP,SIRFSE } // this is the end of the option nose /* end the calculation of asymptotic SE */ /* send MATA matrices to STATA */ if (nose == "") st_matrix("SIRFSE", SIRFSE) st_strscalar("head", invtokens(head, " ")) st_matrix("SIRF", SIRF) } end *********************************** end MATA ***********************************