Matrix is singular

[oriti]

Dear all,

I have a simulation which in some point it stopped because I get singular matrix.

I read Rick wonderful post about the rareness of singularity ( What is the chance that a random matrix is singular? - The DO Loop)

so it really puzzled me that I get singular matrix.

After discussion in the forum

(checking Rick assumption about the var-covar matrix in the program which has one small eigenvalue)

and thinking again about the model I conclude that my problem is technical.

So..

How can I force SAS to inverse a matrix?  The chance that the det=0 is zero..

How can I continue with the simulation without that singular matrix? and how can I print the singular matrix or X data?

Thank you all in advance,

I enclose the code

proc iml;

n=10000;

b=j(6,n);

     do j=1 to n;

       xtl=j(240, 1);

       xsl=j(240, 1);

       y=j(240, 1);

            xtl[1,1]=0.000844;

            xsl[1,1] =2.78896;

            y[1,1] =-0.064516129;

                 a=0.01095 ;*first regression*;

                  k=0.41450;

                 ka=0.01349;

                 kb=0.05634;

                 kc=-0.00228 ;

                 kd=-40.63877;

                 c=-0.00678; *second regression*;

                 e=0.07520;

                 d=0.94798;

                 h=-1.73592;

                 m=-0.00013407; *third regression*;

                 l=0.00144;

                 q=-0.00046335;

                 r=0.99109;

  /** specify the mean and covariance of the population errors **/

  Mean = {0, 0, 0};

        Cov = {0.0047067768 0.0008647583 0.0000199767,

    0.0008647583 0.0081750097 0.0000391468,

   0.0000199767 0.0000391468 0.0000022146};

  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*xsl[(i-1),1]+kb*xtl[(i-1),1]+kc*xsl[(i-1),1]##2+kd*xtl[(i-1),1]##2+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];

    end;

    ylag=lag(y);

    xsl2=xsl##2;

    xtl2=xtl##2;

     /* Step 1: Compute X`X and X`y */

    x = j(nrow(y), 1, 1) || ylag || xsl|| xtl || xsl2||xtl2;            /* 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;

xpxi = inv(xpx); /* form inverse crossproducts */

b[,j] = xpxi * xpy; /* solve for parameter estimates */

  *b[,j]= solve(xpx, xpy);/* another way which is more efficient*/;

  end;

  bt=b`;

  varNames = "b0":"b5";

  experiment="Sim1":"Sim10000";

  print bt[colname=varNames

                rowname=experiment

                label=" Model"];

  quit;

[PaigeMiller]

Sure would be helpful to know which step of this long program produced the error message that a matrix was singular.

Your covariance matrix has a determinant that is close to zero, that may be part of the problem.

[Rick_SAS · Accepted solution]

You are misreading my blog post.  In my post I generated matrices with elements drawn from a uniform or normal distribution, and showed that the (mathematical) probability of getting a singular matrix is zero (although in finite-precision arithmetic the probability in nonzero).

That is not what you are doing here. You are simulating a timne series model and using a (nearly singular) covariance matrix to generate correlated random errors for the sinulation.

Your problem is blowing up because your time series model is diverging to infinity. IFirst, use a nonzero seed so that you can reproduce the results consistently:

call randseed(1);

Next, insert the following code before you compute the inverse:

if abs(det(xpx))<1e-6 then do;
   z = ynew || xnew;
   create SINGULAR from z[c={Y Int Ylag xsl xtl xsl2 xtl2}];
   append from z;
   close SINGULAR;

   create U from U[c={U1 U2 U3}];
   append from U;
   close U;
end;

If you look at the SINGULAR data set that is created, you will see that the Y, YLag, XSL, and XTL variables are all marching off to infinity.

I don't know how you want to deal with this problem in your model. If you want to ignore the trajectories that diverge, you can use IF(abs(DET(xpx)))>0  to protect the code that computes the parameter estimates. You will have missing values for the bad trajectories.

[oriti]

Thank you Rick for your detailed answer!

I really appreciate your help!

It's really interesting and helpful what you wrote.

When I looked at 'SINGULAR' with your help..I saw that the last 20 simulated observations were marching off to infinity,

So..when I simulated 220 instead of 240 obs everything went well..I got 10000 results as I wanted.

Well, I have to think why and how just the last 20 out of 240 simulated results caused the problem?

Is it technical problem?

I'll appreciate your opinion.

Thank you,

Orit

[Rick_SAS]

It's not a technical problem, it is a mathematical problem. An iterated dynamical system will have a deterministic long-term behavior (bounded or diverging). When you add a stochastic component, you introduce a probability that the system will diverge, even if the deterministic system remains bounded.  You can read the main ideas here: Convergence or divergence? A simple iteration with a random component - The DO Loop

A simpler model is a random walk game. Let's play a game where you flip a coin 20 times. You move 1 meter forward whenever the coin is heads; you move 1 meter backwards whenever the coin is tails. Now suppose you begin each game 10 meters aways from the edge of a cliff which has a deadly 1000-foot drop.  You will probably be able to play the game for a long time, but the more you play the more you risk encountering an unlucky string of heads that plunge you off the cliff!

[oriti]

Really interesting!

Thank you Rick!