%Final Model with three shocks and X=0.1 and pnynyt=2 
% Basic Stochastic Growth Model
% Lecture: Solving Nonlinear Dynamic Stochastic Models


%----------------------------------------------------------------
% 0. Housekeeping (close all graphic windows)
%----------------------------------------------------------------

close all;

%----------------------------------------------------------------
% 1. Defining variables
%----------------------------------------------------------------

% list of endogenous variables
var  ty y yt yn k kt kn lt l_n c ct cn i x r pn pc z1 z2 z3 cty ity kty yty yny ct_y iy ky ktkn ltln xy cy ynyt py pnyny;

% the following exogenous variables are shocked
varexo e1 e2 e3 ;

parameters beta alpha delta omega eta mu nu  sigma  psi11 psi12 psi13 psi21 psi22 psi23 psi3  sigma1 sigma2 sigma3 z3_bar  A An X ;

%----------------------------------------------------------------
% 2. Calibration
%----------------------------------------------------------------

alpha   = 0.30;
beta    = 0.95;
delta   = 0.05;
omega   = 0.26;
eta     = 0.5;
nu      = -10;
mu      = 0.316;
psi11   = 0.95;
psi12   = 0;
psi13   = 0;
psi21   = 0.95;
psi22   = 0;
psi23   = 0;
psi3   = 0.95;
%psi32   =0;
%psi33   =0;

z3_bar=1;


sigma   = 5;  
sigma1  = 0.04;
sigma2  = 0.04;
sigma3  =0.05;
X=0.1;
A=1;
An=1.52;
pnynyt=1.3;


kyt=alpha/(1/beta-1+delta);  

kyn=eta/(1/beta-1+delta); 

                                                                             

%----------------------------------------------------------------
% 3. Model
%----------------------------------------------------------------

model; 
  
  c^(-sigma)/pc = beta*((c(+1)^(-sigma))/pc(+1))*r(+1);       % Euler eqution
  r = alpha*A*z1*(kt(-1)/lt)^(alpha-1) + 1-delta;                               % return
  pn = ((1-omega)/omega)*(cn/ct)^(-(1+mu));                              % Substitutions between the consumption of tradable and non-tradable goods
  pc = (omega^(1/(1+mu))+((1-omega)^(1/(1+mu)))*pn^(mu/(1+mu)))^((1+mu)/mu);
  A*z1*alpha*(kt/lt)^(alpha-1)= pn*An*z2*eta*(kn/l_n)^(eta-1) ;               % The cost of capital
  A*z1*(1-alpha)*(kt/lt)^alpha = pn*An*z2*(1-eta)*(kn/l_n)^eta;               % The cost of labor
                                          
   
  ct+i =yt+x ;                                                                  % aggregate resource constraint
  cn = yn;
  c  = (omega*(ct)^(-mu)+ (1-omega)*(cn^(-mu)))^(-1/mu);                        %overall capital
  yt = z1*A*kt(-1)^alpha*(lt)^(1-alpha);                                        % tradable production function
  yn = z2*An*kn(-1)^eta*(l_n)^(1-eta);                                          % non-tradable sector production funstion
  y= yt+pn*yn;                                                                  % GDP
  ty= y+x;
  lt+l_n=1;
  k = i + (1-delta)*k;                                                          % overall labor
  
  k = (kt^(-nu)+kn^(-nu))^(-1/nu);                                              % capital accumulation
  x = X*z3;
  
  cty=c/ty;
  ity=i/ty;  
 
  
 kty=k/ty; 
yty=yt/y;
yny=pn*yn/y;
ct_y=ct/y;
ktkn=kt/kn;

ltln=lt/l_n;
iy=i/y;
ky=k/y;
xy=x/y;
cy=c/y;
ynyt=yn/yt;  
py=pn*yn;
pnyny=pn*yn/y;                                                   
  
  log(z1) =psi11*log(z1(-1))+psi12*log(z2(-1))+psi13*log(x(-1)) +e1;         % technology
  log(z2) =psi21*log(z2(-1))+psi22*log(z1(-1))+psi23*log(x(-1)) +e2;

  log(z3) =(1-psi3)*log(z3_bar) + psi3*log(z3(-1))+e3;  
 

end;

%----------------------------------------------------------------
% 4. Computation
%----------------------------------------------------------------

initval; 



lt=1/(1+((1-eta)/(1-alpha))*pnynyt);
l_n=1-lt;
yt=(A^(1/(1-alpha)))*((kyt)^(alpha/(1-alpha)))*lt;
kt=kyt*yt;
kn=kyn*pnynyt*yt;
k = (kt^(-nu)+(kn)^(-nu))^(-1/nu);
i=delta*k;
yn = An*kn^eta*(l_n)^(1-eta);

y=yt+pnynyt*yt;
x=X;
ct=yt+x-i;
cn=yn;
c  = (omega*(ct)^(-mu)+ (1-omega)*(cn^(-mu)))^(-1/mu);
pn=((1-omega)/omega)*(cn/ct)^(-(1+mu));                              
pc = (omega^(1/(1+mu))+((1-omega)^(1/(1+mu)))*pn^(mu/(1+mu)))^((1+mu)/mu);
ty=y+x;
  cty=c/ty;
  ity=i/ty;  
 
  
 kty=k/ty; 
yty=yt/y;
yny=pn*yn/y;
ct_y=ct/y;
ktkn=kt/kn;

ltln=lt/l_n;
iy=i/y;
ky=k/y;
xy=x/y;
cy=c/y;
ynyt=yn/yt; 
py=pn*yn;  
pnyny=pn*yn/y;   
e1 = 0;
e2 = 0;
e3 = 0;
z1 = 1;
z2 = 1;
z3=1;
end;
steady;
check;


shocks;
var e1=sigma1^2;
var e2=sigma2^2;
var e3=sigma3^2;
var e1,e2=0;
var e1,e3=0.4*sigma1*sigma3;
var e3,e2=0.4*sigma3*sigma2;

end;

resid(1);

steady;

stoch_simul(hp_filter = 6.5, order = 1);