%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%  r=interest rate
%%%%%%  k=capital stock
%%%%%%% pi=inflation
%%%%%%% lambda=shadow price
%%%%%%% d=real profits
%%%%%%% c=consunption
%%%%%%% w=real wage
%%%%%%% n=employment
%%%%%%% rk=capital costs
%%%%%%% y=output
%%%%%%% E_z=technology shock
%%%%%%% E_v=monetary policy shock
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


var  r k pi lambda varphi d c w n rk y q x  E_z  E_v;      %endogenous Variables
										

varexo z v ;   							
                                                                      
parameters beta sigma chi eta mu theta omega rho_r phi_pi phi_x rho_z rho_v alpha delta tau_K;


beta=0.99;
sigma=1;
eta=1;
theta=6;
omega=0.75;
rho_r=0.5;
phi_pi=1.5;
phi_x=0.125;
alpha=0.3;
delta=0.025;
tau_K=20;
mu=(theta/(theta-1));
rho_z=0.9;
rho_v=0.5;


chi=(lambda*w)/((1/3)^(eta));


model(linear);

x=y-STEADY_STATE(y);

-sigma*c=lambda;

lambda+w=eta*n;

lambda=lambda(+1)+r-pi(+1);

lambda-(beta*(1-delta+STEADY_STATE(rk))*lambda(+1))=tau_K*k+tau_K*k(+1)+beta*tau_K*(  ((STEADY_STATE(n)*STEADY_STATE(w)*w(+1))/STEADY_STATE(k))  + 

                                            ((STEADY_STATE(n)*STEADY_STATE(w)*n(+1))/STEADY_STATE(k))+ (STEADY_STATE(rk)*rk(+1))  +  

                                            ((STEADY_STATE(d)*d(+1))/STEADY_STATE(k)) +  ((1-delta+STEADY_STATE(rk))*k(+1))        

                                           - ((STEADY_STATE(c)*c(+1))/STEADY_STATE(k)) ) + beta*STEADY_STATE(rk)*rk(+1);



STEADY_STATE(c)*c+STEADY_STATE(k)*k(+1)-STEADY_STATE(k)*(1-delta+STEADY_STATE(rk)) = STEADY_STATE(n)*STEADY_STATE(w)*w+STEADY_STATE(n)*STEADY_STATE(w)*n+
                    (STEADY_STATE(k)*STEADY_STATE(rk)*rk)+STEADY_STATE(d)*d;

y=E_z+alpha*k+(1-alpha)*n;

w=varphi + E_z + (alpha*k)-(alpha*n);

rk=varphi+E_z+((1-alpha)*n)-((1-alpha)*k);

STEADY_STATE(y)*y=STEADY_STATE(c)*c +STEADY_STATE(k)*k(+1)-(1-delta)*STEADY_STATE(k)*k;


pi*(omega/(1-omega))=(1-(omega*beta))*varphi+(omega*beta)*((omega/(1-omega))*pi(+1)+pi(+1));

q=(omega/(1-omega))*pi;

r=rho_r*r(-1)+(1-rho_r)*(phi_pi*pi+phi_x*x)+E_v; 

E_z=rho_z*E_z(-1)+ z;

E_v=rho_v*E_v(-1)+ v;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Steady State values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



 initval;
     r = 1/beta;

   rk = (1/beta)-1+delta;
 varphi = ((theta-1)/theta);

      n=1/3;

   k= n/((rk/(varphi*alpha))^(1/(alpha-1)));

    w  = varphi*(1-alpha)*( (k/n)^(alpha) );

   y=(k^(alpha))*(n^(1-alpha));	

   c=y-delta*k;

 lambda=c^(-sigma);

d=c+(delta*k)-(w*n)-(rk*k);
 
      pi=1;

       q=1;
   pi=1;
  x=0;
   E_z=1;
  E_v=1;

  end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%steady; 


check;

shocks;
       var z=sigma^2;
  
       var v=sigma^2;
end;

stoch_simul(periods=2100) c w n r pi d x E_v y k rk E_z;