stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Enable filter_covariance option in posterior sampling

time-shift
Johannes Pfeifer 2016-06-02 18:28:36 +02:00 committed by Stéphane Adjemian (Lupi)
parent 65ae97d37d
commit 5d6d1336ef
5 changed files with 116 additions and 38 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
doc/dynare.texi
View File

@ -5809,8 +5809,8 @@ recursions as described by @cite{Herbst, 2015}. This setting is only used with
@item filter_covariance
@anchor{filter_covariance} Saves the series of one step ahead error of
forecast covariance matrices in @ref{oo_.Smoother.Variance}. Not available
with Metropolis.
forecast covariance matrices. With Metropolis, they are saved in @ref{oo_.FilterCovariance},
otherwise in @ref{oo_.Smoother.Variance}.
@item filter_step_ahead = [@var{INTEGER1}:@var{INTEGER2}]
See below.
@ -6287,6 +6287,23 @@ After an estimation with Metropolis, fields are of the form:
@end example
@end defvr
@defvr {MATLAB/Octave variable} oo_.FilterCovariance
@anchor{oo_.FilterCovariance}
Three-dimensional array set by the @code{estimation} command if used with the
@code{smoother} and Metropolis, if the @code{filter_covariance} option
has been requested.
Contains the series of one-step ahead forecast error covariance matrices
from the Kalman smoother. The @code{M_.endo_nbr} times @code{M_.endo_nbr} times
@code{T+1} array contains the variables in declaration order along the first
two dimensions. The third dimension of the array provides the
observation for which the forecast has been made.
Fields are of the form:
@example
@code{oo_.FilterCovariance.@var{MOMENT_NAME}}
@end example
Note that density estimation is not supported.
@end defvr
@defvr {MATLAB/Octave variable} oo_.Smoother.Variance
@anchor{oo_.Smoother.Variance}
Three-dimensional array set by the @code{estimation} command (if used with the

42
matlab/pm3.m
View File

@ -86,6 +86,8 @@ fprintf(['Estimation::mcmc: ' tit1 '\n']);
stock1 = zeros(n1,n2,B);
k = 0;
filter_step_ahead_indicator=0;
filter_covar_indicator=0;
for file = 1:ifil
load([DirectoryName '/' M_.fname var_type int2str(file)]);
if strcmp(var_type,'_filter_step_ahead')
@ -100,10 +102,22 @@ for file = 1:ifil
end
stock = squeeze(stock(1,:,1+1:1+n2,:)); %1 step ahead starts at entry 2
end
k = k(end)+(1:size(stock,3));
stock1(:,:,k) = stock;
if strcmp(var_type,'_filter_covar')
if file==1 %on first run, initialize variable for storing filter_step_ahead
stock1_filter_covar=NaN(n1,n2,size(stock,3),B);
end
filter_covar_indicator=1;
end
if ~strcmp(var_type,'_filter_covar')
k = k(end)+(1:size(stock,3));
stock1(:,:,k) = stock;
else
k = k(end)+(1:size(stock,4));
end
if filter_step_ahead_indicator
stock1_filter_step_ahead(:,:,k,:) = stock_filter_step_ahead;
elseif filter_covar_indicator
stock1_filter_covar(:,:,:,k) = stock;
end
end
clear stock
@ -119,8 +133,30 @@ if filter_step_ahead_indicator
Density_filter_step_ahead = zeros(options_.estimation.moments_posterior_density.gridpoints,2,filter_steps,nvar,n2);
end
end
if filter_covar_indicator
draw_dimension=4;
oo_.FilterCovariance.Mean = squeeze(mean(stock1_filter_covar,draw_dimension));
oo_.FilterCovariance.Median = squeeze(median(stock1_filter_covar,draw_dimension));
oo_.FilterCovariance.var = squeeze(var(stock1_filter_covar,0,draw_dimension));
if size(stock1_filter_covar,draw_dimension)>2
hpd_interval = quantile(stock1_filter_covar,[(1-options_.mh_conf_sig)/2 (1-options_.mh_conf_sig)/2+options_.mh_conf_sig],draw_dimension);
else
size_matrix=size(stock1_filter_covar);
hpd_interval=NaN([size_matrix(1:3),2]);
end
if size(stock1_filter_covar,draw_dimension)>9
post_deciles =quantile(stock1_filter_covar,[0.1:0.1:0.9],draw_dimension);
else
size_matrix=size(stock1_filter_covar);
post_deciles=NaN([size_matrix(1:3),9]);
end
oo_.FilterCovariance.post_deciles=post_deciles;
oo_.FilterCovariance.HPDinf=squeeze(hpd_interval(:,:,:,1));
oo_.FilterCovariance.HPDsup=squeeze(hpd_interval(:,:,:,2));
fprintf(['Estimation::mcmc: ' tit1 ', done!\n']);
return
end
tmp =zeros(B,1);
for i = 1:nvar
for j = 1:n2
if options_.estimation.moments_posterior_density.indicator

57
matlab/prior_posterior_statistics.m
View File

@ -53,17 +53,14 @@ ncx = estim_params_.ncx;
ncn = estim_params_.ncn;
np = estim_params_.np ;
npar = nvx+nvn+ncx+ncn+np;
offset = npar-np;
naK = length(options_.filter_step_ahead);
MaxNumberOfBytes=options_.MaxNumberOfBytes;
endo_nbr=M_.endo_nbr;
exo_nbr=M_.exo_nbr;
meas_err_nbr=length(M_.Correlation_matrix_ME);
nvobs = length(options_.varobs);
iendo = 1:endo_nbr;
horizon = options_.forecast;
filtered_vars = options_.filtered_vars;
IdObs = bayestopt_.mfys;
if horizon
i_last_obs = gend+(1-M_.maximum_endo_lag:0);
@ -106,11 +103,14 @@ end
if horizon
MAX_nforc1 = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*(horizon+maxlag))/8));
MAX_nforc2 = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*(horizon+maxlag))/ ...
8));
MAX_nforc2 = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*(horizon+maxlag))/8));
end
MAX_momentsno = min(B,ceil(MaxNumberOfBytes/(get_moments_size(options_)*8)));
if options_.filter_covariance
MAX_filter_covariance = min(B,ceil(MaxNumberOfBytes/(endo_nbr^2*(gend+1))/8));
end
varlist = options_.varlist;
if isempty(varlist)
varlist = M_.endo_names(1:M_.orig_endo_nbr, :);
@ -123,39 +123,26 @@ for i=1:nvar
end
end
irun = ones(7,1);
ifil = zeros(7,1);
n_variables_to_fill=8;
irun = ones(n_variables_to_fill,1);
ifil = zeros(n_variables_to_fill,1);
stock_param = zeros(MAX_nruns, npar);
stock_logpo = zeros(MAX_nruns,1);
stock_ys = zeros(MAX_nruns,endo_nbr);
run_smoother = 0;
if options_.smoother
stock_smooth = zeros(endo_nbr,gend,MAX_nsmoo);
stock_innov = zeros(exo_nbr,gend,B);
stock_error = zeros(nvobs,gend,MAX_nerro);
stock_update = zeros(endo_nbr,gend,MAX_nsmoo);
if options_.smoother || options_.forecast || options_.filter_step_ahead
run_smoother = 1;
end
if options_.filter_step_ahead
stock_filter_step_ahead = zeros(naK,endo_nbr,gend+ ...
options_.filter_step_ahead(end),MAX_naK);
run_smoother = 1;
filter_covariance=0;
if options_.filter_covariance
filter_covariance=1;
end
if options_.forecast
stock_forcst_mean = zeros(endo_nbr,horizon,MAX_nforc1);
stock_forcst_point = zeros(endo_nbr,horizon,MAX_nforc2);
run_smoother = 1;
end
% Store the variable mandatory for local/remote parallel computing.
localVars.type=type;
localVars.run_smoother=run_smoother;
localVars.filter_covariance=filter_covariance;
localVars.gend=gend;
localVars.Y=Y;
localVars.data_index=data_index;
@ -182,6 +169,9 @@ localVars.MAX_nerro = MAX_nerro;
if naK
localVars.MAX_naK=MAX_naK;
end
if options_.filter_covariance
localVars.MAX_filter_covariance = MAX_filter_covariance;
end
localVars.MAX_nruns=MAX_nruns;
localVars.MAX_momentsno = MAX_momentsno;
localVars.ifil=ifil;
@ -189,6 +179,7 @@ localVars.DirectoryName = DirectoryName;
if strcmpi(type,'posterior'),
b=0;
logpost=NaN(B,1);
while b<=B
b = b + 1;
[x(b,:), logpost(b,1)] = GetOneDraw(type);
@ -206,7 +197,7 @@ if isnumeric(options_.parallel),
% Parallel execution!
else
[nCPU, totCPU, nBlockPerCPU] = distributeJobs(options_.parallel, 1, B);
ifil=zeros(7,totCPU);
ifil=zeros(n_variables_to_fill,totCPU);
for j=1:totCPU-1,
if run_smoother
nfiles = ceil(nBlockPerCPU(j)/MAX_nsmoo);
@ -228,6 +219,10 @@ else
nfiles = ceil(nBlockPerCPU(j)/MAX_nforc2);
ifil(7,j+1) =ifil(7,j)+nfiles;
end
if options_.filter_covariance
nfiles = ceil(nBlockPerCPU(j)/MAX_filter_covariance);
ifil(8,j+1) =ifil(8,j)+nfiles;
end
end
localVars.ifil = ifil;
globalVars = struct('M_',M_, ...
@ -296,11 +291,17 @@ if options_.forecast
pm3(endo_nbr,horizon,ifil(6),B,'Forecasted variables (mean)',...
'',varlist,M_.endo_names_tex,M_.endo_names,...
varlist,'MeanForecast',DirectoryName,'_forc_mean');
pm3(endo_nbr,horizon,ifil(6),B,'Forecasted variables (point)',...
pm3(endo_nbr,horizon,ifil(7),B,'Forecasted variables (point)',...
'',varlist,M_.endo_names_tex,M_.endo_names,...
varlist,'PointForecast',DirectoryName,'_forc_point');
end
if options_.filter_covariance
pm3(endo_nbr,endo_nbr,ifil(8),B,'Filtered covariances',...
'',varlist,M_.endo_names_tex,M_.endo_names,...
varlist,'FilterCovariance',DirectoryName,'_filter_covar');
end
if ~isnumeric(options_.parallel),
options_.parallel_info.leaveSlaveOpen = leaveSlaveOpen;

32
matlab/prior_posterior_statistics_core.m
View File

@ -52,6 +52,7 @@ end
type=myinputs.type;
run_smoother=myinputs.run_smoother;
filter_covariance=myinputs.filter_covariance;
gend=myinputs.gend;
Y=myinputs.Y;
data_index=myinputs.data_index;
@ -73,6 +74,9 @@ end
if naK
MAX_naK=myinputs.MAX_naK;
end
if filter_covariance
MAX_filter_covariance=myinputs.MAX_filter_covariance;
end
exo_nbr=myinputs.exo_nbr;
maxlag=myinputs.maxlag;
@ -127,6 +131,7 @@ if RemoteFlag==1,
OutputFileName_param = {};
OutputFileName_forc_mean = {};
OutputFileName_forc_point = {};
OutputFileName_filter_covar ={};
% OutputFileName_moments = {};
end
@ -149,6 +154,9 @@ end
stock_param = NaN(MAX_nruns,size(myinputs.x,2));
stock_logpo = NaN(MAX_nruns,1);
stock_ys = NaN(MAX_nruns,endo_nbr);
if filter_covariance
stock_filter_covariance = zeros(endo_nbr,endo_nbr,gend+1,MAX_filter_covariance);
end
for b=fpar:B
if strcmpi(type,'prior')
@ -167,9 +175,9 @@ for b=fpar:B
if run_smoother
[dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
[alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK,junk1,junk2,junk3,junk4,junk5,trend_addition] = ...
[alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK,junk1,junk2,P,junk4,junk5,trend_addition] = ...
DsgeSmoother(deep,gend,Y,data_index,missing_value);
if options_.loglinear %reads values from smoother results, which are in dr-order and put them into declaration order
stock_smooth(dr.order_var,:,irun(1)) = alphahat(1:endo_nbr,:)+ ...
repmat(log(SteadyState(dr.order_var)),1,gend);
@ -263,6 +271,9 @@ for b=fpar:B
stock_forcst_mean(:,:,irun(6)) = yf(maxlag+1:end,:)';
stock_forcst_point(:,:,irun(7)) = yf1(maxlag+1:end,:)';
end
if filter_covariance
stock_filter_covariance(dr.order_var,dr.order_var,:,irun(8)) = P;
end
else
[T,R,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
end
@ -270,7 +281,8 @@ for b=fpar:B
stock_logpo(irun(5),1) = logpo;
stock_ys(irun(5),:) = SteadyState';
irun = irun + ones(7,1);
irun = irun + ones(8,1);
if run_smoother && (irun(1) > MAX_nsmoo || b == B),
@ -348,6 +360,17 @@ for b=fpar:B
end
irun(7) = 1;
end
if run_smoother && filter_covariance && (irun(8) > MAX_filter_covariance || b == B)
stock = stock_filter_covariance(:,:,:,1:irun(8)-1);
ifil(8) = ifil(8) + 1;
save([DirectoryName '/' M_.fname '_filter_covar' int2str(ifil(8)) '.mat'],'stock');
if RemoteFlag==1,
OutputFileName_filter_covar = [OutputFileName_filter_covar; {[DirectoryName filesep], [M_.fname '_filter_covar' int2str(ifil(8)) '.mat']}];
end
irun(8) = 1;
end
dyn_waitbar((b-fpar+1)/(B-fpar+1),h);
end
@ -360,7 +383,8 @@ if RemoteFlag==1,
OutputFileName_filter_step_ahead;
OutputFileName_param;
OutputFileName_forc_mean;
OutputFileName_forc_point];
OutputFileName_forc_point
OutputFileName_filter_covar];
end
dyn_waitbar_close(h);

2
tests/TeX/fs2000_corr_ME.mod
View File

@ -144,7 +144,7 @@ stderr gy_obs, 1;
corr gp_obs, gy_obs,0;
end;
estimation(order=1,datafile='../fs2000/fsdat_simul',mode_check,smoother,filter_decomposition,forecast = 8,filtered_vars,filter_step_ahead=[1,3],irf=20,contemporaneous_correlation) m P c e W R k d y gy_obs;
estimation(order=1,datafile='../fs2000/fsdat_simul',mode_check,smoother,filter_covariance,filter_decomposition,forecast = 8,filtered_vars,filter_step_ahead=[1,3],irf=20,contemporaneous_correlation) m P c e W R k d y gy_obs;