%==========================================================
%DSGE model with exog. govt. spending & fixed nominal interest rate
%in a monetary and fiscal union - Nakamura and Steinsson's 2014.
%Regional capital markets specification analogous to the specification of
%Christiano, Eichenbaum 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 
c c_abr c_h c_f c_h_abr c_f_abr y y_abr y_h y_f i_h i_f i_h_abr i_f_abr g 
g_h g_f q l l_abr k k_abr p_f p_h d d_abr r s s_abr s_h s_f pi_h pi_f pi 
pi_abr r_k r_k_abr i p tau i_abr u k_bar u_abr k_bar_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 L_bar u_bar A_c a_c a_pi psi_k psi_c sigma_a kappa_i 
A_double A_prime A Chi Phi I_bar zeta rho_r sigma_c
zeta_c zeta_g mu K_bar;

nu = 1;
sigma = 1;
beta = .99;
theta = 7;
mu = 1;
eta = 2;
delta = .025;
a = 2/3;
alpha = .75;
phi_h = .69;
phi_f = 1 - phi_h;
n = .1; 
phi_h_abr = n/(1-n)*phi_f;
phi_f_abr = 1- phi_h_abr;
kappa = 1;
L_bar = 0.788;
K_bar = 1;
Y_bar = L_bar^a*K_bar^(1-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 = .4;
Chi = 1;
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);
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));
u_bar = 1;
A = 1;
A_double = 1;
A_prime = 1;
sigma_a = A_double/A_prime;
Phi = 1;
kappa_i = 0.25;
sigma_c = (-(1/sigma)*(1 - a/mu*Y_bar/C_bar*1/(1+1/nu))^(-1))^(-1);

model;
k(+1) = (1 - delta)*k + i;
1 = i*(i(-1) - i)/i(-1);
0 = c - c(+1) + sigma*r - pi(+1);
0 = c - c_abr - sigma*q;
c_h = c - eta*p_h;
c_f = c - eta*p_f;

s = psi_nu*y + sigma^(-1)*psi_c*c + psi_k*r_k + psi_c*tau_bar/(1 - tau_bar)*tau;
s_abr = psi_nu*y_abr + sigma^(-1)*psi_c*c_abr + psi_k*r_k_abr + psi_c*q + psi_c*tau_bar/(1-tau_bar)*tau;

s_h = p_h(-1) - alpha*beta/(1-alpha*beta)*pi;
s_f = p_f(-1) - alpha*beta/(1-alpha*beta)*pi;

y_h = phi_h*(C_bar/Y_bar)*c + (1 - n)/n*(C_bar/Y_bar)*phi_h_abr*c_abr + (I_bar/Y_bar)*phi_h*i + (1 - n)/n*(I_bar/Y_bar)*phi_h_abr*i_abr - eta*(C_bar + I_bar)/Y_bar*(phi_h - (1 - n)/n*phi_h_abr)*p_h + eta*(C_bar + I_bar)/Y_bar*(1 - n)/n*phi_h_abr*q + g_h;
y_f = phi_f_abr*(C_bar/Y_bar)*c_abr + n/(1 - n)*(C_bar/Y_bar)*phi_f*c + (I_bar/Y_bar)*phi_f_abr*i_abr + n/(1 - n)*(I_bar/Y_bar)*phi_f*i - eta*(C_bar + I_bar)/Y_bar*(n/(1 - n)*phi_f + phi_f_abr)*p_h + eta*(C_bar + I_bar)/Y_bar*phi_h_abr*q + g_f;
y(+1) = y_h(+1) - theta*(p - p_h(+1) - pi(+1));
y_abr(+1) = y_f(+1) - theta*(p_f - p_f(+1) - pi(+1));

0 = c_h_abr - c_abr + eta*(p_h - q);
0 = c_f_abr - c_abr + eta*(p_f - q);

sigma_a*u = r_k;
sigma_a*u_abr = r_k_abr - q;

pi_h = beta*pi_h(+1) + kappa*zeta*psi_nu*y_h - kappa*zeta*p_h + kappa*zeta*1/sigma*c + kappa*zeta*psi_c*(tau_bar/(1-tau_bar))*tau;
pi_f = beta*pi_f(+1) + kappa*zeta*psi_nu*y_f - kappa*zeta*p_f + kappa*zeta*1/sigma*c + kappa*zeta*psi_c*(tau_bar/(1 - tau_bar))*tau;
pi = phi_h*pi_h + phi_f*pi_f;
pi_abr = phi_h_abr*pi_h + phi_f_abr*pi_f;

l = y + ((1 - a)/a)*r_k - ((1 - a)/a)*s;
k = y - r_k + s;

l_abr = y_abr + (1 - a)/a*r_k_abr - (1 - a)/a*s_abr;
k_abr = y_abr - r_k_abr + s_abr;

g = rho*g(-1) + e;
g_h = phi_f*g;
g_f = phi_f*g;

i_h = phi_h*i;
i_f = phi_f*i;
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);
i_h_abr = i_abr - eta*(p_h - q);
i_f_abr = i_abr - eta*(p_f - q);

k = u + k_bar;
k_abr = u_abr + k_bar_abr;

k_bar(+1) = (1 - delta)*k_bar + (1 - Phi*(i - i(+1))) + i;

d + kappa_i*(beta*(i(+1) - i) - (i - i(-1))) + (1/sigma_c)*c;
d_abr + kappa_i*(beta*(i_abr(+1) - i_abr) - (i_abr - i_abr(-1))) + (1/sigma_c)*c_abr;
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));
r = rho_r*r(-1) + (1-rho)*(phi_pi*(n*pi_h + (1-n)*pi_f) + phi_y*(n*y + (1 - n)*y_abr) + phi_g*(n*g_h + (1 - n)*g_f));

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;
s_abr=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;
q=0;
g_h=0;
g_f=0;
r=0;
i_h_abr=0;
i_f_abr=0;

end;  

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;