%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%  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
%%%%%%% yflex=output with fully flexible prices
%%%%%%% yflex1=output with fully flexible prices(steady state value)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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


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

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=3;
chi= 102.89;
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=5;
mu=(theta/(theta-1));
rho_z=0.9;
rho_v=0.5;

model(linear);

yflex1=(chi*n^(-eta)*varphi^(1/(1-alpha))*(1-alpha)*E_z*(alpha*E_z)^(alpha/(1-alpha)))^(1/sigma)+k(+1)-(1-delta)*k;

yflex=-(1/sigma)*eta*n+(1/sigma)*(1/(1-alpha))*varphi+(1/sigma)*(1/(1-alpha))*E_z+(STEADY_STATE(k)/yflex1)*k(+1)-((1-delta)/yflex1)*STEADY_STATE(k)*k;

x=y-yflex;

-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;



STEADY_STATE(c)*c+STEADY_STATE(k)*k(+1)-STEADY_STATE(k)*(1-delta)*STEADY_STATE(rk)*k = 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+z+((1-alpha)*n)-((1-alpha)*k);

yflex1*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;

varphi = ((theta-1)/theta);

    n  = (((1-beta*(1-delta)) /beta)  * (chi/(lambda *varphi*(1-alpha)))^((1-alpha)/alpha)  * (1/(varphi*alpha)) )^ (-alpha/(eta*(1-alpha)));

    k  = ((n^(eta+alpha)*chi)/ (lambda*varphi*(1-alpha))) ^(1/alpha);

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

lambda = (n^(eta)*chi)/w;

     c = lambda^(-1/sigma);

    rk = (1/beta)-(1-delta);
 
     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);