@ ------------------------------- Test_2.g ------------------------------------
"
    28 December 2010
    Alfred Maußner

    Solve the benchmark business cycle model of Heer and Maußner (2009),
    Section 1.5.1 with the help of the SDRAM program.
    
    Test_2.g uses the model's equations.
"
------------------------------------------------------------------------------- @

@ -------------------- Do not change the next two lines of code ---------------- @
new;                  @ clear memory @
library user, pgraph; @ load standard libraries @
@ ------------------------------------------------------------------------------ @

/* Choice of solution method */
_equations=1;   @ "if =1 you must input the equations of your model, if =0 you must input the coefficients of the matrices of your log-linearized model" @
_GSchur=0;      @ "if =1 the generalized Schur factorization will be used to compute the linear part of the solution" @
_linear=2;      @ "if =1 only linear approximations of the models equations are computed"    @

/* Choice of output */
_Moments=1;    @ "if _moments=0, no moments will be computed"               @
_IR=1;         @ "if _IR=0, no impulse responses will be computed"          @
_confidence=1; @ "if =1 confidence intervalls for the simulated moments will be printed " @

_HPL=1600;     @ "the filter weight of the Hodrick-Prescott filter"         @
               @ "if _HPL=0, the filter will not be applied"                @

nobs=80;       @ "the number of elements in each simulated time series"     @
nobs1=20;      @ "the length of the impulse reponsees"                      @
nofs=300;      @ "the number of simulations from which moment averages are to be computed" @

_scale=0;      @ "if _scale=1 all plots will have the same scale"           @
_nocolor=0;    @ "if=0, lines will be colored"                              @

/* Choice of random numbers */
_makeShocks=0; @ "if=1, an array of shocks will be created that can be used in repeated simulations"   @
_loadShocks=0; @ "if=1, the program looks for the file epsilon_array and loads shocks from this array" @
               @ "if=0, the program draws new random numbers for each simulation"                      @

/* Choice of output file */
outfile="Benchmark_1.txt"; @ "the name of the file to which the program writes its output" @

/* Determine the type and number of the model's variables */
_nx=1;  @ "the number of states, here:  x:=[k]                          "@
_nz=1;  @ "the number of shocks, here:  z:=[Z]                          "@
_ny=7;  @ "the number of non predetermined variables, here: y:=[y, c, i, r, w, n, lambda]  "@
_nu=6;  @ "if _GSchur=0 you must tell the program the number of static equations _nu. If _equations=0 the program assumes that the last _ny-_nu variables in y are costates" @
_ns=0;  @ "if _equations=0 you may include futher variables s in the vector y and provide the respective matrices Lsx and Lsz. See the manual" @

@ ---------------------------------- Do not change the next line of code --------------------- @
#include CoRRAM_1.src;
@ -------------------------------------------------------------------------------------------- @

/*
**
** You must put in the names of all variables in your model
** for which you either want the program to compute and display
** impulse respones or compute and print second moments.
** 
** Thus, variables that are in x, y, or z, for which you do
** not provide informaion, are ignored when plotting impulse
** responses or when printing second moments.
**
**
** The order in which the variables appear below is also the
** order in which results will be printed.
**
*/

Var[1].varname="Output";
Var[1].vartype="y";
Var[1].varpos=1;
Var[1].varprint=1;
Var[1].crosscorr=1;
Var[1].varplot=1;
Var[1].plotno=1;
Var[1].relsx=1;

Var[2].varname="Consumption";
Var[2].vartype="y";
Var[2].varpos=2;
Var[2].varprint=1;
Var[2].crosscorr=0;
Var[2].varplot=1;
Var[2].plotno=2;

Var[3].varname="Investment";
Var[3].vartype="y";
Var[3].varpos=3;
Var[3].varprint=1;
Var[3].crosscorr=0;
Var[3].varplot=1;
Var[3].plotno=2;

Var[4].varname="Hours";
Var[4].vartype="y"; 
Var[4].varpos=6;  
Var[4].varprint=1;  
Var[4].crosscorr=0; 
Var[4].varplot=1;   
Var[4].plotno=4;    

Var[5].varname="Real Wage";
Var[5].vartype="y"; 
Var[5].varpos=5;   
Var[5].varprint=1;  
Var[5].crosscorr=1; 
Var[5].varplot=1;   
Var[5].plotno=3;    

Var[6].varname="User Costs of Capital";
Var[6].vartype="y";
Var[6].varpos=4;
Var[6].varprint=1;
Var[6].crosscorr=0;
Var[6].varplot=1;
Var[6].plotno=3;
    
Var[7].varname="Capital Stock";
Var[7].vartype="x";
Var[7].varpos=1;
Var[7].varprint=0;
Var[7].crosscorr=0;
Var[7].varplot=1;
Var[7].plotno=4;

Var[8].varname="TFP Shock";
Var[8].vartype="z";
Var[8].varpos=1;
Var[8].varprint=0;
Var[8].crosscorr=0;
Var[8].varplot=1;
Var[8].plotno=1;

/* Set legend positions individually */
_plegend[1,1,1]=10;
_plegend[1,1,2]=0.9;

_plegend[1,2,1]=10;
_plegend[1,2,2]=3.5;

_plegend[1,4,1]=6;
_plegend[1,4,2]=0.5;    
  
/*
**
** First, you must insert the parameters of your model.
**
** Second, if _equation=1, you must compute the stationary solution
** of all variables in your model otherwise you must compute
** - if necessary - those quantities that appear the in matrices.
**
** Third, you must put the equations of your model in the
** proc SYS
**
** Fourth, you must stet up the matrices _Rho and _Omega and
** asign a value to _sigma (the parameters of the transition
** equation for the shocks).
*/

@ First Part @
a=1.005;
alpha=0.27;
beta=0.994;
eta=2.0;
delta=0.011;
rhoZ=0.90;
sigmaz=0.0072;
nstar=0.13;

@ Second Part @
/* Compute the solution step by step */
yk=(a^eta-beta*(1-delta))/(beta*alpha);
ck=yk+(1-a-delta);
theta = (1-alpha)*(yk/ck)*(1-nstar)*(1/nstar);
kn=yk^(1/(alpha-1));
kstar=kn*nstar;
ystar=yk*kstar;
cstar=ck*kstar;
istar=ystar-cstar;
wstar=(1-alpha)*(ystar/nstar);
rstar=alpha*yk;
lstar=(cstar^(-eta))*((1-nstar)^(theta*(1-eta)));

/* Solve the respective non-linear system. */
/* Compute an initial vector x0=[k,y,c,i,t,w,r,l] by drawing from a uniform distribution on [0,1]
** Note, I do not compute nstar but fix nstar and compute theta (t in the code below) so that it implies nstar */
x0=rndu(8,1);
{ x1, retcode } = eqSolve(&Sys0,x0); 

proc(1)=Sys0(x); @ the sytem of equations at the stationary solution. I had to cancel lambda in the last equation for EqSolve to work @

    local k, y, c, i, t, w, r, l, fx;

    fx=zeros(8,1);
    k=x[1];
    y=x[2];
    c=x[3];
    i=x[4];
    t=x[5];
    w=x[6];
    r=x[7];
    l=x[8];

    fx[1]=(c^(-eta))*((1-nstar)^(t*(1-eta)))-l;
    fx[2]=t*(c^(1-eta))*((1-nstar)^(t*(1-eta)-1))-w*l;
    fx[3]=w-(1-alpha)*(nstar^(-alpha))*(k^alpha);
    fx[4]=r-alpha*(nstar^(1-alpha))*(k^(alpha-1));
    fx[5]=y-(nstar^(1-alpha))*(k^alpha);
    fx[6]=y-c-i;
    fx[7]=a*k-(1-delta)*k - i;
    fx[8]=l*(beta*(a^(-eta))*(1-delta+r))-l;

retp(fx);

endp;

@ Third Part @

/* You must gather the stationary solutions of the variables
** in the following two vectors: (the _zstar vector is always the zero vector and the program adds this vector for you)
*/
@ "the vector _xstar collects the stationary values of the state variables" @
_xstar=kstar; 
@ "the vector _ystar collects the stationary values of the variables in the vector y" @
_ystar=ystar|cstar|istar|rstar|wstar|nstar|lstar;

/* Since the procedure Sys receives _wstar|_wstar, with
** _wstar=_xstar|_ystar|_zstar, where
** the second entry of _wstar indicates next period values,
** you must use this ordering of the variables in Sys. Also note that
** the line before the endp statement of Sys must contain the
** command:
** if _eqno==0; retp(fx); else; retp(fx[_eqno]); endif;
*/


@ Third Part: in Sys: specify the equations of your model @

proc(1)=Sys(w);

  local fx, c1, y1, i1, r1,w1,n1, k1, l1, z1, c2, y2, i2, r2,w2,n2, k2, l2, z2;
  k1=w[1];
  y1=w[2];
  c1=w[3];
  i1=w[4];
  r1=w[5];
  w1=w[6];
  n1=w[7];
  l1=w[8];
  z1=exp(w[9]);

  k2=w[10];
  y2=w[11];
  c2=w[12];
  i2=w[13];
  r2=w[14];
  w2=w[15];
  n2=w[16];
  l2=w[17];
  z2=exp(w[18]);

  fx=zeros(_nx+_ny,1);

  fx[1]=(c1^(-eta))*((1-n1)^(theta*(1-eta)))-l1;
  fx[2]=theta*(c1^(1-eta))*((1-n1)^(theta*(1-eta)-1))-w1*l1;
  fx[3]=w1-(1-alpha)*z1*(n1^(-alpha))*(k1^alpha);
  fx[4]=r1-alpha*z1*(n1^(1-alpha))*(k1^(alpha-1));
  fx[5]=y1-z1*(n1^(1-alpha))*(k1^alpha);
  fx[6]=y1-c1-i1;

  fx[7]=a*k2-(1-delta)*k1 - i1;
  fx[8]=beta*(a^(-eta))*l2*(1-delta+r2)-l1;

  if _eqno==0; retp(fx); else; retp(fx[_eqno]); endif;
  retp(fx);

endp;


/* Transition equation for shock */
_Rho[1,1]=RhoZ;
_Sigma=SigmaZ;
_Omega=eye(1);


@ -------------- Do not change the next line of code ------------- @
#include CoRRAM_2.src;
@ ----------------------------------------------------------------------- @

/* if you model includes additional variables, _ns>0, you must input Lsx and Lsz here: */
Lsx=0;
Lsz=0;

@ -------------- Do not change the remaining lines of code ------------ @
#include CoRRAM_3.src;
ende:
end;
@ ----------------------------------------------------------------------- @


