% program3.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
% The Neoclassical Growth Model


% variables order

%1  k
%2  z
%3  c
%4  r

clear;

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

% steady state
Rbar=1/betta;
Kbar=(rho*Zbar/(Rbar-1+delta))^(1/(1-rho));
Ybar=Zbar*(Kbar^(rho));
Cbar=Ybar-delta*Kbar;

% eqs
nvar=4;
nk=2;

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

A(1,1)=Kbar/Cbar;
B(1,1)=Kbar/(betta*Cbar);
B(1,2)=Ybar/Cbar;
B(1,3)=-1;

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

B(3,4)=-1;
B(3,2)=(1-betta*(1-delta));
B(3,1)=-(1-betta*(1-delta))*(1-rho);

A(4,3)=eta;
A(4,4)=-1;
B(4,3)=eta;

[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,32);

t=(0:size(IRF,1)-1)';
c=IRF(:,1);
r=IRF(:,2);
k=IRF(:,3);
z=IRF(:,4);

% plot(t,c,t,r,t,k,t,z);

subplot(2,2,1)
plot(t,c)
title('response of consumption to technological shock')

subplot(2,2,2)
plot(t,r)
title('response of interest rate to technological shock')

subplot(2,2,3)
plot(t,k)
title('response of capital to technological shock')

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

shg

%{
This program is checked, it can successfully replicate the results of Uhlig' paper.
Before using this program, please read my lecture note:
"Solving linear rational expectation models".

Some important points are:
1. Put the state variables on the first ordering
2. Add the exogenous shocks, which is treated as state variables
3. The exogenous shocks should be expressed in innovation of period t+1 when entering the equation
4. The number nk is the number of state and exogenous variables
5. (Preferred) When entering equations, the equations for exogenous shocks should be at the same position as the variables ordering
6. Actually, the entering of variables are not relevant
%}


%{
* == ordering does not matter == *
% eqs
nvar=4;
nk=2;

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

A(1,1)=Kbar/Cbar;
B(1,1)=Kbar/(betta*Cbar);
B(1,2)=Ybar/Cbar;
B(1,3)=-1;

B(2,4)=-1;
B(2,2)=(1-betta*(1-delta));
B(2,1)=-(1-betta*(1-delta))*(1-rho);

A(3,3)=eta;
A(3,4)=-1;
B(3,3)=eta;

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

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

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