From 3b09ab5424676631fb094d7f2479d6f18fc8a8b5 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Fri, 3 Jun 2016 17:09:09 +0200 Subject: [PATCH] Save Smoother.Trend and Smoother.Constant after MCMC --- matlab/pm3.m | 139 +++++++++++++++-------- matlab/prior_posterior_statistics.m | 26 ++++- matlab/prior_posterior_statistics_core.m | 74 ++++++++++-- 3 files changed, 178 insertions(+), 61 deletions(-) diff --git a/matlab/pm3.m b/matlab/pm3.m index 7beaab1e2..7a185a042 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') - 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. diff --git a/matlab/prior_posterior_statistics.m b/matlab/prior_posterior_statistics.m index 37e00ff62..31e2bb359 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 a7e34f7d2..8ecf33409 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);