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
Can you
Center and scale the data ?
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.
Hi,
You'd better post it as IML forum since your code is IML code , @Rick is there .
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 = {...};
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
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.
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.