clear all

%%%%% Problem Set #2, International Monetary Economics

%%%%% Declare parameter values

beta  = 0.99;
theta = 0.33;
delta = 0.04;
rho   = 0.95;

%%%%% Solve for non-stochastic steady state

kn = ((theta*beta)/(1-beta+delta*beta))^(1/(1-theta)); %% k/n ratio %

cl = (1-theta)*(kn^theta);                             %% c/(1-n) ratio %

cn = (kn^theta)-(delta*kn);                            %% c/n ratio

nbar  = cl/(cl+cn);
cbar  = cn*nbar;
kbar  = kn*nbar;

%%%%% Other variables we need for coefficients

ibar  = delta*kbar;
ybar  = cbar+ibar;

Fk    = theta*(kn^(theta-1));
Rbar  = 1+Fk-delta; %% check this equals 1/beta!!

%%%%% Coefficients of log-linear model

AA = [0;0;0;0;-kbar;0];
BB = [theta;1;-(1-theta)*Fk;theta*ybar;(1-delta)*kbar;0];

CC = [-1,-(theta+(nbar/(1-nbar))),0,0,0,0;
      0,0,-1,1,0,0;
      0,(1-theta)*Fk,(delta-Fk)*delta,0,-Rbar,0;
      -cbar,(1-theta)*ybar,0,-ibar,0,0;
      0,0,0,ibar,0,0;
      1,0,-delta,0,0,1];

DD = [1;0;Fk;1;0;0];

FF = [0];
GG = [0];
HH = [0];
JJ = [0,0,0,0,1,1];
KK = [0,0,0,0,0,-1];
LL = [0];
MM = [0];

NN = rho;

%%%%% Call Uhli's toolkit

[l_equ,m_states] = size(AA);
[l_equ,n_endog ] = size(CC);
[l_equ,k_exog  ] = size(DD);

message = '                                                                       ';
warnings = [];
OPTIONS
SOLVE


%%%%% recover some plots

T  = 100;

time = [0;cumsum(ones(T,1))];

ZZ = zeros(1,T+1);
XX = zeros(1,T+1);
YY = zeros(6,T+1);

e0 = 1;

for t = 1:T,
    ZZ(:,1)   = e0;
    XX(:,1)   = 0;
    YY(:,1)   = zeros(6,1);
    
    ZZ(:,t+1) = NN*ZZ(:,t);
    
    XX(:,t+1) = PP*XX(:,t)+QQ*ZZ(:,t+1);
    YY(:,t+1) = RR*XX(:,t)+SS*ZZ(:,t+1);
end

figure(1)
plot(time,ZZ,time,YY(1,:),time,YY(2,:),time,YY(4,:),time,YY(5,:))
legend('productivity','consumption','employment','investment','returns')