325 lines
10 KiB
Matlab
325 lines
10 KiB
Matlab
function prior_posterior_statistics(type,dataset)
|
|
|
|
% function prior_posterior_statistics(type,dataset)
|
|
% Computes Monte Carlo filter smoother and forecasts
|
|
%
|
|
% INPUTS
|
|
% type: posterior
|
|
% prior
|
|
% gsa
|
|
% dataset: data structure
|
|
%
|
|
% OUTPUTS
|
|
% none
|
|
%
|
|
% SPECIAL REQUIREMENTS
|
|
% none
|
|
|
|
% PARALLEL CONTEXT
|
|
% See the comments random_walk_metropolis_hastings.m funtion.
|
|
|
|
|
|
% Copyright (C) 2005-2012 Dynare Team
|
|
%
|
|
% This file is part of Dynare.
|
|
%
|
|
% Dynare is free software: you can redistribute it and/or modify
|
|
% it under the terms of the GNU General Public License as published by
|
|
% the Free Software Foundation, either version 3 of the License, or
|
|
% (at your option) any later version.
|
|
%
|
|
% Dynare is distributed in the hope that it will be useful,
|
|
% but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
% GNU General Public License for more details.
|
|
%
|
|
% You should have received a copy of the GNU General Public License
|
|
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
global options_ estim_params_ oo_ M_ bayestopt_
|
|
|
|
localVars=[];
|
|
|
|
Y = dataset.data;
|
|
gend = dataset.info.ntobs;
|
|
data_index = dataset.missing.aindex;
|
|
missing_value = dataset.missing.state;
|
|
bayestopt_.mean_varobs = dataset.descriptive.mean';
|
|
|
|
nvx = estim_params_.nvx;
|
|
nvn = estim_params_.nvn;
|
|
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;
|
|
nvobs = size(options_.varobs,1);
|
|
iendo = 1:endo_nbr;
|
|
horizon = options_.forecast;
|
|
% moments_varendo = options_.moments_varendo;
|
|
filtered_vars = options_.filtered_vars;
|
|
if horizon
|
|
i_last_obs = gend+(1-M_.maximum_endo_lag:0);
|
|
end
|
|
maxlag = M_.maximum_endo_lag;
|
|
%%
|
|
if strcmpi(type,'posterior')
|
|
DirectoryName = CheckPath('metropolis',M_.dname);
|
|
load([ DirectoryName '/' M_.fname '_mh_history'])
|
|
FirstMhFile = record.KeepedDraws.FirstMhFile;
|
|
FirstLine = record.KeepedDraws.FirstLine;
|
|
TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); LastMhFile = TotalNumberOfMhFiles;
|
|
TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
|
|
NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
|
|
clear record;
|
|
if ~isempty(options_.sub_draws)
|
|
B = options_.sub_draws;
|
|
if B > NumberOfDraws
|
|
B = NumberOfDraws;
|
|
end
|
|
else
|
|
B = min(1200, round(0.25*NumberOfDraws));
|
|
end
|
|
elseif strcmpi(type,'gsa')
|
|
RootDirectoryName = CheckPath('gsa',M_.dname);
|
|
if options_.opt_gsa.pprior
|
|
DirectoryName = CheckPath(['gsa',filesep,'prior'],M_.dname);
|
|
load([ RootDirectoryName filesep M_.fname '_prior.mat'],'lpmat0','lpmat','istable')
|
|
else
|
|
DirectoryName = CheckPath(['gsa',filesep,'mc'],M_.dname);
|
|
load([ RootDirectoryName filesep M_.fname '_mc.mat'],'lpmat0','lpmat','istable')
|
|
end
|
|
x=[lpmat0(istable,:) lpmat(istable,:)];
|
|
clear lpmat istable
|
|
NumberOfDraws=size(x,1);
|
|
B=NumberOfDraws;
|
|
elseif strcmpi(type,'prior')
|
|
DirectoryName = CheckPath('prior',M_.dname);
|
|
if ~isempty(options_.subdraws)
|
|
B = options_.subdraws;
|
|
else
|
|
B = 1200;
|
|
end
|
|
end
|
|
%%
|
|
MAX_nruns = min(B,ceil(MaxNumberOfBytes/(npar+2)/8));
|
|
MAX_nsmoo = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*gend)/8));
|
|
MAX_ninno = min(B,ceil(MaxNumberOfBytes/(exo_nbr*gend)/8));
|
|
MAX_nerro = min(B,ceil(MaxNumberOfBytes/(size(options_.varobs,1)*gend)/8));
|
|
if naK
|
|
MAX_naK = min(B,ceil(MaxNumberOfBytes/(size(options_.varobs,1)* ...
|
|
length(options_.filter_step_ahead)*gend)/8));
|
|
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));
|
|
IdObs = bayestopt_.mfys;
|
|
|
|
end
|
|
MAX_momentsno = min(B,ceil(MaxNumberOfBytes/(get_moments_size(options_)*8)));
|
|
%%
|
|
varlist = options_.varlist;
|
|
if isempty(varlist)
|
|
varlist = M_.endo_names(1:M_.orig_endo_nbr, :);
|
|
end
|
|
nvar = size(varlist,1);
|
|
SelecVariables = [];
|
|
for i=1:nvar
|
|
if ~isempty(strmatch(varlist(i,:),M_.endo_names,'exact'))
|
|
SelecVariables = [SelecVariables;strmatch(varlist(i,:),M_.endo_names,'exact')];
|
|
end
|
|
end
|
|
|
|
irun = ones(7,1);
|
|
ifil = zeros(7,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);
|
|
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+maxlag,MAX_nforc1);
|
|
stock_forcst_point = zeros(endo_nbr,horizon+maxlag,MAX_nforc2);
|
|
run_smoother = 1;
|
|
end
|
|
%if moments_varendo
|
|
% stock_moments = cell(MAX_momentsno,1);
|
|
%end
|
|
|
|
|
|
|
|
% Store the variable mandatory for local/remote parallel computing.
|
|
|
|
localVars.type=type;
|
|
localVars.run_smoother=run_smoother;
|
|
localVars.gend=gend;
|
|
localVars.Y=Y;
|
|
localVars.data_index=data_index;
|
|
localVars.missing_value=missing_value;
|
|
localVars.varobs=options_.varobs;
|
|
localVars.irun=irun;
|
|
localVars.endo_nbr=endo_nbr;
|
|
localVars.nvn=nvn;
|
|
localVars.naK=naK;
|
|
localVars.horizon=horizon;
|
|
localVars.iendo=iendo;
|
|
if horizon
|
|
localVars.i_last_obs=i_last_obs;
|
|
localVars.IdObs=IdObs;
|
|
localVars.MAX_nforc1=MAX_nforc1;
|
|
localVars.MAX_nforc2=MAX_nforc2;
|
|
end
|
|
localVars.exo_nbr=exo_nbr;
|
|
localVars.maxlag=maxlag;
|
|
localVars.MAX_nsmoo=MAX_nsmoo;
|
|
localVars.MAX_ninno=MAX_ninno;
|
|
localVars.MAX_nerro = MAX_nerro;
|
|
if naK
|
|
localVars.MAX_naK=MAX_naK;
|
|
end
|
|
localVars.MAX_nruns=MAX_nruns;
|
|
localVars.MAX_momentsno = MAX_momentsno;
|
|
localVars.ifil=ifil;
|
|
localVars.DirectoryName = DirectoryName;
|
|
|
|
if strcmpi(type,'posterior'),
|
|
b=0;
|
|
while b<=B
|
|
b = b + 1;
|
|
[x(b,:), logpost(b,1)] = GetOneDraw(type);
|
|
end
|
|
localVars.logpost=logpost;
|
|
end
|
|
|
|
if ~strcmpi(type,'prior'),
|
|
localVars.x=x;
|
|
end
|
|
|
|
% Like sequential execution!
|
|
if isnumeric(options_.parallel),
|
|
[fout] = prior_posterior_statistics_core(localVars,1,B,0);
|
|
% Parallel execution!
|
|
else
|
|
[nCPU, totCPU, nBlockPerCPU] = distributeJobs(options_.parallel, 1, B);
|
|
ifil=zeros(7,totCPU);
|
|
for j=1:totCPU-1,
|
|
if run_smoother
|
|
nfiles = ceil(nBlockPerCPU(j)/MAX_nsmoo);
|
|
ifil(1,j+1) =ifil(1,j)+nfiles;
|
|
nfiles = ceil(nBlockPerCPU(j)/MAX_ninno);
|
|
ifil(2,j+1) =ifil(2,j)+nfiles;
|
|
nfiles = ceil(nBlockPerCPU(j)/MAX_nerro);
|
|
ifil(3,j+1) =ifil(3,j)+nfiles;
|
|
end
|
|
if naK
|
|
nfiles = ceil(nBlockPerCPU(j)/MAX_naK);
|
|
ifil(4,j+1) =ifil(4,j)+nfiles;
|
|
end
|
|
nfiles = ceil(nBlockPerCPU(j)/MAX_nruns);
|
|
ifil(5,j+1) =ifil(5,j)+nfiles;
|
|
if horizon
|
|
nfiles = ceil(nBlockPerCPU(j)/MAX_nforc1);
|
|
ifil(6,j+1) =ifil(6,j)+nfiles;
|
|
nfiles = ceil(nBlockPerCPU(j)/MAX_nforc2);
|
|
ifil(7,j+1) =ifil(7,j)+nfiles;
|
|
end
|
|
% nfiles = ceil(nBlockPerCPU(j)/MAX_momentsno);
|
|
% ifil(8,j+1) =ifil(8,j)+nfiles;
|
|
end
|
|
localVars.ifil = ifil;
|
|
globalVars = struct('M_',M_, ...
|
|
'options_', options_, ...
|
|
'bayestopt_', bayestopt_, ...
|
|
'estim_params_', estim_params_, ...
|
|
'oo_', oo_);
|
|
|
|
% which files have to be copied to run remotely
|
|
NamFileInput(1,:) = {'',[M_.fname '_static.m']};
|
|
NamFileInput(2,:) = {'',[M_.fname '_dynamic.m']};
|
|
if options_.steadystate_flag,
|
|
NamFileInput(length(NamFileInput)+1,:)={'',[M_.fname '_steadystate.m']};
|
|
end
|
|
[fout] = masterParallel(options_.parallel, 1, B,NamFileInput,'prior_posterior_statistics_core', localVars,globalVars, options_.parallel_info);
|
|
|
|
end
|
|
ifil = fout(end).ifil;
|
|
|
|
|
|
stock_gend=gend;
|
|
stock_data=Y;
|
|
save([DirectoryName '/' M_.fname '_data.mat'],'stock_gend','stock_data');
|
|
|
|
if strcmpi(type,'gsa')
|
|
return
|
|
end
|
|
|
|
if ~isnumeric(options_.parallel),
|
|
leaveSlaveOpen = options_.parallel_info.leaveSlaveOpen;
|
|
if options_.parallel_info.leaveSlaveOpen == 0,
|
|
% Commenting for testing!!!
|
|
options_.parallel_info.leaveSlaveOpen = 1; % Force locally to leave open remote matlab sessions (repeated pm3 calls)
|
|
end
|
|
end
|
|
|
|
if options_.smoother
|
|
pm3(endo_nbr,gend,ifil(1),B,'Smoothed variables',...
|
|
'',M_.endo_names(1:M_.orig_endo_nbr, :),M_.endo_names_tex,M_.endo_names,...
|
|
varlist,'SmoothedVariables',DirectoryName,'_smooth');
|
|
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');
|
|
if nvn
|
|
% needs to be fixed
|
|
% pm3(endo_nbr,gend,ifil(3),B,'Smoothed measurement errors',...
|
|
% M_.endo_names(SelecVariables),M_.endo_names,'tit_tex',M_.endo_names,...
|
|
% 'names2','smooth_errors',[M_.fname '/metropolis'],'_error')
|
|
end
|
|
end
|
|
|
|
if options_.filtered_vars
|
|
pm3(endo_nbr,gend,ifil(1),B,'Updated Variables',...
|
|
'',varlist,'tit_tex',M_.endo_names,...
|
|
varlist,'UpdatedVariables',DirectoryName, ...
|
|
'_update');
|
|
pm3(endo_nbr,gend+1,ifil(4),B,'One step ahead forecast',...
|
|
'',varlist,'tit_tex',M_.endo_names,...
|
|
varlist,'FilteredVariables',DirectoryName,'_filter_step_ahead');
|
|
end
|
|
|
|
if options_.forecast
|
|
pm3(endo_nbr,horizon+maxlag,ifil(6),B,'Forecasted variables (mean)',...
|
|
'',varlist,'tit_tex',M_.endo_names,...
|
|
varlist,'MeanForecast',DirectoryName,'_forc_mean');
|
|
pm3(endo_nbr,horizon+maxlag,ifil(6),B,'Forecasted variables (point)',...
|
|
'',varlist,'tit_tex',M_.endo_names,...
|
|
varlist,'PointForecast',DirectoryName,'_forc_point');
|
|
end
|
|
|
|
|
|
if ~isnumeric(options_.parallel),
|
|
options_.parallel_info.leaveSlaveOpen = leaveSlaveOpen;
|
|
if leaveSlaveOpen == 0,
|
|
closeSlave(options_.parallel,options_.parallel_info.RemoteTmpFolder),
|
|
end
|
|
end
|