%%%%% Scrap of Matlab code for stochastic growth model

clear all
close all

%%%%% parameters  

alpha = 0.3;          %% capital share in production function
beta  = 1/1.05;       %% time discount factor
delta = 0.05;         %% depreciation rate
sigma = 1;            %% coefficient of relative risk aversion
phi   = 0.95;         %% AR(1) coefficient for log productivity

sigeps= 0.015;        %% std dev of innovations to productivity


rho   = (1/beta)-1;   %% rate of time preferences 

%%%%% steady state

rbar  = rho + delta;                       %% steady state marginal product of capital

kbar  = (alpha/(rho+delta))^(1/(1-alpha)); %% steady state capital stock (per worker)
ybar  = kbar^alpha;                        %% steady state output        (per worker)
cbar  = ybar-delta*kbar;                   %% steady state investment

%%%%% We need to determine four coefficients
%
%      (psi_kk, psi_ck, psi_kz, psi_cz)
%

%%%%% As shown in the notes, psi_kk solves a quadratic of the form
%
%      psi_kk^2 - trace*psi_kk + determinant = 0
% 

trace       = 1+(1/beta)-(beta*rbar/sigma)*(alpha-1)*(cbar/kbar);
determinant = 1/beta;

lambdas = roots([1,-trace,determinant]); %% calculate the roots of a polynomial with these coefficients
lambda  = min(abs(lambdas));             %% choose stable root


psi_kk = lambda;
psi_ck = ((1/beta)-psi_kk)*(kbar/cbar);

%%%%% With these coefficients solved for, can recover psi_kz, psi_cz by
%%%%% solving two linear equations in two unknowns

temp1  = beta*rbar*phi + sigma*(1-phi)*(ybar/cbar);
temp2  = beta*rbar*(1-alpha) + sigma*psi_ck + sigma*(1-phi)*(kbar/cbar); 

psi_kz = temp1/temp2;
psi_cz = (ybar/cbar) - (kbar/cbar)*psi_kz;

display('Psi = [psi_kk,psi_kz;psi_ck,psi_cz]')

Psi    = [psi_kk,psi_kz;psi_ck,psi_cz]


%%%%% do a simple impulse response to a unit productivity shock

T = 150;          %% calculate for T periods after shock

time = (-1:1:T)'; %% steady state at t=-1, shock at t=0

shock    = zeros(size(time));
shock(2) = 0.01;     %% set the impulse 

zt   = zeros(size(time));
kt   = zeros(size(time));
ct   = zeros(size(time));
yt   = zeros(size(time));
it   = zeros(size(time));

for t=2:length(time),
    
    zt(t) = phi*zt(t-1) + shock(t);                %% productivity
    
    kt(t) = psi_kk*kt(t-1) + psi_kz*zt(t);         %% capital stock
    ct(t) = psi_ck*kt(t-1) + psi_cz*zt(t);         %% consumption
    
    it(t) = (kt(t) - (1-delta)*kt(t-1))/delta;     %% investment
    
    yt(t) = zt(t) + alpha*kt(t-1);                 %% output 
    
end    

figure(1)
plot(time,zt,'b-',time,ct,'g-',time,yt,'r-',time,kt,'--');
legend('productivity','consumption','output','capital','Location','Best')
xlabel('periods after shock')
ylabel('log-deviation from steady state') 
axis([-5, 150, 0 0.0125])


%%%%% do a simple simulation

randn('state',0); %% set seed of random number generator

T = 250;          %% calculate for T periods after shock

time = (-1:1:T)'; %% steady state at t=-1, shock at t=0

shock    = randn(size(time))*sigeps;

%shock(2) = 1;     %% set the impulse 

zt   = zeros(size(time));
kt   = zeros(size(time));
ct   = zeros(size(time));
yt   = zeros(size(time));
it   = zeros(size(time));

for t=2:length(time),
    
    zt(t) = phi*zt(t-1) + shock(t);                %% productivity
    
    kt(t) = psi_kk*kt(t-1) + psi_kz*zt(t);         %% capital stock
    ct(t) = psi_ck*kt(t-1) + psi_cz*zt(t);         %% consumption
    
    it(t) = (kt(t) - (1-delta)*kt(t-1))/delta;     %% investment
    
    yt(t) = zt(t) + alpha*kt(t-1);                 %% output 
    
end    

data = [zt,ct,yt];                 

figure(2)
plot(time,[zt,ct,yt])
plot(time,zt,'b-',time,ct,'g-',time,yt,'r-');
legend('productivity','consumption','output','Location','Best')
xlabel('time')
ylabel('log deviation from steady state')
axis([-5, T, 1.2*min(min(data)), 1.2*max(max(data))])


% with investment

data = [zt,ct,yt,it];                 

figure(3)
plot(time,[zt,ct,yt])
plot(time,zt,'b-',time,ct,'g-',time,yt,'r-',time,it,'c-');
legend('productivity','consumption','output','investment','Location','Best')
xlabel('time')
ylabel('log deviation from steady state')
axis([-5, T, 1.2*min(min(data)), 1.2*max(max(data))])


