Dear all, I wrote a simulation that about ten percent of the time end with a matrix data that contains one 'zero' column. I identify those cases with if abs(det(xpx))=0. I clean the zero column and I try to calculate the OLS coefficients. But unfortunately I get a matrix without the OLS coefficients. The program works but I don't get the 10 percent results matrix with data, even though I record the data files and there should be coefficients for the OLS. (data with the zero columns). I don't understand why. I enclose the program. You can 'run' it and see the results. 1. Can someone understand why I get matrix without data for the 10 percent cases of the simulation? 2. Any suggestions how to correct that? 3. Maybe there is more efficient way to write this program that at the end I can get one matrix with also the 10 percent result, Ideas are welcome.. Thanks in advance! Orit proc iml; n=100; b=j(6,n); bpsl=j(5,n); bnsl=j(5,n); bptl=j(5,n); bntl=j(5,n); do j=1 to n; xtl=j(240, 1); xsl=j(240, 1); ptl=j(240, 1); psl=j(240, 1); ntl=j(240, 1); nsl=j(240, 1); y=j(240, 1); xtl[1,1]=-1.73754E-16;*/SAS avg obs*/; xsl[1,1] =-1.12557E-15; ptl[1,1]=0.0051111; psl[1,1] =0.2604981; ntl[1,1]=-0.0051111; nsl[1,1] =-0.2604981; y[1,1] =0.9783607; a=0.18935 ;*first regression*; k=0.83843; ka=0.17226; kb=-0.39250; kc=-0.21518; kd=5.16021; c=-0.03537; *second regression*; e=0.02979; d=0.94756; h=-2.04415; m=-0.00034939; *third regression*; l=0.00023025; q=-0.00045077; r=0.98906; /** specify the mean and covariance of the population **/ Mean = {0, 0, 0}; Cov = {0.0467930604 0.0024632274 0.0000520215, 0.0024632274 0.0080371426 0.0000384574, 0.0000520215 0.0000384574 0.0000022165}; NumSamples = 240; call randseed(0); /** set seed for the RandNormal module **/ U = RandNormal(NumSamples, Mean, Cov); do i=2 to 240; Y[i,1]=a+k*y[(i-1),1]+ka*psl[(i-1),1]+kb*ptl[(i-1),1]+kc*nsl[(i-1),1]+kd*ntl[(i-1),1]+u[i,1]; xsl[i,1]=c+e*y[(i-1),1]+d*xsl[(i-1),1]+h*xtl[(i-1),1]+u[i,2]; xtl[i,1]=m+l*y[(i-1),1]+q*xsl[(i-1),1]+r*xtl[(i-1),1]+u[i,3]; if xsl[i,1]>=0 then do; psl[i,1]=xsl[i,1]; nsl[i,1]=0; end; if xsl[i,1]<0 then do; psl[i,1]=0; nsl[i,1]=xsl[i,1]; end; if xtl[i,1]>=0 then do; ptl[i,1]=xtl[i,1]; ntl[i,1]=0; end; if xtl[i,1]<0 then do; ptl[i,1]=0; ntl[i,1]=xtl[i,1]; end; end; ylag=lag(y); /* Step 1: Compute X`X and X`y */ x = j(nrow(y), 1, 1) || ylag || psl|| ptl ||nsl||ntl; /* add intercept column */ v = {1}; /* specify rows to exclude */ idx = setdif(1:nrow(x), v); /* start with 1:240, exclude values in v */ xnew = x[idx, ]; /* extract submatrix */ xpx = xnew`* xnew; /* cross-products */ v = {1}; /* specify rows to exclude */ idy = setdif(1:nrow(y), v); /* start with 1:5, exclude values in v */ ynew = y[idx, ]; /* extract submatrix*/ xpy = xnew` * ynew; if abs(det(xpx))=0 then do; z1 = ynew || xnew; create SINGULAR from z1[c={Y Ylag psl ptl nsl ntl}]; append from z1; close SINGULAR; end; if abs(det(xpx))=0 then do; if ssq(psl)=0 then do; z = j(nrow(y), 1, 1) || ylag || ptl ||nsl||ntl; /* add intercept column */ v = {1}; /* specify rows to exclude */ idx = setdif(1:nrow(z), v); /* start with 1:240, exclude values in v */ znew = z[idx, ]; /* extract submatrix */ zpz = znew`* znew; /* cross-products */ v = {1}; /* specify rows to exclude */ idy = setdif(1:nrow(y), v); /* start with 1:5, exclude values in v */ ynew = y[idx, ]; /* extract submatrix*/ zpy = znew` * ynew; bpsl[,j]= solve(zpz, xpy); end; if ssq(ptl)=0 then do; z = j(nrow(y), 1, 1) || ylag || psl||nsl||ntl; /* add intercept column */ v = {1}; /* specify rows to exclude */ idx = setdif(1:nrow(z), v); /* start with 1:240, exclude values in v */ znew = z[idx, ]; /* extract submatrix */ zpz = znew`* znew; /* cross-products */ v = {1}; /* specify rows to exclude */ idy = setdif(1:nrow(y), v); /* start with 1:5, exclude values in v */ ynew = y[idx, ]; /* extract submatrix*/ zpy = znew` * ynew; bptl[,j]= solve(zpz, xpy); end; if ssq(nsl)=0 then do; z = j(nrow(y), 1, 1) || ylag || psl|| ptl ||ntl; /* add intercept column */ v = {1}; /* specify rows to exclude */ idx = setdif(1:nrow(z), v); /* start with 1:240, exclude values in v */ znew = x[idx, ]; /* extract submatrix */ zpz = znew`* xnew; /* cross-products */ v = {1}; /* specify rows to exclude */ idy = setdif(1:nrow(y), v); /* start with 1:5, exclude values in v */ ynew = y[idx, ]; /* extract submatrix*/ zpy = znew` * ynew; bnsl[,j]= solve(zpz, xpy); end; if ssq(ntl)=0 then do; z = j(nrow(y), 1, 1) || ylag || psl|| ptl ||nsl; /* add intercept column */ v = {1}; /* specify rows to exclude */ idx = setdif(1:nrow(x), v); /* start with 1:240, exclude values in v */ znew = x[idx, ]; /* extract submatrix */ zpz = znew`* znew; /* cross-products */ v = {1}; /* specify rows to exclude */ idy = setdif(1:nrow(y), v); /* start with 1:5, exclude values in v */ ynew = y[idx, ]; /* extract submatrix*/ zpy = znew` * ynew; bntl[,j]= solve(zpz, xpy); end; end; IF(abs(DET(xpx)))>0 then do; *xpxi = inv(xpx); /* form inverse crossproducts */ *b[,j] = xpxi * xpy; /* solve for parameter estimates */ *b= solve(xpx, xpy); *print b; b[,j]= solve(xpx, xpy); end; end; bt=b`; varNames = "b0":"b5"; experiment="Sim1":"Sim100"; print bt[colname=varNames rowname=experiment label="Expanded Model"]; bpslt=bpsl`; varNames = "b0":"b4"; experiment="Sim1":"Sim100"; print bpslt[colname=varNames rowname=experiment label="Expanded Model"]; bnslt=bnsl`; varNames = "b0":"b4"; experiment="Sim1":"Sim100"; print bnslt[colname=varNames rowname=experiment label="Expanded Model"]; bptlt=bptl`; varNames = "b0":"b4"; experiment="Sim1":"Sim100"; print bptlt[colname=varNames rowname=experiment label="Expanded Model"]; bntlt=bntl`; varNames = "b0":"b4"; experiment="Sim1":"Sim100"; print bntlt[colname=varNames rowname=experiment label="Expanded Model"]; quit; proc print data=singular; run;
... View more