clear all

%%%%% Problem Set #3

%%%%% Declare parameter values

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

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

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

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

cn = (kn^alpha)-(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;

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



%%%%% Coefficients of log-linear model

AA = [0;-kbar];
BB = [alpha;alpha*ybar+(1-delta)*kbar];
CC = [-1,-(alpha+(nbar/(1-nbar))); -cbar, (1-alpha)*ybar];
DD = [1;1];

FF = [0];
GG = [-beta*rbar*(1-alpha)];
HH = [0];
JJ = [-1,beta*rbar*(1-alpha)];
KK = [1,0];
LL = [beta*rbar];
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(2,T+1);

e0 = 1;

for t = 1:T,
    ZZ(:,1)   = e0;
    XX(:,1)   = 0;
    YY(:,1)   = zeros(2,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

productivity = ZZ;
consumption  = YY(1,:);
employment   = YY(2,:);
capital      = XX;

output       = (productivity+alpha*ybar*capital+(1-alpha)*ybar*employment)/ybar;

investment   = (ybar*output-cbar*consumption)/ibar;

figure(1)
plot(time,productivity,time,consumption,time,employment,time,investment)
legend('productivity','consumption','employment','investment')
    
%%%%% SIMULATIONS

T = 1000;

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

%%%%% Begin with draw of innovations to technology

EE = randn(T,1);

%%%%% Consrtruct paths for ZZ and XX,YY

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

%%%%% Drop first D observations

D  = 500 %% D<T

ZZ = ZZ(:,D+1:T+1);
XX = XX(:,D+1:T+1);
YY = YY(:,D+1:T+1);

productivity = ZZ;
consumption  = YY(1,:);
employment   = YY(2,:);
capital      = XX;

output       = (productivity+alpha*ybar*capital+(1-alpha)*ybar*employment)/ybar;

investment   = (ybar*output-cbar*consumption)/ibar;

time = cumsum(ones((T-D)+1,1));

figure(2)
plot(time,productivity,time,consumption,time,employment,time,investment)
legend('productivity','consumption','employment','investment')

std_z = std(productivity)
std_c = std(consumption)
std_n = std(employment)
std_i = std(investment)
std_y = std(output)