BookmarkSubscribeRSS Feed
manto
Fluorite | Level 6

Hi all,

I am doing some simulation work using SAS/IML. Some of my numbers are very small (elements of a Park's type VCV matrix constructed from regression residuals). I use the ROOT function to get the Cholesky decomposition of a matrix that is positive definite by construction (VCV matrix), but it returns the error message that the matrix provided was not positive definite. I suspect a rounding error but I do not know how to solve it. I got the same problem before with Stata (Mata), and SAS/IML helped counter it. But now I am using a new data set, and facing again the problem with SAS.

The portion of the code up to the Cholesky decomposition is provided below:

proc iml;

start program;

use Newdata.AD;

read all var {Germany_aid} into y;

read all var {GDP_per_capita} into x;

read all var ('c1':'c77') into s;

read all var ('t1':'t25') into years;

***THIS SECTION CREATES THE OMEGA MATRIX AND;

***INITIALIZES THE MONTE CARLO PARAMETERS;

n=50;

periods=20;

t=periods;

nt=n*periods;

yi=j(nt,1,0);

xi=j(nt,1,0);

si=j(nt,n,0);

yearsi=j(nt,(periods-1),0);

omega=j(nt,nt,0);

ones=j(nt,1,1);

r2=0;

yy=j(nt,1,0);

xx=j(nt,1,0);

phitotal=j(n,n,0);

rhohattotal=0;

yy=0;

xx=0;

r2=0;

do k=1 to (25-periods+1);

   do j=1 to n;

      yi[1+(j-1)*periods:j*periods]=y[k+(j-1)*25:k+(periods-1)+(j-1)*25];

      xi[1+(j-1)*periods:j*periods]=x[k+(j-1)*25:k+(periods-1)+(j-1)*25];

      si[1+(j-1)*periods:j*periods,]=s[k+(j-1)*25:k+(periods-1)+(j-1)*25,1:n];

      yearsi[1+(j-1)*periods:j*periods,1:(periods-1)]=years[k+(j-1)*25:k+(periods-1)+(j-1)*25,k+1:k+(periods-1)];

   end;

zi=si||yearsi||xi;

ei=yi-zi*inv(zi`*zi)*zi`*yi;

a=i(nt)-(1/nt)*ones*ones`;

r2i=1-((ei`*ei)/(yi`*a*yi));

eols1=j(n*(periods-1),1,0);

eols2=j(n*(periods-1),1,0);

pp=j(periods,periods,0);

bols=inv(zi`*zi)*zi`*yi;

do i=1 to n;

   eols1[1+((i-1)*(periods-1)):(i*(periods-1))]=ei[(2+((i-1)*periods)):i*periods];

   eols2[1+((i-1)*(periods-1)):(i*(periods-1))]=ei[(1+(i-1)*periods):((i*periods)-1)];

end;

***NOTE THAT I AM USING A DIFFERENT FORMULA;

***THE GREENE FORMULA;

***TO CALCULATE RHOHAT;

rhohati=(eols1`*eols2)/(ei`*ei);

*rhohat=0;

pp[1,1]=sqrt(1-(rhohati**2));

do i=2 to periods;

pp[i,i]=1;

pp[i,i-1]=-rhohati;

end;

p=i(n)@pp;

bstep2=inv(zi`*p`*p*zi)*zi`*p`*p*yi;

estep2=p*yi-(p*zi*bstep2);

ee=j(periods,n,0);

do i=1 to n;

   ee[,i]=estep2[1+(i-1)*periods:i*periods];

end;

phii=(ee`*ee)/periods;

*sigmai=phi/(1-(rhohati**2));

****THIS SECTION IS USED TO CALCULATE MEAN OMEGA;

****AND Y AND X VALUES OVER ALL SUBSAMPLES; 

phitotal=phitotal+phii;

rhohattotal=rhohattotal+rhohati;

end;

phibar=phitotal/(25-periods+1);

rhohatbar=rhohattotal/(25-periods+1);

omegabar=j(periods,periods,0);

do i=1 to periods;

   do j=1 to periods;

      omegabar[i,j]=rhohatbar**(abs(i-j));

   end;

end;

sigmabar=phibar/(1-(rhohatbar**2));

vbar=sigmabar@omegabar;

Q=root(vbar);

print Q;

finish;

run program;

quit;

run;

Your help will be highly appreciated.

Mantobaye

5 REPLIES 5
Ksharp
Super User

Can you

Center and scale the data  ?

manto
Fluorite | Level 6

Hi Xia Keshan,

Thanks for getting back to me with a suggestion. I have not done so yet. Actually, there are many N - T combinations, and the works for most of them expect for few cases.

Thanks again.

Ksharp
Super User

Hi,

You'd better post it as IML forum since your code is IML code , @Rick is there .

Rick_SAS
SAS Super FREQ

We don't have your data, so seeing your program is not helpful. Could you show us the final

sigmabar and omegabar matrices ?

proc iml;

sigmabar = {...};

omegabar = {...};


manto
Fluorite | Level 6

Hi Rick,

I am willing to provide the data set as well, but if failed when I tried to upload it with the initial post. Could you possibly contact me at moundigbaye.mantobaye@pg.canterbury.ac.nz so I can send it to you? The data is not confidential so I can share it.

Mantobaye

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!

How to Concatenate Values

Learn how use the CAT functions in SAS to join values from multiple variables into a single value.

Find more tutorials on the SAS Users YouTube channel.

Click image to register for webinarClick image to register for webinar

Classroom Training Available!

Select SAS Training centers are offering in-person courses. View upcoming courses for:

View all other training opportunities.

Discussion stats
  • 5 replies
  • 468 views
  • 0 likes
  • 3 in conversation