## Rounding error suspected. Could you help with the solution?

Occasional Contributor
Posts: 6

# Rounding error suspected. Could you help with the solution?

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;

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,1periods-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

Super User
Posts: 10,770

## Re: Rounding error suspected. Could you help with the solution?

Can you

Center and scale the data  ?

Occasional Contributor
Posts: 6

## Re: Rounding error suspected. Could you help with the solution?

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.

Super User
Posts: 10,770

## Re: Rounding error suspected. Could you help with the solution?

Hi,

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

SAS Super FREQ
Posts: 4,240

## Re: Rounding error suspected. Could you help with the solution?

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 = {...};

Occasional Contributor
Posts: 6

## Re: Rounding error suspected. Could you help with the solution?

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

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