diff --git a/matlab/pm3.m b/matlab/pm3.m
index 7beaab1e266462023d86f76f1cdd3d075d825a5d..7a185a04270754fac012eec38f24463b163baa81 100644
--- a/matlab/pm3.m
+++ b/matlab/pm3.m
@@ -83,7 +83,6 @@ if options_.estimation.moments_posterior_density.indicator
Density = zeros(options_.estimation.moments_posterior_density.gridpoints,2,n2,nvar);
end
fprintf(['Estimation::mcmc: ' tit1 '\n']);
-stock1 = zeros(n1,n2,B);
k = 0;
filter_step_ahead_indicator=0;
filter_covar_indicator=0;
@@ -93,6 +92,7 @@ for file = 1:ifil
if strcmp(var_type,'_filter_step_ahead')
if file==1 %on first run, initialize variable for storing filter_step_ahead
stock1_filter_step_ahead=NaN(n1,n2,B,length(options_.filter_step_ahead));
+ stock1 = zeros(n1,n2,B);
end
filter_step_ahead_indicator=1;
stock_filter_step_ahead=zeros(n1,n2,size(stock,4),length(options_.filter_step_ahead));
@@ -101,23 +101,29 @@ for file = 1:ifil
stock_filter_step_ahead(:,:,:,ii)=stock(ii,:,1+K_step_ahead:n2+K_step_ahead,:);
end
stock = squeeze(stock(1,:,1+1:1+n2,:)); %1 step ahead starts at entry 2
- end
- if strcmp(var_type,'_filter_covar')
+ k = k(end)+(1:size(stock,3));
+ stock1(:,:,k) = stock;
+ stock1_filter_step_ahead(:,:,k,:) = stock_filter_step_ahead;
+ elseif 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);
+ 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;
+ elseif strcmp(var_type,'_trend_coeff')
+ if file==1 %on first run, initialize variable for storing filter_step_ahead
+ stock1_filter_step_ahead=NaN(n1,n2,B,length(options_.filter_step_ahead));
+ stock1 = zeros(n1,B);
+ end
+ k = k(end)+(1:size(stock,2));
+ stock1(:,k) = stock;
+ else
+ if file==1 %on first run, initialize variable for storing filter_step_ahead
+ stock1 = zeros(n1,n2,B);
+ end
+ k = k(end)+(1:size(stock,3));
+ stock1(:,:,k) = stock;
end
end
clear stock
@@ -157,30 +163,43 @@ if filter_covar_indicator
return
end
-for i = 1:nvar
- for j = 1:n2
+if strcmp(var_type,'_trend_coeff') %two dimensional arrays
+ for i = 1:nvar
if options_.estimation.moments_posterior_density.indicator
- [Mean(j,i),Median(j,i),Var(j,i),HPD(:,j,i),Distrib(:,j,i),Density(:,:,j,i)] = ...
- posterior_moments(squeeze(stock1(SelecVariables(i),j,:)),1,options_.mh_conf_sig,options_.estimation.moments_posterior_density);
+ [Mean(1,i),Median(1,i),Var(1,i),HPD(:,1,i),Distrib(:,1,i),Density(:,:,1,i)] = ...
+ posterior_moments(squeeze(stock1(SelecVariables(i),:)),1,options_.mh_conf_sig,options_.estimation.moments_posterior_density);
else
- [Mean(j,i),Median(j,i),Var(j,i),HPD(:,j,i),Distrib(:,j,i)] = ...
- posterior_moments(squeeze(stock1(SelecVariables(i),j,:)),0,options_.mh_conf_sig);
+ [Mean(1,i),Median(1,i),Var(1,i),HPD(:,1,i),Distrib(:,1,i)] = ...
+ posterior_moments(squeeze(stock1(SelecVariables(i),:)),0,options_.mh_conf_sig);
end
- if filter_step_ahead_indicator
+ end
+else %three dimensional arrays
+ for i = 1:nvar
+ for j = 1:n2
if options_.estimation.moments_posterior_density.indicator
- for K_step = 1:length(options_.filter_step_ahead)
- [Mean_filter_step_ahead(K_step,i,j),Median_filter_step_ahead(K_step,i,j),Var_filter_step_ahead(K_step,i,j),HPD_filter_step_ahead(:,K_step,i,j),Distrib_filter_step_ahead(:,K_step,i,j),Density_filter_step_ahead(:,:,K_step,i,j) ] = ...
- posterior_moments(squeeze(stock1_filter_step_ahead(SelecVariables(i),j,:,K_step)),1,options_.mh_conf_sig,options_.estimation.moments_posterior_density);
- end
+ [Mean(j,i),Median(j,i),Var(j,i),HPD(:,j,i),Distrib(:,j,i),Density(:,:,j,i)] = ...
+ posterior_moments(squeeze(stock1(SelecVariables(i),j,:)),1,options_.mh_conf_sig,options_.estimation.moments_posterior_density);
else
- for K_step = 1:length(options_.filter_step_ahead)
- [Mean_filter_step_ahead(K_step,i,j),Median_filter_step_ahead(K_step,i,j),Var_filter_step_ahead(K_step,i,j),HPD_filter_step_ahead(:,K_step,i,j),Distrib_filter_step_ahead(:,K_step,i,j)] = ...
- posterior_moments(squeeze(stock1_filter_step_ahead(SelecVariables(i),j,:,K_step)),0,options_.mh_conf_sig);
+ [Mean(j,i),Median(j,i),Var(j,i),HPD(:,j,i),Distrib(:,j,i)] = ...
+ posterior_moments(squeeze(stock1(SelecVariables(i),j,:)),0,options_.mh_conf_sig);
+ end
+ if filter_step_ahead_indicator
+ if options_.estimation.moments_posterior_density.indicator
+ for K_step = 1:length(options_.filter_step_ahead)
+ [Mean_filter_step_ahead(K_step,i,j),Median_filter_step_ahead(K_step,i,j),Var_filter_step_ahead(K_step,i,j),HPD_filter_step_ahead(:,K_step,i,j),Distrib_filter_step_ahead(:,K_step,i,j),Density_filter_step_ahead(:,:,K_step,i,j) ] = ...
+ posterior_moments(squeeze(stock1_filter_step_ahead(SelecVariables(i),j,:,K_step)),1,options_.mh_conf_sig,options_.estimation.moments_posterior_density);
+ end
+ else
+ for K_step = 1:length(options_.filter_step_ahead)
+ [Mean_filter_step_ahead(K_step,i,j),Median_filter_step_ahead(K_step,i,j),Var_filter_step_ahead(K_step,i,j),HPD_filter_step_ahead(:,K_step,i,j),Distrib_filter_step_ahead(:,K_step,i,j)] = ...
+ posterior_moments(squeeze(stock1_filter_step_ahead(SelecVariables(i),j,:,K_step)),0,options_.mh_conf_sig);
+ end
end
end
end
end
end
+
clear stock1
if filter_step_ahead_indicator %write matrices corresponding to ML
clear stock1_filter_step_ahead
@@ -194,31 +213,51 @@ if filter_step_ahead_indicator %write matrices corresponding to ML
oo_.FilteredVariablesKStepAheadVariances=FilteredVariablesKStepAheadVariances;
end
-for i = 1:nvar
- name = deblank(names1(SelecVariables(i),:));
- oo_.(name3).Mean.(name) = Mean(:,i);
- oo_.(name3).Median.(name) = Median(:,i);
- oo_.(name3).Var.(name) = Var(:,i);
- oo_.(name3).deciles.(name) = Distrib(:,:,i);
- oo_.(name3).HPDinf.(name) = HPD(1,:,i)';
- oo_.(name3).HPDsup.(name) = HPD(2,:,i)';
- if options_.estimation.moments_posterior_density.indicator
- oo_.(name3).density.(name) = Density(:,:,:,i);
+if strcmp(var_type,'_trend_coeff') || strcmp(var_type,'_smoothed_trend') || strcmp(var_type,'_smoothed_trend')
+ for i = 1:nvar
+ name = deblank(names1(SelecVariables(i),:));
+ oo_.Smoother.(name3).Mean.(name) = Mean(:,i);
+ oo_.Smoother.(name3).Median.(name) = Median(:,i);
+ oo_.Smoother.(name3).Var.(name) = Var(:,i);
+ oo_.Smoother.(name3).deciles.(name) = Distrib(:,:,i);
+ oo_.Smoother.(name3).HPDinf.(name) = HPD(1,:,i)';
+ oo_.Smoother.(name3).HPDsup.(name) = HPD(2,:,i)';
+ if options_.estimation.moments_posterior_density.indicator
+ oo_.Smoother.(name3).density.(name) = Density(:,:,:,i);
+ end
end
- if filter_step_ahead_indicator
- for K_step = 1:length(options_.filter_step_ahead)
- name4=['Filtered_Variables_',num2str(K_step),'_step_ahead'];
- oo_.(name4).Mean.(name) = squeeze(Mean_filter_step_ahead(K_step,i,:));
- oo_.(name4).Median.(name) = squeeze(Median_filter_step_ahead(K_step,i,:));
- oo_.(name4).Var.(name) = squeeze(Var_filter_step_ahead(K_step,i,:));
- oo_.(name4).deciles.(name) = squeeze(Distrib_filter_step_ahead(:,K_step,i,:));
- oo_.(name4).HPDinf.(name) = squeeze(HPD_filter_step_ahead(1,K_step,i,:));
- oo_.(name4).HPDsup.(name) = squeeze(HPD_filter_step_ahead(2,K_step,i,:));
- if options_.estimation.moments_posterior_density.indicator
- oo_.(name4).density.(name) = squeeze(Density_filter_step_ahead(:,:,K_step,i,:));
+else
+ for i = 1:nvar
+ name = deblank(names1(SelecVariables(i),:));
+ oo_.(name3).Mean.(name) = Mean(:,i);
+ oo_.(name3).Median.(name) = Median(:,i);
+ oo_.(name3).Var.(name) = Var(:,i);
+ oo_.(name3).deciles.(name) = Distrib(:,:,i);
+ oo_.(name3).HPDinf.(name) = HPD(1,:,i)';
+ oo_.(name3).HPDsup.(name) = HPD(2,:,i)';
+ if options_.estimation.moments_posterior_density.indicator
+ oo_.(name3).density.(name) = Density(:,:,:,i);
+ end
+ if filter_step_ahead_indicator
+ for K_step = 1:length(options_.filter_step_ahead)
+ name4=['Filtered_Variables_',num2str(K_step),'_step_ahead'];
+ oo_.(name4).Mean.(name) = squeeze(Mean_filter_step_ahead(K_step,i,:));
+ oo_.(name4).Median.(name) = squeeze(Median_filter_step_ahead(K_step,i,:));
+ oo_.(name4).Var.(name) = squeeze(Var_filter_step_ahead(K_step,i,:));
+ oo_.(name4).deciles.(name) = squeeze(Distrib_filter_step_ahead(:,K_step,i,:));
+ oo_.(name4).HPDinf.(name) = squeeze(HPD_filter_step_ahead(1,K_step,i,:));
+ oo_.(name4).HPDsup.(name) = squeeze(HPD_filter_step_ahead(2,K_step,i,:));
+ if options_.estimation.moments_posterior_density.indicator
+ oo_.(name4).density.(name) = squeeze(Density_filter_step_ahead(:,:,K_step,i,:));
+ end
end
end
- end
+ end
+end
+
+if strcmp(var_type,'_trend_coeff') || max(max(abs(Mean(:,:))))<=10^(-6) || all(all(isnan(Mean)))
+ fprintf(['Estimation::mcmc: ' tit1 ', done!\n']);
+ return %not do plots
end
%%
%% Finally I build the plots.
diff --git a/matlab/prior_posterior_statistics.m b/matlab/prior_posterior_statistics.m
index 37e00ff62d6acd04b1530a312e1e6a2af3aac68d..31e2bb35962868d7598f70a4ca6c946be17440df 100644
--- a/matlab/prior_posterior_statistics.m
+++ b/matlab/prior_posterior_statistics.m
@@ -93,6 +93,9 @@ end
MAX_nruns = min(B,ceil(MaxNumberOfBytes/(npar+2)/8));
MAX_nsmoo = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*gend)/8));
+MAX_n_smoothed_constant = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*gend)/8));
+MAX_n_smoothed_trend = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*gend)/8));
+MAX_n_trend_coeff = min(B,ceil(MaxNumberOfBytes/endo_nbr/8));
MAX_ninno = min(B,ceil(MaxNumberOfBytes/(exo_nbr*gend)/8));
MAX_nerro = min(B,ceil(MaxNumberOfBytes/(size(options_.varobs,1)*gend)/8));
@@ -123,7 +126,7 @@ for i=1:nvar
end
end
-n_variables_to_fill=8;
+n_variables_to_fill=11;
irun = ones(n_variables_to_fill,1);
ifil = zeros(n_variables_to_fill,1);
@@ -172,6 +175,9 @@ end
if options_.filter_covariance
localVars.MAX_filter_covariance = MAX_filter_covariance;
end
+localVars.MAX_n_smoothed_constant=MAX_n_smoothed_constant;
+localVars.MAX_n_smoothed_trend=MAX_n_smoothed_trend;
+localVars.MAX_n_trend_coeff=MAX_n_trend_coeff;
localVars.MAX_nruns=MAX_nruns;
localVars.MAX_momentsno = MAX_momentsno;
localVars.ifil=ifil;
@@ -223,6 +229,14 @@ else
nfiles = ceil(nBlockPerCPU(j)/MAX_filter_covariance);
ifil(8,j+1) =ifil(8,j)+nfiles;
end
+ if run_smoother
+ nfiles = ceil(nBlockPerCPU(j)/MAX_n_trend_coeff);
+ ifil(9,j+1) =ifil(9,j)+nfiles;
+ nfiles = ceil(nBlockPerCPU(j)/MAX_n_smoothed_constant);
+ ifil(10,j+1) =ifil(10,j)+nfiles;
+ nfiles = ceil(nBlockPerCPU(j)/MAX_n_smoothed_trend);
+ ifil(11,j+1) =ifil(11,j)+nfiles;
+ end
end
localVars.ifil = ifil;
globalVars = struct('M_',M_, ...
@@ -264,6 +278,16 @@ if options_.smoother
pm3(exo_nbr,gend,ifil(2),B,'Smoothed shocks',...
'',M_.exo_names,M_.exo_names_tex,M_.exo_names,...
M_.exo_names,'SmoothedShocks',DirectoryName,'_inno');
+ pm3(endo_nbr,1,ifil(9),B,'Trend_coefficients',...
+ '',M_.endo_names(1:M_.orig_endo_nbr,:),M_.endo_names_tex,M_.endo_names,...
+ varlist,'TrendCoeff',DirectoryName,'_trend_coeff');
+ pm3(endo_nbr,gend,ifil(10),B,'Smoothed constant',...
+ '',M_.endo_names(1:M_.orig_endo_nbr, :),M_.endo_names_tex,M_.endo_names,...
+ varlist,'Constant',DirectoryName,'_smoothed_constant');
+ pm3(endo_nbr,gend,ifil(11),B,'Smoothed trend',...
+ '',M_.endo_names(1:M_.orig_endo_nbr, :),M_.endo_names_tex,M_.endo_names,...
+ varlist,'Trend',DirectoryName,'_smoothed_trend');
+
if nvn
for obs_iter=1:length(options_.varobs)
meas_error_names{obs_iter,1}=['SE_EOBS_' M_.endo_names(strmatch(options_.varobs{obs_iter},M_.endo_names,'exact'),:)];
diff --git a/matlab/prior_posterior_statistics_core.m b/matlab/prior_posterior_statistics_core.m
index a7e34f7d2ddbcec60a913610a4bab8db784d1ba5..8ecf334096ad934962d40b049014655b8410ae03 100644
--- a/matlab/prior_posterior_statistics_core.m
+++ b/matlab/prior_posterior_statistics_core.m
@@ -15,7 +15,11 @@ function myoutput=prior_posterior_statistics_core(myinputs,fpar,B,whoiam, ThisMa
% _filter_step_ahead;
% _param;
% _forc_mean;
-% _forc_point
+% _forc_point;
+% _filter_covar;
+% _trend_coeff;
+% _smoothed_trend;
+% _smoothed_constant;
%
% ALGORITHM
% Portion of prior_posterior.m function.
@@ -77,11 +81,13 @@ end
if filter_covariance
MAX_filter_covariance=myinputs.MAX_filter_covariance;
end
-
exo_nbr=myinputs.exo_nbr;
maxlag=myinputs.maxlag;
MAX_nsmoo=myinputs.MAX_nsmoo;
MAX_ninno=myinputs.MAX_ninno;
+MAX_n_smoothed_constant=myinputs.MAX_n_smoothed_constant;
+MAX_n_smoothed_trend=myinputs.MAX_n_smoothed_trend;
+MAX_n_trend_coeff=myinputs.MAX_n_trend_coeff;
MAX_nerro = myinputs.MAX_nerro;
MAX_nruns=myinputs.MAX_nruns;
MAX_momentsno = myinputs.MAX_momentsno;
@@ -132,6 +138,9 @@ if RemoteFlag==1,
OutputFileName_forc_mean = {};
OutputFileName_forc_point = {};
OutputFileName_filter_covar ={};
+ OutputFileName_trend_coeff = {};
+ OutputFileName_smoothed_trend = {};
+ OutputFileName_smoothed_constant = {};
% OutputFileName_moments = {};
end
@@ -140,6 +149,9 @@ if run_smoother
stock_smooth=NaN(endo_nbr,gend,MAX_nsmoo);
stock_update=NaN(endo_nbr,gend,MAX_nsmoo);
stock_innov=NaN(M_.exo_nbr,gend,MAX_ninno);
+ stock_smoothed_constant=NaN(endo_nbr,gend,MAX_n_smoothed_constant);
+ stock_smoothed_trend=NaN(endo_nbr,gend,MAX_n_smoothed_trend);
+ stock_trend_coeff = zeros(endo_nbr,MAX_n_trend_coeff);
if horizon
stock_forcst_mean= NaN(endo_nbr,horizon,MAX_nforc1);
stock_forcst_point = NaN(endo_nbr,horizon,MAX_nforc2);
@@ -177,18 +189,23 @@ for b=fpar:B
[dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
[alphahat,etahat,epsilonhat,alphatilde,SteadyState,trend_coeff,aK,junk1,junk2,P,junk4,junk5,trend_addition] = ...
DsgeSmoother(deep,gend,Y,data_index,missing_value);
-
+
+ stock_trend_coeff(options_.varobs_id,irun(9))=trend_coeff;
+ stock_smoothed_trend(IdObs,:,irun(11))=trend_addition;
if options_.loglinear %reads values from smoother results, which are in dr-order and put them into declaration order
+ constant_part=repmat(log(SteadyState(dr.order_var)),1,gend);
stock_smooth(dr.order_var,:,irun(1)) = alphahat(1:endo_nbr,:)+ ...
- repmat(log(SteadyState(dr.order_var)),1,gend);
+ constant_part;
stock_update(dr.order_var,:,irun(1)) = alphatilde(1:endo_nbr,:)+ ...
- repmat(log(SteadyState(dr.order_var)),1,gend);
+ constant_part;
else
+ constant_part=repmat(SteadyState(dr.order_var),1,gend);
stock_smooth(dr.order_var,:,irun(1)) = alphahat(1:endo_nbr,:)+ ...
- repmat(SteadyState(dr.order_var),1,gend);
+ constant_part;
stock_update(dr.order_var,:,irun(1)) = alphatilde(1:endo_nbr,:)+ ...
- repmat(SteadyState(dr.order_var),1,gend);
- end
+ constant_part;
+ end
+ stock_smoothed_constant(dr.order_var,:,irun(10))=constant_part;
%% Compute constant for observables
if options_.prefilter == 1 %as mean is taken after log transformation, no distinction is needed here
constant_part=repmat(mean_varobs',1,gend);
@@ -205,6 +222,7 @@ for b=fpar:B
else
mean_correction=-repmat(SteadyState(IdObs),1,gend)+constant_part;
end
+ stock_smoothed_constant(IdObs,:,irun(10))=stock_smoothed_constant(IdObs,:,irun(10))+mean_correction;
%smoothed variables are E_T(y_t) so no trend shift is required
stock_smooth(IdObs,:,irun(1))=stock_smooth(IdObs,:,irun(1))+trend_addition+mean_correction;
%updated variables are E_t(y_t) so no trend shift is required
@@ -282,7 +300,7 @@ for b=fpar:B
stock_ys(irun(5),:) = SteadyState';
- irun = irun + ones(8,1);
+ irun = irun + ones(11,1);
if run_smoother && (irun(1) > MAX_nsmoo || b == B),
@@ -370,6 +388,39 @@ for b=fpar:B
end
irun(8) = 1;
end
+
+ irun_index=9;
+ if run_smoother && (irun(irun_index) > MAX_n_trend_coeff || b == B)
+ stock = stock_trend_coeff(:,1:irun(irun_index)-1);
+ ifil(irun_index) = ifil(irun_index) + 1;
+ save([DirectoryName '/' M_.fname '_trend_coeff' int2str(ifil(irun_index)) '.mat'],'stock');
+ if RemoteFlag==1,
+ OutputFileName_trend_coeff = [OutputFileName_trend_coeff; {[DirectoryName filesep], [M_.fname '_trend_coeff' int2str(ifil(irun_index)) '.mat']}];
+ end
+ irun(irun_index) = 1;
+ end
+
+ irun_index=10;
+ if run_smoother && (irun(irun_index) > MAX_n_smoothed_constant || b == B)
+ stock = stock_smoothed_constant(:,:,1:irun(irun_index)-1);
+ ifil(irun_index) = ifil(irun_index) + 1;
+ save([DirectoryName '/' M_.fname '_smoothed_constant' int2str(ifil(irun_index)) '.mat'],'stock');
+ if RemoteFlag==1,
+ OutputFileName_smoothed_constant = [OutputFileName_smoothed_constant; {[DirectoryName filesep], [M_.fname '_smoothed_constant' int2str(ifil(irun_index)) '.mat']}];
+ end
+ irun(irun_index) = 1;
+ end
+
+ irun_index=11;
+ if run_smoother && (irun(irun_index) > MAX_n_smoothed_trend || b == B)
+ stock = stock_smoothed_trend(:,:,1:irun(irun_index)-1);
+ ifil(irun_index) = ifil(irun_index) + 1;
+ save([DirectoryName '/' M_.fname '_smoothed_trend' int2str(ifil(irun_index)) '.mat'],'stock');
+ if RemoteFlag==1,
+ OutputFileName_smoothed_trend = [OutputFileName_smoothed_trend; {[DirectoryName filesep], [M_.fname '_smoothed_trend' int2str(ifil(irun_index)) '.mat']}];
+ end
+ irun(irun_index) = 1;
+ end
dyn_waitbar((b-fpar+1)/(B-fpar+1),h);
end
@@ -384,7 +435,10 @@ if RemoteFlag==1,
OutputFileName_param;
OutputFileName_forc_mean;
OutputFileName_forc_point
- OutputFileName_filter_covar];
+ OutputFileName_filter_covar
+ OutputFileName_trend_coeff
+ OutputFileName_smoothed_trend
+ OutputFileName_smoothed_constant];
end
dyn_waitbar_close(h);