/**********************************************************
                              Problem  7.4
                 Finding a bifurcation point
                           t.h.zang, april 1999
***********************************************************/

new;
library pgraph;

/* starting value for theta */
theta=-1.91;

/* ending value for theta in the search */
endtheta=0.92;

/*step length */
step=0.001;

/* specification of parameters */
alpha=0.67;
beta=0.97;
delta=0.1;

/* solve for steady state values */
ystar=(1/alpha)*((1/beta)-1+delta);
cstar=ystar-delta;
lstar=(1-alpha)*ystar/(cstar+(1-alpha)*ystar);

/* initializing the output */
out1=1;
out2=0;
out3=0;
histtheta=theta;
out=out1~out2~out3;

DO WHILE theta<=endtheta;

/* define the matrices of the linearized system */
a1=zeros(2,2);
a1[1,1]=alpha/(1-theta);
a1[2,2]=1;

a2=zeros(2,2);
a2[1,1]=-1;
a2[1,2]=(1-alpha)/(1-theta);
a2[2,1]=-1;
a2[2,2]=1/(1-lstar);
                                                         
a3=zeros(2,2);
a3[1,1]=1-delta;
a3[1,2]=-cstar;
a3[2,2]=1;

a4=zeros(2,2);
a4[1,1]=ystar;

a5=zeros(2,2);
a5[1,1]=-1;
a5[2,1]=-alpha*beta*ystar;
a5[2,2]=-1;

a6=zeros(2,2);
a6[2,1]=alpha*beta*ystar;


/* matrix in the reduced form system
  of stochastic difference equations for the state variables */
j=-inv(a3-a4*inv(a2)*a1)*(a5-a6*inv(a2)*a1);

/* diagonalize the system */
{ eval,evec }=eigv(j);       @find the eigenvalues and eigenvectors@

length=abs(eval);                   @taking the modulus@
sorted=sortc(length,1);

out=out|(1~(sorted'));
histtheta=histtheta|theta;        @the list of values for the plot@

theta=theta+step;

ENDO;


xy(histtheta,out);                @graphical output@