%  
% Simulation program for Chapter 18 
%
% Figure 18.4: Impulse-response functions
%    Transitory technology shock
%
% B.J. Heijdra 
%
% 20 March 2017
% 
% The controlled parameters are:
%
xi_Z      = 0;
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)
%
%  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
%
%  Impulse-response functions
%
%  Note that X plays the role of xi_Z
%
tlin     =linspace(0,100,101);
xlin     =linspace(0.5,1,51);
[X,Y]    =meshgrid(xlin,tlin);
Aterm    = 1 - (1-lambda_1).^Y ;
Tterm    = ( X.^Y - (1-lambda_1).^Y ) ./ (X - (1-lambda_1));
gamma_K  = phi*ystar;
gamma_C  = phi*zeta*(X - ystar*(1-phi*alpha));
%mesh(X,Y,Tterm)
%mesh(X,Y,Aterm)


figure(1)
Zpath = eta_0*X.^Y;
mesh(X,Y,Zpath)
view([-150,40])
xlabel('$\xi_Z$','interpreter','latex','FontSize',16)
ylabel('$t$','interpreter','latex','FontSize',16)
zlabel('$\tilde{Z}_{t} \times 100\%$','interpreter','latex','FontSize',16)
title('(a) Productivity','FontSize',16)
print('-deps2','-loose','F18_04_Z')
pause(1)


figure(2)

Kpath = (eta_0*(gamma_K*(Delta22+1-X) - gamma_C*Delta12)./(lambda_2+1-X)) .*Tterm ;
mesh(X,Y,Kpath)
%view([-150,40])
xlabel('$\xi_Z$','interpreter','latex','FontSize',16)
ylabel('$t$','interpreter','latex','FontSize',16)
zlabel('$\tilde{K}_{t} \times 100\%$','interpreter','latex','FontSize',16)
title('(b) Capital stock','FontSize',16)
print('-deps2','-loose','F18_04_K')
pause(1)

figure(3)
Cnul = - eta_0*(gamma_C*Delta12 + gamma_K*(lambda_2-Delta22)) ...
   ./ ((lambda_2 + (1-X) )*Delta12);
Cpath  = (1-Aterm).*Cnul + ...
  eta_0*((gamma_C.*(Delta11+1-X) - gamma_K*Delta21)./(lambda_2+1-X)).*Tterm;
mesh(X,Y,Cpath)
%view([-150,40])
xlabel('$\xi_Z$','interpreter','latex','FontSize',16)
ylabel('$t$','interpreter','latex','FontSize',16)
zlabel('$\tilde{C}_{t} \times 100\%$','interpreter','latex','FontSize',16)
title('(c) Consumption','FontSize',16)
print('-deps2','-loose','F18_04_C')
pause(1)

figure(4)
Ypath  = phi*(Zpath +alpha*Kpath) - (phi-1)*Cpath ;
mesh(X,Y,Ypath)
view([-150,40])
xlabel('$\xi_Z$','interpreter','latex','FontSize',16)
ylabel('$t$','interpreter','latex','FontSize',16)
ylabel('$t$','interpreter','latex','FontSize',16)
zlabel('$\tilde{Y}_{t} \times 100\%$','interpreter','latex','FontSize',16)
title('(e) Output','FontSize',16)
print('-deps2','-loose','F18_04_Y')
pause(1)

figure(5)
Lpath  = ( Ypath - Zpath - alpha*Kpath )/(1-alpha);
mesh(X,Y,Lpath)
view([150,40])
xlabel('$\xi_Z$','interpreter','latex','FontSize',16)
ylabel('$t$','interpreter','latex','FontSize',16)
zlabel('$\tilde{L}_{t} \times 100\%$','interpreter','latex','FontSize',16)
title('(f) Employment','FontSize',16)
print('-deps2','-loose','F18_04_L')
pause(1)

figure(6)
Wpath  = Ypath - Lpath;
mesh(X,Y,Wpath)
%view([-150,40])
xlabel('$\xi_Z$','interpreter','latex','FontSize',16)
ylabel('$t$','interpreter','latex','FontSize',16)
zlabel('$\tilde{w}_{t} \times 100\%$','interpreter','latex','FontSize',16)
title('(g) Wage','FontSize',16)
print('-deps2','-loose','F18_04_W')
pause(1)

figure(7)
Rpath  =((R0+delta)/R0)*(Ypath - Kpath);
mesh(X,Y,Rpath)
view([-150,40])
xlabel('$\xi_Z$','interpreter','latex','FontSize',16)
ylabel('$t$','interpreter','latex','FontSize',16)
zlabel('$\tilde{r}_{t} \times 100$','interpreter','latex','FontSize',16)
title('(h) Interest rate','FontSize',16)
print('-deps2','-loose','F18_04_R')
pause(1)

figure(8)
Ipath    = (1/omega_I)*(Ypath - omega_C*Cpath);
mesh(X,Y,Ipath)
view([150,40])
xlabel('$\xi_Z$','interpreter','latex','FontSize',16)
ylabel('$t$','interpreter','latex','FontSize',16)
zlabel('$\tilde{I}_{t} \times 100\%$','interpreter','latex','FontSize',16)
title('(d) Investment','FontSize',16)
print('-deps2','-loose','F18_04_I')
pause(1)







