var y_hat k_hat n_hat c_hat i_hat g_hat V_hat x_hat m_hat u_hat v_hat q_hat s_hat
    kp_hat ez ei lam_hat r_hat pi_hat lam eb  qk_hat rk_hat pw_hat a_hat khi_hat
mu_hat epsilon_hat w_hat wo_hat b_hat ew eeta er eg theta_hat ep;

varexo tieta tib tii tiz tig tir tip tiw;

parameters alpha  yV yx sigma n_ss u_ss xi delta gammaz h_til beta etaV kawV etak
beta_til xia xiw xilam  xii khi rho lambda gamma mu phia phis phix phib phikhi s b a_ss
w_ss pw_ss H_ss kappa y_ss  k_ss x_ss rhoeta phieta gammab gammao gammaf
tau2 phi sig tau tau1 TAU lb lo lf gammap phip sigp taup lambdap epsilonp rpi ry rhos 
 rhob rhoi rhoz rhog rhor rhop sigmaz sigmai sigmag sigmab sigmar sigmaw sigmap h
eta epsilon rk_ss v_ss sigmam c_y_ss i_y_ss g_y_ss psi ;

%shocks
rhob = 0.5;
rhoi = 0.5;
rhoz = 0.5;
rhog = 0.5;
rhor=0.5;
rhop=0.5;
rhow=0.5;
rhoeta=0.5;
sigmab = 0.15;
sigmai = 0.15;
sigmaz = 0.15;
sigmag = 0.15;
sigmap=0.15;
sigmaw=0.15;
sigmar=0.15;

%calibrated
alpha=0.33;
sigma=0.5;
rho=0.895;
delta=0.025;
xi=10;
beta=0.99;
s=0.95;
sigmam=1;
g_y_ss=0.2;

%estimated
kawV=0.5;
etak=4;
h=0.5;
eta=0.5;
b_til=0.5;
lambda=0.75;
lambdap=0.66;
gamma=0;
epsilonp=1.15;
rpi=1.7;
ry=0.125;
rhos=0.75;
gammaz=1;
b=0.5;%guess
etaV=1;
rhoeta=0.5;%guess
gammap=1;%guess
H_ss=0.39;%guess
kappa=3;%guess

%steady 
rk_ss=(gammaz/beta)-1+delta;
x_ss=1-rho;
u_ss=x_ss/(s+x_ss);
v_ss=u_ss*((s/sigmam)^(1/(1-sigma)));
n_ss=1-u_ss;
qk_ss=1;
k_ss=n_ss*((rk_ss/alpha)^(1/(alpha-1)));
mu=1/(1-lambda*beta);
epsilon=1/(1-rho*lambda*beta);
khi=eta/(eta+(1-eta)*mu/epsilon);
a_ss=(1-alpha)*(k_ss/n_ss)^(alpha);
w_ss= khi*(a_ss+beta*(kappa/2)*(x_ss^2)+beta*kappa*s*x_ss)+(1-khi)*b;

n_y_ss=(k_ss/n_ss)^(-alpha);
i_y_ss=(1-((1-delta)/gammaz))*gammaz*(k_ss/n_ss)^(1-alpha);
c_y_ss=1-g_y_ss-i_y_ss-((kappa/2)*(x_ss^2)*n_y_ss);
pw_ss=(1/a_ss)*(kappa*x_ss+w_ss-(beta*(kappa/2)*x_ss^2)-beta*rho*kappa*x_ss);
y_ss = (k_ss^alpha)*(n_ss^(1-alpha));

h_til=h/gammaz;
xii=(kappa*x_ss)^(-1);
xiw=xii*w_ss;

psi = khi*mu*beta*lambda*mu + (1-khi)*rho*beta*lambda*epsilon;%
tau =psi/(1+psi);



TAU=(1-eta*x_ss*beta*lambda*mu)*(1/eta)*mu*xiw;
phix=khi*beta*kappa*(x_ss^(2))*w_ss^(-1);
phikhi=khi*(1-khi)^(-1)*kappa*x_ss*w_ss^(-1);
phis=(1-khi)*s*beta*H_ss*w_ss^(-1);
tau1=(xiw*mu*phix+phikhi*(1-khi)*(x_ss*beta*lambda)*(xiw*mu)*mu*(rho*beta)+phis*TAU)*(1-tau);

tau2=(-(xiw*mu)*phikhi*(1-khi)*(x_ss*beta*lambda)*mu*(1-tau));
sig=(1-lambda)*(1-tau)*(1/lambda);
phi=(1+tau2)+sig+(tau*lambda^(-1)-tau1);
phip=1+beta*gammap;
lb=gammap*(1/phip);
sigp=(1-lambdap)*(1-lambdap*beta)/lambdap;
taup=1+(epsilonp-1)*xi;
lo=(sigp/taup)*(1/phip);
lf=beta/phip;
gammao=sig/phi;
gammaf= (tau*lambda^(-1)-tau1)*phi^(-1);
gammab= (1+tau2)*(1/phi);


phia= khi*pw_ss*a_ss*w_ss^(-1);
phib=(1-khi)*b*w_ss^(-1);

yV=rk_ss*k_ss/y_ss;
yx=(kappa/2)*(x_ss^2*n_ss/y_ss);


beta_til= beta/gammaz;
xia=xii*pw_ss*a_ss;

xilam=beta*(1+rho)/2;
phieta=phikhi*(1-khi)*(1-eta)^(-1);






model(linear);


%technology(1)
y_hat=alpha*k_hat+(1-alpha)*n_hat;
%ressource constraint(2)
y_hat=c_y_ss*c_hat+i_y_ss*i_hat+g_y_ss*g_hat+yV*V_hat+yx*(2*x_hat+n_hat(-1));
%Matching(3)
m_hat=sigma*u_hat+(1-sigma)*v_hat;
%employment dynamics(4)
n_hat=n_hat(-1)+(1-rho)*x_hat;
%transition probabilities (5-6)
q_hat=m_hat-v_hat;
s_hat=m_hat-u_hat;
%Unemployment(7)
u_hat=-(n_ss/u_ss)*n_hat(-1);
%effective capital(8)
k_hat=V_hat+kp_hat(-1)-ez;
%physical capital dynamics (9)
kp_hat=xi*(kp_hat(-1)-ez)+(1-xi)*(i_hat+ei);
%aggregate vacancies(10)
x_hat=q_hat+v_hat-n_hat(-1);
%consumpion saving(11)
0=lam_hat(+1)+(r_hat-pi_hat(+1))-ez(+1);
%marginal utility(12)
lam=(1/((1-h_til)*(1-beta*h_til)))*(h_til*(c_hat(-1)-ez)-(1+beta*h_til^2)*c_hat
+beta*h_til*(c_hat(+1)+ez(+1))+(1-h_til)*(eb-beta*h_til*eb(+1)));
%capital utilization(13)
V_hat=etaV*rk_hat;
%investment(14)
i_hat=1/(1+beta)*(i_hat(-1)-ez)+(1/(etak*((gammaz)^2)))/(1+beta)*(qk_hat+ei)+
(beta/1+beta)*(i_hat(+1)+ez(+1));
%Capital renting(15)
rk_hat=pw_hat+y_hat-k_hat;
%Tobin's q (16)
qk_hat=beta_til*(1-delta)*qk_hat(+1)+(1-beta_til*(1-delta))*rk_hat(+1)-(r_hat-pi_hat(+1));
%Aggregate hiring(17)
x_hat=xia*(pw_hat+a_hat)-xiw*w_hat+xilam*lam_hat+beta*x_hat(+1);
%marginal product of labor (18)
a_hat=y_hat-n_hat;
%weight in Nash bargaining(19)
khi_hat=-(1-khi)*(mu_hat-epsilon_hat);
%(20)
epsilon_hat=(rho*lambda*beta)*(lam_hat-pi_hat(+1)+gamma*pi_hat+epsilon_hat-ez(+1));
%21
mu_hat= (x_ss*lambda*beta)*(x_hat(+1))-(x_ss*lambda*beta)*(xiw*mu)*mu*(w_hat+gamma*pi_hat
-pi_hat(+1)-ez(+1)-w_hat(+1))+(lambda*beta)*(mu_hat(+1)+lam_hat+gamma*pi_hat-pi_hat(+1)
-ez(+1));
%spillover free target wage (22)

wo_hat=phia*(pw_hat+a_hat)+(phis+phix)*x_hat(+1)+phis*s_hat(+1)+phib*b_hat+(phis+(phix/2))*lam_hat
+phikhi*(khi_hat-(rho-s)*beta*khi_hat(+1))+ew;
%agregate wage (23)
w_hat= gammab*(w_hat(-1)-pi_hat+gamma*pi_hat(-1)-ez)+gammao*wo_hat+gammaf*(w_hat(+1)+pi_hat(+1)
-gamma*pi_hat+ez(+1));

%phillips curve (24)
pi_hat= lb*pi_hat(-1)+lo*(pw_hat+ep)+lf*pi_hat(+1);
%taylor rule(25)

r_hat= rhos*r_hat(-1)+(1-rhos)*(rpi*pi_hat+ry*(y_hat-y_hat(-1)))+er;
//gov(26)
g_hat=y_hat+((1-g_y_ss)/g_y_ss)*eg;
%market tightness (27)

theta_hat=v_hat-u_hat;
%benefits (28)
b_hat=kp_hat;
%shock

ew=phieta*(1-(rho-s)*beta*rhoeta)*eeta;

eeta=rhoeta*eeta(-1)+tieta;
% preferences
eb=rhob*eb(-1)+tib;
% investment 
ei=rhoi*ei(-1)+tii;

% labour productivity    
ez=(1-rhoz)*gammaz+rhoz*ez(-1)+tiz;

% govenrment spendings    
eg=rhog*eg(-1)+tig;
%policy shock
er=rhor*er(-1)+tir;
ep=rhop*ep(-1)+tip;

lam_hat=lam(+1)-lam;

end;
initval;
y_hat=0;
k_hat=0;
 n_hat=0;
 c_hat=0;
 i_hat=0;
 g_hat=0;
 V_hat=0;
 x_hat=0;
 m_hat=0;
 u_hat=0;
 v_hat=0;
 q_hat=0;
 s_hat=0;
    kp_hat=0;
 ez=0;
 ei=0;
 lam_hat=0;
 r_hat=0;
 pi_hat=0;
 lam=0;
 eb =0;
 qk_hat=0;
 rk_hat=0;
 pw_hat=0;
 a_hat=0;
 khi_hat=0;

mu_hat=0;
 epsilon_hat=0;
 w_hat=0;
 wo_hat=0;
 b_hat=0;
 ew=0;
 eeta=0;
 er=0;
 eg=0;
 theta_hat=0;
 ep=0;

end;
steady;

shocks;

var tiz = sigmaz^2 ;     
var tii = sigmai^2;
var tig = sigmag^2 ;
var tib = sigmab^2 ;
var tir = sigmar^2;
var tiw= sigmaw^2;
var tip= sigmap^2;
end;

stoch_simul(order=1,IRF=25);