%=======================================================================
%    DSGE model with exog. govt. spending & fixed nominal interest rate
%    in a monetary and fiscal union - Nakamura and Steinsson's 2014 base
%    line model with "Volcker Greenspan" monetary policy rule.
%=======================================================================

% Specification à la Woodford (2003, 2005). Firms own their own capital
% stock and face adjustment costs at the firm level in adjusting their
% capital stock. Households' beaviour is governed by the same equations
% as in our baseline GHH model of section E. Firms belong to industries
% x, which make use of a specific type of labour.

var c c_h c_f c_abr 
p p_abr p_h p_f pi_h pi_f pi pi_abr s s_h k k_tilde
q r g_h g_f L L_abr tau y y_abr k_h i_h u_c p_tilde;

parameters a n nu  mu beta alpha sigma eta theta phi_h phi_f phi_h_abr 
phi_f_abr rho rho_pi rho_y rho_g phi_pi phi_y phi_g  psi_nu zeta zeta_c 
chi zeta_g kappa C_bar Y_bar G_bar tau_bar A_c a_c a_pi L_bar sigma_c 
sigma_l xi_l xi_y r_bar epsilon_psi omega_bar nu_bar Theta L Q(L) A phi;

nu = 1;
epsilon_psi = 3;
sigma = 0.75;
beta = 0.99;
theta = 7;
eta = 2;
a = 2/3;
delta = 0.025;
alpha = 0.75;
phi_h = 0.69;
phi_f = 0.31;
n = 0.1;
phi_h_abr = n/(1-n)*phi_f;
phi_f_abr = 1- phi_h_abr;
kappa =(1-alpha)*(1-alpha*beta)/alpha;
Y_bar = 1;
G_bar = 0.2;
C_bar=Y_bar-G_bar; 

tau_bar = 0.35;  
rho = 0.8;
rho_pi = 0.75;
rho_y = 0.85;
rho_g = 0.933;
psi_nu = (1+nu*(1-a))/(nu*a);
zeta = 1/(1+psi_nu*theta);
phi_pi = 1.5;
phi_y = 0.5;
phi_g = 0.2;

sigma_c = (1/sigma)*(1-a*((theta-1)/theta)*(1/C_bar)*(1+1/nu)^(-1))^(-1);
xi_y = (1/C_bar)*((theta-1)/theta);

omega_bar = (1/nu + 1 -a)/a;
nu_bar = 1/nu;

Theta = (1 - beta*(1 - delta))*omega_bar*(y - y_h)*theta*epsilon_psi^(-1);
Q(L) = beta - (1 - beta + (1 - beta*(1 - delta))*(omega_bar)*(k - k_h)*epsilon_psi^(-1))*L + L^2;

A = 1 + beta + (1 - beta*(1 - delta))*rho_k*epsilon_psi^(-1);
phi = 1 + omega_bar*theta - (omega_bar - nu_bar)*(gamma*alpha*beta)/(1 - alpha*beta*lambda);

model

k(+1) = (1 - delta)*k + delta*i;
k_h(+1) = (1 - delta)*k_h + delta*i_h;
s = (1/nu + 1 - a)*l - (1 - a)*k - p_h + tau_bar/(1 - tau_bar)*tau;
r_k = s + a*l - a*k;
u_c + epsilon_psi*(k(+1) - k) + u_c(+1) + beta*(k(+2) - k(+1)) + (1 - beta*(1 - delta))*(r_k(+1) + p_h(+1));
y = a*l + (1 - a)*k;
s_h = omega_bar*y_h - (omega_bar - nu_bar)*k_h - p_h + tau_bar/(1 - tau_bar)*tau;
s = s_h + omega_bar*(y - y_h) - (omega_bar - nu_bar)*(k - k_h);
1/(1 - alpha*beta)*(p - pi_h(+1) - s(+1)) = 0;
1/(1 - alpha*beta)*(p - pi_h(+1) - s_h(+1) + omega_bar*theta*(p - pi(+1)) + (omega_bar - nu_bar)*(k(+1) - k)) = 0;
Theta*p(+1) = (Q(L)*(k(+2) - k(+1)));
p = p_h - psi*k;
pi_h = (1 - alpha)/alpha*p_h;
p_tilde(+1) = alpha*(p_tilde - pi_h(+1)) + (1 - alpha)*p(+1) ;
p_tilde(+1) = alpha*(p_tilde - pi_h(+1)) + (1 - alpha)*(p_h(+1) - psi*k_tilde(+1));
k_tilde(+1) = lambda*k_bar - gamma*p_tilde;
k_tilde(+2) = lambda*k_bar(+1) - Tau*p_tilde(+1) = (lambda + (1 - alpha)*gamma*psi)*k_tilde(+1) - alpha*gamma*p_tilde(+2);
Theta*alpha*p - Theta*(1 - alpha)*psi*k_tilde(+1) = 1/lambda*k_tilde(+1) + 1/lambda*gamma*p - A*k_tilde(+1) + beta*(lambda + (1 - alpha)*gamma*psi)*k_tilde(+1) - alpha*beta*theta*p_tilde;
(1 - alpha*beta*lambda)*gamma = Theta*alpha*lambda;
(1/beta - alpha*lambda)*Q(beta*lambda) + (1 - alpha)*Theta*psi*lambda;
k_tilde(+1) = lambda*k_tilde(+1) - gamma*(p_tilde - pi_h(+1));
(k(+1) - k_h(+1)) = lambda*k_tilde - gamma*p_tilde;
1/(1 - alpha*beta)*k_tilde(+1) = k_tilde/(1 - alpha*beta*lambda) - gamma*alpha*beta/(1 - alpha*beta*lambda)*1/(1 - alpha*beta)*p_tilde(+1);
1/(1 - alpha*beta)*p_tilde(+1) = 1/(1 - alpha*beta)*pi_h(+1)*(p_tilde - pi_h(1));
(1 + omega_bar*theta)*p = s_h(+1) + (1 - alpha*beta)*(1 + omega_bar*theta);
1/(1 - alpha*beta)*pi_h(+1) = (1 - alpha*beta)(omega_bar - nu_bar)/(1 - alpha*beta*lambda)*k_tilde + (gamma*alpha*beta*(omega_bar - nu_bar))/(1 - alpha*beta*lambda)*1/(1 - alpha*beta)*pi(+1);
phi*p = s_h(+1) + 1/(1 - alpha*beta)*pi_h(+1) - (omega_bar - nu_bar)*(1 - alpha*beta)/(1 - alpha*beta*lambda)*k_tilde;
phi*p_h = (1 - alpha*beta)/(1 - alpha*beta)*s_h(+1) + phi/(1 - alpha*beta)*pi_h(+1);
phi(lambda) = (omega_bar - nu_bar)*(1 - alpha*beta)/(1 - alpha*beta*lambda);
psi(lambda) = (1/beta - alpha*lambda)*Q(beta*lambda)/(1 - alpha)*Theta*lambda;
gamma(lambda) = Theta*alpha*lambda/(1 - alpha*beta*lambda);
V(lambda) = ((1 + omega_ber*theta)*(1 - alpha*beta*lambda)^2 - alpha^2*beta*(omega_bar - nu_bar)*Theta*lambda);
Q(beta*lambda) = beta*(1 - alpha)*(1 - alpha*beta)*(omega_bar - nu_bar)*Theta*lambda = 0;
p_h - alpha*beta*p_h(+1) = (1 - alpha)*1/phi*s_h + alpha*beta*pi_h(+1);
alpha/(1 - alpha)*pi_h - alpha^2*beta/(1 - alpha)*pi_h(+1) = (1 - alpha*beta)*1/(phi)*s_h + alpha*beta;
pi_h = kappa*phi^(-1)*s_h + beta*pi_h(+1);

initval; 
c=0;
c_abr=0;
p=0;
d=0;
d_abr=0;
r_k_abr=0;
r=0;
r_k=0;
i=0;
i_abr=0;
i_f=0;
i_h=0;
u=0;
u_abr=0;
k=0;
k_abr=0;
k_bar=0;
k_bar_abr=0;
s_h=0;
s_f=0;
s=0;
pi=0;
pi_h=0;
pi_f=0;
pi_abr=0;
l=(y/(k^(1-a)))^(1/a);
l_abr=0;
tau=0;
c_h=0;
c_f=0;
p_h=0;
p_f=0;
tau=0;
y=0;
y_h=0;
y_f=0;
y_abr=0;
q=0;
g_h=0;
g_f=0;
r=0;
i_h_abr=0;
i_f_abr=0;

end;  

steady(maxit=1000);

check;// (qz_zero_threshold=1e-20);

shocks;
var e; 
stderr 0.01;
end;

stoch_simul(order=1,IRF=25) c, g_h, pi, r, y, y_abr, c_h, c_f, p_h, p_f k_h k_f;