var pi r n v x theta a y w lambda u G jdr jcr lambdan RM phi;
varexo e eA;

parameters beta lambdad rho nss mss uss vss 
xi q yss chi eta sigma Hss lambdanss Hauc 
at mu nu sig lambdass b wss xss thetass Ass 
rhoA Fss m D fss Rss psi e_Ha;
                                                                                                                                                       
beta=0.99;                                 % Discounting factor (equivalent to 1% nominal interest rate in ss)  
mu=0;                                      % Mean distribution parameter
sig=0.12;                                  % Variance distribution parameter                
lambdad=0.068;                             % Job destruction Rate
rho=0.5;                                   % AR(1) process schock parameter
nss=0.88;                                  % Employment rate steady state
lambdass=0.1;                              % Separations steady state
mss=(lambdass/(1-lambdass))*nss;           % Matches steady state
lambdanss=0.034;                           % Endogenous separations steady state
uss=1-nss;                                 % Unemployment steady state
q=0.7;                                     % job filling rate
vss=mss/q;                                 % Vacancies steady state
xi=0.4;                                    % search elasticity of matches
chi=0.07;                                  % Vacancy posting costs
nu=11;                                     % Elasticity of substitution
eta=0.5;                                   % Worker's bargaining power
b=0.6;                                     % Unemployment benefit
xss=(nu-1)/nu;                             % steady state markup over marginal costs
sigma=2;                                   % Risk aversion parameter
at=logninv(lambdanss,mu,sig);              % Productivity threshold steady state
Ass=1;                                     % Overall productivity steady state
rhoA=0.95;                                 % AR(1) parameter overall productivity
Hauc=quad(@afa,at,8);                      % Integral of productivity threshold
Hss=Hauc/(1-logncdf(at,mu,sig));           % Conditional expectation of productivity
Fss=logncdf(at,mu,sig);                    % Threshold steady State in cdf
fss=lognpdf(at,mu,sig);                    % Threshold steady State in pdf
yss=nss*Hss*Ass-chi*vss;                   % production steady state
thetass=vss/uss;                           % Steady state labor market tightness
wss=eta*thetass*chi+(1-eta)*b+eta*((nu-1)/nu)*Hss; % steady state wage
m=q*thetass^(xi);                          % search efficiency
D=((1-eta)*((nu-1)/nu)*Ass*Hss)-chi*eta*thetass-(1-eta)*b+((chi/m)*(thetass^(xi)));
Rss=1/beta;
psi=40;
e_Ha=((at*fss*(Hss-at))/(Hss*(1-Fss)));


% Production function
% Overall separation dynamics
% Labor market tightness
% Employment dynamics
% Job creation condition
% Productivity threshold
% Identity
% Euler consumption equation
% Phillips Curve
% Taylor Rule
% Unemployment dynamics
% Shock

model(linear);
y=((Ass*nss*Hss)/yss)*(G+n+((e_Ha)/Hss)*a)-((chi*vss)/yss)*v;
lambda=(((1-lambdad)*lambdanss)/(lambdass))*lambdan;
lambdan=at*(fss/Fss)*a;
theta=v-u;
n(+1)+(lambdass/(1-lambdass))*lambda(+1)=(1-lambdass)*n+lambdass*(1-xi)*v+lambdass*xi*u;
xi*theta-sigma*y=(chi*(xi/m)*thetass^(xi)-eta*chi*thetass)*(1/D)*theta(+1)
+(((1-eta)*Ass)*((nu-1)/nu)*(e_Ha)*(1/D))*a(+1)+((((1-eta)*Ass*Hss)*((nu-1)/nu))*(1/D))*G(+1)
+((Ass*Hss*(1-eta))*(1/(nu*D)))*x(+1)-(lambdass/(1-lambdass))*lambda(+1)-sigma*y(+1);
a=-(1/(nu-1))*x-G+(((chi/(1-eta))*(eta*thetass-(xi/m)*(thetass)^(xi)))/(b+(chi/(1-eta))*(eta*thetass-(1/m)*thetass^(xi))))*theta;
sigma*y(+1)=sigma*y+r-pi(+1);
RM=-(1/(Rss-1))*r+sigma*y;
RM=RM(-1)+phi-pi;                                                                                                                                                                                           
psi*pi=x+beta*psi*pi(+1);
wss*w=eta*chi*thetass*theta+eta*((nu-1)/nu)*Ass*Hss*G
+((eta*Ass*Hss)/nu)*x+eta*((nu-1)/nu)*Ass*(e_Ha)*a;
uss*u=-nss*n;
jdr=(lambdass/(lambdass-lambdad))*lambda; 
((lambdass-lambdad)/(lambdass))*jcr+(lambdass/(1-lambdass))*lambda=n(-1)+xi*u(-1)+(1-xi)*v(-1); 
G=rhoA*G(-1)+eA;
phi=rho*phi(-1)+e; 
end; 
                                               
shocks;
var e=0.0064;
end;

stoch_simul(irf=5);