%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Jeff Thurk %
% Example of RBC Model from Corbae Notes %
% Macro II %
% Corbae %
% 2/10/2006 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following program computes decision rules, an impulse response
% function for capital, and simulated moments for the RBC model presented
% in class. The results differ from those calculated in class since this
% program calculates the coefficients using absolute deviations from SS
% rather than percentage deviations.
% Let me know if you have any questions, or find any errors.
clear all; % Clears all variables
clc % Clears the command window
% Parameters
gamma = .016;
theta = .36;
delta = .06;
rho = .95;
sigma = .007;
beta = .9691;
% Step One- Find the FOC (done by hand)
% Step Two- Find the SS
guess = [1,1,1]; % steady-state values for c, k, z
options = []; % placeholder for the fsolve command
steady = fsolve(@ss,guess,options,beta,theta,delta,rho);
% Name the solutions from the fsolve command
cbar = steady(1); kbar = steady(2); zbar = steady(3);
% Step Three- Linearize the FOC
syms k0 k1 k2 c0 c1 c2 z0 z1 % These variables are now symbolic (can take derivatives of)
% In order to save space and complication in the FOC, define consumption as
% a function of the other symbolic variables (Could also do for Y and I)
c0 = z0*k0^theta + (1-delta)*k0 - k1;
c1 = z1*k1^theta + (1-delta)*k1 - k2;
% The lone FOC set equal to zero. When more than 1, than each is a row
FOC = [c1/c0-beta*(theta*z1*k1^(theta-1)+(1-delta))];
% The V matrix will tell matlab what order to take the partial derivatives
V = [k2 k1 k0 z1 z0];
% Take partial derivatives of the FOC (We're linearizing the FOC)
J = jacobian(FOC,V);
% Substitute in the Steady-States (i.e., evaluate Taylor expansion at SS)
k0=kbar; k1=kbar; k2=kbar; c0=cbar; c1=cbar; c2=cbar; z0=zbar; z1=zbar;
J = subs(J);
% Transform the zt [from (2.2)] part into a first-order difference equation
% Note that V = [k2 k1 k0 z1 z0]
% A0 A1 A2 B1 B2
% You can follow christiano's paper, or just think about matching
% (k(t+2)-kbar) with the FOC partial derivative with respect to k(t+2). If
% we had another component in the Zt vector (say labor (n)) then A0 would
% be a 2x2 matrix covering the variables with t+2 subscripts.
alpha0 = J(:,1);
alpha1 = J(:,2);
alpha2 = J(:,3);
beta0 = J(:,4);
beta1 = J(:,5);
% Create the a and b martrices (Just a matrix algebra transformation of
% Christiano's equation 2.2)
a = [alpha0 0;0 1];
b = [alpha1 alpha2;-1 0];
% Here a is invertible. Otherwise you would have to parse out the a,b
% matrices using a QZ decomposition.
[E,D]=eig(-inv(a)*b);
% The eigenvalue in the bottom right is nonexplosive therefore
D1 = 1; % In effect, negate the explosive root by setting equal to one
D2 = -D(2,2);
A = -inv(D1)*D2;
% Display the results in the command window
disp('*********************')
disp(['A matrix ' num2str(A)])
disp('*********************')
B = inv(alpha0*rho + A*alpha0 + alpha1)*(-beta0*rho-beta1);
disp('*********************')
disp(['B matrix ' num2str(B)])
disp('*********************')
% We now have the decision rules so we can compute a simulated economy
% Note- Given the decision Rules for capital, we could calculate series of
% consumption, output, and investment.
% IMPULSE RESPONSE FUNCTION
N = 100; % Length of IRF and Simulation
for t = 2:N;
% Create the shock matrix given an initial epsilon shock
zhat(1,1) = 1;
zhat(t,1) = rho*zhat(t-1,1);
end
for t = 2:N
% Use the decision rule and the TFP vector to create the IRF
IRF(:,1) = B*zhat(1,1);
IRF(:,t) = A*IRF(:,t-1,1) + B*zhat(t,1);
end
% SIMULATION
S = 100; % Number of simulations
T = 200;
for s=1:S
eps = randn(T)*sigma;
for t = 2:T;
zhat(1,1) = 0;
zhat(t,1) = rho*zhat(t-1,1) + eps(t);
end
for t = 2:T;
SIM(1,1,s) = B*zhat(1);
SIM(1,t,s) = A*SIM(1,t-1) + B*zhat(t);
end
end
% Plot the Results
figure (1)
subplot(2,1,1)
plot(1:N,IRF,'-.r')
xlabel('Time Period')
ylabel('Capital Deviation from SS')
title('IMPULSE RESPONSE FUNCTION')
subplot(2,1,2)
plot(1:T,SIM(1,:,1),'--k')
xlabel('Time Period')
ylabel('Capital Deviation from SS')
title('SIMULATED ECONOMY')
% CALCULATE MOMENTS
% Calculate the std of each simulated time series (along 2nd dimension)
STD = std(SIM,0,2);
% Average the standard deviations across the simulations
STD = mean(STD,3);
disp('*********************')
disp(' Model Moments ')
disp('*********************')
disp(['Var K ' num2str(STD^2)])
disp('*********************')