%Taking ratio pnyn/yt as a parameter and adding X=(0;0.01) as a parameter.In between these numbers the model works
% 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  z1 z2 z3 i  kt kn k yt yn y ct cn c pn pc lt l_n r x ;

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

parameters beta alpha delta omega eta mu nu  sigma  psi1 psi2 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.3157;
psi1   = 0.95;
psi2   = 0.95;
psi3  =0.95;
z3_bar=1;


sigma   = 5;  
sigma1  = 0.11;
sigma2  = 0.11;
sigma3  =0.74;
A=1;
An=1.52;
X=0.1;
pnynyt=2;

klt=((alpha*A)/(1/beta-1+delta))^(1/(1-alpha));                                    %Kt/Lt
kln=((1-alpha)/(1-eta))*(eta/alpha)*(klt);                                          %Kn/Ln


                                                                             

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

model; 
  
  c^(-sigma)/pc = beta*((c(+1)^(-sigma))/pc(+1))*r(+1);       % Euler eqution
  r = alpha*A*z1*(kt/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^alpha*(lt)^(1-alpha);                                            % tradable production function
  yn = z2*An*kn^eta*(l_n)^(1-eta);                                                % non-tradable sector production funstion
  y= yt+pn*yn;                                                                      % GDP
  lt+l_n=1;
  k = i + (1-delta)*k(-1);                                                          % overall labor
  
  k = (kt^(-nu)+(kn)^(-nu))^(-1/nu);                                     % capital accumulation
  x=z3*X;
  
  
                                                   
  
  log(z1) =psi1*log(z1(-1)) +e1;                                                    % technology
  log(z2) =psi2*log(z2(-1)) +e2;
  
  log(z3) =(1-psi3)*log(z3_bar)+psi3*(log(z3(-1)))+e3;   

end;

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

initval; 
%pnynyt=1.75;

lt=1/(1+((1-eta)/(1-alpha))*pnynyt);
l_n=1-lt;
kt=klt*lt;
kn=kln*l_n;
k = (kt^(-nu)+(kn)^(-nu))^(-1/nu);
i=delta*k;
yt=A*lt*(klt)^alpha;
yn=An*l_n*(kln)^alpha;
x=X;
cn=yn;
ct=yt-i+x;
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);
y=yt+pn*yn;
  
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 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 = 1600, order = 1);