-
Notifications
You must be signed in to change notification settings - Fork 13
/
gk.m
139 lines (108 loc) · 4.27 KB
/
gk.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
clear;
addpath('../functions/');
% Empirical example
% Based on Gertler & Karadi (AEJ:Macro 2015), DOI: 10.1257/mac.20130329
% MPM 2020-07-11
%% Settings
% Data files from Gertler & Karadi (2015)
file_var = '../data/gk/VAR_data.csv'; % VAR series
file_iv = '../data/gk/factor_data.csv'; % External instrument series
% Specification
vars = {'logcpi', 'logip', 'ff', 'ebp', 'ff4_tc'}; % Variables in VAR
p = 12; % Lag length
% Reporting
resp_var = 4; % Index of response variable
innov = 5; % Index of innovation
maxhorz = 48; % Maximum horizon to plot
alpha = 0.1; % Significance level
plot_ylim = [-4 4]; % y-axis limits
plot_xtick = [1 6:6:maxhorz]; % x-axis ticks
% Bootstrap
boot_num = 500; % # of repetitions
poolobj = gcp;
boot_workers = poolobj.NumWorkers; % # of parallel workers
rng(20200711, 'twister'); % Set random number seed
%% Load data
dat = innerjoin(readtable(file_var), readtable(file_iv));
Y = dat{:,vars}; % Select variables
Y = Y(~any(isnan(Y),2),:); % Remove missing
Y = detrend(Y,0); % Remove mean
horzs = 1:maxhorz;
%% Lag-augmented local projection
disp('Lag-augmented LP: bootstrapping...');
[irs_lp_la, ~, cis_dm_lp_la, cis_boot_lp_la] = ...
ir_estim(Y, p, horzs, ...
'estimator', 'lp', 'lag_aug', true, ...
'resp_var', resp_var, 'innov', innov, 'alpha', alpha, ...
'bootstrap', 'var', 'boot_num', boot_num, ...
'boot_workers', boot_workers, 'verbose', true);
figure;
subplot(2,1,1);
plot_band(horzs, irs_lp_la, ...
cis_dm_lp_la(1,:), cis_dm_lp_la(2,:), ...
'Lag-augmented LP, delta method interval', ...
plot_ylim, plot_xtick);
subplot(2,1,2);
plot_band(horzs, irs_lp_la, ...
cis_boot_lp_la(1,:,3), cis_boot_lp_la(2,:,3), ...
'Lag-augmented LP, percentile-t interval', ...
plot_ylim, plot_xtick);
drawnow;
%% Non-augmented local projection
disp('Non-augmented LP: bootstrapping...');
[irs_lp_na, ~, cis_dm_lp_na, cis_boot_lp_na] = ...
ir_estim(Y, p, horzs, ...
'estimator', 'lp', 'lag_aug', false, ...
'resp_var', resp_var, 'innov', innov, 'alpha', alpha, ...
'har', @(Y,bw) ewc(Y,bw), ... % HAR estimator
'har_bw', @(T) round(0.4*T.^(2/3)), ... % HAR bandwidth
'har_cv', @(bw) tinv(1-alpha/2,bw), ... % HAR CV
'bootstrap', 'var', 'boot_num', boot_num, ...
'boot_workers', boot_workers, 'verbose', true);
figure;
subplot(2,1,1);
plot_band(horzs, irs_lp_na, ...
cis_dm_lp_na(1,:), cis_dm_lp_na(2,:), ...
'Non-augmented LP, delta method interval', ...
plot_ylim, plot_xtick);
subplot(2,1,2);
plot_band(horzs, irs_lp_na, ...
cis_boot_lp_na(1,:,3), cis_boot_lp_na(2,:,3), ...
'Non-augmented LP, percentile-t interval', ...
plot_ylim, plot_xtick);
drawnow;
%% Non-augmented VAR
disp('Non-augmented VAR: bootstrapping...');
[irs_var_na, ~, cis_dm_var_na, cis_boot_var_na] = ...
ir_estim(Y, p, horzs, ...
'estimator', 'var', 'lag_aug', false, ...
'resp_var', resp_var, 'innov', innov, 'alpha', alpha, ...
'bootstrap', 'var', 'boot_num', boot_num, ...
'boot_workers', boot_workers, 'verbose', true);
figure;
subplot(2,1,1);
plot_band(horzs, irs_var_na, ...
cis_dm_var_na(1,:), cis_dm_var_na(2,:), ...
'Non-augmented VAR, delta method interval', ...
plot_ylim, plot_xtick);
subplot(2,1,2);
plot_band(horzs, irs_var_na, ...
cis_boot_var_na(1,:,3), cis_boot_var_na(2,:,3), ...
'Non-augmented VAR, percentile-t interval', ...
plot_ylim, plot_xtick);
drawnow;
delete(poolobj);
%% Auxiliary plot function
function plot_band(x, y, lower, upper, plot_title, plot_ylim, plot_xtick)
% Plot IRF and error bands
plot(x, y, '-k', 'LineWidth', 2);
hold on;
plot(x, [lower; upper], 'LineStyle', '--', 'Color', 'k');
line([min(x) max(x)], [0 0], 'Color', 'k');
hold off;
xlim([min(x) max(x)]);
ylim(plot_ylim);
set(gca, 'XTick', plot_xtick);
set(gca, 'FontSize', 12);
title(plot_title, 'FontSize', 14);
end