Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

🔒 This topic is **solved** and **locked**.
Need further help from the community? Please
sign in and ask a **new** question.

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

Posted 02-03-2014 04:49 AM
(3960 views)

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

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.

5 REPLIES 5

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

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

Paige Miller

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

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!

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

Really interesting!

Thank you Rick!

**Don't miss out on SAS Innovate - Register now for the FREE Livestream!**

Can't make it to Vegas? No problem! Watch our general sessions LIVE or on-demand starting April 17th. Hear from SAS execs, best-selling author Adam Grant, Hot Ones host Sean Evans, top tech journalist Kara Swisher, AI expert Cassie Kozyrkov, and the mind-blowing dance crew iLuminate! Plus, get access to over 20 breakout sessions.

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.