*! Mata code for splagvar *! By P. Wilner Jeanty *! Date: March 2009 version 10.1 mata: mata set matastrict on void splagvar_lagmyvar(string scalar xvars1, string scalar xvars2 , string scalar ftouse) { transmorphic fh, wyindxes, nwy, stubwy, wy, wxindxes, stubwx, w2xindxes, stubw2x, w3xindxes, stubw3x, nwx, nw2x, nw3x real scalar n1, n2, m1, m2, i, j real matrix ytolag, xtolag, wx, w2x, w3x external real scalar splagvar_Ord external string scalar splagvar_TOMAT1 external real matrix splagvar_w // To make the weights matrix global st_view(ytolag, ., tokens(xvars1), ftouse) st_view(xtolag, ., tokens(xvars2), ftouse) if (splagvar_TOMAT1=="Mata") { if (!fileexists(st_local("wname"))) { "" errprintf("File %s not found\n", st_local("wname")) exit(601) } else { fh = fopen(st_local("wname"), "r") splagvar_w=fgetmatrix(fh) fclose(fh) } } else { splagvar_w=st_matrix(st_local("wname")) if (cols(splagvar_w)==0 | rows(splagvar_w)==0) { "" errprintf("Matrix %s not found\n", st_local("wname")) exit(601) } } n1=cols(ytolag); m1=rows(ytolag) n2=cols(xtolag); m2=rows(xtolag) if (n1!=0) { if (rows(splagvar_w)!=m1) { "" errprintf("Number of observations (%f) not conformable with weights matrix size (%f by %f)\n", m1, rows(splagvar_w), rows(splagvar_w)) exit(3200) } wyindxes=J(1,n1,.) stubwy="weird_wy" for (i=1; i<=n1; i++) { nwy=stubwy+strofreal(i) wy=splagvar_w*ytolag[.,i] wyindxes[i]=st_addvar("double", nwy) st_store(., nwy, wy) } } if (n2!=0) { if (rows(splagvar_w)!=m2) { "" errprintf("Number of observations (%f) not conformable with weights matrix size (%f by %f)\n", m2, rows(splagvar_w), rows(splagvar_w)) exit(3200) } wxindxes=J(1,n2,.) stubwx="weird_wx" if (splagvar_Ord==2 | splagvar_Ord==3) { w2xindxes=J(1,n2,.) stubw2x="weird_w2x" if (splagvar_Ord==3) { w3xindxes=J(1,n2,.) stubw3x="weird_w3x" } } for (j=1; j<=n2; j++) { nwx=stubwx+strofreal(j) wx=splagvar_w*xtolag[.,j] wxindxes[j]=st_addvar("double", nwx) st_store(., nwx, wx) if (splagvar_Ord==2 | splagvar_Ord==3) { nw2x=stubw2x+strofreal(j) w2x=splagvar_w*wx w2xindxes[j]=st_addvar("double", nw2x) st_store(., nw2x, w2x) if (splagvar_Ord==3) { nw3x=stubw3x+strofreal(j) w3x=splagvar_w*w2x w3xindxes[j]=st_addvar("double", nw3x) st_store(., nw3x, w3x) } } } } } void splagvar_CalcMoran(string scalar zvar, string scalar ztouse) { real scalar vz2, I, n, C, C2, swyy, ei_n, ei_r, A, B, D1, D2, vn, vr, ZI_n, ZI_r, p_n, p_r, E, i, MI real colvector z, zd, zd2, zd4, B1 real matrix sww external real matrix splagvar_w // To use the weights matrix retrieved by the previous function st_view(z=., ., zvar, ztouse) zd=z:-mean(z,1) zd2=zd:^2 vz2=sum(zd2) // vz=crossdev(z,meanz, z,meanz) would do as well swyy=sum(zd'*splagvar_w*zd) C=sum(splagvar_w) C2=C^2 n=rows(z) I=(n/C)*(swyy/vz2) // Moran's I st_numscalar("r(I)", I) ei_n=ei_r=-1/(n-1) // Mean // Calculating the variances zd4=zd2:^2 sww=(splagvar_w+splagvar_w'):^2 A=0.5*sum(sww) // A B1=J(n,1,.) for (i=1; i<=n; i++) { B1[i]=(rowsum(splagvar_w[i,.])+ colsum(splagvar_w[.,i]))^2 } B=sum(B1) D1=((n^2)-3*n+3)*A - n*B + 3*C2 D2=((n^2)-n)*A - 2*n*B + 6*C2 E=mean(zd4,1)/(mean(zd2,1))^2 // Variance under normal approximation assumption vn=(((n^2)*A - n*B + 3*C2)/(((n^2)-1)*C2))-(ei_n)^2 ZI_n=(I-ei_n)/sqrt(vn) p_n=(1-normal(abs(ZI_n)))*2 // Variance under randomization assumption vr=((n*D1-E*D2)/((n-1)*(n-2)*(n-3)*C2))-(ei_r)^2 ZI_r=(I-ei_r)/sqrt(vr) p_r=(1-normal(abs(ZI_r)))*2 MI=(I,I\ei_n,ei_r\sqrt(vn),sqrt(vr)\ZI_n,ZI_r\p_n,p_r) st_matrix("morstat", MI) } mata mlib create lsplagvar, dir(PLUS) replace mata mlib add lsplagvar splagvar_*(), dir(PLUS) mata mlib index end