OPTIONS MEMORY=32; set s = 4; ? サンプル・サイズを変えるときはここをいじる set n = 2*s-1; set L = n*n; ? total number of grid smpl 1,L; ? numbering grid for 1 to L mform(type=gen, nrow=L,ncol=2)xm = 0; do j = 1 to n by 1 ; do i = 1 to n by 1 ; set k = i+n*(j-1); set xm(k,1) = j; set xm(k,2) = i; enddo; enddo; ? make a Rook type contiguity matrix mform(type=gen, nrow=L,ncol=L)D = 0; mform(type=gen, nrow=L,ncol=L)CM = 0; do q = 1 to L by 1 ; do p = 1 to L by 1 ; set D(p,q) = abs(xm(p,1) - xm(q,1)) + abs(xm(p,2) - xm(q,2)); if D(p,q) = 1; then; set CM(p,q) = 1; enddo; enddo; ? row-standardized spatial weight matrix mat link = CM*C; unmake link links; mat W = (diag(CM*C)"CM); set trial = 1000; ? the number of trial mmake rho -0.9 -0.7 -0.5 -0.3 -0.1 0.1 0.3 0.5 0.7 0.9; unmake rho rho1-rho10; mform(type=gen,nrow=trial,ncol=10)ZMIMAT = 0; set beta1 = 1; set beta2 = 1; set sigma2 = 1; x1 = C; trend x2; mmake X x1 x2; mmake beta beta1 beta2; mat coefn=ncol(X); dot(index=i) 1-10; mat Amat = (ident(L) - rho.*W)"; dot(index=j) 1-1000; ? trial rand(stdev=sigma2)epsilon; mat ymat = Amat*X*beta + Amat*epsilon; unmake ymat y; olsq(silent) y X1 X2; mat e = @res; mat eWe = e'W*e; mat MI = eWe/e'e ; ? Moran's I mat M = ident(L) - x*(x'x)"x'; ? Projection Matrix mat EI = tr(M*W)/(L-coefn); set df = (L-coefn)*(L-coefn+2); mat tr1 = tr(M*W*M*W'); mat tr2 = tr((M*W)*(M*W)); mat tr3 = (tr(M*W))^2; mat VI = ( tr1 + tr2 + tr3 )/df - EI^2 ; set zMI. = (MI - EI)/sqrt(VI); set ZMIMAT(j,i) = zmi.; enddot; enddot; smpl 1,trial; unmake ZMIMAT ZMIS1-ZMIS10; msd ZMIS1 -ZMIS10; dot 1-10; print rho.; hist ZMIS.; enddot; ? write(file="ZMI.xls")ZMIS1-ZMIS10; date; end;