forked from dake-li/lp_var_simul
-
Notifications
You must be signed in to change notification settings - Fork 0
/
run_plot_dgp.m
307 lines (238 loc) · 14.2 KB
/
run_plot_dgp.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
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
%% DFM SIMULATION STUDY: TABLES AND PLOTS FOR DGP FEATURES
% Dake Li, Mikkel Plagborg-Møller and Christian Wolf
% This version: 02/23/2021
clc
close all
clear all;
addpath('Plotting_Functions');
addpath(genpath(fullfile('..', 'Subroutines')))
%% SETTINGS
% select robustness check mode
mode_select = 1; % options: 1 (baseline), 2 (small sample), 3 (large sample),
% 4 (salient series), 5 (more observables), 6 (first diff)
% select lag length specifications
lags_select = 2; % options: 1 (AIC), 2 (4 lags), 3 (8 lags), 4 (12 lags)
% select and group experiments
exper_select_group = {[2,5], [3,6], [1,4]}; % combine G and MP for observed shock, recursive, and IV
% select estimation methods for each experiment
methods_iv_select = [1 2 3 4 5 6 7 8];
methods_obsshock_select = [1 2 3 4 5 6 7];
methods_recursive_select = [1 2 3 4 5 6 7];
% select a subset of DGPs
DGP_select = 0; % options: 0 (all DGPs), 1 (specifications with asset price & sentiment),
% 2 (low degree of invertibility), 3 (high degree of invertibility)
% Apply shared settings
settings_shared;
% Summary statistics for table
tab_stat = {'LRV_Cov_tr_ratio', 'dLRV_dCov_tr_ratio', 'VAR_largest_root', 'VAR_quant_root', 'frac_coef_for_large_lags', 'R0_sq'}; % Summary stats to copy (in addition to IRF stats defined below)
tab_quants = [0.1 0.25 0.5 0.75 0.9]; % Quantiles to report across specifications
% Table with model specification tests
spec_lags_cutoffs = [0 1 2]; % Report how often quantile of lag length strictly exceeds each of these numbers
spec_lags_quant = 0.9; % Quantile of lag length to report
spec_lm_powers = [0.1 0.25 0.5]; % Report fraction of DGPs with at least this power of LM test
spec_lm_signif = 0.1; % Significance level for LM test
% IRF examples to plot
spec_select = [10 20 30 3010 3020 3030];
linestyles = {'-', '--', ':', '-o', '--o', ':o'};
colors = lines(6);
%% CREATE TABLES AND PLOTS
for n_mode=1:length(mode_folders) % For each robustness check mode...
for nf=1:length(lags_folders) % For each lag-order folder...
for ne=1:length(exper_files) % For each experiment in folder...
%----------------------------------------------------------------
% Load Results
%----------------------------------------------------------------
load_results;
% see if ready to plot for this group of experiments
if exper_group_end(ne) == 0
continue;
end
% keep only the selected subset of DGPs
if DGP_select > 0
DGP_selected = arrayfun(@(x) select_DGP_fn(x,res), 1:res.settings.specifications.n_spec)'; % binary DGP selection label
res = combine_struct(res,[],[],DGP_selected);
end
% Record the index of median across MCs
median_idx = stat_index(0.5, res.settings); % index of median number in the quantile list (including mean and std)
%----------------------------------------------------------------
% Tables of Summary Statistics
%----------------------------------------------------------------
tab = table;
% Basic DGP summary statistics
for is=1:length(tab_stat)
tab.(tab_stat{is}) = res.DF_model.(tab_stat{is});
end
% IV F-stat
if isfield(res.DF_model, 'IV_strength')
tab.iv_F = res.settings.simul.T./(1./res.DF_model.IV_strength-1);
end
% IRF summary statistics
tab.irf_num_local_extrema = sum(diff(sign(diff(res.DF_model.target_irf',1,2)),1,2)~=0,2);
[~,I] = max(abs(res.DF_model.target_irf)',[],2);
tab.irf_maxabs_h = res.settings.est.IRF_select(I)'-1;
tab.irf_mean_max_ratio = mean(res.DF_model.target_irf',2)./max(abs(res.DF_model.target_irf)',[],2);
% R-squared from regressing IRFs on quadratic
X = [ones(res.settings.est.IRF_hor,1) (res.settings.est.IRF_select').^(1:2)];
betas = X\res.DF_model.target_irf;
resid = res.DF_model.target_irf-X*betas;
tab.irf_R2s = 1-(var(resid)./var(res.DF_model.target_irf))';
% Create table of quantiles
tab_summ = varfun(@(x) [min(x) quantile(x,tab_quants) max(x)]', tab);
tab_summ.Properties.VariableNames=regexprep(tab_summ.Properties.VariableNames, 'Fun_', '');
tab_summ2 = table;
tab_summ2.quantile = [0; tab_quants(:); 1];
tab_summ = [tab_summ2 tab_summ];
% Write to file
writetable(tab_summ, fullfile(output_folder, 'dgp_summ.csv'));
clearvars I X betas resid tab tab_summ tab_summ2;
%----------------------------------------------------------------
% Tables of Specification Tests
%----------------------------------------------------------------
tab_spec = table;
% Number of lags
for ii=1:length(spec_lags_cutoffs)
tab_spec.(sprintf('%s%d', 'nlag_exceed_', spec_lags_cutoffs(ii))) ...
= mean(res.results.n_lags.svar(stat_index(spec_lags_quant, res.settings),:)>spec_lags_cutoffs(ii));
end
% Power of LM test
for ii=1:length(spec_lm_powers)
tab_spec.(sprintf('%s%02d', 'lm_power_', round(100*spec_lm_powers(ii)))) ...
= mean(res.results.LM_pvalue.svar(stat_index(spec_lm_powers(ii), res.settings),:)<spec_lm_signif);
end
% Write to file
writetable(tab_spec, fullfile(output_folder, 'dgp_spec.csv'));
clearvars tab_spec;
%----------------------------------------------------------------
% Plots to summarize DGP Features
%----------------------------------------------------------------
% Degree of invertibility
figure;
histogram(res.DF_model.R0_sq, 'Normalization', 'probability');
title(strjoin({exper_plotname, ': degree of invertibility'}), 'Interpreter', 'none');
plot_save(fullfile(output_folder, 'dgp_R0sq'), output_suffix);
% Persistence
figure;
histogram(res.DF_model.LRV_Cov_tr_ratio, 'Normalization', 'probability');
title(strjoin({exper_plotname, ': LRV to Var ratio'}), 'Interpreter', 'none');
plot_save(fullfile(output_folder, 'dgp_LRVVar'), output_suffix);
figure;
histogram(res.DF_model.dLRV_dCov_tr_ratio, 'Normalization', 'probability');
title(strjoin({exper_plotname, ': LRV to Var ratio, first differences'}), 'Interpreter', 'none');
plot_save(fullfile(output_folder, 'dgp_dLRVdVar'), output_suffix);
figure;
histogram(res.DF_model.VAR_largest_root, 'Normalization', 'probability');
title(strjoin({exper_plotname, ': largest VAR root'}), 'Interpreter', 'none');
plot_save(fullfile(output_folder, 'dgp_largroot'), output_suffix);
figure;
histogram(res.DF_model.VAR_quant_root, 'Normalization', 'probability');
title(strjoin({exper_plotname, ': quantile of VAR roots'}), 'Interpreter', 'none');
plot_save(fullfile(output_folder, 'dgp_quantroot'), output_suffix);
figure;
histogram(res.DF_model.frac_coef_for_large_lags, 'Normalization', 'probability');
title(strjoin({exper_plotname, ': fraction of long-lag VAR coefs'}), 'Interpreter', 'none');
plot_save(fullfile(output_folder, 'dgp_longlag'), output_suffix);
% IV strength
if isfield(res.DF_model, 'IV_strength')
figure;
histogram(res.DF_model.IV_strength, 'Normalization', 'probability');
title(strjoin({exper_plotname, ': IV strength'}), 'Interpreter', 'none');
plot_save(fullfile(output_folder, 'dgp_IVstrength'), output_suffix);
end
% Some IRFs
if DGP_select == 0
figure('Units', 'inches', 'Position', [0 0 8 4]);
hold on;
norm_irf = @(x) x/max(abs(x));
for i_spec_indx = 1:length(spec_select)
i_spec = spec_select(i_spec_indx);
plot(res.settings.est.IRF_select-1, norm_irf(res.DF_model.target_irf(:,i_spec)), ...
linestyles{i_spec_indx}, 'Color', colors(i_spec_indx,:), 'Linewidth', 2, ...
'MarkerSize', 4, 'MarkerFaceColor', colors(i_spec_indx,:));
end
hold off;
xlabel('Horizon','interpreter','latex','FontSize',12);
set(gca,'XTick',[min(res.settings.est.IRF_select-1) 2:2:max(res.settings.est.IRF_select-1)]);
xlim([min(res.settings.est.IRF_select-1) max(res.settings.est.IRF_select-1)]);
grid on;
set(gca,'TickLabelInterpreter','latex');
set(gca,'FontSize',12);
plot_save(fullfile(output_folder, 'dgp_irfs'), output_suffix);
end
%----------------------------------------------------------------
% Tuning Parameters for Estimation Methods
%----------------------------------------------------------------
% Number of lags
the_nlags = res.results.n_lags.svar;
figure;
histogram(the_nlags(median_idx,:), 'Normalization', 'probability');
title(strjoin({exper_plotname, ': median number of lags (across specs)'}), 'Interpreter', 'none');
plot_save(fullfile(output_folder, 'dgp_nlags'), output_suffix);
figure;
histogram(the_nlags(2,:), 'Normalization', 'probability');
title(strjoin({exper_plotname, ': std (across sims) of number of lags'}), 'Interpreter', 'none');
plot_save(fullfile(output_folder, 'dgp_nlags_std'), output_suffix);
% Shrinkage penalty
the_lambda = res.results.lambda.lp_penalize;
figure;
[~,the_edges] = histcounts(log10(the_lambda(median_idx,:)));
histogram(the_lambda(median_idx,:),10.^the_edges);
set(gca, 'xscale','log'); % Log scale for x axis
title(strjoin({exper_plotname, ': median shrinkage penalty (across specs)'}), 'Interpreter', 'none');
plot_save(fullfile(output_folder, 'dgp_lambda'), output_suffix);
% Giannone, Lenza & Primiceri BVAR hyper-parameters
if isfield(res.results,'GLP_hyper')
the_GLP = res.results.GLP_hyper;
the_fields = fieldnames(the_GLP);
figure;
for nj = 1:length(the_fields)
subplot(1,3,nj);
histogram(the_GLP.(the_fields{nj})(median_idx,:), 'Normalization', 'probability');
title(the_fields{nj}, 'Interpreter', 'none');
end
sgtitle(strjoin({exper_plotname, ': median GLP hyper-param (across specs)'}), 'FontSize', 11, 'FontWeight', 'bold', 'Interpreter', 'none');
plot_save(fullfile(output_folder, 'dgp_glphyper'), output_suffix);
end
% Model-averaging weights
if isfield(res.results, 'weight')
the_weights = res.results.weight.var_avg;
the_store_weights = res.settings.est.average_store_weight;
the_maxlag = res.settings.est.n_lags_max;
figure;
for j=1:length(the_store_weights) % For each horizon where weights are stored...
subplot(1,length(the_store_weights),j);
plot(1:the_maxlag, mean(reshape(the_weights(1:the_maxlag,j,1,:), the_maxlag, []), 2)); % AR weights
hold on;
plot(1:the_maxlag, mean(reshape(the_weights(the_maxlag+1:end,j,1,:), the_maxlag, []), 2)); % VAR weights
hold off;
title(sprintf('%s%d', 'h = ', the_store_weights(j)));
xlabel('no. of lags');
legend('AR', 'VAR', 'Location', 'NorthEast');
end
sgtitle(strjoin({exper_plotname, ': average model-avg weight (across specs+sims)'}), 'FontSize', 11, 'FontWeight', 'bold', 'Interpreter', 'none');
plot_save(fullfile(output_folder, 'dgp_weight'), output_suffix);
end
%----------------------------------------------------------------
% IV F-Statistics
%----------------------------------------------------------------
if isfield(res.results, 'F_stat')
the_Fstats = res.results.F_stat.svar_iv;
% Average F stat
figure;
histogram(the_Fstats(1,:), 'Normalization', 'probability');
title(strjoin({exper_plotname, ': average F stat (across specs)'}), 'Interpreter', 'none');
plot_save(fullfile(output_folder, 'dgp_F'), output_suffix);
end
%----------------------------------------------------------------
% Granger Causality
%----------------------------------------------------------------
if isfield(res.results, 'Granger_stat')
the_Grangerstats = res.results.Granger_stat.svar;
% Average Granger causality Wald-stat
figure;
histogram(the_Grangerstats(1,:), 'Normalization', 'probability');
title(strjoin({exper_plotname, ': average Granger causality Wald stat (across specs)'}), 'Interpreter', 'none');
plot_save(fullfile(output_folder, 'dgp_Granger'), output_suffix);
end
end
end
end