%% RBC_Model Chapter 4, Calibration for US

%% define parameters and do model calibration
sigma=1;
alpha   = 0.33;
i_y=0.25;
k_y=10.4;
x=0.0055;
n=0.0027;
gammax=(1+n)*(1+x);
delta=i_y/k_y-x-n-n*x;
beta=(1+x)*(1+n)/(alpha/k_y+(1-delta));
rhoz     = 0.97;  
rhog    = 0.98;
gshare     = 0.2038;

%% compute steady states
l_ss=0.33;
k_ss = ((1/beta*(1+n)*(1+x)-(1-delta))/alpha)^(1/(alpha-1))*l_ss; 
i_ss = (x+n+delta+n*x)*k_ss;
y_ss=k_ss^alpha*l_ss^(1-alpha);
g_ss=gshare*y_ss;
c_ss = (1-gshare)*k_ss^(alpha)*l_ss^(1-alpha)-i_ss;
psi=(1-alpha)*(k_ss/l_ss)^alpha*(1-l_ss)/c_ss^sigma;

%% define coefficients and auxiliary variables
gamma_l=(l_ss/(1-l_ss)+alpha)^-1;

alpha1=1/gammax*(alpha*y_ss/k_ss*(1+(1-alpha)*gamma_l)+(1-delta))
alpha2=-(c_ss/(gammax*k_ss)+y_ss/(gammax*k_ss)*(1-alpha)*gamma_l*sigma)
alpha3=y_ss/(gammax*k_ss)*(1+(1-alpha)*gamma_l)
alpha4=-g_ss/(gammax*k_ss)

alpha5=-beta/(gammax*sigma)*alpha*(k_ss/l_ss)^(alpha-1)*(alpha-1)*(1-alpha*gamma_l)
alpha6=1-beta/gammax*alpha*(k_ss/l_ss)^(alpha-1)*(alpha-1)*gamma_l
alpha7=-beta/(gammax*sigma)*alpha*(k_ss/l_ss)^(alpha-1)*(1-(alpha-1)*gamma_l)

phi_kk=(1-alpha5*alpha2+alpha6*alpha1)/(2*alpha6)-sqrt(((1-alpha5*alpha2+alpha6*alpha1)/(2*alpha6))^2-alpha1/alpha6)
phi_ck=phi_kk/alpha2-alpha1/alpha2
phi_cz=(rhoz*alpha7+alpha6*phi_ck*alpha3+alpha3*alpha5)/(1-alpha6*rhoz-alpha6*phi_ck*alpha2-alpha5*alpha2)
phi_kz=alpha3+alpha2*phi_cz
phi_cg=(alpha6*phi_ck+alpha5)*alpha4/(1-(alpha6*phi_ck+alpha5)*alpha2-alpha6*rhog)
phi_kg=alpha4+alpha2*phi_cg



%% Solve Model using eigenvalue decomposition
%% matrices of linear RE model: A*E_t(x_{t+1})=B*(x_{t})
% here; x=[k_t z_t g_t c_t]
A= [1       0           0           0; 
    alpha5  alpha7      0           alpha6;
    0       0           1           0 
    0       1           0           0
   ];
    
B= [alpha1  alpha3  alpha4  alpha2 ;
    0       0           0       1
    0       0           rhog    0;
    0       rhoz        0       0];


%% multiplying with inv(A) 
W=A\B;

%% partition matricies (the first three elements in x are states)
%% we need w11 and w12 for M
n_x=3;
n_x1=1; %endogenous states
n_x2=2; %exogenous states
n_u=1;

w11=W(1:n_x,1:n_x);
w12=W(1:n_x,n_x+1:end);

%% jordan decomposition: W=D*\Lambda*D^{-1}
[evec,eval]=eig(W);
[mu,ind]=sort(abs(diag(eval)));% sort eigenvalues in increasing order
dinv=inv(evec(:,ind));         % take inverse of matrix D

if sum(mu<1) ~= n_x % check BK condition
    alarm_flag=1;
    F = NaN;
    H = NaN;
    return
end

%% partition D^{-1} 
q21=dinv(n_x+1:end,1:n_x);               
q22=dinv(n_x+1:end,n_x+1:end);                 

%% Model solution 
M= w11-w12*q21/q22;            % transition under optimal policy 
PI = -q22\q21;      % policy function mapping state into control

disp('The model solution for the state vector M is:');
disp(M)
disp('The policy function PI is:');
disp(PI)

n_irf=40;
x=zeros(n_x,n_irf);
cons(1,1)=PI*x(:,1);
eta_mat=[zeros(n_x1,n_x2);eye(n_x2,n_x2)];
eps_vec=[1;0]; %eps_vec=[eps_z_t; eps_g_t]

x(:,2)=M*x(:,1)+eta_mat*eps_vec;
cons(1,2)=PI*x(:,2);

for ii=3:n_irf
    x(:,ii)=M*x(:,ii-1);
    cons(1,ii)=PI*x(:,ii);
end
capital=x(1,:);
technology=x(2,:);
g_spend=x(3,:);

figure

subplot(4,1,1)
plot(-1:n_irf-2,technology)
title('z')
axis tight

subplot(4,1,2)
plot(-1:n_irf-2,capital)
title('k')
axis tight

subplot(4,1,3)
plot(-1:n_irf-2,cons)
title('c')
axis tight

subplot(4,1,4)
plot(-1:n_irf-2,g_spend)
title('g')
axis tight
