add utilities for forecast, irf, plot_irf, and their relevant options

covariance-quadratic-approximation
Marco Ratto 2022-01-31 14:15:16 +01:00
parent b1aa88e8da
commit 91a2cd2496
4 changed files with 445 additions and 4 deletions

142
matlab/+occbin/forecast.m Normal file
View File

@ -0,0 +1,142 @@
function [oo_, error_flag] = forecast(options_,M_,oo_,forecast) %,hist_period)
%function oo_ = forecast(options_,M_,oo_,forecast)
% forecast
opts = options_.occbin.forecast;
options_.occbin.simul.maxit = opts.maxit;
options_.occbin.simul.check_ahead_periods = opts.check_ahead_periods;
options_.occbin.simul.periods = forecast;
SHOCKS0 = opts.SHOCKS0;
if ~isempty(SHOCKS0)
if iscell(SHOCKS0)
for j=1:length(SHOCKS0)
sname = SHOCKS0{j}{1};
inds(j)=strmatch(sname,M_.exo_names);
SHOCKS1(j,:)=SHOCKS0{j}{2};
end
elseif isreal(SHOCKS0)
SHOCKS1=SHOCKS0;
inds = 1:M_.exo_nbr;
end
end
if opts.replic
h = dyn_waitbar(0,'Please wait occbin forecast replic ...');
ishock = find(sqrt(diag((M_.Sigma_e))));
effective_exo_nbr= length(ishock);
effective_exo_names = M_.exo_names(ishock);
effective_Sigma_e = M_.Sigma_e(ishock,ishock);
[U,S,V] = svd(effective_Sigma_e);
if opts.qmc
opts.replic =2^(round(log2(opts.replic+1)))-1;
SHOCKS_ant = qmc_sequence(forecast*effective_exo_nbr, int64(1), 1, opts.replic)';
else
SHOCKS_ant = randn(forecast*effective_exo_nbr,opts.replic)';
end
zlin0=zeros(forecast,M_.endo_nbr,opts.replic);
zpiece0=zeros(forecast,M_.endo_nbr,opts.replic);
for iter=1:opts.replic
if ~isempty(SHOCKS0)
for j=1:length(SHOCKS0)
SHOCKS(:,inds(j))=SHOCKS1(j,:);
end
end
error_flagx=1;
% while error_flagx,
% SHOCKS=transpose(sqrt(diag(diag(effective_Sigma_e)))*(reshape(SHOCKS_ant(iter,:),forecast,effective_exo_nbr))');
SHOCKS=transpose(U*sqrt(S)*(reshape(SHOCKS_ant(iter,:),forecast,effective_exo_nbr))');
% SHOCKS=transpose(U*sqrt(S)*randn(forecast,M_.exo_nbr)'); %realized shocks
options_.occbin.simul.endo_init = M_.endo_histval(:,1)-oo_.steady_state;
options_.occbin.simul.init_regime = opts.frcst_regimes;
options_.occbin.simul.init_binding_indicator = [];
options_.occbin.simul.exo_pos=ishock;
options_.occbin.simul.SHOCKS = SHOCKS;
options_.occbin.simul.waitbar=0;
[~, out] = occbin.solver(M_,oo_,options_);
zlin0(:,:,iter)=out.linear;
zpiece0(:,:,iter)=out.piecewise;
ys=out.ys;
frcst_regime_history(iter,:)=out.regime_history;
error_flag(iter)=out.error_flag;
error_flagx = error_flag(iter);
% end
simul_SHOCKS(:,:,iter) = SHOCKS;
if error_flag(iter) && debug_flag
% display('no solution')
% keyboard;
save no_solution SHOCKS zlin0 zpiece0 iter frcst_regime_history
end
dyn_waitbar(iter/opts.replic,h,['OccBin MC forecast replic ',int2str(iter),'/',int2str(opts.replic)])
end
dyn_waitbar_close(h);
save temp zlin0 zpiece0 simul_SHOCKS error_flag
inx=find(error_flag==0);
zlin0=zlin0(:,:,inx);
zpiece0=zpiece0(:,:,inx);
zlin = mean(zlin0,3);
zpiece = mean(zpiece0,3);
zpiecemin = min(zpiece0,[],3);
zpiecemax = max(zpiece0,[],3);
zlinmin = min(zlin0,[],3);
zlinmax = max(zlin0,[],3);
for i=1:M_.endo_nbr
for j=1:forecast
[post_mean(j,1), post_median(j,1), post_var(j,1), hpd_interval(j,:), post_deciles(j,:)] = posterior_moments(squeeze(zlin0(j,i,:)),1,options_.forecasts.conf_sig);
end
oo_.occbin.linear_forecast.Mean.(M_.endo_names{i})=post_mean;
oo_.occbin.linear_forecast.Median.(M_.endo_names{i})=post_median;
oo_.occbin.linear_forecast.Var.(M_.endo_names{i})=post_var;
oo_.occbin.linear_forecast.HPDinf.(M_.endo_names{i})=hpd_interval(:,1);
oo_.occbin.linear_forecast.HPDsup.(M_.endo_names{i})=hpd_interval(:,2);
oo_.occbin.linear_forecast.Deciles.(M_.endo_names{i})=post_deciles;
oo_.occbin.linear_forecast.Min.(M_.endo_names{i})=zlinmin(:,i);
oo_.occbin.linear_forecast.Max.(M_.endo_names{i})=zlinmax(:,i);
for j=1:forecast
[post_mean(j,1), post_median(j,1), post_var(j,1), hpd_interval(j,:), post_deciles(j,:)] = posterior_moments(squeeze(zpiece0(j,i,:)),1,options_.forecasts.conf_sig);
end
oo_.occbin.forecast.Mean.(M_.endo_names{i})=post_mean;
oo_.occbin.forecast.Median.(M_.endo_names{i})=post_median;
oo_.occbin.forecast.Var.(M_.endo_names{i})=post_var;
oo_.occbin.forecast.HPDinf.(M_.endo_names{i})=hpd_interval(:,1);
oo_.occbin.forecast.HPDsup.(M_.endo_names{i})=hpd_interval(:,2);
oo_.occbin.forecast.Deciles.(M_.endo_names{i})=post_deciles;
oo_.occbin.forecast.Min.(M_.endo_names{i})=zpiecemin(:,i);
oo_.occbin.forecast.Max.(M_.endo_names{i})=zpiecemax(:,i);
% eval([M_.endo_names{i},'_ss=zdatass(i);']);
end
else
SHOCKS = zeros(forecast,M_.exo_nbr);
if ~isempty(SHOCKS0)
for j=1:length(SHOCKS0)
SHOCKS(:,inds(j))=SHOCKS1(j,:);
end
end
options_.occbin.simul.endo_init = M_.endo_histval(:,1)-oo_.steady_state;
options_.occbin.simul.init_regime = opts.frcst_regimes;
options_.occbin.simul.init_violvecbool = [];
options_.occbin.simul.irfshock = M_.exo_names;
options_.occbin.simul.SHOCKS = SHOCKS;
[~, out] = occbin.solver(M_,oo_,options_);
zlin=out.linear;
zpiece=out.piecewise;
frcst_regime_history=out.regime_history;
error_flag=out.error_flag;
for i=1:M_.endo_nbr
oo_.occbin.linear_forecast.Mean.(M_.endo_names{i})= zlin(:,i);
oo_.occbin.forecast.Mean.(M_.endo_names{i})= zpiece(:,i);
oo_.occbin.forecast.HPDinf.(M_.endo_names{i})= nan;
oo_.occbin.forecast.HPDsup.(M_.endo_names{i})= nan;
end
end
oo_.occbin.forecast.regimes=frcst_regime_history;

140
matlab/+occbin/irf.m Normal file
View File

@ -0,0 +1,140 @@
function [oo_] = irf(M_,oo_,options_)
shocknames = options_.occbin.irf.exo_names;
shocksigns = options_.occbin.irf.shocksigns;
shocksize = options_.occbin.irf.shocksize;
t0 = options_.occbin.irf.t0;
options_.occbin.simul.init_regime = options_.occbin.irf.init_regime;
options_.occbin.simul.check_ahead_periods = options_.occbin.irf.check_ahead_periods;
options_.occbin.simul.maxit = options_.occbin.irf.maxit;
options_.occbin.simul.periods = options_.irf;
% Run inital conditions + other shocks
if t0 == 0
shocks0= zeros(options_.occbin.simul.periods+1,M_.exo_nbr);
options_.occbin.simul.endo_init = [];
else
% girf conditional to smoothed states in t0 and shocks in t0+1
shocks0= [oo_.occbin.smoother.etahat(:,t0+1)'; zeros(options_.occbin.simul.periods,M_.exo_nbr)];
options_.occbin.simul.SHOCKS=shocks0;
options_.occbin.simul.endo_init = oo_.occbin.smoother.alphahat(oo_.dr.inv_order_var,t0);
end
options_.occbin.simul.SHOCKS=shocks0;
[~, out0] = occbin.solver(M_,oo_,options_);
zlin0 = out0.linear;
zpiece0 = out0.piecewise;
% Select shocks of interest
jexo_all =zeros(size(shocknames,1),1);
for i=1:length(shocknames)
jexo_all(i) = strmatch(shocknames{i},M_.exo_names,'exact');
end
oo_.occbin.linear_irfs = struct();
oo_.occbin.irfs = struct();
% Set shock size
if isempty(shocksize)
% if isfield(oo_.posterior_mode.shocks_std,M_.exo_names{jexo})
shocksize = sqrt(diag(M_.Sigma_e(jexo_all,jexo_all))); %oo_.posterior_mode.shocks_std.(M_.exo_names{jexo});
if any(shocksize < 1.e-9)
shocksize(shocksize < 1.e-9) = 0.01;
end
end
if numel(shocksize)==1
shocksize=repmat(shocksize,[length(shocknames),1]);
end
% Run IRFs
for counter = 1:length(jexo_all)
jexo = jexo_all(counter);
if ~options_.noprint
fprintf('Producing GIRFs for shock %s. Simulation %d out of %d. \n',M_.exo_names{jexo},counter,size(jexo_all,1));
end
if ismember('pos',shocksigns)
% (+) shock
shocks1=shocks0;
shocks1(1,jexo)=shocks0(1,jexo)+shocksize(counter);
if t0 == 0
options_.occbin.simul.SHOCKS=shocks1;
options_.occbin.simul.endo_init = [];
[~, out_pos] = occbin.solver(M_,oo_,options_);
else
options_.occbin.simul.SHOCKS=shocks1;
options_.occbin.simul.endo_init = oo_.occbin.smoother.alphahat(oo_.dr.inv_order_var,t0);
[~, out_pos] = occbin.solver(M_,oo_,options_);
end
if out_pos.error_flag
warning('Occbin error.')
return
end
zlin_pos = out_pos.linear;
zpiece_pos = out_pos.piecewise;
% Substract inital conditions + other shocks
zlin_pos_diff = zlin_pos-zlin0;
zpiece_pos_diff = zpiece_pos-zpiece0;
end
if ismember('neg',shocksigns)
% (-) shock
shocks_1=shocks0;
shocks_1(1,jexo)=shocks0(1,jexo)-shocksize(counter);
if t0 == 0
options_.occbin.simul.SHOCKS=shocks_1;
options_.occbin.simul.endo_init = [];
[~, out_neg] = occbin.solver(M_,oo_,options_);
else
options_.occbin.simul.SHOCKS=shocks_1;
options_.occbin.simul.endo_init = oo_.occbin.smoother.alphahat(oo_.dr.inv_order_var,t0);
[~, out_neg] = occbin.solver(M_,oo_,options_);
end
if out_neg.error_flag
warning('Occbin error.')
return
end
zlin_neg = out_neg.linear;
zpiece_neg = out_neg.piecewise;
zlin_neg_diff = zlin_neg-zlin0;
zpiece_neg_diff = zpiece_neg-zpiece0;
end
% Save
if ~isfield(oo_.occbin,'linear_irfs')
oo_.occbin.linear_irfs = struct();
end
if ~isfield(oo_.occbin,'irfs')
oo_.occbin.irfs = struct();
end
for jendo=1:M_.endo_nbr
% oo_.occbin.irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '1']) = zpiece_pos (:,jendo);
% oo_.occbin.irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_1']) = zpiece_neg (:,jendo);
% oo_.occbin.linear_irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '1']) = zlin_pos (:,jendo);
% oo_.occbin.linear_irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_1']) = zlin_neg(:,jendo);
if ismember('pos',shocksigns)
oo_.occbin.irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_pos']) = zpiece_pos_diff (:,jendo);
oo_.occbin.linear_irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_pos']) = zlin_pos_diff (:,jendo);
end
if ismember('neg',shocksigns)
oo_.occbin.irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_neg']) = zpiece_neg_diff (:,jendo);
oo_.occbin.linear_irfs.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_neg']) = zlin_neg_diff (:,jendo);
end
% %
% oo_.occbin.irfs0.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '1']) = zpiece0(:,jendo);
% oo_.occbin.linear_irfs0.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '1']) = zlin0(:,jendo);
% oo_.occbin.irfs0.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_1']) = zpiece0(:,jendo);
% oo_.occbin.linear_irfs0.([M_.endo_names{jendo} '_' M_.exo_names{jexo} '_1']) = zlin0(:,jendo);
end
end
end

141
matlab/+occbin/plot_irfs.m Normal file
View File

@ -0,0 +1,141 @@
function plot_irfs(M_,oo_,options_,irf3,irf4)
shocknames = options_.occbin.plot_irf.exo_names;
simulname = options_.occbin.plot_irf.simulname;
if isempty(simulname)
simulname_ = simulname;
else
simulname_ = [ simulname '_' ];
end
vars_irf = options_.occbin.plot_irf.endo_names;
endo_names_long = options_.occbin.plot_irf.endo_names_long;
endo_scaling_factor = options_.occbin.plot_irf.endo_scaling_factor;
length_irf = options_.occbin.plot_irf.tplot;
if isempty(length_irf)
length_irf = options_.irf;
end
irflocation_lin = oo_.occbin.linear_irfs;
irflocation_piece = oo_.occbin.irfs;
steps_irf = 1;
warning('off','all')
DirectoryName = CheckPath('graphs',M_.dname);
iexo=[];
for i=1:size(shocknames,1)
itemp = strmatch(shocknames{i},M_.exo_names,'exact');
if isempty(itemp)
error(['Shock ',shocknames{i},' is not defined!'])
else
iexo=[iexo, itemp];
end
end
ncols = options_.occbin.plot_irf.ncols;
nrows = options_.occbin.plot_irf.nrows;
npan = ncols*nrows;
plot_grid = options_.occbin.plot_irf.grid;
shocksigns = options_.occbin.plot_irf.shocksigns;
threshold = options_.occbin.plot_irf.threshold;
% Add steady_state
if options_.occbin.plot_irf.add_steadystate
add_stst = options_.occbin.plot_irf.add_steadystate;
else
add_stst = 0;
end
for sss = 1:numel(shocksigns)
shocksign = shocksigns{sss};
for j=1:size(shocknames,1)
%shocknames = M_.exo_names{j};
j1 = 0;
isub = 0;
ifig = 0;
% Variables
% ----------------------
for i = 1:length(vars_irf)
j1=j1+1;
if mod(j1,npan)==1
% vector corresponds to [left bottom width height]. 680 and 678 for the left and bottom elements correspond to the default values used by MATLAB while creating a figure and width, .
hfig = dyn_figure(options_.nodisplay,'name',['OccbinIRFs ' shocknames{j} ' ' simulname],'PaperPositionMode', 'auto','PaperType','A4','PaperOrientation','portrait','renderermode','auto','position',[10 10 950 650]);
ifig=ifig+1;
isub=0;
end
isub=isub+1;
if isempty(endo_scaling_factor)
exofactor = 1;
else
exofactor = endo_scaling_factor{i};
end
subplot(nrows,ncols,isub)
irf_field = strcat(vars_irf{i,1},'_',shocknames{j},'_',shocksign);
irfvalues = irflocation_lin.(irf_field);
if add_stst
irfvalues = irfvalues + get_mean(vars_irf{i,1});
end
irfvalues(abs(irfvalues) <threshold) = 0;
plot(irfvalues(1:steps_irf:length_irf)*exofactor,'linewidth',2);
hold on
irfvalues = irflocation_piece.(irf_field);
if add_stst
irfvalues = irfvalues + get_mean(vars_irf{i,1});
end
irfvalues(abs(irfvalues) <threshold) = 0;
plot(irfvalues(1:steps_irf:length_irf)*exofactor,'r--','linewidth',2);
hold on
plot(irfvalues(1:steps_irf:length_irf)*0,'k-','linewidth',1.5);
% Optional additional IRFs
if nargin > 10
irfvalues = irf3.(irf_field) ;
irfvalues(abs(irfvalues) <threshold) = 0;
plot(irfvalues(1:steps_irf:length_irf)*exofactor,'k:','linewidth',2);
end
if nargin > 11
irfvalues = irf4.(irf_field) ;
irfvalues(abs(irfvalues) <threshold) = 0;
plot(irfvalues(1:steps_irf:length_irf)*exofactor,'g-.','linewidth',2);
end
if plot_grid
grid on
end
xlim([1 (length_irf/steps_irf)]);
% title
if isempty(endo_names_long)
title(regexprep(vars_irf{i},'_',' '))
else
title(endo_names_long{i})
end
% Annotation Box + save figure
% ----------------------
if mod(j1,npan)==0 || (mod(j1,npan)~=0 && i==length(vars_irf))
annotation('textbox', [0.1,0,0.35,0.05],'String', 'Linear','Color','Blue','horizontalalignment','center','interpreter','none');
annotation('textbox', [0.55,0,0.35,0.05],'String', 'Piecewise','Color','Red','horizontalalignment','center','interpreter','none');
dyn_saveas(hfig,[DirectoryName,filesep,M_.fname,'_irf_occbin_',simulname_,shocknames{j},'_',shocksign,'_',int2str(ifig)],options_.nodisplay,options_.graph_format);
end
end
end
end
warning('on','all')
end

View File

@ -46,23 +46,26 @@ if ismember(flag,{'filter','all'})
end
if ismember(flag,{'forecast','all'})
options_occbin_.forecast.check_ahead_periods=30;
options_occbin_.forecast.debug_flag=false;
options_occbin_.forecast.frcst_regimes=[];
options_occbin_.forecast.maxit=30;
options_occbin_.forecast.periods=30;
options_occbin_.forecast.qmc=0;
options_occbin_.forecast.replic=0;
options_occbin_.forecast.sepath=0;
options_occbin_.forecast.SHOCKS0=[];
options_occbin_.forecast.treepath=1; % number of branches
end
if ismember(flag,{'irf','all'})
if ismember(flag,{'irf','all'})
options_occbin_.irf.check_ahead_periods=30;
options_occbin_.irf.exo_names=M_.exo_names;
options_occbin_.irf.init_regime=[];
options_occbin_.irf.maxit=30;
options_occbin_.irf.threshold = 10^-6;
% options_occbin_.irf.periods=options_.irf;
options_occbin_.irf.shocksize=[];
options_occbin_.irf.shocksigns = {'1','_1'};
options_occbin_.irf.shocksigns = {'pos','neg'};
options_occbin_.irf.t0=0;
end
if ismember(flag,{'likelihood','all'})
@ -87,6 +90,21 @@ if ismember(flag,{'likelihood','all'})
options_occbin_.likelihood.waitbar=false;
end
if ismember(flag,{'plot_irf','all'})
options_occbin_.plot_irf.add_steadystate = 0;
options_occbin_.plot_irf.exo_names = [];
options_occbin_.plot_irf.endo_names = M_.endo_names;
options_occbin_.plot_irf.endo_names_long = [];
options_occbin_.plot_irf.endo_scaling_factor = [];
options_occbin_.plot_irf.grid = true;
options_occbin_.plot_irf.ncols = 3;
options_occbin_.plot_irf.nrows = 3;
options_occbin_.plot_irf.shocksigns = {'pos','neg'};
options_occbin_.plot_irf.simulname='';
options_occbin_.plot_irf.threshold = 10^-6;
options_occbin_.plot_irf.tplot = [];
end
if ismember(flag,{'plot_shock_decomp','all'})
options_occbin_.plot_shock_decomp.add_steadystate = false;
options_occbin_.plot_shock_decomp.add_zero_line = false;