/**********************************************************
                      Problem  3.4
              Simulation of an RBC economy
       thomas hintermaier, bauersknecht, march 1999
***********************************************************/

new;

/* !!!! specify question a., b., or c. of problem 3.4  !!!!*/
question=3; @1=a., 2=b., any other number=c.@

/* number of periods for the simulation */
periods=40;

/* specification of parameters */
rho=0.0163;
gam=0;
alpha=0.4;
lam=0.98;
varu=0.07;
theta=0.9;
ub=0.9957;
vare=0.01;
delta=0.0272;

/* solve for steady state values */
gss=ub^(1/(1-theta));
m1=((rho+delta)/alpha)^((alpha+gam)/(1-alpha)); @defined to reduce the mess 
                                                 in the computation of ss@
m2=(rho+delta-alpha*delta)/alpha;

eqsolveset;                            @ solve for steady state capital @
fn f(k)=(1-alpha)*k^(-gam)-(m1*(m2*k-gss));
__output=0;
{ kss,s }=eqsolve(&f,gss);

css=m2*kss-gss;
lss=((rho+delta)/alpha)^(1/(1-alpha))*kss;
yss=((rho+delta)/alpha)*kss;

/* define the matrices of the linearized system */
a1=zeros(2,4);
a1[1,2]=alpha;
a1[1,3]=1;
a1[2,1]=-1;
a2=zeros(2,2);
a2[1,1]=-1;
a2[1,2]=1-alpha;
a2[2,1]=1;
a2[2,2]=-(1+gam);
 
b1=1-delta;       @some coefficients needed in the following matrices@
b2=yss/kss;
b3=-css/kss;
b4=-gss/kss;
r=alpha/(1+rho)*yss/kss;

a3=zeros(4,4);
a3[1,1]=b3;
a3[1,2]=b1;
a3[1,4]=b4;
a3[2,1]=1;
a3[3,3]=lam;
a3[4,4]=theta;
a4=zeros(4,2);
a4[1,1]=b2;
a5=zeros(4,4);
a5[1,2]=-1;
a5[2,1]=-1;
a5[2,2]=-r;
a5[3,3]=-1;
a5[4,4]=-1;
a6=zeros(4,2);
a6[2,1]=r;
a7=zeros(4,3);
a7[2,3]=1;
a7[3,1]=1;
a7[4,2]=1;

/* matrices in  the system of stochastic difference equations for the state variables */
j1=-inv(a3-a4*inv(a2)*a1)*(a5-a6*inv(a2)*a1);
j2=-inv(a3-a4*inv(a2)*a1)*a7;

/* diagonalize the system */
{ evalu,evecu }=eigv(j1);          @find the eigenvalues and eigenvectors@
 
modev=abs(evalu);                                @and order them@
temp=sortc(modev~evalu~evecu',1);
eval=temp[.,2];
evec=(temp[.,3:cols(temp)])';

/* decide on existence and uniqueness of equilibrium */
count=0;
ind=1;
DO WHILE ind<=rows(eval);
      IF abs(eval[ind])<1;
         count=count+1;
      ENDIF;
      ind=ind+1;
ENDO;

/* terminate program in case of no solution or multiple equilibria */
IF count<1;
    print "No eigenvalue inside unit circle: multiple equilibria ";
    waitc;
    end;
ELSEIF count>1;
    print "More than one eigenvalue inside unit circle: no solution ";
    waitc;
    end;
ENDIF;

/* generate the sequences of shocks as specified in questions a., b., and c., respectively */
shocks=zeros(2,periods);
IF  question==1;
   shocks[1,1]=sqrt(varu);
ELSEIF question==2;
   shocks[2,1]=sqrt(vare);
ELSE;
   seed=151271;           @nota bene: this is NOT MY birthdate@
   rndseed seed;
   omega=varu~0|0~vare;
   shocks=chol(omega)*rndn(cols(shocks),rows(shocks))';
ENDIF;

/* initializing a matrix to store the sequences of relative deviations from  ss
    for c, k, a, g, y, l (this ordering) */
devstore=zeros(6,periods+1);

/**** simulating the economy for the specified number of periods ****/
z=zeros(4,periods+1);   @initializes matrix of obsvervations for the
                                                    diagonalized system@
@!!!! z[1,.]=0, that is, the first row of z is zero for all periods, 
    since it corresponds to the eigenvalue inside the unit circle !!!!@

combi=real(inv(evec)*j2);   @by this matrix a linear combination of the
                              fundamental and the expectational errors
                              is added to each element of the 
                              diagonalized system@
                                                               
rela=combi[1,.];   @since, for all t, z[1,.] is zero, this row of the matrix
                     tells us how the expectational error is a function 
                     of the fundamental errors, ...@
w=zeros(1,periods);
t=1;
DO WHILE t<=periods;      @...which is used here to ...@
    w[t]=-1/rela[3]*rela[1]*shocks[1,t]-1/rela[3]*rela[2]*shocks[2,t];
                         @ ...generate the sequence of expectational errors@
     t=t+1;
ENDO;
errors=shocks|w; @concatenation of all error sequences into one matrix@

/* generating the sequences for the observations in the diagonalized system */
k=2;
DO WHILE k<=4;
    t=1;
    DO WHILE t<=periods;
        z[k,t+1]=1/eval[k]*z[k,t]-1/eval[k]*combi[k,.]*errors[.,t];
        t=t+1;
    ENDO;
   k=k+1;
ENDO;

/* recovering the proportional deviations of the state variables from the diagonalized system */
j=1;
DO WHILE j<=periods+1;
   devstore[1:4,j]=evec*z[.,j];
    j=j+1;
ENDO;

/* get the proportional deviations of the other variables from their 
   contemporaneous relationship with the state variables */
j=1;
DO WHILE j<=periods+1;
   devstore[5:6,j]=-inv(a2)*a1*devstore[1:4,j];
    j=j+1;
ENDO;

/* convert relative deviations from ss to levels of variables 
   by multiplying by the steady state values and 
   adding the steady state values */
cstore=devstore[1,.]*css+css;
kstore=devstore[2,.]*kss+kss;
astore=devstore[3,.]+1;
gstore=devstore[4,.]*gss+gss;
ystore=devstore[5,.]*yss+yss;
lstore=devstore[6,.]*lss+lss;

/* the last section generates the graphical output */
library pgraph;
graphset;
_ptitlht=0.5;
xax=seqa(1,1,cols(devstore));      @an additive sequence defining the x-axis@
begwind;
window(2,3,0);
setwind(1);
  _pcolor=1;
  title("Consumption");
  xy(xax,cstore');
nextwind;
  _pcolor=2;
  title("Capital");
  xy(xax,kstore');
nextwind;
  _pcolor=3;
  title("Prod. disturb.");
  xy(xax,astore');
nextwind;
  _pcolor=4;
  title("Government");
   xy(xax,gstore');
nextwind;
  _pcolor=5;
  title("Output");
  xy(xax,ystore');
nextwind;
  _pcolor=6;
  title("Labour");
  xy(xax,lstore');
endwind;

end;