The left hand side is referring only to the diagonal elements, while the right hand side contains n^2 elements this is the reason for the non-conforming error message.  The solution is to extract only the diagonal element of AR using the vecdiag function.   y[do(1,n#n,n+1)] = 1 - vecdiag(AR); ... View more

I think the following is what you need:   y[do(1,n#n,n+1)] = vecdiag(y)^=.; which sets the diagonal to 1 wherever it is non-missing after the division, otherwise it is set to zero.   ... View more

I think it would be best to create a single matrix large enough to hold all 500 results.  Then the matrix is filled up one row at a time for each iteration of the loop.   Afterwards results can be extracted by row sub setting.  The syntax would look like this:   n = sqrt(ncol(x)); want = j(nobs, n); do i = 1 to nobs; temp = shape(x[i,], n); want[i,] = y[i,]*temp; end; w2 = want[2,]; print w2; ... View more

If you move PRINT LEFT inside the DO loop then it will execute 3 times, once for each group. ... View more

want = y[i,]*temp; would be more simple.   ... View more

Just count the number of rows that have been read into the matrix x.   To store the matrices you could use indirect assignment with VALSET and VALUE.   proc iml; use have(keep = pr_st:); read all var _num_ into x [colname=vn]; rcidx = t( num(scan(vn,3,'t_')) // num(scan(vn,4,'t_'))); n = nrow(x); pn = 'p1' : 'p'+strip(char(n)); do i = 1 to n; call valset( pn[i], full ( t( x[i,]) || rcidx ) ); end; q = value(pn[1]); print q; quit; ... View more

I am not sure that you can use the SHAPE function if the data is like the example you give above, since the 9 variables are not in row major order.  For example the 2nd variable gives the value for the 1st row and 3rd column. ... View more

Have a look at the set functions XSECT and ELEMENT.   x = xsect(A, B);   will give a row vector with the values that are common to both A and B.   idx = loc(element(B, x));   will give a vector with the values of the loop variable i that correspond to a "yes". ... View more

I don't think that there is any way to recover from the singular matrix error with the INV function, so the initial test must be done if you go this route.   Please check the 3rd code example above as this is not doing what you think.  Specifically, the statement  "if rank(MATRIX)=&SIZE." will always evaluate as false, so GINV is used for all the simulations, and the timings are more a reflection of how much slower GINV is when compared to INV.  The RANK function in IML returns the rank order of the elements within a matrix.  To get the matrix rank, you will have to use one of the methods in Rick's article,  see link in the original topic which contains code for a RankMatrix function, however this will be expensive as the preferred methods involve the calculation of an SVD or a generalized inverse. ... View more