BookmarkSubscribeRSS Feed
🔒 This topic is solved and locked. Need further help from the community? Please sign in and ask a new question.
oriti
Fluorite | Level 6

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;

1 ACCEPTED SOLUTION

Accepted Solutions
Rick_SAS
SAS Super FREQ

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.

View solution in original post

5 REPLIES 5
PaigeMiller
Diamond | Level 26

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.

--
Paige Miller
Rick_SAS
SAS Super FREQ

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
Fluorite | Level 6

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
SAS Super FREQ

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
Fluorite | Level 6

Really interesting!

Thank you Rick!

sas-innovate-2024.png

Join us for SAS Innovate April 16-19 at the Aria in Las Vegas. Bring the team and save big with our group pricing for a limited time only.

Pre-conference courses and tutorials are filling up fast and are always a sellout. Register today to reserve your seat.

 

Register now!

Multiple Linear Regression in SAS

Learn how to run multiple linear regression models with and without interactions, presented by SAS user Alex Chaplin.

Find more tutorials on the SAS Users YouTube channel.

From The DO Loop
Want more? Visit our blog for more articles like these.
Discussion stats
  • 5 replies
  • 3903 views
  • 7 likes
  • 3 in conversation