/******************************************************************
                             Problem 8.4
                  Roots of a difference equation
                  Thomas Hintermaier, Sept. 1999
******************************************************************/
new;

format 8,3;

/* defining the points for the externality 
for which the system is evaluated */
points=seqa(0,0.01,100);

outlist=zeros(rows(points),7);        /* initialization of output */

z=1;
DO WHILE z<=rows(points);

/* specification of the parameters */
sigma=0.75;
beta=0.95;
a=0.33;
delta=0.1;
gam=0;
lam=0.007;
      
theta=points[z];


/* computing some steady state values and elasticitiy parameters */
yokss=(1/beta+delta-1)*1/a;
cokss=yokss-delta;
lss=(1-a)*yokss/cokss;
pa1=1-delta;
pa2=yokss;
pa3=-cokss;
pr=beta*a*yokss;
pchi=-(1-sigma)*(lss^(gam+1))/sigma;


/* definition of the matrices in the linearized dynamic system */
a1=zeros(3,3);
a1[1,2]=1/sigma;
a1[2,1]=-a*(1+theta);
a1[2,3]=-1;

a2=zeros(3,3);
a2[1,2]=1;
a2[1,3]=-pchi;
a2[2,1]=1;
a2[2,3]=-(1-a)*(1+theta);
a2[3,1]=-1;
a2[3,2]=1;
a2[3,3]=1+gam;

a3=zeros(3,3);
a3[1,1]=pa1;
a3[2,2]=-1;
a3[3,3]=lam;

a4=zeros(3,3);
a4[1,1]=pa2;
a4[1,2]=pa3;

a5=zeros(3,3);
a5[1,1]=-1;
a5[2,1]=-pr;
a5[2,2]=1;
a5[3,3]=-1;

a6=zeros(3,3);
a6[2,1]=pr;


/* defining the matrix j1 in the reduced dynamic system */
j1=-inv(a3-a4*inv(a2)*a1)*(a5-a6*inv(a2)*a1);

/* finding the eigenvalues of j1 and their modulus */
eval=eig(j1);
length=abs(eval);
@ print eval~length; @

IF z==1;
"theta           eval1                   eval2                 eval3     	 modulus1      modulus2       modulus3";
ENDIF;

outlist[z,1]=theta;
outlist[z,2]=eval[1]; outlist[z,3]=eval[2]; outlist[z,4]=eval[3];
outlist[z,5]=length[1]; outlist[z,6]=length[2]; outlist[z,7]=length[3];

z=z+1;
ENDO;

outlist;