1function prior_posterior_statistics(type,dataset,dataset_info) 2 3% function prior_posterior_statistics(type,dataset) 4% Computes Monte Carlo filter smoother and forecasts 5% 6% INPUTS 7% type: posterior 8% prior 9% gsa 10% dataset: data structure 11% 12% OUTPUTS 13% none 14% 15% SPECIAL REQUIREMENTS 16% none 17% 18% PARALLEL CONTEXT 19% See the comments in the posterior_sampler.m funtion. 20 21 22% Copyright (C) 2005-2020 Dynare Team 23% 24% This file is part of Dynare. 25% 26% Dynare is free software: you can redistribute it and/or modify 27% it under the terms of the GNU General Public License as published by 28% the Free Software Foundation, either version 3 of the License, or 29% (at your option) any later version. 30% 31% Dynare is distributed in the hope that it will be useful, 32% but WITHOUT ANY WARRANTY; without even the implied warranty of 33% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 34% GNU General Public License for more details. 35% 36% You should have received a copy of the GNU General Public License 37% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 38 39global options_ estim_params_ oo_ M_ bayestopt_ 40 41localVars=[]; 42 43Y = transpose(dataset.data); 44gend = dataset.nobs; 45data_index = dataset_info.missing.aindex; 46missing_value = dataset_info.missing.state; 47mean_varobs = dataset_info.descriptive.mean; 48 49 50nvx = estim_params_.nvx; 51nvn = estim_params_.nvn; 52ncx = estim_params_.ncx; 53ncn = estim_params_.ncn; 54np = estim_params_.np ; 55npar = nvx+nvn+ncx+ncn+np; 56naK = length(options_.filter_step_ahead); 57 58MaxNumberOfBytes=options_.MaxNumberOfBytes; 59endo_nbr=M_.endo_nbr; 60exo_nbr=M_.exo_nbr; 61meas_err_nbr=length(M_.Correlation_matrix_ME); 62iendo = 1:endo_nbr; 63horizon = options_.forecast; 64IdObs = bayestopt_.mfys; 65if horizon 66 i_last_obs = gend+(1-M_.maximum_endo_lag:0); 67end 68maxlag = M_.maximum_endo_lag; 69 70if strcmpi(type,'posterior') 71 DirectoryName = CheckPath('metropolis',M_.dname); 72 B = options_.sub_draws; 73elseif strcmpi(type,'gsa') 74 RootDirectoryName = CheckPath('gsa',M_.dname); 75 if options_.opt_gsa.pprior 76 DirectoryName = CheckPath(['gsa',filesep,'prior'],M_.dname); 77 load([ RootDirectoryName filesep M_.fname '_prior.mat'],'lpmat0','lpmat','istable') 78 else 79 DirectoryName = CheckPath(['gsa',filesep,'mc'],M_.dname); 80 load([ RootDirectoryName filesep M_.fname '_mc.mat'],'lpmat0','lpmat','istable') 81 end 82 if ~isempty(lpmat0) 83 x=[lpmat0(istable,:) lpmat(istable,:)]; 84 else 85 x=lpmat(istable,:); 86 end 87 clear lpmat lpmat0 istable 88 B = size(x,1); 89elseif strcmpi(type,'prior') 90 DirectoryName = CheckPath('prior',M_.dname); 91 B = options_.prior_draws; 92end 93 94MAX_nruns = min(B,ceil(MaxNumberOfBytes/(npar+2)/8)); 95MAX_nsmoo = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*gend)/8)); 96MAX_n_smoothed_constant = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*gend)/8)); 97MAX_n_smoothed_trend = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*gend)/8)); 98MAX_n_trend_coeff = min(B,ceil(MaxNumberOfBytes/endo_nbr/8)); 99MAX_ninno = min(B,ceil(MaxNumberOfBytes/(exo_nbr*gend)/8)); 100MAX_nerro = min(B,ceil(MaxNumberOfBytes/(size(options_.varobs,1)*gend)/8)); 101 102if naK 103 MAX_naK = min(B,ceil(MaxNumberOfBytes/(endo_nbr* ... 104 length(options_.filter_step_ahead)*(gend+max(options_.filter_step_ahead)))/8)); 105end 106 107if horizon 108 MAX_nforc1 = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*(horizon+maxlag))/8)); 109 MAX_nforc2 = min(B,ceil(MaxNumberOfBytes/((endo_nbr)*(horizon+maxlag))/8)); 110 if ~isequal(M_.H,0) 111 MAX_nforc_ME = min(B,ceil(MaxNumberOfBytes/((size(options_.varobs,1))*(horizon+maxlag))/8)); 112 end 113end 114MAX_momentsno = min(B,ceil(MaxNumberOfBytes/(get_moments_size(options_)*8))); 115 116if options_.filter_covariance 117 MAX_filter_covariance = min(B,ceil(MaxNumberOfBytes/(endo_nbr^2*(gend+1))/8)); 118end 119 120if options_.smoothed_state_uncertainty 121 MAX_n_smoothed_state_uncertainty = min(B,ceil(MaxNumberOfBytes/((endo_nbr*endo_nbr)*gend)/8)); 122end 123 124varlist = options_.varlist; 125if isempty(varlist) 126 varlist = sort(M_.endo_names(1:M_.orig_endo_nbr)); 127end 128nvar = length(varlist); 129SelecVariables = []; 130for i=1:nvar 131 if ~isempty(strmatch(varlist{i}, M_.endo_names, 'exact')) 132 SelecVariables = [SelecVariables; strmatch(varlist{i}, M_.endo_names, 'exact')]; 133 end 134end 135 136n_variables_to_fill=13; 137 138irun = ones(n_variables_to_fill,1); 139ifil = zeros(n_variables_to_fill,1); 140 141run_smoother = 0; 142if options_.smoother || options_.forecast || options_.filter_step_ahead || options_.smoothed_state_uncertainty 143 run_smoother = 1; 144 if options_.loglinear 145 oo_.Smoother.loglinear = true; 146 else 147 oo_.Smoother.loglinear = false; 148 end 149end 150 151filter_covariance=0; 152if options_.filter_covariance 153 filter_covariance=1; 154end 155 156smoothed_state_uncertainty=0; 157if options_.smoothed_state_uncertainty 158 smoothed_state_uncertainty=1; 159end 160 161% Store the variable mandatory for local/remote parallel computing. 162 163localVars.type=type; 164localVars.run_smoother=run_smoother; 165localVars.filter_covariance=filter_covariance; 166localVars.smoothed_state_uncertainty=smoothed_state_uncertainty; 167localVars.gend=gend; 168localVars.Y=Y; 169localVars.data_index=data_index; 170localVars.missing_value=missing_value; 171localVars.varobs=options_.varobs; 172localVars.mean_varobs=mean_varobs; 173localVars.irun=irun; 174localVars.endo_nbr=endo_nbr; 175localVars.nvn=nvn; 176localVars.naK=naK; 177localVars.horizon=horizon; 178localVars.iendo=iendo; 179localVars.IdObs=IdObs; 180if horizon 181 localVars.i_last_obs=i_last_obs; 182 localVars.MAX_nforc1=MAX_nforc1; 183 localVars.MAX_nforc2=MAX_nforc2; 184 if ~isequal(M_.H,0) 185 localVars.MAX_nforc_ME = MAX_nforc_ME; 186 end 187end 188localVars.exo_nbr=exo_nbr; 189localVars.maxlag=maxlag; 190localVars.MAX_nsmoo=MAX_nsmoo; 191localVars.MAX_ninno=MAX_ninno; 192localVars.MAX_nerro = MAX_nerro; 193if naK 194 localVars.MAX_naK=MAX_naK; 195end 196if options_.filter_covariance 197 localVars.MAX_filter_covariance = MAX_filter_covariance; 198end 199if options_.smoothed_state_uncertainty 200 localVars.MAX_n_smoothed_state_uncertainty = MAX_n_smoothed_state_uncertainty ; 201end 202localVars.MAX_n_smoothed_constant=MAX_n_smoothed_constant; 203localVars.MAX_n_smoothed_trend=MAX_n_smoothed_trend; 204localVars.MAX_n_trend_coeff=MAX_n_trend_coeff; 205localVars.MAX_nruns=MAX_nruns; 206localVars.MAX_momentsno = MAX_momentsno; 207localVars.ifil=ifil; 208localVars.DirectoryName = DirectoryName; 209 210if strcmpi(type,'posterior') 211 BaseName = [DirectoryName filesep M_.fname]; 212 load_last_mh_history_file(DirectoryName, M_.fname); 213 FirstMhFile = record.KeepedDraws.FirstMhFile; 214 FirstLine = record.KeepedDraws.FirstLine; 215 TotalNumberOfMhFiles = sum(record.MhDraws(:,2)); 216 LastMhFile = TotalNumberOfMhFiles; 217 TotalNumberOfMhDraws = sum(record.MhDraws(:,1)); 218 NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws); 219 mh_nblck = options_.mh_nblck; 220 if B==NumberOfDraws*mh_nblck 221 % we load all retained MH runs ! 222 logpost=GetAllPosteriorDraws(0, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws); 223 for column=1:npar 224 x(:,column) = GetAllPosteriorDraws(column, FirstMhFile, FirstLine, TotalNumberOfMhFiles, NumberOfDraws); 225 end 226 else 227 logpost=NaN(B,1); 228 for b=1:B 229 [x(b,:), logpost(b)] = GetOneDraw(type,M_,estim_params_,oo_,options_,bayestopt_); 230 end 231 end 232 localVars.logpost=logpost; 233end 234 235if ~strcmpi(type,'prior') 236 localVars.x=x; 237end 238 239% Like sequential execution! 240if isnumeric(options_.parallel) 241 [fout] = prior_posterior_statistics_core(localVars,1,B,0); 242 % Parallel execution! 243else 244 [nCPU, totCPU, nBlockPerCPU] = distributeJobs(options_.parallel, 1, B); 245 ifil=zeros(n_variables_to_fill,totCPU); 246 for j=1:totCPU-1 247 if run_smoother 248 nfiles = ceil(nBlockPerCPU(j)/MAX_nsmoo); 249 ifil(1,j+1) =ifil(1,j)+nfiles; 250 nfiles = ceil(nBlockPerCPU(j)/MAX_ninno); 251 ifil(2,j+1) =ifil(2,j)+nfiles; 252 nfiles = ceil(nBlockPerCPU(j)/MAX_nerro); 253 ifil(3,j+1) =ifil(3,j)+nfiles; 254 end 255 if naK 256 nfiles = ceil(nBlockPerCPU(j)/MAX_naK); 257 ifil(4,j+1) =ifil(4,j)+nfiles; 258 end 259 nfiles = ceil(nBlockPerCPU(j)/MAX_nruns); 260 ifil(5,j+1) =ifil(5,j)+nfiles; 261 if horizon 262 nfiles = ceil(nBlockPerCPU(j)/MAX_nforc1); 263 ifil(6,j+1) =ifil(6,j)+nfiles; 264 nfiles = ceil(nBlockPerCPU(j)/MAX_nforc2); 265 ifil(7,j+1) =ifil(7,j)+nfiles; 266 if ~isequal(M_.H,0) 267 nfiles = ceil(nBlockPerCPU(j)/MAX_nforc_ME); 268 ifil(12,j+1) =ifil(12,j)+nfiles; 269 end 270 end 271 if options_.filter_covariance 272 nfiles = ceil(nBlockPerCPU(j)/MAX_filter_covariance); 273 ifil(8,j+1) =ifil(8,j)+nfiles; 274 end 275 if run_smoother 276 nfiles = ceil(nBlockPerCPU(j)/MAX_n_trend_coeff); 277 ifil(9,j+1) =ifil(9,j)+nfiles; 278 nfiles = ceil(nBlockPerCPU(j)/MAX_n_smoothed_constant); 279 ifil(10,j+1) =ifil(10,j)+nfiles; 280 nfiles = ceil(nBlockPerCPU(j)/MAX_n_smoothed_trend); 281 ifil(11,j+1) =ifil(11,j)+nfiles; 282 if smoothed_state_uncertainty 283 nfiles = ceil(nBlockPerCPU(j)/MAX_n_smoothed_state_uncertainty); 284 ifil(13,j+1) =ifil(13,j)+nfiles; 285 end 286 end 287 end 288 localVars.ifil = ifil; 289 globalVars = struct('M_',M_, ... 290 'options_', options_, ... 291 'bayestopt_', bayestopt_, ... 292 'estim_params_', estim_params_, ... 293 'oo_', oo_); 294 % which files have to be copied to run remotely 295 NamFileInput(1,:) = {'',[M_.fname '.static.m']}; 296 NamFileInput(2,:) = {'',[M_.fname '.dynamic.m']}; 297 if M_.set_auxiliary_variables 298 NamFileInput(3,:) = {'',[M_.fname '.set_auxiliary_variables.m']}; 299 end 300 if options_.steadystate_flag 301 if options_.steadystate_flag == 1 302 NamFileInput(length(NamFileInput)+1,:)={'',[M_.fname '_steadystate.m']}; 303 else 304 NamFileInput(length(NamFileInput)+1,:)={'',[M_.fname '.steadystate.m']}; 305 end 306 end 307 [fout] = masterParallel(options_.parallel, 1, B,NamFileInput,'prior_posterior_statistics_core', localVars,globalVars, options_.parallel_info); 308 309end 310ifil = fout(end).ifil; 311 312stock_gend=gend; 313stock_data=Y; 314save([DirectoryName '/' M_.fname '_data.mat'],'stock_gend','stock_data'); 315 316if strcmpi(type,'gsa') 317 return 318end 319 320if ~isnumeric(options_.parallel), 321 leaveSlaveOpen = options_.parallel_info.leaveSlaveOpen; 322 if options_.parallel_info.leaveSlaveOpen == 0, 323 % Commenting for testing!!! 324 options_.parallel_info.leaveSlaveOpen = 1; % Force locally to leave open remote matlab sessions (repeated pm3 calls) 325 end 326end 327 328if options_.smoother 329 pm3(endo_nbr,gend,ifil(1),B,'Smoothed variables',... 330 '',varlist, M_.endo_names_tex,M_.endo_names,... 331 varlist,'SmoothedVariables',DirectoryName,'_smooth'); 332 pm3(exo_nbr,gend,ifil(2),B,'Smoothed shocks',... 333 '',M_.exo_names,M_.exo_names_tex,M_.exo_names,... 334 M_.exo_names,'SmoothedShocks',DirectoryName,'_inno'); 335 pm3(endo_nbr,1,ifil(9),B,'Trend_coefficients',... 336 '',varlist,M_.endo_names_tex,M_.endo_names,... 337 varlist,'TrendCoeff',DirectoryName,'_trend_coeff'); 338 pm3(endo_nbr,gend,ifil(10),B,'Smoothed constant',... 339 '',varlist,M_.endo_names_tex,M_.endo_names,... 340 varlist,'Constant',DirectoryName,'_smoothed_constant'); 341 pm3(endo_nbr,gend,ifil(11),B,'Smoothed trend',... 342 '',varlist,M_.endo_names_tex,M_.endo_names,... 343 varlist,'Trend',DirectoryName,'_smoothed_trend'); 344 pm3(endo_nbr,gend,ifil(1),B,'Updated Variables',... 345 '',varlist,M_.endo_names_tex,M_.endo_names,... 346 varlist,'UpdatedVariables',DirectoryName, ... 347 '_update'); 348 if smoothed_state_uncertainty 349 pm3(endo_nbr,endo_nbr,ifil(13),B,'State Uncertainty',... 350 '',varlist,M_.endo_names_tex,M_.endo_names,... 351 varlist,'StateUncertainty',DirectoryName,'_state_uncert'); 352 end 353 354 if nvn 355 for obs_iter=1:length(options_.varobs) 356 meas_error_names{obs_iter,1}=['SE_EOBS_' M_.endo_names{strmatch(options_.varobs{obs_iter},M_.endo_names,'exact')}]; 357 texnames{obs_iter,1}=['SE_EOBS_' M_.endo_names_tex{strmatch(options_.varobs{obs_iter},M_.endo_names,'exact')}]; 358 end 359 pm3(meas_err_nbr,gend,ifil(3),B,'Smoothed measurement errors',... 360 '',meas_error_names,texnames,meas_error_names,... 361 meas_error_names,'SmoothedMeasurementErrors',DirectoryName,'_error') 362 end 363end 364 365if options_.filtered_vars 366 pm3(endo_nbr,gend,ifil(4),B,'One step ahead forecast (filtered variables)',... 367 '',varlist,M_.endo_names_tex,M_.endo_names,... 368 varlist,'FilteredVariables',DirectoryName,'_filter_step_ahead'); 369end 370 371if options_.forecast 372 pm3(endo_nbr,horizon,ifil(6),B,'Forecasted variables (mean)',... 373 '',varlist,M_.endo_names_tex,M_.endo_names,... 374 varlist,'MeanForecast',DirectoryName,'_forc_mean'); 375 pm3(endo_nbr,horizon,ifil(7),B,'Forecasted variables (point)',... 376 '',varlist,M_.endo_names_tex,M_.endo_names,... 377 varlist,'PointForecast',DirectoryName,'_forc_point'); 378 if ~isequal(M_.H,0) && ~isempty(intersect(options_.varobs,varlist)) 379 texnames = cell(length(options_.varobs), 1); 380 obs_names = cell(length(options_.varobs), 1); 381 for obs_iter=1:length(options_.varobs) 382 obs_names{obs_iter}=M_.endo_names{strmatch(options_.varobs{obs_iter},M_.endo_names,'exact')}; 383 texnames{obs_iter}=M_.endo_names_tex{strmatch(options_.varobs{obs_iter},M_.endo_names,'exact')}; 384 end 385 varlist_forecast_ME=intersect(options_.varobs,varlist); 386 pm3(meas_err_nbr,horizon,ifil(12),B,'Forecasted variables (point) with ME',... 387 '',varlist_forecast_ME,texnames,obs_names,... 388 varlist_forecast_ME,'PointForecastME',DirectoryName,'_forc_point_ME') 389 end 390end 391 392if options_.filter_covariance 393 pm3(endo_nbr,endo_nbr,ifil(8),B,'Filtered covariances',... 394 '',varlist,M_.endo_names_tex,M_.endo_names,... 395 varlist,'FilterCovariance',DirectoryName,'_filter_covar'); 396end 397 398 399if ~isnumeric(options_.parallel) 400 options_.parallel_info.leaveSlaveOpen = leaveSlaveOpen; 401 if leaveSlaveOpen == 0 402 closeSlave(options_.parallel,options_.parallel_info.RemoteTmpFolder), 403 end 404end