Rounding error suspected. Could you help with the solution?

[manto]

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

[Ksharp]

Can you

Center and scale the data  ?

[manto]

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]

Hi,

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

[Rick_SAS]

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]

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