%  
% Simulation program for Chapter 18 
%
% Figure 18.2: Impulse-response functions
%    Permanent technology shock (xi_Z=1)
%
% B.J. Heijdra 
%
% 20 March 2017
% 
% The controlled parameters are:
%
xi_Z      = 1;
eta_0     = 1;
rho       = (1.065)^(0.25) - 1;
R0        = rho;
delta     = (1.1)^(0.25) - 1;
alpha     = 1/3;
ystar     = (rho+delta)/alpha;
omega_I   = delta/ystar;
omega_G   = 0.2;
omega_C   = 1 - omega_I - omega_G;
L0        = 0.2;   
omega_LL  = (1-L0)/L0;
epsilon   = omega_C / (omega_C + omega_LL*(1-alpha)); 
GAMY      = L0^(alpha-1) * ystar^(alpha);
phi       = (1+omega_LL) / (1 + alpha * omega_LL);

%
%  Jacobian matrix
%

zeta       = (rho+delta) / (1-delta+phi*(rho+delta));
theta_21   = zeta*( 1 - alpha * phi);
Delta11    =  ystar*( phi*alpha - omega_I );
Delta12    = -ystar*( omega_C + phi - 1 );
Delta21    = -theta_21*(1+Delta11); 
Delta22    = -theta_21*Delta12-zeta*(phi-1) ;
%
%  Matrix of shock vectors (system is y_t+1 - y_t = Delta*y_t + Gam_t)
%
gamma_K    = phi*ystar;
gamma_C    = phi*zeta*(xi_Z - ystar*(1-phi*alpha)); 
%
%  Determinant and characteristic roots
%
DDelta     =  Delta11*Delta22 - Delta12*Delta21;
TDelta     =  Delta11 + Delta22;
XX3        =  sqrt(TDelta*TDelta - 4*DDelta);
lambda_2   =  ( TDelta + XX3)/2;
lambda_1   =  (-TDelta + XX3)/2;
%
%  Results based on linearisation
%
%
%  Laplace transform of shocks
%
%  If xi_Z = 1 the technology shock is permanent
%  If xi_Z > 0 the technology shock is temporary
%
LAP        = eta_0/(lambda_2 + (1-xi_Z) );
%
%  Short-run results
%
ZZERO      =  eta_0 ;
GZERO      =  0;
KZERO      =  0;
CZERO      = - (gamma_C*Delta12 + gamma_K*(lambda_2-Delta22))*LAP/Delta12 ;
YZERO      =  phi*ZZERO -(phi-1)*CZERO ;
LZERO      =  (YZERO -  ZZERO)/(1-alpha) ;
WZERO      =  YZERO - LZERO ;
RZERO      =  ((R0+delta)/R0)*YZERO ;
DRZERO     =  R0*RZERO ;
IZERO      =  (1/omega_I)*(YZERO - omega_C*CZERO - omega_G*GZERO) ;
%
%  Long-run resuts
%
TCAP       =  10e6 ;
ZINF       =  exp((xi_Z-1)*TCAP)*eta_0 ;
GINF       =  0;
KINF       =  ( Delta22*gamma_K - Delta12*gamma_C)*ZINF/(lambda_2*lambda_1) ;
CINF       =  (-Delta21*gamma_K + Delta11*gamma_C)*ZINF/(lambda_2*lambda_1) ;
YINF       =  phi*(ZINF+alpha*KINF) - (phi-1)*CINF ;
LINF       =  ( YINF - ZINF - alpha*KINF ) / (1-alpha) ;
WINF       =  YINF - LINF ;
RINF       =  ((R0+delta)/R0)*(YINF - KINF ) ;
IINF       =  (1/omega_I)*(YINF - omega_C*CINF - omega_G*GINF);
%
%  Impulse-response functions
%
time=linspace(0,100,101);
Aterm= 1 - (1-lambda_1).^time ;
Tterm= ( xi_Z.^time - (1-lambda_1).^time ) ./ (xi_Z - (1-lambda_1));
%mesh(X,Y,Tterm)
%mesh(X,Y,Aterm)
%

figure(1)
Zpath = eta_0 * xi_Z.^time;
plot(time,Zpath,'b.','MarkerSize',20,'LineWidth',1.5)
xlim ([0 40])
ylim ([0 1])
box off
xlabel('$t$','interpreter','latex','FontSize',16)
ylabel('$\tilde{Z}_{t} \times 100\%$','interpreter','latex','FontSize',16)
title('(a) Productivity','FontSize',16)
print('-deps2','-loose','F18_03_Z')
pause(1)


figure(2)
gamma_C  = phi*zeta*(xi_Z - ystar*(1-phi*alpha));
Kpath = (eta_0*(gamma_K*(Delta22+1-xi_Z) - gamma_C*Delta12)./(lambda_2+1-xi_Z)) .*Tterm ;
plot(time,Kpath,'b.','MarkerSize',20,'LineWidth',1.5)
xlim ([0 40])
box off
xlabel('$t$','interpreter','latex','FontSize',16)
ylabel('$\tilde{K}_{t} \times 100\%$','interpreter','latex','FontSize',16)
title('(b) Capital stock','FontSize',16)
print('-deps2','-loose','F18_03_K')
pause(1)

figure(3)
Cnul = - eta_0*(gamma_C*Delta12 + gamma_K*(lambda_2-Delta22)) ...
   ./ ((lambda_2 + (1-xi_Z) )*Delta12);
Cpath  = (1-Aterm).*Cnul + ...
  eta_0*((gamma_C.*(Delta11+1-xi_Z) - gamma_K*Delta21)./(lambda_2+1-xi_Z)).*Tterm;
plot(time,Cpath,'b.','MarkerSize',20,'LineWidth',1.5)
xlim ([0 40])
box off
xlabel('$t$','interpreter','latex','FontSize',16)
ylabel('$\tilde{C}_{t} \times 100\%$','interpreter','latex','FontSize',16)
title('(c) Consumption','FontSize',16)
print('-deps2','-loose','F18_03_C')
pause(1)

figure(4)
Ypath  = phi*(Zpath +alpha*Kpath) - (phi-1)*Cpath ;
plot(time,Ypath,'b.','MarkerSize',20,'LineWidth',1.5)
xlim ([0 40])
box off
xlabel('$t$','interpreter','latex','FontSize',16)
ylabel('$\tilde{Y}_{t} \times 100\%$','interpreter','latex','FontSize',16)
title('(e) Output','FontSize',16)
print('-deps2','-loose','F18_03_Y')
pause(1)

figure(5)
Lpath  = ( Ypath - Zpath - alpha*Kpath )/(1-alpha);
plot(time,Lpath,'b.','MarkerSize',20,'LineWidth',1.5)
xlim ([0 40])
box off
xlabel('$t$','interpreter','latex','FontSize',16)
ylabel('$\tilde{L}_{t} \times 100\%$','interpreter','latex','FontSize',16)
title('(f) Employment','FontSize',16)
print('-deps2','-loose','F18_03_L')
pause(1)

figure(6)
Wpath  = Ypath - Lpath;
plot(time,Wpath,'b.','MarkerSize',20,'LineWidth',1.5)
xlim ([0 40])
box off
xlabel('$t$','interpreter','latex','FontSize',16)
ylabel('$\tilde{w}_{t} \times 100\%$','interpreter','latex','FontSize',16)
title('(g) Wage','FontSize',16)
print('-deps2','-loose','F18_03_W')
pause(1)

figure(7)
Rpath  =((R0+delta)/R0)*(Ypath - Kpath);
plot(time,Rpath,'b.','MarkerSize',20,'LineWidth',1.5)
xlim ([0 40])
box off
xlabel('$t$','interpreter','latex','FontSize',16)
ylabel('$\tilde{r}_{t} \times 100$','interpreter','latex','FontSize',16)
title('(h) Interest rate','FontSize',16)
print('-deps2','-loose','F18_03_R')
pause(1)

figure(8)
Ipath    = (1/omega_I)*(Ypath - omega_C*Cpath);
plot(time,Ipath,'b.','MarkerSize',20,'LineWidth',1.5)
xlim ([0 40])
box off
xlabel('$t$','interpreter','latex','FontSize',16)
ylabel('$\tilde{I}_{t} \times 100\%$','interpreter','latex','FontSize',16)
title('(d) Investment','FontSize',16)
print('-deps2','-loose','F18_03_I')
pause(1)







