Save Smoother.Trend and Smoother.Constant after MCMC

time-shift
Johannes Pfeifer 2016-06-03 17:09:09 +02:00 committed by Stéphane Adjemian (Lupi)
parent 1c84ccdc34
commit 3b09ab5424
3 changed files with 178 additions and 61 deletions

View File

@ -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')
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
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);
end
filter_covar_indicator=1;
k = k(end)+(1:size(stock,4));
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.

View File

@ -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'),:)];

View File

@ -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);