1function dynare_estimation_1(var_list_,dname) 2% function dynare_estimation_1(var_list_,dname) 3% runs the estimation of the model 4% 5% INPUTS 6% var_list_: selected endogenous variables vector 7% dname: alternative directory name 8% 9% OUTPUTS 10% none 11% 12% SPECIAL REQUIREMENTS 13% none 14 15% Copyright (C) 2003-2018 Dynare Team 16% 17% This file is part of Dynare. 18% 19% Dynare is free software: you can redistribute it and/or modify 20% it under the terms of the GNU General Public License as published by 21% the Free Software Foundation, either version 3 of the License, or 22% (at your option) any later version. 23% 24% Dynare is distributed in the hope that it will be useful, 25% but WITHOUT ANY WARRANTY; without even the implied warranty of 26% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 27% GNU General Public License for more details. 28% 29% You should have received a copy of the GNU General Public License 30% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 31 32global M_ options_ oo_ estim_params_ bayestopt_ dataset_ dataset_info 33 34if isempty(estim_params_) 35 mode_compute_o = options_.mode_compute; 36 mh_replic_o = options_.mh_replic; 37 options_.mode_compute = 0; 38 options_.mh_replic = 0; 39 reset_options_related_to_estimation = true; 40else 41 reset_options_related_to_estimation = false; 42end 43 44 45%store qz_criterium 46qz_criterium_old=options_.qz_criterium; 47if isnan(options_.first_obs) 48 first_obs_nan_indicator=true; 49else 50 first_obs_nan_indicator=false; 51end 52 53% Set particle filter flag. 54if options_.order > 1 55 if options_.particle.status 56 skipline() 57 disp('Estimation using a non linear filter!') 58 skipline() 59 if ~options_.nointeractive && ismember(options_.mode_compute,[1,3,4]) && ~strcmpi(options_.particle.filter_algorithm,'gf')% Known gradient-based optimizers 60 disp('You are using a gradient-based mode-finder. Particle filtering introduces discontinuities in the') 61 disp('objective function w.r.t the parameters. Thus, should use a non-gradient based optimizer.') 62 fprintf('\nPlease choose a mode-finder:\n') 63 fprintf('\t 0 - Continue using gradient-based method (it is most likely that you will no get any sensible result).\n') 64 fprintf('\t 6 - Monte Carlo based algorithm\n') 65 fprintf('\t 7 - Nelder-Mead simplex based optimization routine (Matlab optimization toolbox required)\n') 66 fprintf('\t 8 - Nelder-Mead simplex based optimization routine (Dynare''s implementation)\n') 67 fprintf('\t 9 - CMA-ES (Covariance Matrix Adaptation Evolution Strategy) algorithm\n') 68 choice = []; 69 while isempty(choice) 70 choice = input('Please enter your choice: '); 71 if isnumeric(choice) && isint(choice) && ismember(choice,[0 6 7 8 9]) 72 if choice 73 options_.mode_compute = choice; 74 end 75 else 76 fprintf('\nThis is an invalid choice (you have to choose between 0, 6, 7, 8 and 9).\n') 77 choice = []; 78 end 79 end 80 end 81 else 82 error('For estimating the model with a second order approximation using a non linear filter, one should have options_.particle.status=true;') 83 end 84end 85 86if ~options_.dsge_var 87 if options_.particle.status 88 objective_function = str2func('non_linear_dsge_likelihood'); 89 if strcmpi(options_.particle.filter_algorithm, 'sis') 90 options_.particle.algorithm = 'sequential_importance_particle_filter'; 91 elseif strcmpi(options_.particle.filter_algorithm, 'apf') 92 options_.particle.algorithm = 'auxiliary_particle_filter'; 93 elseif strcmpi(options_.particle.filter_algorithm, 'gf') 94 options_.particle.algorithm = 'gaussian_filter'; 95 elseif strcmpi(options_.particle.filter_algorithm, 'gmf') 96 options_.particle.algorithm = 'gaussian_mixture_filter'; 97 elseif strcmpi(options_.particle.filter_algorithm, 'cpf') 98 options_.particle.algorithm = 'conditional_particle_filter'; 99 elseif strcmpi(options_.particle.filter_algorithm, 'nlkf') 100 options_.particle.algorithm = 'nonlinear_kalman_filter'; 101 else 102 error(['Estimation: Unknown filter ' options_.particle.filter_algorithm]) 103 end 104 else 105 objective_function = str2func('dsge_likelihood'); 106 end 107else 108 objective_function = str2func('dsge_var_likelihood'); 109end 110 111[dataset_, dataset_info, xparam1, hh, M_, options_, oo_, estim_params_, bayestopt_, bounds] = ... 112 dynare_estimation_init(var_list_, dname, [], M_, options_, oo_, estim_params_, bayestopt_); 113 114if options_.dsge_var 115 check_dsge_var_model(M_, estim_params_, bayestopt_); 116 if dataset_info.missing.state 117 error('Estimation::DsgeVarLikelihood: I cannot estimate a DSGE-VAR model with missing observations!') 118 end 119 if options_.noconstant 120 var_sample_moments(options_.dsge_varlag, -1, dataset_); 121 else 122 % The steady state is non zero ==> a constant in the VAR is needed! 123 var_sample_moments(options_.dsge_varlag, 0, dataset_); 124 end 125end 126 127% Set sigma_e_is_diagonal flag (needed if the shocks block is not declared in the mod file). 128M_.sigma_e_is_diagonal = true; 129if estim_params_.ncx || any(nnz(tril(M_.Correlation_matrix,-1))) || isfield(estim_params_,'calibrated_covariances') 130 M_.sigma_e_is_diagonal = false; 131end 132 133data = dataset_.data; 134rawdata = dataset_info.rawdata; 135data_index = dataset_info.missing.aindex; 136missing_value = dataset_info.missing.state; 137 138% Set number of observations 139gend = dataset_.nobs; 140 141% Set the number of observed variables. 142n_varobs = length(options_.varobs); 143 144% Get the number of parameters to be estimated. 145nvx = estim_params_.nvx; % Variance of the structural innovations (number of parameters). 146nvn = estim_params_.nvn; % Variance of the measurement innovations (number of parameters). 147ncx = estim_params_.ncx; % Covariance of the structural innovations (number of parameters). 148ncn = estim_params_.ncn; % Covariance of the measurement innovations (number of parameters). 149np = estim_params_.np ; % Number of deep parameters. 150nx = nvx+nvn+ncx+ncn+np; % Total number of parameters to be estimated. 151 152% Set the names of the priors. 153pnames = {''; 'beta'; 'gamm'; 'norm'; 'invg'; 'unif'; 'invg2'; ''; 'weibl'}; 154 155dr = oo_.dr; 156 157if ~isempty(estim_params_) 158 M_ = set_all_parameters(xparam1,estim_params_,M_); 159end 160 161 162%% perform initial estimation checks; 163try 164 oo_ = initial_estimation_checks(objective_function,xparam1,dataset_,dataset_info,M_,estim_params_,options_,bayestopt_,bounds,oo_); 165catch % if check fails, provide info on using calibration if present 166 e = lasterror(); 167 if estim_params_.full_calibration_detected %calibrated model present and no explicit starting values 168 skipline(1); 169 fprintf('ESTIMATION_CHECKS: There was an error in computing the likelihood for initial parameter values.\n') 170 fprintf('ESTIMATION_CHECKS: If this is not a problem with the setting of options (check the error message below),\n') 171 fprintf('ESTIMATION_CHECKS: you should try using the calibrated version of the model as starting values. To do\n') 172 fprintf('ESTIMATION_CHECKS: this, add an empty estimated_params_init-block with use_calibration option immediately before the estimation\n') 173 fprintf('ESTIMATION_CHECKS: command (and after the estimated_params-block so that it does not get overwritten):\n'); 174 skipline(2); 175 end 176 rethrow(e); 177end 178 179if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && options_.mh_posterior_mode_estimation==0 180 if options_.order==1 && ~options_.particle.status 181 if options_.smoother 182 [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,options_,bayestopt_] = DsgeSmoother(xparam1,gend,transpose(data),data_index,missing_value,M_,oo_,options_,bayestopt_,estim_params_); 183 [oo_]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty); 184 if options_.forecast > 0 185 oo_.forecast = dyn_forecast(var_list_,M_,options_,oo_,'smoother',dataset_info); 186 end 187 end 188 %reset qz_criterium 189 options_.qz_criterium=qz_criterium_old; 190 return 191 else %allow to continue, e.g. with MCMC_jumping_covariance 192 if options_.smoother 193 error('Estimation:: Particle Smoothers are not yet implemented.') 194 end 195 end 196end 197 198%% Estimation of the posterior mode or likelihood mode 199 200if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation 201 %prepare settings for newrat 202 if options_.mode_compute==5 203 %get whether outer product Hessian is requested 204 newratflag=[]; 205 if ~isempty(options_.optim_opt) 206 options_list = read_key_value_string(options_.optim_opt); 207 for i=1:rows(options_list) 208 if strcmp(options_list{i,1},'Hessian') 209 newratflag=options_list{i,2}; 210 end 211 end 212 end 213 if options_.analytic_derivation 214 options_analytic_derivation_old = options_.analytic_derivation; 215 options_.analytic_derivation = -1; 216 if ~isempty(newratflag) && newratflag~=0 %numerical hessian explicitly specified 217 error('newrat: analytic_derivation is incompatible with numerical Hessian.') 218 else %use default 219 newratflag=0; %exclude DYNARE numerical hessian 220 end 221 elseif ~options_.analytic_derivation 222 if isempty(newratflag) 223 newratflag=options_.newrat.hess; %use default numerical dynare hessian 224 end 225 end 226 end 227 228 [xparam1, fval, exitflag, hh, options_, Scale, new_rat_hess_info] = dynare_minimize_objective(objective_function,xparam1,options_.mode_compute,options_,[bounds.lb bounds.ub],bayestopt_.name,bayestopt_,hh,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_); 229 fprintf('\nFinal value of minus the log posterior (or likelihood):%f \n', fval); 230 231 if isnumeric(options_.mode_compute) && options_.mode_compute==5 && options_.analytic_derivation==-1 %reset options changed by newrat 232 options_.analytic_derivation = options_analytic_derivation_old; %reset 233 elseif isnumeric(options_.mode_compute) && options_.mode_compute==6 %save scaling factor 234 save([M_.fname '_optimal_mh_scale_parameter.mat'],'Scale'); 235 options_.mh_jscale = Scale; 236 bayestopt_.jscale(:) = options_.mh_jscale; 237 end 238 if ~isnumeric(options_.mode_compute) || ~isequal(options_.mode_compute,6) %always already computes covariance matrix 239 if options_.cova_compute == 1 %user did not request covariance not to be computed 240 if options_.analytic_derivation && strcmp(func2str(objective_function),'dsge_likelihood') 241 ana_deriv_old = options_.analytic_derivation; 242 options_.analytic_derivation = 2; 243 [~,~,~,~,hh] = feval(objective_function,xparam1, ... 244 dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_); 245 options_.analytic_derivation = ana_deriv_old; 246 elseif ~isnumeric(options_.mode_compute) || ~(isequal(options_.mode_compute,5) && newratflag~=1 && strcmp(func2str(objective_function),'dsge_likelihood')) 247 % with flag==0, we force to use the hessian from outer product gradient of optimizer 5 248 if options_.hessian.use_penalized_objective 249 penalized_objective_function = str2func('penalty_objective_function'); 250 hh = hessian(penalized_objective_function, xparam1, options_.gstep, objective_function, fval, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_); 251 else 252 hh = hessian(objective_function, xparam1, options_.gstep, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_); 253 end 254 hh = reshape(hh, nx, nx); 255 elseif isnumeric(options_.mode_compute) && isequal(options_.mode_compute,5) 256 % other numerical hessian options available with optimizer 5 257 % 258 % if newratflag == 0 259 % compute outer product gradient of optimizer 5 260 % 261 % if newratflag == 2 262 % compute 'mixed' outer product gradient of optimizer 5 263 % with diagonal elements computed with numerical second order derivatives 264 % 265 % uses univariate filters, so to get max # of available 266 % densitities for outer product gradient 267 kalman_algo0 = options_.kalman_algo; 268 compute_hessian = 1; 269 if ~((options_.kalman_algo == 2) || (options_.kalman_algo == 4)) 270 options_.kalman_algo=2; 271 if options_.lik_init == 3 272 options_.kalman_algo=4; 273 end 274 elseif newratflag==0 % hh already contains outer product gradient with univariate filter 275 compute_hessian = 0; 276 end 277 if compute_hessian 278 crit = options_.newrat.tolerance.f; 279 newratflag = newratflag>0; 280 hh = reshape(mr_hessian(xparam1,objective_function,fval,newratflag,crit,new_rat_hess_info,dataset_, dataset_info, options_,M_,estim_params_,bayestopt_,bounds,oo_), nx, nx); 281 end 282 options_.kalman_algo = kalman_algo0; 283 end 284 end 285 end 286 parameter_names = bayestopt_.name; 287 if options_.cova_compute || options_.mode_compute==5 || options_.mode_compute==6 288 save([M_.fname '_mode.mat'],'xparam1','hh','parameter_names','fval'); 289 else 290 save([M_.fname '_mode.mat'],'xparam1','parameter_names','fval'); 291 end 292end 293 294if ~options_.mh_posterior_mode_estimation && options_.cova_compute 295 try 296 chol(hh); 297 catch 298 skipline() 299 disp('POSTERIOR KERNEL OPTIMIZATION PROBLEM!') 300 disp(' (minus) the hessian matrix at the "mode" is not positive definite!') 301 disp('=> posterior variance of the estimated parameters are not positive.') 302 disp('You should try to change the initial values of the parameters using') 303 disp('the estimated_params_init block, or use another optimization routine.') 304 params_at_bound=find(abs(xparam1-bounds.ub)<1.e-10 | abs(xparam1-bounds.lb)<1.e-10); 305 if ~isempty(params_at_bound) 306 for ii=1:length(params_at_bound) 307 params_at_bound_name{ii,1}=get_the_name(params_at_bound(ii),0,M_,estim_params_,options_); 308 end 309 disp_string=[params_at_bound_name{1,:}]; 310 for ii=2:size(params_at_bound_name,1) 311 disp_string=[disp_string,', ',params_at_bound_name{ii,:}]; 312 end 313 fprintf('\nThe following parameters are at the prior bound: %s\n', disp_string) 314 fprintf('Some potential solutions are:\n') 315 fprintf(' - Check your model for mistakes.\n') 316 fprintf(' - Check whether model and data are consistent (correct observation equation).\n') 317 fprintf(' - Shut off prior_trunc.\n') 318 fprintf(' - Change the optimization bounds.\n') 319 fprintf(' - Use a different mode_compute like 6 or 9.\n') 320 fprintf(' - Check whether the parameters estimated are identified.\n') 321 fprintf(' - Check prior shape (e.g. Inf density at bound(s)).\n') 322 fprintf(' - Increase the informativeness of the prior.\n') 323 end 324 warning('The results below are most likely wrong!'); 325 end 326end 327 328if options_.mode_check.status && ~options_.mh_posterior_mode_estimation 329 ana_deriv_old = options_.analytic_derivation; 330 options_.analytic_derivation = 0; 331 mode_check(objective_function,xparam1,hh,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_); 332 options_.analytic_derivation = ana_deriv_old; 333end 334 335oo_.posterior.optimization.mode = []; 336oo_.posterior.optimization.Variance = []; 337oo_.posterior.optimization.log_density=[]; 338 339invhess=[]; 340if ~options_.mh_posterior_mode_estimation 341 oo_.posterior.optimization.mode = xparam1; 342 if exist('fval','var') 343 oo_.posterior.optimization.log_density=-fval; 344 end 345 if options_.cova_compute 346 hsd = sqrt(diag(hh)); 347 invhess = inv(hh./(hsd*hsd'))./(hsd*hsd'); 348 stdh = sqrt(diag(invhess)); 349 oo_.posterior.optimization.Variance = invhess; 350 end 351else 352 variances = bayestopt_.p2.*bayestopt_.p2; 353 idInf = isinf(variances); 354 variances(idInf) = 1; 355 invhess = options_.mh_posterior_mode_estimation*diag(variances); 356 xparam1 = bayestopt_.p5; 357 idNaN = isnan(xparam1); 358 xparam1(idNaN) = bayestopt_.p1(idNaN); 359 outside_bound_pars=find(xparam1 < bounds.lb | xparam1 > bounds.ub); 360 xparam1(outside_bound_pars) = bayestopt_.p1(outside_bound_pars); 361end 362 363if ~options_.cova_compute 364 stdh = NaN(length(xparam1),1); 365end 366 367if any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation 368 % display results table and store parameter estimates and standard errors in results 369 oo_ = display_estimation_results_table(xparam1, stdh, M_, options_, estim_params_, bayestopt_, oo_, pnames, 'Posterior', 'posterior'); 370 % Laplace approximation to the marginal log density: 371 if options_.cova_compute 372 estim_params_nbr = size(xparam1,1); 373 if ispd(invhess) 374 log_det_invhess = log(det(invhess./(stdh*stdh')))+2*sum(log(stdh)); 375 likelihood = feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_); 376 oo_.MarginalDensity.LaplaceApproximation = .5*estim_params_nbr*log(2*pi) + .5*log_det_invhess - likelihood; 377 else 378 oo_.MarginalDensity.LaplaceApproximation = NaN; 379 end 380 skipline() 381 disp(sprintf('Log data density [Laplace approximation] is %f.',oo_.MarginalDensity.LaplaceApproximation)) 382 skipline() 383 end 384 if options_.dsge_var 385 [~,~,~,~,~,~,~,oo_.dsge_var.posterior_mode.PHI_tilde,oo_.dsge_var.posterior_mode.SIGMA_u_tilde,oo_.dsge_var.posterior_mode.iXX,oo_.dsge_var.posterior_mode.prior] =... 386 feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_); 387 end 388 389elseif ~any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation 390 oo_=display_estimation_results_table(xparam1, stdh, M_, options_, estim_params_, bayestopt_, oo_, pnames, 'Maximum Likelihood', 'mle'); 391end 392 393if np > 0 394 pindx = estim_params_.param_vals(:,1); 395 save([M_.fname '_params.mat'],'pindx'); 396end 397 398switch options_.MCMC_jumping_covariance 399 case 'hessian' %Baseline 400 %do nothing and use hessian from mode_compute 401 case 'prior_variance' %Use prior variance 402 if any(isinf(bayestopt_.p2)) 403 error('Infinite prior variances detected. You cannot use the prior variances as the proposal density, if some variances are Inf.') 404 else 405 hh = diag(1./(bayestopt_.p2.^2)); 406 end 407 hsd = sqrt(diag(hh)); 408 invhess = inv(hh./(hsd*hsd'))./(hsd*hsd'); 409 case 'identity_matrix' %Use identity 410 invhess = eye(nx); 411 otherwise %user specified matrix in file 412 try 413 load(options_.MCMC_jumping_covariance,'jumping_covariance') 414 hh=jumping_covariance; 415 catch 416 error(['No matrix named ''jumping_covariance'' could be found in ',options_.MCMC_jumping_covariance,'.mat']) 417 end 418 [nrow, ncol]=size(hh); 419 if ~isequal(nrow,ncol) && ~isequal(nrow,nx) %check if square and right size 420 error(['jumping_covariance matrix must be square and have ',num2str(nx),' rows and columns']) 421 end 422 try %check for positive definiteness 423 chol(hh); 424 hsd = sqrt(diag(hh)); 425 invhess = inv(hh./(hsd*hsd'))./(hsd*hsd'); 426 catch 427 error(['Specified jumping_covariance is not positive definite']) 428 end 429end 430 431if (any(bayestopt_.pshape >0 ) && options_.mh_replic) || ... 432 (any(bayestopt_.pshape >0 ) && options_.load_mh_file) %% not ML estimation 433 bounds = prior_bounds(bayestopt_, options_.prior_trunc); %reset bounds as lb and ub must only be operational during mode-finding 434 outside_bound_pars=find(xparam1 < bounds.lb | xparam1 > bounds.ub); 435 if ~isempty(outside_bound_pars) 436 for ii=1:length(outside_bound_pars) 437 outside_bound_par_names{ii,1}=get_the_name(ii,0,M_,estim_params_,options_); 438 end 439 disp_string=[outside_bound_par_names{1,:}]; 440 for ii=2:size(outside_bound_par_names,1) 441 disp_string=[disp_string,', ',outside_bound_par_names{ii,:}]; 442 end 443 if options_.prior_trunc>0 444 error(['Estimation:: Mode value(s) of ', disp_string ,' are outside parameter bounds. Potentially, you should set prior_trunc=0.']) 445 else 446 error(['Estimation:: Mode value(s) of ', disp_string ,' are outside parameter bounds.']) 447 end 448 end 449 % Tunes the jumping distribution's scale parameter 450 if options_.mh_tune_jscale.status 451 if strcmp(options_.posterior_sampler_options.posterior_sampling_method, 'random_walk_metropolis_hastings') 452 options = options_.mh_tune_jscale; 453 options.rwmh = options_.posterior_sampler_options.rwmh; 454 options_.mh_jscale = calibrate_mh_scale_parameter(objective_function, ... 455 invhess, xparam1, [bounds.lb,bounds.ub], ... 456 options, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds, oo_); 457 bayestopt_.jscale(:) = options_.mh_jscale; 458 disp(sprintf('mh_jscale has been set equal to %s', num2str(options_.mh_jscale))) 459 skipline() 460 else 461 warning('mh_tune_jscale is only available with Random Walk Metropolis Hastings!') 462 end 463 end 464 % runs MCMC 465 if options_.mh_replic || options_.load_mh_file 466 posterior_sampler_options = options_.posterior_sampler_options.current_options; 467 posterior_sampler_options.invhess = invhess; 468 [posterior_sampler_options, options_] = check_posterior_sampler_options(posterior_sampler_options, options_); 469 % store current options in global 470 options_.posterior_sampler_options.current_options = posterior_sampler_options; 471 if options_.mh_replic 472 ana_deriv_old = options_.analytic_derivation; 473 options_.analytic_derivation = 0; 474 posterior_sampler(objective_function,posterior_sampler_options.proposal_distribution,xparam1,posterior_sampler_options,bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_); 475 options_.analytic_derivation = ana_deriv_old; 476 end 477 end 478 %% Here I discard first mh_drop percent of the draws: 479 CutSample(M_, options_, estim_params_); 480 if options_.mh_posterior_mode_estimation 481 %reset qz_criterium 482 options_.qz_criterium=qz_criterium_old; 483 return 484 else 485 %get stored results if required 486 if options_.load_mh_file && options_.load_results_after_load_mh 487 oo_load_mh=load([M_.fname '_results'],'oo_'); 488 end 489 if ~options_.nodiagnostic 490 if (options_.mh_replic>0 || (options_.load_mh_file && ~options_.load_results_after_load_mh)) 491 oo_= McMCDiagnostics(options_, estim_params_, M_,oo_); 492 elseif options_.load_mh_file && options_.load_results_after_load_mh 493 if isfield(oo_load_mh.oo_,'convergence') 494 oo_.convergence=oo_load_mh.oo_.convergence; 495 end 496 end 497 end 498 %% Estimation of the marginal density from the Mh draws: 499 if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh) 500 [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_); 501 % Store posterior statistics by parameter name 502 oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_, pnames); 503 if ~options_.nograph 504 oo_ = PlotPosteriorDistributions(estim_params_, M_, options_, bayestopt_, oo_); 505 end 506 % Store posterior mean in a vector and posterior variance in 507 % a matrix 508 [oo_.posterior.metropolis.mean,oo_.posterior.metropolis.Variance] ... 509 = GetPosteriorMeanVariance(M_,options_.mh_drop); 510 elseif options_.load_mh_file && options_.load_results_after_load_mh 511 %% load fields from previous MCMC run stored in results-file 512 field_names={'posterior_mode','posterior_std_at_mode',...% fields set by marginal_density 513 'posterior_mean','posterior_hpdinf','posterior_hpdsup','posterior_median','posterior_variance','posterior_std','posterior_deciles','posterior_density',...% fields set by GetPosteriorParametersStatistics 514 'prior_density',...%fields set by PlotPosteriorDistributions 515 }; 516 for field_iter=1:size(field_names,2) 517 if isfield(oo_load_mh.oo_,field_names{1,field_iter}) 518 oo_.(field_names{1,field_iter})=oo_load_mh.oo_.(field_names{1,field_iter}); 519 end 520 end 521 % field set by marginal_density 522 if isfield(oo_load_mh.oo_,'MarginalDensity') && isfield(oo_load_mh.oo_.MarginalDensity,'ModifiedHarmonicMean') 523 oo_.MarginalDensity.ModifiedHarmonicMean=oo_load_mh.oo_.MarginalDensity.ModifiedHarmonicMean; 524 end 525 % field set by GetPosteriorMeanVariance 526 if isfield(oo_load_mh.oo_,'posterior') && isfield(oo_load_mh.oo_.posterior,'metropolis') 527 oo_.posterior.metropolis=oo_load_mh.oo_.posterior.metropolis; 528 end 529 end 530 [error_flag,~,options_]= metropolis_draw(1,options_,estim_params_,M_); 531 if ~(~isempty(options_.sub_draws) && options_.sub_draws==0) 532 if options_.bayesian_irf 533 if error_flag 534 error('Estimation::mcmc: I cannot compute the posterior IRFs!') 535 end 536 PosteriorIRF('posterior'); 537 end 538 if options_.moments_varendo 539 if error_flag 540 error('Estimation::mcmc: I cannot compute the posterior moments for the endogenous variables!') 541 end 542 if options_.load_mh_file && options_.mh_replic==0 %user wants to recompute results 543 [MetropolisFolder, info] = CheckPath('metropolis',M_.dname); 544 if ~info 545 generic_post_data_file_name={'Posterior2ndOrderMoments','decomposition','PosteriorVarianceDecomposition','correlation','PosteriorCorrelations','conditional decomposition','PosteriorConditionalVarianceDecomposition'}; 546 for ii=1:length(generic_post_data_file_name) 547 delete_stale_file([MetropolisFolder filesep M_.fname '_' generic_post_data_file_name{1,ii} '*']); 548 end 549 end 550 end 551 oo_ = compute_moments_varendo('posterior',options_,M_,oo_,var_list_); 552 end 553 if options_.smoother || ~isempty(options_.filter_step_ahead) || options_.forecast 554 if error_flag 555 error('Estimation::mcmc: I cannot compute the posterior statistics!') 556 end 557 if options_.order==1 && ~options_.particle.status 558 prior_posterior_statistics('posterior',dataset_,dataset_info); %get smoothed and filtered objects and forecasts 559 else 560 error('Estimation::mcmc: Particle Smoothers are not yet implemented.') 561 end 562 end 563 else 564 fprintf('Estimation:mcmc: sub_draws was set to 0. Skipping posterior computations.') 565 end 566 xparam1 = get_posterior_parameters('mean',M_,estim_params_,oo_,options_); 567 M_ = set_all_parameters(xparam1,estim_params_,M_); 568 end 569end 570 571if options_.particle.status 572 %reset qz_criterium 573 options_.qz_criterium=qz_criterium_old; 574 return 575end 576 577if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.pshape> 0) && options_.load_mh_file)) ... 578 || ~options_.smoother ) && ~options_.partial_information % to be fixed 579 %% ML estimation, or posterior mode without Metropolis-Hastings or Metropolis without Bayesian smoothes variables 580 [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp,Trend,state_uncertainty,M_,oo_,options_,bayestopt_] = DsgeSmoother(xparam1,dataset_.nobs,transpose(dataset_.data),dataset_info.missing.aindex,dataset_info.missing.state,M_,oo_,options_,bayestopt_,estim_params_); 581 [oo_,yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,dataset_info,atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,P,PK,decomp,Trend,state_uncertainty); 582 583 if ~options_.nograph 584 [nbplt,nr,nc,lr,lc,nstar] = pltorg(M_.exo_nbr); 585 if ~exist([M_.fname '/graphs'],'dir') 586 mkdir(M_.fname,'graphs'); 587 end 588 if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format))) 589 fidTeX = fopen([M_.fname, '/graphs/' M_.fname '_SmoothedShocks.tex'],'w'); 590 fprintf(fidTeX,'%% TeX eps-loader file generated by dynare_estimation_1.m (Dynare).\n'); 591 fprintf(fidTeX,['%% ' datestr(now,0) '\n']); 592 fprintf(fidTeX,' \n'); 593 end 594 for plt = 1:nbplt 595 fh = dyn_figure(options_.nodisplay,'Name','Smoothed shocks'); 596 NAMES = []; 597 if options_.TeX, TeXNAMES = []; end 598 nstar0=min(nstar,M_.exo_nbr-(plt-1)*nstar); 599 if gend==1 600 marker_string{1,1}='-ro'; 601 marker_string{2,1}='-ko'; 602 else 603 marker_string{1,1}='-r'; 604 marker_string{2,1}='-k'; 605 end 606 for i=1:nstar0 607 k = (plt-1)*nstar+i; 608 subplot(nr,nc,i); 609 plot([1 gend],[0 0],marker_string{1,1},'linewidth',.5) 610 hold on 611 plot(1:gend,innov(k,:),marker_string{2,1},'linewidth',1) 612 hold off 613 name = M_.exo_names{k}; 614 if isempty(NAMES) 615 NAMES = name; 616 else 617 NAMES = char(NAMES,name); 618 end 619 if ~isempty(options_.XTick) 620 set(gca,'XTick',options_.XTick) 621 set(gca,'XTickLabel',options_.XTickLabel) 622 end 623 if gend>1 624 xlim([1 gend]) 625 end 626 if options_.TeX 627 texname = M_.exo_names_tex{k}; 628 if isempty(TeXNAMES) 629 TeXNAMES = ['$ ' deblank(texname) ' $']; 630 else 631 TeXNAMES = char(TeXNAMES,['$ ' deblank(texname) ' $']); 632 end 633 end 634 title(name,'Interpreter','none') 635 end 636 dyn_saveas(fh,[M_.fname, '/graphs/' M_.fname '_SmoothedShocks' int2str(plt)],options_.nodisplay,options_.graph_format); 637 if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format))) 638 fprintf(fidTeX,'\\begin{figure}[H]\n'); 639 for jj = 1:nstar0 640 fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:))); 641 end 642 fprintf(fidTeX,'\\centering \n'); 643 fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_SmoothedShocks%s}\n',options_.figures.textwidth*min(i/nc,1),[M_.fname, '/graphs/' M_.fname],int2str(plt)); 644 fprintf(fidTeX,'\\caption{Smoothed shocks.}'); 645 fprintf(fidTeX,'\\label{Fig:SmoothedShocks:%s}\n',int2str(plt)); 646 fprintf(fidTeX,'\\end{figure}\n'); 647 fprintf(fidTeX,'\n'); 648 end 649 end 650 if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format))) 651 fprintf(fidTeX,'\n'); 652 fprintf(fidTeX,'%% End of TeX file.\n'); 653 fclose(fidTeX); 654 end 655 end 656 if nvn 657 number_of_plots_to_draw = 0; 658 index = []; 659 for obs_iter=1:n_varobs 660 if max(abs(measurement_error(obs_iter,:))) > options_.ME_plot_tol; 661 number_of_plots_to_draw = number_of_plots_to_draw + 1; 662 index = cat(1,index,obs_iter); 663 end 664 end 665 if ~options_.nograph 666 [nbplt,nr,nc,lr,lc,nstar] = pltorg(number_of_plots_to_draw); 667 if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format))) 668 fidTeX = fopen([M_.fname, '/graphs/' M_.fname '_SmoothedObservationErrors.tex'],'w'); 669 fprintf(fidTeX,'%% TeX eps-loader file generated by dynare_estimation_1.m (Dynare).\n'); 670 fprintf(fidTeX,['%% ' datestr(now,0) '\n']); 671 fprintf(fidTeX,' \n'); 672 end 673 for plt = 1:nbplt 674 fh = dyn_figure(options_.nodisplay,'Name','Smoothed observation errors'); 675 NAMES = []; 676 if options_.TeX, TeXNAMES = []; end 677 nstar0=min(nstar,number_of_plots_to_draw-(plt-1)*nstar); 678 if gend==1 679 marker_string{1,1}='-ro'; 680 marker_string{2,1}='-ko'; 681 else 682 marker_string{1,1}='-r'; 683 marker_string{2,1}='-k'; 684 end 685 for i=1:nstar0 686 k = (plt-1)*nstar+i; 687 subplot(nr,nc,i); 688 plot([1 gend],[0 0],marker_string{1,1},'linewidth',.5) 689 hold on 690 plot(1:gend,measurement_error(index(k),:),marker_string{2,1},'linewidth',1) 691 hold off 692 name = options_.varobs{index(k)}; 693 if gend>1 694 xlim([1 gend]) 695 end 696 if isempty(NAMES) 697 NAMES = name; 698 else 699 NAMES = char(NAMES,name); 700 end 701 if ~isempty(options_.XTick) 702 set(gca,'XTick',options_.XTick) 703 set(gca,'XTickLabel',options_.XTickLabel) 704 end 705 if options_.TeX 706 idx = strmatch(options_.varobs{index(k)}, M_.endo_names, 'exact'); 707 texname = M_.endo_names_tex{idx}; 708 if isempty(TeXNAMES) 709 TeXNAMES = ['$ ' texname ' $']; 710 else 711 TeXNAMES = char(TeXNAMES,['$ ' texname ' $']); 712 end 713 end 714 title(name,'Interpreter','none') 715 end 716 dyn_saveas(fh,[M_.fname, '/graphs/' M_.fname '_SmoothedObservationErrors' int2str(plt)],options_.nodisplay,options_.graph_format); 717 if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format))) 718 fprintf(fidTeX,'\\begin{figure}[H]\n'); 719 for jj = 1:nstar0 720 fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:))); 721 end 722 fprintf(fidTeX,'\\centering \n'); 723 fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_SmoothedObservationErrors%s}\n',options_.figures.textwidth*min(i/nc,1),[M_.fname, '/graphs/' M_.fname],int2str(plt)); 724 fprintf(fidTeX,'\\caption{Smoothed observation errors.}'); 725 fprintf(fidTeX,'\\label{Fig:SmoothedObservationErrors:%s}\n',int2str(plt)); 726 fprintf(fidTeX,'\\end{figure}\n'); 727 fprintf(fidTeX,'\n'); 728 end 729 end 730 if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format))) 731 fprintf(fidTeX,'\n'); 732 fprintf(fidTeX,'%% End of TeX file.\n'); 733 fclose(fidTeX); 734 end 735 end 736 end 737 %% 738 %% Historical and smoothed variabes 739 %% 740 if ~options_.nograph 741 [nbplt,nr,nc,lr,lc,nstar] = pltorg(n_varobs); 742 if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format))) 743 fidTeX = fopen([M_.fname, '/graphs/' M_.fname '_HistoricalAndSmoothedVariables.tex'],'w'); 744 fprintf(fidTeX,'%% TeX eps-loader file generated by dynare_estimation_1.m (Dynare).\n'); 745 fprintf(fidTeX,['%% ' datestr(now,0) '\n']); 746 fprintf(fidTeX,' \n'); 747 end 748 for plt = 1:nbplt 749 fh = dyn_figure(options_.nodisplay,'Name','Historical and smoothed variables'); 750 NAMES = []; 751 if options_.TeX, TeXNAMES = []; end 752 nstar0=min(nstar,n_varobs-(plt-1)*nstar); 753 if gend==1 754 marker_string{1,1}='-ro'; 755 marker_string{2,1}='--ko'; 756 else 757 marker_string{1,1}='-r'; 758 marker_string{2,1}='--k'; 759 end 760 for i=1:nstar0 761 k = (plt-1)*nstar+i; 762 subplot(nr,nc,i); 763 plot(1:gend,yf(k,:),marker_string{1,1},'linewidth',1) 764 hold on 765 plot(1:gend,rawdata(:,k),marker_string{2,1},'linewidth',1) 766 hold off 767 name = options_.varobs{k}; 768 if isempty(NAMES) 769 NAMES = name; 770 else 771 NAMES = char(NAMES,name); 772 end 773 if ~isempty(options_.XTick) 774 set(gca,'XTick',options_.XTick) 775 set(gca,'XTickLabel',options_.XTickLabel) 776 end 777 if gend>1 778 xlim([1 gend]) 779 end 780 if options_.TeX 781 idx = strmatch(options_.varobs{k}, M_.endo_names,'exact'); 782 texname = M_.endo_names_tex{idx}; 783 if isempty(TeXNAMES) 784 TeXNAMES = ['$ ' texname ' $']; 785 else 786 TeXNAMES = char(TeXNAMES,['$ ' texname ' $']); 787 end 788 end 789 title(name,'Interpreter','none') 790 end 791 dyn_saveas(fh,[M_.fname, '/graphs/' M_.fname '_HistoricalAndSmoothedVariables' int2str(plt)],options_.nodisplay,options_.graph_format); 792 if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format))) 793 fprintf(fidTeX,'\\begin{figure}[H]\n'); 794 for jj = 1:nstar0 795 fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:))); 796 end 797 fprintf(fidTeX,'\\centering \n'); 798 fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_HistoricalAndSmoothedVariables%s}\n',options_.figures.textwidth*min(i/nc,1),[M_.fname, '/graphs/' M_.fname],int2str(plt)); 799 fprintf(fidTeX,'\\caption{Historical and smoothed variables.}'); 800 fprintf(fidTeX,'\\label{Fig:HistoricalAndSmoothedVariables:%s}\n',int2str(plt)); 801 fprintf(fidTeX,'\\end{figure}\n'); 802 fprintf(fidTeX,'\n'); 803 end 804 end 805 if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format))) 806 fprintf(fidTeX,'\n'); 807 fprintf(fidTeX,'%% End of TeX file.\n'); 808 fclose(fidTeX); 809 end 810 end 811end 812 813if options_.forecast > 0 && options_.mh_replic == 0 && ~options_.load_mh_file 814 oo_.forecast = dyn_forecast(var_list_,M_,options_,oo_,'smoother',dataset_info); 815end 816 817if np > 0 818 pindx = estim_params_.param_vals(:,1); 819 save([M_.fname '_pindx.mat'] ,'pindx'); 820end 821 822%reset qz_criterium 823options_.qz_criterium=qz_criterium_old; 824 825if reset_options_related_to_estimation 826 options_.mode_compute = mode_compute_o; 827 options_.mh_replic = mh_replic_o; 828end 829if first_obs_nan_indicator 830 options_.first_obs=NaN; 831end 832