prior_posterior_statistics

 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

0001 function prior_posterior_statistics(type,dataset)
0002 
0003 % function prior_posterior_statistics(type,dataset)
0004 % Computes Monte Carlo filter smoother and forecasts
0005 %
0006 % INPUTS
0007 %    type:         posterior
0008 %                  prior
0009 %                  gsa
0010 %    dataset:      data structure
0011 %
0012 % OUTPUTS
0013 %    none
0014 %
0015 % SPECIAL REQUIREMENTS
0016 %    none
0017 
0018 % PARALLEL CONTEXT
0019 % See the comments random_walk_metropolis_hastings.m funtion.
0020 
0021 
0022 % Copyright (C) 2005-2011 Dynare Team
0023 %
0024 % This file is part of Dynare.
0025 %
0026 % Dynare is free software: you can redistribute it and/or modify
0027 % it under the terms of the GNU General Public License as published by
0028 % the Free Software Foundation, either version 3 of the License, or
0029 % (at your option) any later version.
0030 %
0031 % Dynare is distributed in the hope that it will be useful,
0032 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0033 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0034 % GNU General Public License for more details.
0035 %
0036 % You should have received a copy of the GNU General Public License
0037 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0038 
0039 global options_ estim_params_ oo_ M_ bayestopt_
0040 
0041 localVars=[];
0042 
0043 Y = dataset.data;
0044 gend = dataset.info.ntobs;
0045 data_index = dataset.missing.aindex;
0046 missing_value = dataset.missing.state;
0047 bayestopt_.mean_varobs = dataset.descriptive.mean';
0048 
0049 nvx  = estim_params_.nvx;
0050 nvn  = estim_params_.nvn;
0051 ncx  = estim_params_.ncx;
0052 ncn  = estim_params_.ncn;
0053 np   = estim_params_.np ;
0054 npar = nvx+nvn+ncx+ncn+np;
0055 offset = npar-np;
0056 naK = length(options_.filter_step_ahead);
0057 %%
0058 MaxNumberOfBytes=options_.MaxNumberOfBytes;
0059 endo_nbr=M_.endo_nbr;
0060 exo_nbr=M_.exo_nbr;
0061 nvobs     = size(options_.varobs,1);
0062 iendo = 1:endo_nbr;
0063 horizon = options_.forecast;
0064 % moments_varendo = options_.moments_varendo;
0065 filtered_vars = options_.filtered_vars;
0066 if horizon
0067     i_last_obs = gend+(1-M_.maximum_endo_lag:0);
0068 end
0069 maxlag = M_.maximum_endo_lag;
0070 %%
0071 if strcmpi(type,'posterior')
0072     DirectoryName = CheckPath('metropolis',M_.dname);
0073     load([ DirectoryName '/'  M_.fname '_mh_history'])
0074     FirstMhFile = record.KeepedDraws.FirstMhFile;
0075     FirstLine = record.KeepedDraws.FirstLine;
0076     TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); LastMhFile = TotalNumberOfMhFiles;
0077     TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
0078     NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
0079     clear record;
0080     if ~isempty(options_.sub_draws)
0081         B = options_.sub_draws;
0082         if B > NumberOfDraws
0083             B = NumberOfDraws;
0084         end
0085     else
0086         B = min(1200, round(0.25*NumberOfDraws));
0087     end
0088 elseif strcmpi(type,'gsa')
0089     RootDirectoryName = CheckPath('gsa',M_.dname);
0090     if options_.opt_gsa.pprior
0091         DirectoryName = CheckPath(['gsa',filesep,'prior'],M_.dname);
0092         load([ RootDirectoryName filesep  M_.fname '_prior.mat'],'lpmat0','lpmat','istable')
0093     else
0094         DirectoryName = CheckPath(['gsa',filesep,'mc'],M_.dname);
0095         load([ RootDirectoryName filesep  M_.fname '_mc.mat'],'lpmat0','lpmat','istable')
0096     end
0097     x=[lpmat0(istable,:) lpmat(istable,:)];
0098     clear lpmat istable
0099     NumberOfDraws=size(x,1);
0100     B=NumberOfDraws; 
0101 elseif strcmpi(type,'prior')
0102     DirectoryName = CheckPath('prior',M_.dname);
0103     if ~isempty(options_.subdraws)
0104         B = options_.subdraws;
0105     else
0106         B = 1200;
0107     end
0108 end
0109 %%
0110 MAX_nruns = min(B,ceil(MaxNumberOfBytes/(npar+2)/8));
0111 MAX_nsmoo = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*gend)/8));
0112 MAX_ninno = min(B,ceil(MaxNumberOfBytes/(exo_nbr*gend)/8));
0113 MAX_nerro = min(B,ceil(MaxNumberOfBytes/(size(options_.varobs,1)*gend)/8));
0114 if naK
0115     MAX_naK   = min(B,ceil(MaxNumberOfBytes/(size(options_.varobs,1)* ...
0116                                              length(options_.filter_step_ahead)*gend)/8));
0117 end
0118 if horizon
0119     MAX_nforc1 = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*(horizon+maxlag))/8));
0120     MAX_nforc2 = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*(horizon+maxlag))/ ...
0121                             8));
0122     IdObs    = bayestopt_.mfys;
0123     
0124 end
0125 MAX_momentsno = min(B,ceil(MaxNumberOfBytes/(get_moments_size(options_)*8)));
0126 %%
0127 varlist = options_.varlist;
0128 if isempty(varlist)
0129     varlist = M_.endo_names(1:M_.orig_endo_nbr, :);
0130 end
0131 nvar = size(varlist,1);
0132 SelecVariables = [];
0133 for i=1:nvar
0134     if ~isempty(strmatch(varlist(i,:),M_.endo_names,'exact'))
0135         SelecVariables = [SelecVariables;strmatch(varlist(i,:),M_.endo_names,'exact')];
0136     end
0137 end
0138 
0139 irun = ones(7,1);
0140 ifil = zeros(7,1);
0141 
0142 
0143 stock_param = zeros(MAX_nruns, npar);
0144 stock_logpo = zeros(MAX_nruns,1);
0145 stock_ys = zeros(MAX_nruns,endo_nbr);
0146 run_smoother = 0;
0147 if options_.smoother
0148     stock_smooth = zeros(endo_nbr,gend,MAX_nsmoo);
0149     stock_innov  = zeros(exo_nbr,gend,B);
0150     stock_error = zeros(nvobs,gend,MAX_nerro);
0151     stock_update = zeros(endo_nbr,gend,MAX_nsmoo);
0152     run_smoother = 1;
0153 end
0154 
0155 if options_.filter_step_ahead
0156     stock_filter_step_ahead = zeros(naK,endo_nbr,gend+ ...
0157                                     options_.filter_step_ahead(end),MAX_naK);
0158     run_smoother = 1;
0159 end
0160 if options_.forecast
0161     stock_forcst_mean = zeros(endo_nbr,horizon+maxlag,MAX_nforc1);
0162     stock_forcst_point = zeros(endo_nbr,horizon+maxlag,MAX_nforc2);
0163     run_smoother = 1;
0164 end
0165 %if moments_varendo
0166 %    stock_moments = cell(MAX_momentsno,1);
0167 %end
0168 
0169 
0170 
0171 % Store the variable mandatory for local/remote parallel computing.
0172 
0173 localVars.type=type;
0174 localVars.run_smoother=run_smoother;
0175 localVars.gend=gend;
0176 localVars.Y=Y;
0177 localVars.data_index=data_index;
0178 localVars.missing_value=missing_value;
0179 localVars.varobs=options_.varobs;
0180 localVars.irun=irun;
0181 localVars.endo_nbr=endo_nbr;
0182 localVars.nvn=nvn;
0183 localVars.naK=naK;
0184 localVars.horizon=horizon;
0185 localVars.iendo=iendo;
0186 if horizon
0187     localVars.i_last_obs=i_last_obs;
0188     localVars.IdObs=IdObs;
0189     localVars.MAX_nforc1=MAX_nforc1;
0190     localVars.MAX_nforc2=MAX_nforc2;
0191 end
0192 localVars.exo_nbr=exo_nbr;
0193 localVars.maxlag=maxlag;
0194 localVars.MAX_nsmoo=MAX_nsmoo;
0195 localVars.MAX_ninno=MAX_ninno;
0196 localVars.MAX_nerro = MAX_nerro;
0197 if naK
0198     localVars.MAX_naK=MAX_naK;
0199 end
0200 localVars.MAX_nruns=MAX_nruns;
0201 localVars.MAX_momentsno = MAX_momentsno;
0202 localVars.ifil=ifil;
0203 localVars.DirectoryName = DirectoryName;
0204 
0205 if strcmpi(type,'posterior'),
0206     b=0;
0207     while b<=B
0208         b = b + 1;
0209         [x(b,:), logpost(b,1)] = GetOneDraw(type);
0210     end
0211     localVars.logpost=logpost;
0212 end
0213 
0214 if ~strcmpi(type,'prior'),
0215     localVars.x=x;
0216 end
0217 
0218 % Like sequential execution!
0219 if isnumeric(options_.parallel),
0220     [fout] = prior_posterior_statistics_core(localVars,1,B,0);
0221     % Parallel execution!
0222 else
0223     [nCPU, totCPU, nBlockPerCPU] = distributeJobs(options_.parallel, 1, B);
0224     for j=1:totCPU-1,
0225         nfiles = ceil(nBlockPerCPU(j)/MAX_nsmoo);
0226         ifil(1,j+1) =ifil(1,j)+nfiles;
0227         nfiles = ceil(nBlockPerCPU(j)/MAX_ninno);
0228         ifil(2,j+1) =ifil(2,j)+nfiles;
0229         nfiles = ceil(nBlockPerCPU(j)/MAX_nerro);
0230         ifil(3,j+1) =ifil(3,j)+nfiles;
0231         if naK
0232             nfiles = ceil(nBlockPerCPU(j)/MAX_naK);
0233             ifil(4,j+1) =ifil(4,j)+nfiles;
0234         end
0235         nfiles = ceil(nBlockPerCPU(j)/MAX_nruns);
0236         ifil(5,j+1) =ifil(5,j)+nfiles;
0237         if horizon
0238             nfiles = ceil(nBlockPerCPU(j)/MAX_nforc1);
0239             ifil(6,j+1) =ifil(6,j)+nfiles;
0240             nfiles = ceil(nBlockPerCPU(j)/MAX_nforc2);
0241             ifil(7,j+1) =ifil(7,j)+nfiles;
0242         end
0243         %       nfiles = ceil(nBlockPerCPU(j)/MAX_momentsno);
0244         %       ifil(8,j+1) =ifil(8,j)+nfiles;
0245     end
0246     localVars.ifil = ifil;
0247     globalVars = struct('M_',M_, ...
0248                         'options_', options_, ...
0249                         'bayestopt_', bayestopt_, ...
0250                         'estim_params_', estim_params_, ...
0251                         'oo_', oo_);
0252     
0253     % which files have to be copied to run remotely
0254     NamFileInput(1,:) = {'',[M_.fname '_static.m']};
0255     NamFileInput(2,:) = {'',[M_.fname '_dynamic.m']};
0256     if options_.steadystate_flag,
0257         NamFileInput(length(NamFileInput)+1,:)={'',[M_.fname '_steadystate.m']};
0258     end
0259     [fout] = masterParallel(options_.parallel, 1, B,NamFileInput,'prior_posterior_statistics_core', localVars,globalVars, options_.parallel_info);
0260     
0261 end
0262 ifil = fout(end).ifil;
0263 
0264 
0265 stock_gend=gend;
0266 stock_data=Y;
0267 save([DirectoryName '/' M_.fname '_data.mat'],'stock_gend','stock_data');
0268 
0269 if strcmpi(type,'gsa')
0270     return
0271 end
0272 
0273 if ~isnumeric(options_.parallel),
0274     leaveSlaveOpen = options_.parallel_info.leaveSlaveOpen;
0275     if options_.parallel_info.leaveSlaveOpen == 0,
0276         % Commenting for testing!!!
0277         % options_.parallel_info.leaveSlaveOpen = 1; % Force locally to leave open remote matlab sessions (repeated pm3 calls)
0278     end
0279 end
0280 
0281 if options_.smoother
0282     pm3(endo_nbr,gend,ifil(1),B,'Smoothed variables',...
0283         '',M_.endo_names(1:M_.orig_endo_nbr, :),'tit_tex',M_.endo_names,...
0284         varlist,'SmoothedVariables',DirectoryName,'_smooth');
0285     pm3(exo_nbr,gend,ifil(2),B,'Smoothed shocks',...
0286         '',M_.exo_names,'tit_tex',M_.exo_names,...
0287         M_.exo_names,'SmoothedShocks',DirectoryName,'_inno');
0288     if nvn
0289         % needs to  be fixed
0290         %        pm3(endo_nbr,gend,ifil(3),B,'Smoothed measurement errors',...
0291         %            M_.endo_names(SelecVariables),M_.endo_names,'tit_tex',M_.endo_names,...
0292         %            'names2','smooth_errors',[M_.fname '/metropolis'],'_error')
0293     end
0294 end
0295 
0296 if options_.filtered_vars
0297     pm3(endo_nbr,gend,ifil(1),B,'Updated Variables',...
0298         '',varlist,'tit_tex',M_.endo_names,...
0299         varlist,'UpdatedVariables',DirectoryName, ...
0300         '_update');
0301     pm3(endo_nbr,gend+1,ifil(4),B,'One step ahead forecast',...
0302         '',varlist,'tit_tex',M_.endo_names,...
0303         varlist,'FilteredVariables',DirectoryName,'_filter_step_ahead');
0304 end
0305 
0306 if options_.forecast
0307     pm3(endo_nbr,horizon+maxlag,ifil(6),B,'Forecasted variables (mean)',...
0308         '',varlist,'tit_tex',M_.endo_names,...
0309         varlist,'MeanForecast',DirectoryName,'_forc_mean');
0310     pm3(endo_nbr,horizon+maxlag,ifil(6),B,'Forecasted variables (point)',...
0311         '',varlist,'tit_tex',M_.endo_names,...
0312         varlist,'PointForecast',DirectoryName,'_forc_point');
0313 end
0314 
0315 
0316 if ~isnumeric(options_.parallel),
0317     options_.parallel_info.leaveSlaveOpen = leaveSlaveOpen;
0318     if leaveSlaveOpen == 0,
0319         closeSlave(options_.parallel,options_.parallel_info.RemoteTmpFolder),
0320     end
0321 end