% program6.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 Otsu (2007)
% The 4 wedges are taken as given

% Otsu'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 Chen-chiao Kuo's note "Derivation and Log-linearization of Otsu (2007)'s Small Open Economy Model"
% Otsu, Keisuke (2007), IMES Discussion Paper No. 2007-E-16
% Use Paul Klein's code to solve the model


% variables order

%1  k
%2  d
%3  z
%4  taul
%5  taux
%6  taud
%7  c
%8  x
%9  y
%10 l

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

disp('clear');
clear;

% Parameters Values (page 26 of Otsu (2007), Table 1, Singapore)
% structural parameters
betta=0.911; % discount factor
theta=1/3;   % capital share
delta=0.018; % depreciation rate
chi=1.089;   % level parameter of GHH preference
nu=1.409;    % curvature parameter of GHH prefernce
Gamma=1.052; % Growth trend = (1+gamma)(1+r)
psi=0.00001;

% VAR parameters
p11=0.9;
p12=0.0;
p13=0.0;
p14=0.0;
p21=0.0;
p22=0.9;
p23=0.0;
p24=0.0;
p31=0.0;
p32=0.0;
p33=0.9;
p34=0.0;
p41=0.0;
p42=0.0;
p43=0.0;
p44=0.9;

% steady state
z=1;
taul=1;
taux=1;
taud=1;

% steady state
R=Gamma/(betta*taud);
kl=((1/(z*theta))*((taux*Gamma)/betta-taux*(1-delta)))^(1/(theta-1));
l=((1-theta)*(z/(chi*nu))*(kl^theta)*(1/taul))^(1/nu);
k=kl*l;
y=z*(kl^theta)*l;
tby=-0.03;                               % key, calibrated from observed data
tb=tby*y;
d=tb/(1-betta);
c=z*(kl^theta)*l-(Gamma-(1-delta))*(kl)*l+d*(betta-1);
x=(Gamma-(1-delta))*k;

% checked!
disp('y');  disp(y);
disp('x');  disp(x);
disp('l');  disp(l);
disp('c');  disp(c);
disp('tb'); disp(tb);

% eqs
nvar=10;
nk=6;

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

A(1,2)=(Gamma*d)/(R*taud*y);
B(1,9)=-1;
B(1,7)=c/y;
B(1,8)=x/y;
B(1,2)=d/y;
B(1,6)=(Gamma*d)/(R*taud*y);

B(2,9)=-1;
B(2,3)=1;
B(2,1)=theta;
B(2,10)=1-theta;

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

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

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

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

B(7,4)=-1;
B(7,3)=1;
B(7,1)=theta;
B(7,10)=(1-theta-nu);

A(8,9)=-1*(betta*theta*y)/(Gamma*taux*k);
A(8,1)=(betta*theta*y)/(Gamma*taux*k);
A(8,7)=c/(c-chi*(l^nu));
A(8,10)=-1*(chi*nu*(l^nu))/(c-chi*(l^nu));
A(8,5)=-1*(betta*(1-delta))/Gamma;
B(8,5)=-1;
B(8,7)=c/(c-chi*(l^nu));
B(8,10)=-1*(chi*nu*(l^nu))/(c-chi*(l^nu));

A(9,7)=-1*c/((c-chi*(l^nu)));
A(9,10)=(chi*(l^nu)*nu)/(c-chi*(l^nu));
A(9,2)=(psi*d)/(betta);
B(9,6)=-1;
B(9,7)=-1*c/((c-chi*(l^nu)));
B(9,10)=(chi*(l^nu)*nu)/(c-chi*(l^nu));

A(10,1)=Gamma;
B(10,8)=x/k;
B(10,1)=1-delta;

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

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

