%==========================================================
%    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.
%==========================================================

%% Regional capital markets analogous to the specification of
% Christiano, Eichenbau and Evans (2005). Household preferences are given
% by equations (15) - (17). Households decisions about consumption,
% saving, and labour supply are the same as before. However, in addition to
% these choices, households own the capital stock, they choose how much to
% invest and they choose the rate of utilization of capital stock.

var d d_abr r_k r_k_abr pi_h pi_f pi_abr q l k k_bar l_abr y_abr y k_abr u y_f 
y_h i_f i_h i_f_abr i_h_abr i i_abr s s_abr c_abr p_h p_f p c_h_abr c_f_abr 
b tau w c_h c_f pi c M r g_h g_f u_abr k_bar_abr g s_h s_f g_abr;

varexo e;

parameters a n nu 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 delta 
kappa C_bar Y_bar G_bar tau_bar A_c a_c a_pi L_bar psi_k psi_c 
u_bar sigma_a kappa_i A_double A_prime A Chi Phi I_bar zeta 
rho_r sigma_c zeta_c zeta_g;

nu = 1;
sigma = 1.5;
beta = .99;
theta = 7;
eta = 2;
a = 2/3;
alpha = .75;
phi_h = .69;
phi_f = .31;
n = .1; 
phi_h_abr = n/(1-n)*phi_f;
phi_f_abr = 1- phi_h_abr;
kappa = 1;
L_bar = 0.788; 
Y_bar = L_bar^a;  
G_bar = .2*Y_bar;
C_bar = .2634;
I_bar = n*Y_bar;
tau_bar = .27;
rho = .8;
rho_pi = .75;
rho_y = .85;
rho_g = .933;
rho_r = 1;
phi_pi = 1.5;
phi_y = .5;
phi_g = 0;
Chi = 1;
zeta = 1/(1+psi_nu*theta);
zeta_c = zeta*(sigma^(-1)+(C_bar/Y_bar)*psi_nu);
zeta_g = zeta*psi_nu;
A_c = (1-rho_g)*(1-beta*rho_g)-rho_g*kappa*zeta_c;
a_c = (rho_g*kappa*zeta_g)/A_c;
a_pi = kappa*(zeta_c/(1-beta*rho_g))*a_c+kappa*(zeta_g/(1-beta*rho_g));
psi_nu = (a/nu)/((1-a)/nu+1);
psi_k = (1+1/nu)*(1-a)/((1-a)/nu+1);
psi_c = a/((1-a)/nu+1);
delta = .025;
u_bar = 1;
A = 1;
A_double = 1;
A_prime = 1;
sigma_a = A_double/A_prime;
Phi = 1;
kappa_i = 0.25;


model;

M = 1/(1 + r);
0 = c - c(+1) + sigma*r - pi(+1);
0 = c - c_abr - sigma*q;
r = rho_r*r(-1) + (1-rho)*(phi_pi*pi + phi_y*y + phi_g*g);
p_h = s_h(+1) + alpha*beta/(1-alpha*beta)*pi(+1);
p_f = s_f(+1) + alpha*beta/(1-alpha*beta)*pi(+1);
phi_h*p_h + phi_f*p_f =0;
pi = phi_h*pi_h + phi_f*pi_f;
phi_h_abr*pi_h + phi_f_abr*p_f = q;
pi_abr = phi_h_abr*pi_h + phi_f_abr*pi_f;
p_h = p_h(+1) - pi(-1);
p_f = p_f(+1) - pi(-1);
pi = n*pi_h + (1-n)*pi_abr;
g = n*g + (1-n)*g_abr;
g = rho*g(-1) + e;
c_h = c - eta*p_h;
c_f = c - eta*p_f;
0 = c_h_abr - c_abr + eta*(p_h - q);
0 = c_f_abr - c_abr + eta*(p_f - q);
A_prime*u = r_k - p;
y = (n*c_h + (1-n)*c_h_abr + n*i_h + (1-n)*i_abr + n*g_h)*(p_h/p);
3*p + c + i + A_prime*(u_bar)/A*(u_bar)*u + k_bar + M(+1) + b(+1) = b - 1/(1 - tau_bar)*tau 
+ w + l + r_k + u + Chi;
s = (1/nu + 1 - a)*(y + (1 - a)/a*r_k - (1 - a)/a*s) - (1 - a)*(y - r_k + s) + 1/sigma*c + tau_bar/(1 - tau_bar)*tau;
s_abr = psi_nu*y_abr + 1/sigma*psi_c*c_abr + psi_k*r_k_abr + psi_c*q + psi_c*tau_bar/(1-tau_bar)*tau;
i_h = phi_h*i;
i_f = phi_f*i;
pi_h = beta*pi_h(+1) + kappa*zeta*psi_nu*y_h - kappa*zeta*p_h + kappa*zeta*1/sigma*c + kappa*zeta*(tau_bar/(1-tau_bar))*tau;
p = (phi_h*p_h^(1-eta) + phi_f*p_f^(1-eta))^(eta/(1-eta));
i = phi_h*i_h + phi_f*i_f;
i_abr = phi_h_abr*i_h_abr + phi_f_abr*i_f_abr;
i_h = i - eta*p_h;
i_f = i - eta*(p_f-q);
y_h = phi_h*(C_bar/Y_bar)*c_h + phi_h*(I_bar/Y_bar)*i_h + phi_h_abr*(1-n)/n*(C_bar/Y_bar)*c_h_abr + phi_h_abr*n/(1-n)*(I_bar/Y_bar)*i_abr + g_h;
y_f = phi_f_abr*(C_bar/Y_bar)*c_f_abr + phi_f_abr*(I_bar/Y_bar)*i_f_abr + phi_f*n/(1-n)*(C_bar/Y_bar)*c_f + phi_f_abr*n/(1-n)*(I_bar/Y_bar)*i_f;
k = u + k_bar;
k_abr = u_abr + k_bar_abr;
y = a*l + (1-a)*k;
r_k = a*l - a*k + s;
y_abr = a*l_abr + (1-a)*k_abr;
r_k_abr = a*l_abr - a*k_abr + s_abr;
k_bar(+1) = (1 - delta)*k_bar + (1 - Phi*(i - i(+1))) + i;
k = y - r_k + s;
l = y + (1 - a)/a*s;
pi_h = beta*pi_h(+1) + kappa*zeta*(psi_nu*y_h - p_h + 1/sigma*psi_c*c + psi_c*tau_bar/(1 - tau_bar)*tau + psi_k*r_k);
pi_f = beta*pi_f(+1) + kappa*zeta*(psi_nu*y_f - p_f + 1/sigma*psi_c*c_abr + psi_c*tau_bar/(1 - tau_bar)*tau + psi_k*r_k_abr + psi_c*q);
sigma_a*u = r_k;
sigma_a*u_abr = r_k_abr + q;
d = beta*(1 - delta)*d(+1) + (1 - beta*(1 - delta))*(r_k(+1) - 1/sigma_c*c(+1));
d_abr = beta*(1 - delta)*d_abr(+1) + (1 - beta*(1 - delta))*(r_k(+1) + pi(+1) - q - pi_abr(+1) - 1/sigma_c*c_abr(+1));

end;

initval; 
c=0;
c_abr=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=0;
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;
M=0;
b=0;
q=0;
g_h=0;
g_f=0;
g_abr=0;
w=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;