%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%        Matlab code for solving a stochastic growth model in HW # 6.		
%            Jinhui Bai (jinhui.bai@yale.edu), Feb 17, 2004
% 
% Discription of the model: CES utility, AR(1) shocks
% Algorithm: Linear Quadratic Method from Euler Equation
% Purpose: Calculate Impulse-Response to techonology shocks
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                         Outline of the Code
% Step 1: Calculate Steady State Value		
% Step 2: Form dynamic system and find impulse-reponse function
% Step 3: Graph the result
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
diary brock_hw6.out;
time_start = cputime;
disp('solving a stochastic growth model in HW # 6');
disp('');

%  set parameter values

ALPHA  = 0.36;              % capital share
BETA   = 0.99;              % subjective discount factor 
DELTA  = 0.025;             % capital depreciation rate
RHO    = 0.95;              % persistence of shocks
SIGMA  = 0.007;             % variance of technology shocks
GAMMA  = 1;                 % inverse of Intertemporal Elasticity of Substitution
% or GAMMA = 2;

% Calculate Steady State Value

KSS = ((1/BETA - 1 + DELTA)/ALPHA)^(1/(ALPHA - 1)); % steady-state capital stock
YSS = KSS^ALPHA;                                    % steady-state output
ISS = DELTA * KSS;                                  % steady-state investment
CSS = YSS - ISS;                                    % steady-state consumption

% Form Dynamic System and Find Impulse-Response Function

A1 = -GAMMA * (KSS)^ALPHA;
A2 = -GAMMA/BETA; 
A3 = GAMMA;
B1 = (KSS^ALPHA) * (1- GAMMA - BETA + BETA*(1 - ALPHA)*DELTA);
B2 = -GAMMA/BETA + (1-1/ALPHA) * (1/BETA -1 + DELTA) * (1 - BETA);
B3 = GAMMA;

A = roots([B3 (B2-A3) -A2]);                  % Solve for the quadratic equation
A
CK = A(abs(A) < 1);                             % Get |k coefficient|< 1
CZ = (A1 - B1*RHO)/(B2 + B3*(CK + RHO) - A3);   % coefficient for z(t)
N  = 101;                   % number of periods to simulate

z = RHO(ones(N,1));        % initialize shock matrix
z(1) = SIGMA;              % initial shock
z = cumprod(z);            % the shock matrix in impulse-response function

k = zeros(N,1);            % initialize capital matrix
k = filter([0 CZ],[1 -CK],z); % the capital matrix in reponse of shocks
kdev = k/KSS;               % percentage deviation for k

y = ALPHA * (KSS^(ALPHA - 1)) * k ; % the output matrix in reponse of shocks
ydev = y/YSS;                       % percentage deviation for y

i = (CK - 1 + DELTA)*k + CZ*z;      % the invest matrix in reponse of shocks
idev = i/ISS;                       % percentage deviation for i

c = y - i;                          % the consumption matrix in reponse of shocks
cdev = c/CSS;                       % percentage deviation for c


disp('Time Spent in Computation')
disp('');
time_spent = cputime - time_start;
time_spent
diary off;

% Display and Plot Result

disp('  Absolute Deviation of different variables');
disp('  technology capital output consumption investment');
[z k y c i]

disp('  Parameter Values in recursive law of motion');
disp('k coefficient z coefficient');
[CK CZ]

period = (1:1:N);
figure(1)
plot(period,kdev,'-',period,ydev,':');
title('Impulse Response: Capital Versus Output');
print output.ps

figure(2)
plot(period,kdev,'-',period,cdev,':');
title('Impulse Response: Capital Versus Consumption');
print consumption.ps

figure(3)
plot(period,kdev,'-',period,idev,':');
title('Impulse Response: Capital Versus Investment');
print invest.ps

