%----------------------------------------------------------------
% 0. Housekeeping (close all graphic windows)
%----------------------------------------------------------------
close all;

%----------------------------------------------------------------
% 1. Defining variables
%----------------------------------------------------------------

// Endogenous variables (11)
var  c i l y k b r_w tau_l tau_k tau_b tfp;

// Exogenous variables
varexo u_tfp u_l u_k u_b;

// Parameters
parameters  b_y delta eta rbar alpha nu rho_tfp rho_l rho_k rho_b gamma n psi beta phi B sigma_l sigma_k sigma_b sigma_tfp;

%----------------------------------------------------------------
% 2. Calibration
%----------------------------------------------------------------
n=0.0146;
delta=0.05;
phi=3.674;
gamma=0.0034;
rbar=1.04;
b_y=0.31;
v=0.0001;
alpha=0.3;
nu=1.6;
psi=3.530;
beta=0.954;
rho_tfp=0.987;
rho_l=0.995;
rho_k=0.884;
rho_b=0.994;
sigma_tfp=0.047;
sigma_l=0.041;
sigma_k=0.042;
sigma_b=0.000;

%----------------------------------------------------------------
% 3. Model (11 equations)
%----------------------------------------------------------------
model;
// Consumption leisure decision
psi*exp(c)/(1-exp(l))=exp(tfp)*(1-alpha)*exp(tau_l)*((exp(k(-1))/exp(l))^alpha);

// Capital Euler Equation
(1/exp(c))/(1-phi*(exp(i)/exp(k(-1))-delta-gamma-n-gamma*n))=beta/(1+gamma)*(1/exp(c(+1))*(exp(tau_k)*alpha*((exp(k)/exp(l(+1)))^(alpha-1))+1/(1-phi*(exp(i(+1))/exp(k)-delta-gamma-n-gamma*n))*((1-delta)-(phi/2)*(exp(i(+1))/exp(k)-delta-gamma-n-gamma*n)+phi*(exp(i)/exp(k(-1))-delta-gamma-n-gamma*n)*exp(i(+1))/exp(k))));

// Euler Equation for bond
1/exp(c)=beta/(1+gamma)*1/exp(c(+1))*exp(tau_b)*exp(r_w(+1));

// Accumulation of net foreign assets
(1+gamma)*(1+n)*exp(b(+1))=exp(c)+exp(i)+exp(b)*exp(r_w)-exp(y);

// Law of motion of capital
(1+gamma)*(1+n)*exp(k)=(1-delta)*exp(k(-1))+exp(i)-(phi/2)*(exp(i)/exp(k(-1))-delta-gamma-n-gamma*n)*exp(k(-1));

// Supply of funds
exp(r_w(+1))=rbar*((exp(b(+1))/B)^eta);

// Definition of output
exp(y)=exp(tfp)*(exp(k(-1))^alpha)*(exp(l)^alpha);

// Productivity
log(exp(tfp))=rho_tfp*log(exp(tfp(-1)))+u_tfp;

// Labor wedge
log(exp(tau_l))=rho_l*log(exp(tau_l(-1)))+u_l;

// Capital wedge
log(exp(tau_k))=rho_k*log(exp(tau_k(-1)))+u_k;

// Bond wedge
log(exp(tau_b))=rho_b*log(exp(tau_b(-1)))+u_b;

end;

%----------------------------------------------------------------
% 4. Steady State
%----------------------------------------------------------------
steady_state_model;
tau_l_ss=1;
tau_k_ss=1;
tau_b_ss=1;
tfp_ss=1;
b_to_y_ss=b_y;
r_w_ss=(1+gamma)/(beta*tau_b_ss);
k_to_y_ss=alpha*tau_k_ss/(((1+gamma)/beta)-1+delta);
y_to_k_ss=(((1+gamma)/beta)-1+delta)/(alpha*tau_k_ss);
i_to_y_ss=k_to_y_ss*(gamma+nu+gamma*nu+delta);
c_to_y_ss=1-i_to_y_ss-b_to_y_ss*(gamma+nu+gamma*nu+delta);
l_ss=1/(1+c_to_y_ss*psi/((1-alpha)*tau_l_ss));
k_to_l_ss=y_to_k_ss^(1/(alpha-1));
k_ss=k_to_l_ss*l_ss;
y_ss=(k_ss^alpha)*(l_ss^(1-alpha));
b_ss=b_to_y_ss*y_ss;
B=b_ss;
c_ss=c_to_y_ss*y_ss;
i_ss=i_to_y_ss*y_ss;

tau_l=log(tau_l_ss);
tau_k=log(tau_k_ss);
tau_b=log(tau_b_ss);
tfp=log(tfp_ss);
r_w=log(r_w_ss);
l=log(l_ss);
k=log(k_ss);
y=log(y_ss);
b=log(b_ss);
c=log(c_ss);
i=log(i_ss);

end;

%----------------------------------------------------------------
% 5. Computation
%----------------------------------------------------------------
check;
steady; 


shocks;
var u_tfp=sigma_tfp^2;
var u_l=sigma_l^2;
var u_k=sigmau_k^2;
var u_b=sigmau_b^2;
end;


stoch_simul(order=1, irf=24, nograph);
save final04102.mat;