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
Are you ready for the spotlight? We're accepting content ideas for SAS Innovate 2025 to be held May 6-9 in Orlando, FL. The call is open until September 25. Read more here about why you should contribute and what is in it for you!
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.