%==========================================================================
%           RBC capital utilization. Question 5. Paper parameters
%==========================================================================

clear all
close all
clc

%--------------------------------------------------------------------------
% 1. Parameters (Bergoing & Soto, 2002)
%--------------------------------------------------------------------------

par.T  = 200;                        % length of simulation

par.beta  = 0.97;                       % discount factor
par.alpha = 0.4;                       % curvature prod. function
par.depre = 0.048;                      % depreciation rate
par.rho   = 0.95;                       % persistence parameter on shocks

par.gamma  = 2;                        % coefficient of risk aversion
par.sigshock = 0.0095;             % values for standard deviations of shocks

%Micro calibration
%par.phi = 30;                       %scaling parameter of labor supply
%par.eta = 2;                       %inverse elasticity of labor supply

%Macro calibration
par.phi = 6;                       %scaling parameter of labor supply       
par.eta = 0.5;                       %inverse elasticity of labor supply

par.kappa= ((1-par.beta)/(par.depre*par.beta))+1; %kappa parameter of capital utilization

%steady state
u_ss      = ((1-par.beta)/(par.depre*par.beta*(par.kappa-1)))^(1/par.kappa) ; %steady state capital utilization
delta_ss  = par.depre*(u_ss^par.kappa); %steady state depreciation
k = (((1/par.beta)-1+delta_ss)/(par.alpha*(u_ss^par.alpha)))^(1/(par.alpha-1));  %steady state capital per hour
h_ss = (((1/par.phi)*(1-par.alpha)*(u_ss*k)^par.alpha)/(((u_ss*k)^par.alpha)-delta_ss*k))^(1/par.eta+par.gamma);
k_ss      = k*h_ss;  %steady state capital 
y_ss      = ((u_ss^par.alpha)*(k_ss^par.alpha))*h_ss;  %steady state output 
c_ss      = (((1-delta_ss)*k)+((u_ss^par.alpha)*(k^par.alpha)))*h_ss; %steady state consumption 
i_ss      = (k-(1-delta_ss)*k)*h_ss; %steady state investment


par.k_init  = log(k_ss);                       %initial value used in .mod file
par.y_init  = log(y_ss);                       %initial value used in .mod file
par.c_init  = log(c_ss);                       %initial value used in .mod file
par.i_init  = log(i_ss);                       %initial value used in .mod file
par.h_init  = log(h_ss); 
par.u_init   = log(u_ss); 
par.delta_init   = log(delta_ss); 

%save the values
save parametervalues par

%--------------------------------------------------------------------------
% 2. Call Dynare
%--------------------------------------------------------------------------

addpath(genpath('C:\dynare\5.0'));
dynare rbc_u noclearall
rmpath(genpath('C:\dynare\5.0\matlab'));

%--------------------------------------------------------------------------
%Impulse response functions to preference shocks
%--------------------------------------------------------------------------
%order that variables are declared
p_c = 1;
p_k = 2;
p_a = 3;
p_y = 4; 
p_i = 5;
p_h = 6;
p_u = 7;
p_delta = 8;

%dynamics of states to states
A = [oo_.dr.ghx(oo_.dr.inv_order_var(p_k),:);
    oo_.dr.ghx(oo_.dr.inv_order_var(p_a),:)];

%dynamics of shocks to states
B = [oo_.dr.ghu(oo_.dr.inv_order_var(p_k),:);
    oo_.dr.ghu(oo_.dr.inv_order_var(p_a),:)];

%dynamics of states to controlls
C = [oo_.dr.ghx(oo_.dr.inv_order_var(p_c),:);
    oo_.dr.ghx(oo_.dr.inv_order_var(p_y),:);
    oo_.dr.ghx(oo_.dr.inv_order_var(p_i),:);
    oo_.dr.ghx(oo_.dr.inv_order_var(p_h),:);
    oo_.dr.ghx(oo_.dr.inv_order_var(p_u),:);
    oo_.dr.ghx(oo_.dr.inv_order_var(p_delta),:)];

%dynamics of shocks to controlls
D = [oo_.dr.ghu(oo_.dr.inv_order_var(p_c),:);
    oo_.dr.ghu(oo_.dr.inv_order_var(p_y),:);
    oo_.dr.ghu(oo_.dr.inv_order_var(p_i),:);
    oo_.dr.ghu(oo_.dr.inv_order_var(p_h),:);
    oo_.dr.ghu(oo_.dr.inv_order_var(p_u),:);
    oo_.dr.ghu(oo_.dr.inv_order_var(p_delta),:)];

%steady state
ss = oo_.dr.ys;

T = 30;
Simp = zeros(3,T);
Ximp = zeros(6,T);

Simp(:,1) = B(:,2)*par.sigshock;
Ximp(:,1) = D(:,2)*par.sigshock;

for t = 2:T
    Simp(:,t) = A*Simp(:,t-1);
    Ximp(:,t) = C*Simp(:,t-1);
end


%impulse response functions I
plot(1:T,Simp(3,:),'LineWidth',2.5,'Marker','*')
hold on
plot(1:T,Ximp(4,:),'LineWidth',2.5)
hold on
plot(1:T,Ximp(5,:),'LineWidth',2.5,'LineStyle','--')
hold on
plot(1:T,Ximp(6,:),'LineWidth',2.5,'LineStyle',':')
hold off
axis tight
legend('A','H','W','R')
xlabel('Year','fontsize',13,'FontWeight','bold')
ylabel('% deviation from steady state','fontsize',13,'FontWeight','bold')
axis tight

%impulse response functions II
plot(1:T,Simp(3,:),'LineWidth',2.5,'Marker','o')
hold on
plot(1:T,Ximp(2,:),'LineWidth',2.5)
hold on
plot(1:T,Ximp(1,:),'LineWidth',2.5,'LineStyle','--')
hold on
plot(1:T,Ximp(3,:),'LineWidth',2.5,'LineStyle',':')
hold off
axis tight
legend('A','Y','C','I')
xlabel('Year','fontsize',13,'FontWeight','bold')
ylabel('% deviation from steady state','fontsize',13,'FontWeight','bold')
axis tight

%impulse response functions exam
plot(1:T,Simp(3,:),'LineWidth',1.5,'Marker','o')
hold on
plot(1:T,Ximp(1,:),'LineWidth',1.5)
hold on
plot(1:T,Ximp(4,:),'LineWidth',1.5,'LineStyle','--')
hold off
axis tight
legend('A','C','H')
xlabel('Year','fontsize',10,'FontWeight','bold')
ylabel('% deviation from steady state','fontsize',13,'FontWeight','bold')
axis tight