diff --git a/doc/dynare.texi b/doc/dynare.texi
index 0fbc242b278e43fd66faf4653a975b6f19e5eb01..16528a3e75cca5d7d0242102d8fcb52dca6d3f42 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -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
diff --git a/matlab/pm3.m b/matlab/pm3.m
index 747389d5eee8d421cd7155bb67c3f9783d72b61f..7beaab1e266462023d86f76f1cdd3d075d825a5d 100644
--- a/matlab/pm3.m
+++ b/matlab/pm3.m
@@ -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
diff --git a/matlab/prior_posterior_statistics.m b/matlab/prior_posterior_statistics.m
index bfce7bcb3ac4269ba950d6da347e13278c7bfe44..37e00ff62d6acd04b1530a312e1e6a2af3aac68d 100644
--- a/matlab/prior_posterior_statistics.m
+++ b/matlab/prior_posterior_statistics.m
@@ -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;
-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;
+filter_covariance=0;
+if options_.filter_covariance
+ filter_covariance=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;
diff --git a/matlab/prior_posterior_statistics_core.m b/matlab/prior_posterior_statistics_core.m
index 267029e62e0c0606965554466a4ac87dd1312840..a7e34f7d2ddbcec60a913610a4bab8db784d1ba5 100644
--- a/matlab/prior_posterior_statistics_core.m
+++ b/matlab/prior_posterior_statistics_core.m
@@ -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);
diff --git a/tests/TeX/fs2000_corr_ME.mod b/tests/TeX/fs2000_corr_ME.mod
index a26447fb058b8480deac5ac4dad9072a7de024c4..9413e4f3b641d180085232b3ad66e3361459bd87 100644
--- a/tests/TeX/fs2000_corr_ME.mod
+++ b/tests/TeX/fs2000_corr_ME.mod
@@ -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;