% program4.m
%
% By Dr. Tai-kuang Ho
% Department of Quantitative Finance
% National Tsing Hua University, Taiwan
% ****** Benutzung auf eigene Gefahr! ******
%
% I use Paul Klein's code to solve the RBC example of Uhlig's paper
% Uhlig, A toolkit for analyzing nonlinear dynamic stochastic models easily
% Hansen's Real Business Cycle Model


% variables order

%1  k
%2  z
%3  i
%4  y
%5  n
%6  c
%7  r

disp('clear');
clear;

% parameters values
Nbar=1/3;
betta=1/1.01;
eta=1.0;
delta=0.025;
rho=0.36;
Zbar=1;
phi=0.95;

% steady state

Rbar=1/betta;
YbarKbar=(Rbar-1+delta)/rho;
KbarNbar=(YbarKbar/Zbar)^(1/(rho-1));
Kbar=KbarNbar*Nbar;
Ibar=delta*Kbar;
Ybar=YbarKbar*Kbar;
Cbar=Ybar-delta*Kbar;
Abar=(Cbar^(-1*eta))*(1-rho)*(Ybar/Nbar);

Ybarcheck=Zbar*(Kbar^(rho))*(Nbar^(1-rho));

disp('Cbar='); disp(Cbar);
disp('Ibar='); disp(Ibar);
disp('Ybar='); disp(Ybar);
disp('Kbar='); disp(Kbar);
disp('Nbar='); disp(Nbar);
disp('Rbar='); disp(Rbar);
disp('Zbar='); disp(Zbar);
disp('Abar='); disp(Abar);
disp('Ybarcheck='); disp(Ybarcheck);


% eqs
nvar=7;
nk=2;

A=zeros(nvar,nvar);
B=zeros(nvar,nvar);

B(1,6)=-Cbar;
B(1,3)=-Ibar;
B(1,4)=Ybar;

A(2,2)=1;
B(2,2)=phi;

A(3,1)=Kbar;
B(3,3)=Ibar;
B(3,1)=(1-delta)*Kbar;

B(4,4)=-1;
B(4,2)=1;
B(4,1)=rho;
B(4,5)=(1-rho);

B(5,6)=-eta;
B(5,4)=1;
B(5,5)=-1;

A(6,6)=eta;
A(6,7)=-1;
B(6,6)=eta;

B(7,7)=-Rbar;
B(7,4)=rho*(Ybar/Kbar);
B(7,1)=-rho*(Ybar/Kbar);

[f,p]=solab(A,B,nk);

disp('policy     function'); disp(f);
disp('transition function'); disp(p);


% impulse responses
% use ir.m to compute impulse response function (IRF)
gx=f;
hx=p;

x0=zeros(size(hx,1),1);
x0(end)=1;
IRF=ir(gx,hx,x0,100);

t=(0:size(IRF,1)-1)';
i=IRF(:,1);
y=IRF(:,2);
n=IRF(:,3);
c=IRF(:,4);
r=IRF(:,5);
k=IRF(:,6);
z=IRF(:,7);

subplot(4,2,1)
plot(t,i)
title('response of investment')

subplot(4,2,2)
plot(t,y)
title('response of output')

subplot(4,2,3)
plot(t,n)
title('response of labour hour')

subplot(4,2,4)
plot(t,c)
title('response of consumption')

subplot(4,2,5)
plot(t,r)
title('response of interest rate')

subplot(4,2,6)
plot(t,k)
title('response of capital')

subplot(4,2,7)
plot(t,z)
title('productivity shocks')

shg


%{
This program is checked.
One important notice is that the IRF produced by ir.m of SGU
put the control variables at the first ordering, followed by
the state variables (state and exogenou shocks).
So it is important not so confuse the IRF of one for another.
%}