% this is the .mod file

var x a b c_e c_u c_n l u e s k v fs fv 
    Lambda xi chi MRS_NC Omega Gamma tau pi i Y D dw P p_star p_a p_b q varsigma
    z rho b_bar eps ;

varexo e_z e_r e_b e_i; 

parameters beta alpham alpha vartheta varrho chi_bar xi_bar iota phi eps_p eta theta psi 
           w tau_u z_ss pi_ss i_bar fs_ss fv_ss b_bar_ss B l_ss e_ss u_ss s_ss k_ss v_ss  
           x_ss a_ss b_ss c_e_ss c_u_ss c_n_ss xi_ss chi_ss MRS_NC_ss Omega_ss Gamma_ss 
           P_ss p_a_ss p_b_ss q_ss tau_ss Y_ss D_ss dw_ss
           rho_ss rho_r sigma_r rho_b sigma_b rho_z sigma_z rho_i sigma_i  ; 


load param;
load vecSS;

set_param_value('beta', beta );
set_param_value('alpham', alpham);
set_param_value('alpha', alpha);
set_param_value('vartheta', vartheta);
set_param_value('varrho', varrho);
set_param_value('chi_bar', chi_bar);
set_param_value('xi_bar', xi_bar);
set_param_value('iota',  iota);
set_param_value('phi',   phi);
set_param_value('eps_p', eps_p);
set_param_value('eta',   eta);
set_param_value('theta', theta);
set_param_value('psi',   psi);
set_param_value('tau_u', tau_u);

set_param_value('w', w);
%set_param_value('w_ss', w_ss);
set_param_value('z_ss',  z_ss);
set_param_value('pi_ss', pi_ss);
set_param_value('i_bar', i_bar);
set_param_value('fs_ss', fs_ss);
set_param_value('fv_ss', fv_ss);
set_param_value('b_bar_ss', b_bar_ss);
set_param_value('B', B);
set_param_value('l_ss', l_ss);
set_param_value('u_ss', u_ss);
set_param_value('e_ss', e_ss);
set_param_value('s_ss', s_ss);
set_param_value('k_ss', k_ss);
set_param_value('v_ss',  v_ss);

set_param_value('x_ss', x_ss);
set_param_value('b_ss', b_ss);
set_param_value('a_ss', a_ss);
set_param_value('c_e_ss', c_e_ss);
set_param_value('c_u_ss', c_u_ss);
set_param_value('c_n_ss', c_n_ss);
set_param_value('xi_ss', xi_ss);
set_param_value('chi_ss', chi_ss);
set_param_value('MRS_NC_ss', MRS_NC_ss);
set_param_value('Omega_ss', Omega_ss);
set_param_value('Gamma_ss', Gamma_ss);

set_param_value('P_ss', P_ss);
set_param_value('p_a_ss', p_a_ss);
set_param_value('p_b_ss', p_b_ss);
set_param_value('q_ss', q_ss);
set_param_value('tau_ss', tau_ss);
set_param_value('Y_ss', Y_ss);
set_param_value('D_ss', D_ss);
set_param_value('dw_ss', dw_ss);

set_param_value('rho_ss', rho_ss);
set_param_value('rho_r', rho_r);
set_param_value('sigma_r', sigma_r);
set_param_value('rho_b', rho_b);
set_param_value('sigma_b', sigma_b);
set_param_value('rho_z', rho_z);
set_param_value('sigma_z', sigma_z);
set_param_value('rho_i', rho_i);
set_param_value('sigma_i', sigma_i);

model; 

% Euler
c_e^(-varrho)= beta*((1+i(+1))/pi(+1))*(e(+1)*(c_e(+1))^(-varrho) + u(+1)*(c_u(+1))^(-varrho)+(1-l(+1))*(c_n(+1))^(-varrho));

% BC and liquidity constraints
x=(b(+1)/P)+(1+i)*(a/P)-(1+i)*(b/P);
c_u=x+(tau_u/P);
c_n=x;
a(+1)/P=x+(1-tau)*e*(w/P)+(1-tau)*D+u*(tau_u/P) -(e*c_e+u*c_u+(1-l)*c_n);

% employment
e=rho*e(-1)+fs*s;

% searchers
s=l-rho*e(-1);

% participation condition


MRS_NC-Omega=fs*((1-tau)*(w/P) -chi_bar*(e^(1+iota))*(c_e^(varrho))-Gamma)+fs*beta*rho(+1)*(c_e(+1)^(-varrho))*((1-fs(+1))/fs(+1))*(MRS_NC(+1)-Omega(+1));

% where:

chi=chi_bar*(e^(1+iota))/(1+iota);
xi=xi_bar*log(1-l);
MRS_NC=xi_bar*(c_e)^(varrho);

%xi=xi_bar*((1-l)^(1-phi))/(1-phi);
%MRS_NC=xi_bar*((1-l)^(1-phi))/mu_ce;

Omega=(tau_u/P) + (c_n-c_u)+(c_e^varrho)*(log(c_u) - log(c_n) - xi);
Gamma=(tau_u/P) + (c_e-c_u)+(c_e^varrho)*(log(c_u) - log(c_e) + chi);

% aggregate employment and unemployment 
l=e+u;

% discount factor 
Lambda=beta*((1-tau(+1))/(1-tau))*(c_e(+1)/c_e)^(-varrho);

% hiring costs
k=z*B*(fs)^eta;

% asset market equilibrium
b(+1)/P=b_bar;
a(+1)/P=b_bar;

% Firms:

%optimal hiring
k/fv=(q/P)*z-(w/P)+rho(+1)*Lambda(+1)*k(+1)/fv(+1);

% dividends definition

dw=(q/P)*z*e-(w/P)*e-k*v;

% price
 
p_star=p_a/p_b;
p_a=((eps_p-1)/eps_p)*q*Y + Lambda(+1)*(1-theta)*p_a(+1)*(pi(+1)^(eps_p-1));
p_b=Y+Lambda(+1)*(1-theta)*p_b(+1)*(pi(+1)^(eps_p-1));

% inflation
%pi=((1-theta)/(1-theta*(p_star/P)^(1-eps_p)))^(1/(1-eps_p));

P=(theta*(p_star)^(1-eps_p)+(1-theta)*P(-1)^(1-eps_p))^(1/(1-eps_p));

pi=P/P(-1);

%Output
Y=z*e/varsigma;

% price dispersion
varsigma=theta*(p_star/P)^(-eps_p)+(1-theta)*(pi^(eps_p))*varsigma(-1);

% dividends
D=Y-(q/P)*z*e+dw;

% government
tau*(e*(w/P)+D)=u*tau_u;


% Taylor rule 
1+i(+1) = (1+i_bar)*(pi^(psi))*exp(eps);

% job finding rates
fs=alpham*(v/s)^(1-alpha);
fv=alpham*(v/s)^(-alpha);

% bargained wage - to simplify, suppose there is binding min wage so this condition does not bind
%w= vartheta*(q*z+Lambda(+1)*rho(+1)*(1-(c_e^(-varrho))*(1-fs(+1)))*k(+1)/fv(+1))+(1-vartheta)*((e^(1+iota))*chi_bar*(c_e^varrho)+Gamma)/(1-tau);

% exogenous variables

log(z)=(1-rho_z)*log(z_ss)+rho_z*log(z(-1))+sigma_z*e_z;
log(rho)=(1-rho_r)*log(rho_ss)+rho_r*log(rho(-1))+sigma_r*e_r;
b_bar=(1-rho_b)*b_bar_ss+rho_b*(b_bar(-1))+sigma_b*e_b;
eps=rho_i*eps(-1)+sigma_i*e_i;

end;

initval; 
z=z_ss;
pi=pi_ss;
Lambda=beta;
i=i_bar;
rho=rho_ss;
b_bar=b_bar_ss;
varsigma=1;
fs=fs_ss;
fv=fv_ss;
l=l_ss;
u=u_ss;
e=e_ss;
s=s_ss;
k=k_ss;
v=v_ss;
x=x_ss;
b=b_ss;
a=a_ss;
%w=w_ss;
tau=tau_ss;
c_e=c_e_ss;
c_u=c_u_ss;
c_n=c_n_ss;
xi=xi_ss;
chi=chi_ss;
MRS_NC=MRS_NC_ss;
Omega=Omega_ss;
Gamma=Gamma_ss;
p_star=P_ss;
P=P_ss;
q=q_ss;
Y=Y_ss;
p_a=p_a_ss;
p_b=p_b_ss;
D=D_ss;
dw=dw_ss;
end;


shocks; 
var e_z;  stderr sigma_z; 
var e_r; stderr sigma_r; 
var e_b;  stderr sigma_b;
var e_i;  stderr sigma_i;

end; 


resid(1);
steady;

check;
stoch_simul(order=1, irf=40, nograph) c_e ;