% program5.m
%
% By Dr. Tai-kuang Ho
% Department of Quantitative Finance
% National Tsing Hua University, Taiwan
% ****** Benutzung auf eigene Gefahr! ******
%
% This program solves the log-linearized prototype model of Chari et al. (2007)
% The 4 wedges are taken as given

% Chari et al.'s model is for purpose of business cycle accounting.
% This program simply serves to solve the decision rule of the endogenous variables.
% It makes no sense to compute the impulse reponses to the wedges.

% See Chi-hong Chen's note "Derivation and Log-linearization of Chari, Kehoe and McGrattan(2007)'s Closed Economy Model"
% Chari, Kehoe, and McGrattan (2007), Econometrica
% Use Paul Klein's code to solve the model


% variables order

%1 k
%2 Z
%3 taul
%4 taux
%5 g
%6 c
%7 x
%8 y
%9 l

% observables
%1 y
%2 x
%3 l
%4 g


% Parameters Values
alpha=1/3;
psi  =2.24;
delta=0.0464;
gama =0.016;
gaman=0.015;
betta=0.9722;

% VAR parameters
p11=0.9;
p12=0;
p13=0;
p14=0;

p21=0;
p22=0.9;
p23=0;
p24=0;

p31=0;
p32=0;
p33=0.9;
p34=0;

p41=0;
p42=0;
p43=0;
p44=0.9;


Zbar=1;
taulbar=0;
tauxbar=0;
gbar=1;
% steady state
% as input: Zbar,taulbar,tauxbar,gbar

klbar=((1/(alpha*Zbar^(1-alpha)))*((1+tauxbar)*(1+gama)/betta-(1-delta)*(1+tauxbar)))^(1/(alpha-1));
lbar=(gbar+(1/psi)*(1-taulbar)*Zbar^(1-alpha)*(klbar^alpha))...
    /(Zbar^(1-alpha)*(klbar^alpha)-((1+gaman)*(1+gama)-(1-delta))*klbar+(1/psi)*(1-taulbar)*Zbar^(1-alpha)*(klbar^alpha));
kbar=klbar*lbar;
cbar=((1-lbar)/psi)*(1-taulbar)*Zbar^(1-alpha)*(1-alpha)*(klbar^alpha);
ybar=Zbar^(1-alpha)*(klbar^alpha)*lbar;
xbar=((1+gaman)*(1+gama)-(1-delta))*kbar;


% eqs
nvar=9;
nk=5;

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

B(1,6)=-1*(cbar/ybar);
B(1,7)=-1*(xbar/ybar);
B(1,5)=-1*(gbar/ybar);
B(1,8)=1;

A(2,2)=1;
B(2,2)=p11;
B(2,3)=p12;
B(2,4)=p13;
B(2,5)=p14;

A(3,3)=1;
B(3,2)=p21;
B(3,3)=p22;
B(3,4)=p23;
B(3,5)=p24;

A(4,4)=1;
B(4,2)=p31;
B(4,3)=p32;
B(4,4)=p33;
B(4,5)=p34;

A(5,5)=1;
B(5,2)=p41;
B(5,3)=p42;
B(5,4)=p43;
B(5,5)=p44;

B(6,8)=-1;
B(6,2)=(1-alpha);
B(6,1)=alpha;
B(6,9)=(1-alpha);

B(7,6)=-1;
B(7,2)=(1-alpha);
B(7,1)=alpha;
B(7,9)=-1*(alpha+lbar/(1-lbar));
B(7,3)=-1*(1/(1-taulbar));

A(8,6)=(1+tauxbar)*(1+gama)/betta;
A(8,8)=-1*alpha*(ybar/kbar);
A(8,1)=alpha*(ybar/kbar);
A(8,4)=-1*(1-delta);
B(8,4)=-1*(1+gama)/betta;
B(8,6)=(1+tauxbar)*(1+gama)/betta;

A(9,1)=(1+gaman)*(1+gama);
B(9,1)=(1-delta);
B(9,7)=xbar/kbar;


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

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

gx=f;
hx=p;





