dynare_estimation_1.m 37.3 KB
Newer Older
1
2
3
function dynare_estimation_1(var_list_,dname)
% function dynare_estimation_1(var_list_,dname)
% runs the estimation of the model
4
%
5
6
7
% INPUTS
%   var_list_:  selected endogenous variables vector
%   dname:      alternative directory name
8
%
9
10
11
12
13
14
% OUTPUTS
%   none
%
% SPECIAL REQUIREMENTS
%   none

15
% Copyright (C) 2003-2018 Dynare Team
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
%
% 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/>.

32
global M_ options_ oo_ estim_params_ bayestopt_ dataset_ dataset_info
33

34
35
36
37
38
39
40
41
42
43
44
if isempty(estim_params_)
    mode_compute_o = options_.mode_compute;
    mh_replic_o = options_.mh_replic;
    options_.mode_compute = 0;
    options_.mh_replic = 0;
    reset_options_related_to_estimation = true;
else
    reset_options_related_to_estimation = false;
end


45
46
%store qz_criterium
qz_criterium_old=options_.qz_criterium;
47
48
49
50
51
if isnan(options_.first_obs)
    first_obs_nan_indicator=true;
else
    first_obs_nan_indicator=false;
end
52

53
54
% Set particle filter flag.
if 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')
Stéphane Adjemian's avatar
Stéphane Adjemian committed
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
66
            fprintf('\t 7 - Nelder-Mead simplex based optimization routine (Matlab optimization toolbox required)\n')
            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
69
70
71
72
73
74
75
76
77
78
            choice = [];
            while isempty(choice)
                choice = input('Please enter your choice: ');
                if isnumeric(choice) && isint(choice) && ismember(choice,[0 6 7 8 9])
                    if choice
                        options_.mode_compute = choice;
                    end
                else
                    fprintf('\nThis is an invalid choice (you have to choose between 0, 6, 7, 8 and 9).\n')
                    choice = [];
                end
79
80
            end
        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
84
85
end

86
if ~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
92
93
94
95
96
97
98
        elseif strcmpi(options_.particle.filter_algorithm, 'apf')
            options_.particle.algorithm = 'auxiliary_particle_filter';
        elseif strcmpi(options_.particle.filter_algorithm, 'gf')
            options_.particle.algorithm = 'gaussian_filter';
        elseif strcmpi(options_.particle.filter_algorithm,  'gmf')
            options_.particle.algorithm = 'gaussian_mixture_filter';
        elseif strcmpi(options_.particle.filter_algorithm, 'cpf')
            options_.particle.algorithm = 'conditional_particle_filter';
99
100
        elseif strcmpi(options_.particle.filter_algorithm, 'nlkf')
            options_.particle.algorithm = 'nonlinear_kalman_filter';
101
        else
102
            error(['Estimation: Unknown filter ' options_.particle.filter_algorithm])
103
        end
104
105
106
    else
        objective_function = str2func('dsge_likelihood');
    end
107
else
108
    objective_function = str2func('dsge_var_likelihood');
109
110
end

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

114
115
if options_.dsge_var
    check_dsge_var_model(M_, estim_params_, bayestopt_);
116
117
118
119
120
121
122
123
124
    if dataset_info.missing.state
        error('Estimation::DsgeVarLikelihood: I cannot estimate a DSGE-VAR model with missing observations!')
    end
    if options_.noconstant
        var_sample_moments(options_.dsge_varlag, -1, dataset_);
    else
        % The steady state is non zero ==> a constant in the VAR is needed!
        var_sample_moments(options_.dsge_varlag, 0, dataset_);
    end
125
126
end

127
% Set sigma_e_is_diagonal flag (needed if the shocks block is not declared in the mod file).
128
M_.sigma_e_is_diagonal = true;
129
if estim_params_.ncx || any(nnz(tril(M_.Correlation_matrix,-1))) || isfield(estim_params_,'calibrated_covariances')
130
    M_.sigma_e_is_diagonal = false;
131
132
end

133
data = dataset_.data;
134
rawdata = dataset_info.rawdata;
135
136
data_index = dataset_info.missing.aindex;
missing_value = dataset_info.missing.state;
137

138
% Set number of observations
139
gend = dataset_.nobs;
140

141
% Set the number of observed variables.
142
n_varobs = length(options_.varobs);
143

144
% Get the number of parameters to be estimated.
145
146
147
148
149
150
nvx = estim_params_.nvx;  % Variance of the structural innovations (number of parameters).
nvn = estim_params_.nvn;  % Variance of the measurement innovations (number of parameters).
ncx = estim_params_.ncx;  % Covariance of the structural innovations (number of parameters).
ncn = estim_params_.ncn;  % Covariance of the measurement innovations (number of parameters).
np  = estim_params_.np ;  % Number of deep parameters.
nx  = nvx+nvn+ncx+ncn+np; % Total number of parameters to be estimated.
151

152
dr = oo_.dr;
153
154

if ~isempty(estim_params_)
155
    M_ = set_all_parameters(xparam1,estim_params_,M_);
156
end
157

158

159
160
%% perform initial estimation checks;
try
161
    oo_ = initial_estimation_checks(objective_function,xparam1,dataset_,dataset_info,M_,estim_params_,options_,bayestopt_,bounds,oo_);
162
163
catch % if check fails, provide info on using calibration if present
    e = lasterror();
164
    if estim_params_.full_calibration_detected %calibrated model present and no explicit starting values
165
166
        skipline(1);
        fprintf('ESTIMATION_CHECKS: There was an error in computing the likelihood for initial parameter values.\n')
167
168
        fprintf('ESTIMATION_CHECKS: If this is not a problem with the setting of options (check the error message below),\n')
        fprintf('ESTIMATION_CHECKS: you should try using the calibrated version of the model as starting values. To do\n')
169
        fprintf('ESTIMATION_CHECKS: this, add an empty estimated_params_init-block with use_calibration option immediately before the estimation\n')
170
171
172
        fprintf('ESTIMATION_CHECKS: command (and after the estimated_params-block so that it does not get overwritten):\n');
        skipline(2);
    end
173
    rethrow(e);
174
end
175

176
if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && options_.mh_posterior_mode_estimation==0
177
178
179
180
181
182
183
184
185
186
187
188
189
190
    if options_.order==1 && ~options_.particle.status
        if options_.smoother
            [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_);
            [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);
            if options_.forecast > 0
                oo_.forecast = dyn_forecast(var_list_,M_,options_,oo_,'smoother',dataset_info);
            end
            %reset qz_criterium
            options_.qz_criterium=qz_criterium_old;
            return
        end
    else %allow to continue, e.g. with MCMC_jumping_covariance
        if options_.smoother
            error('Estimation:: Particle Smoothers are not yet implemented.')
191
        end
192
193
194
    end
end

195
196
%% Estimation of the posterior mode or likelihood mode

197
if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
198
199
    %prepare settings for newrat
    if options_.mode_compute==5
200
        %get whether outer product Hessian is requested
201
        newratflag=[];
202
        if ~isempty(options_.optim_opt)
203
204
            options_list = read_key_value_string(options_.optim_opt);
            for i=1:rows(options_list)
205
206
                if strcmp(options_list{i,1},'Hessian')
                    newratflag=options_list{i,2};
207
208
209
                end
            end
        end
210
        if options_.analytic_derivation
211
            options_analytic_derivation_old = options_.analytic_derivation;
212
            options_.analytic_derivation = -1;
213
            if ~isempty(newratflag) && newratflag~=0 %numerical hessian explicitly specified
214
215
                error('newrat: analytic_derivation is incompatible with numerical Hessian.')
            else %use default
216
                newratflag=0; %exclude DYNARE numerical hessian
217
            end
218
219
220
        elseif ~options_.analytic_derivation
            if isempty(newratflag)
                newratflag=options_.newrat.hess; %use default numerical dynare hessian
221
222
            end
        end
223
    end
224

225
    [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_);
226
227
    fprintf('\nFinal value of minus the log posterior (or likelihood):%f \n', fval);

228
    if isnumeric(options_.mode_compute) && options_.mode_compute==5 && options_.analytic_derivation==-1 %reset options changed by newrat
229
        options_.analytic_derivation = options_analytic_derivation_old; %reset
230
    elseif isnumeric(options_.mode_compute) && options_.mode_compute==6 %save scaling factor
231
232
        save([M_.fname '_optimal_mh_scale_parameter.mat'],'Scale');
        options_.mh_jscale = Scale;
233
        bayestopt_.jscale(:) = options_.mh_jscale;
234
    end
235
    if ~isnumeric(options_.mode_compute) || ~isequal(options_.mode_compute,6) %always already computes covariance matrix
236
        if options_.cova_compute == 1 %user did not request covariance not to be computed
237
            if options_.analytic_derivation && strcmp(func2str(objective_function),'dsge_likelihood')
238
                ana_deriv_old = options_.analytic_derivation;
239
                options_.analytic_derivation = 2;
240
241
                [~,~,~,~,hh] = feval(objective_function,xparam1, ...
                                     dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
242
                options_.analytic_derivation = ana_deriv_old;
Stéphane Adjemian's avatar
Stéphane Adjemian committed
243
            elseif ~isnumeric(options_.mode_compute) || ~(isequal(options_.mode_compute,5) && newratflag~=1 && strcmp(func2str(objective_function),'dsge_likelihood'))
244
245
246
247
248
249
250
251
                % with flag==0, we force to use the hessian from outer product gradient of optimizer 5
                if options_.hessian.use_penalized_objective
                    penalized_objective_function = str2func('penalty_objective_function');
                    hh = hessian(penalized_objective_function, xparam1, options_.gstep, objective_function, fval, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_);
                else
                    hh = hessian(objective_function, xparam1, options_.gstep, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_);
                end
                hh = reshape(hh, nx, nx);
252
            elseif isnumeric(options_.mode_compute) && isequal(options_.mode_compute,5)
253
254
255
256
257
258
259
260
261
262
263
264
265
                % other numerical hessian options available with optimizer 5
                %
                % if newratflag == 0
                % compute outer product gradient of optimizer 5
                %
                % if newratflag == 2
                % compute 'mixed' outer product gradient of optimizer 5
                % with diagonal elements computed with numerical second order derivatives
                %
                % uses univariate filters, so to get max # of available
                % densitities for outer product gradient
                kalman_algo0 = options_.kalman_algo;
                compute_hessian = 1;
266
                if ~((options_.kalman_algo == 2) || (options_.kalman_algo == 4))
267
                    options_.kalman_algo=2;
268
                    if options_.lik_init == 3
269
270
                        options_.kalman_algo=4;
                    end
271
                elseif newratflag==0 % hh already contains outer product gradient with univariate filter
272
                    compute_hessian = 0;
273
                end
274
                if compute_hessian
275
276
                    crit = options_.newrat.tolerance.f;
                    newratflag = newratflag>0;
277
                    hh = reshape(mr_hessian(xparam1,objective_function,fval,newratflag,crit,new_rat_hess_info,[bounds.lb bounds.ub],bayestopt_.p2,dataset_, dataset_info, options_,M_,estim_params_,bayestopt_,bounds,oo_), nx, nx);
278
279
                end
                options_.kalman_algo = kalman_algo0;
280
            end
281
        end
282
283
    end
    parameter_names = bayestopt_.name;
284
    if options_.cova_compute || options_.mode_compute==5 || options_.mode_compute==6
285
        save([M_.fname '_mode.mat'],'xparam1','hh','parameter_names','fval');
286
    else
287
        save([M_.fname '_mode.mat'],'xparam1','parameter_names','fval');
288
    end
289
290
end

291
if ~options_.mh_posterior_mode_estimation && options_.cova_compute
292
293
294
    try
        chol(hh);
    catch
295
        skipline()
296
297
298
        disp('POSTERIOR KERNEL OPTIMIZATION PROBLEM!')
        disp(' (minus) the hessian matrix at the "mode" is not positive definite!')
        disp('=> posterior variance of the estimated parameters are not positive.')
299
        disp('You should try to change the initial values of the parameters using')
300
        disp('the estimated_params_init block, or use another optimization routine.')
301
        params_at_bound=find(abs(xparam1-bounds.ub)<1.e-10 | abs(xparam1-bounds.lb)<1.e-10);
302
303
        if ~isempty(params_at_bound)
            for ii=1:length(params_at_bound)
304
                params_at_bound_name{ii,1}=get_the_name(params_at_bound(ii),0,M_,estim_params_,options_);
305
306
307
308
309
310
311
312
313
314
            end
            disp_string=[params_at_bound_name{1,:}];
            for ii=2:size(params_at_bound_name,1)
                disp_string=[disp_string,', ',params_at_bound_name{ii,:}];
            end
            fprintf('\nThe following parameters are at the prior bound: %s\n', disp_string)
            fprintf('Some potential solutions are:\n')
            fprintf('   - Check your model for mistakes.\n')
            fprintf('   - Check whether model and data are consistent (correct observation equation).\n')
            fprintf('   - Shut off prior_trunc.\n')
315
            fprintf('   - Change the optimization bounds.\n')
316
317
            fprintf('   - Use a different mode_compute like 6 or 9.\n')
            fprintf('   - Check whether the parameters estimated are identified.\n')
318
            fprintf('   - Check prior shape (e.g. Inf density at bound(s)).\n')
319
320
            fprintf('   - Increase the informativeness of the prior.\n')
        end
321
322
        warning('The results below are most likely wrong!');
    end
323
324
end

325
if options_.mode_check.status && ~options_.mh_posterior_mode_estimation
326
    ana_deriv_old = options_.analytic_derivation;
327
    options_.analytic_derivation = 0;
328
    mode_check(objective_function,xparam1,hh,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
329
    options_.analytic_derivation = ana_deriv_old;
330
331
end

332
oo_.posterior.optimization.mode = [];
333
oo_.posterior.optimization.Variance = [];
334
335
oo_.posterior.optimization.log_density=[];

336
invhess=[];
337
if ~options_.mh_posterior_mode_estimation
338
    oo_.posterior.optimization.mode = xparam1;
339
340
341
    if exist('fval','var')
        oo_.posterior.optimization.log_density=-fval;
    end
342
    if options_.cova_compute
343
344
        hsd = sqrt(diag(hh));
        invhess = inv(hh./(hsd*hsd'))./(hsd*hsd');
345
        stdh = sqrt(diag(invhess));
346
        oo_.posterior.optimization.Variance = invhess;
347
    end
348
else
349
350
351
352
353
354
355
    variances = bayestopt_.p2.*bayestopt_.p2;
    idInf = isinf(variances);
    variances(idInf) = 1;
    invhess = options_.mh_posterior_mode_estimation*diag(variances);
    xparam1 = bayestopt_.p5;
    idNaN = isnan(xparam1);
    xparam1(idNaN) = bayestopt_.p1(idNaN);
356
357
    outside_bound_pars=find(xparam1 < bounds.lb | xparam1 > bounds.ub);
    xparam1(outside_bound_pars) = bayestopt_.p1(outside_bound_pars);
358
359
end

360
361
362
if ~options_.cova_compute
    stdh = NaN(length(xparam1),1);
end
363

364
if any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
365
    % display results table and store parameter estimates and standard errors in results
366
    oo_ = display_estimation_results_table(xparam1, stdh, M_, options_, estim_params_, bayestopt_, oo_, prior_dist_names, 'Posterior', 'posterior');
367
    % Laplace approximation to the marginal log density:
368
369
    if options_.cova_compute
        estim_params_nbr = size(xparam1,1);
370
371
372
373
374
375
376
        if ispd(invhess)
            log_det_invhess = log(det(invhess./(stdh*stdh')))+2*sum(log(stdh));
            likelihood = feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
            oo_.MarginalDensity.LaplaceApproximation = .5*estim_params_nbr*log(2*pi) + .5*log_det_invhess - likelihood;
        else
            oo_.MarginalDensity.LaplaceApproximation = NaN;
        end
377
        skipline()
378
        disp(sprintf('Log data density [Laplace approximation] is %f.',oo_.MarginalDensity.LaplaceApproximation))
379
        skipline()
380
    end
381
    if options_.dsge_var
382
        [~,~,~,~,~,~,~,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] =...
383
384
385
            feval(objective_function,xparam1,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
    end

386
elseif ~any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
387
    oo_=display_estimation_results_table(xparam1, stdh, M_, options_, estim_params_, bayestopt_, oo_, prior_dist_names, 'Maximum Likelihood', 'mle');
388
389
390
391
392
393
394
end

if np > 0
    pindx = estim_params_.param_vals(:,1);
    save([M_.fname '_params.mat'],'pindx');
end

395
switch options_.MCMC_jumping_covariance
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
  case 'hessian' %Baseline
                 %do nothing and use hessian from mode_compute
  case 'prior_variance' %Use prior variance
    if any(isinf(bayestopt_.p2))
        error('Infinite prior variances detected. You cannot use the prior variances as the proposal density, if some variances are Inf.')
    else
        hh = diag(1./(bayestopt_.p2.^2));
    end
    hsd = sqrt(diag(hh));
    invhess = inv(hh./(hsd*hsd'))./(hsd*hsd');
  case 'identity_matrix' %Use identity
    invhess = eye(nx);
  otherwise %user specified matrix in file
    try
        load(options_.MCMC_jumping_covariance,'jumping_covariance')
        hh=jumping_covariance;
    catch
        error(['No matrix named ''jumping_covariance'' could be found in ',options_.MCMC_jumping_covariance,'.mat'])
    end
    [nrow, ncol]=size(hh);
    if ~isequal(nrow,ncol) && ~isequal(nrow,nx) %check if square and right size
        error(['jumping_covariance matrix must be square and have ',num2str(nx),' rows and columns'])
    end
    try %check for positive definiteness
        chol(hh);
421
422
        hsd = sqrt(diag(hh));
        invhess = inv(hh./(hsd*hsd'))./(hsd*hsd');
423
424
425
    catch
        error(['Specified jumping_covariance is not positive definite'])
    end
426
427
end

428
429
if (any(bayestopt_.pshape  >0 ) && options_.mh_replic) || ...
        (any(bayestopt_.pshape >0 ) && options_.load_mh_file)  %% not ML estimation
430
    bounds = prior_bounds(bayestopt_, options_.prior_trunc); %reset bounds as lb and ub must only be operational during mode-finding
431
    outside_bound_pars=find(xparam1 < bounds.lb | xparam1 > bounds.ub);
432
433
434
435
436
437
438
    if ~isempty(outside_bound_pars)
        for ii=1:length(outside_bound_pars)
            outside_bound_par_names{ii,1}=get_the_name(ii,0,M_,estim_params_,options_);
        end
        disp_string=[outside_bound_par_names{1,:}];
        for ii=2:size(outside_bound_par_names,1)
            disp_string=[disp_string,', ',outside_bound_par_names{ii,:}];
439
        end
440
441
442
443
444
        if options_.prior_trunc>0
            error(['Estimation:: Mode value(s) of ', disp_string ,' are outside parameter bounds. Potentially, you should set prior_trunc=0.'])
        else
            error(['Estimation:: Mode value(s) of ', disp_string ,' are outside parameter bounds.'])
        end
445
    end
446
447
    % Tunes the jumping distribution's scale parameter
    if options_.mh_tune_jscale.status
448
        if strcmp(options_.posterior_sampler_options.posterior_sampling_method, 'random_walk_metropolis_hastings')
449
450
451
            options = options_.mh_tune_jscale;
            options.rwmh = options_.posterior_sampler_options.rwmh;
            options_.mh_jscale = calibrate_mh_scale_parameter(objective_function, ...
452
453
454
455
456
                                                              invhess, xparam1, [bounds.lb,bounds.ub], ...
                                                              options, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds, oo_);
            bayestopt_.jscale(:) = options_.mh_jscale;
            disp(sprintf('mh_jscale has been set equal to %s', num2str(options_.mh_jscale)))
            skipline()
457
458
459
460
        else
            warning('mh_tune_jscale is only available with Random Walk Metropolis Hastings!')
        end
    end
461
    % runs MCMC
462
    if options_.mh_replic || options_.load_mh_file
463
        posterior_sampler_options = options_.posterior_sampler_options.current_options;
464
465
        posterior_sampler_options.invhess = invhess;
        [posterior_sampler_options, options_] = check_posterior_sampler_options(posterior_sampler_options, options_);
466
467
        % store current options in global
        options_.posterior_sampler_options.current_options = posterior_sampler_options;
468
469
470
471
472
473
        if options_.mh_replic
            ana_deriv_old = options_.analytic_derivation;
            options_.analytic_derivation = 0;
            posterior_sampler(objective_function,posterior_sampler_options.proposal_distribution,xparam1,posterior_sampler_options,bounds,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
            options_.analytic_derivation = ana_deriv_old;
        end
474
    end
475
476
    %% Here I discard first mh_drop percent of the draws:
    CutSample(M_, options_, estim_params_);
477
    if options_.mh_posterior_mode_estimation
478
479
        %reset qz_criterium
        options_.qz_criterium=qz_criterium_old;
480
        return
481
    else
482
483
        %get stored results if required
        if options_.load_mh_file && options_.load_results_after_load_mh
484
            oo_load_mh=load([M_.fname '_results'],'oo_');
485
        end
486
        if ~options_.nodiagnostic
487
488
489
490
491
492
493
            if (options_.mh_replic>0 || (options_.load_mh_file && ~options_.load_results_after_load_mh))
                oo_= McMCDiagnostics(options_, estim_params_, M_,oo_);
            elseif options_.load_mh_file && options_.load_results_after_load_mh
                if isfield(oo_load_mh.oo_,'convergence')
                    oo_.convergence=oo_load_mh.oo_.convergence;
                end
            end
494
        end
495
        %% Estimation of the marginal density from the Mh draws:
496
        if options_.mh_replic || (options_.load_mh_file && ~options_.load_results_after_load_mh)
497
            [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_, bayestopt_);
MichelJuillard's avatar
MichelJuillard committed
498
            % Store posterior statistics by parameter name
499
            oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_, prior_dist_names);
500
501
502
            if ~options_.nograph
                oo_ = PlotPosteriorDistributions(estim_params_, M_, options_, bayestopt_, oo_);
            end
MichelJuillard's avatar
MichelJuillard committed
503
504
            % Store posterior mean in a vector and posterior variance in
            % a matrix
505
            [oo_.posterior.metropolis.mean,oo_.posterior.metropolis.Variance] ...
506
                = GetPosteriorMeanVariance(M_,options_.mh_drop);
507
        elseif options_.load_mh_file && options_.load_results_after_load_mh
508
            %% load fields from previous MCMC run stored in results-file
509
            field_names={'posterior_mode','posterior_std_at_mode',...% fields set by marginal_density
510
511
512
                         'posterior_mean','posterior_hpdinf','posterior_hpdsup','posterior_median','posterior_variance','posterior_std','posterior_deciles','posterior_density',...% fields set by GetPosteriorParametersStatistics
                         'prior_density',...%fields set by PlotPosteriorDistributions
                        };
513
514
515
516
            for field_iter=1:size(field_names,2)
                if isfield(oo_load_mh.oo_,field_names{1,field_iter})
                    oo_.(field_names{1,field_iter})=oo_load_mh.oo_.(field_names{1,field_iter});
                end
517
            end
518
519
520
            % field set by marginal_density
            if isfield(oo_load_mh.oo_,'MarginalDensity') && isfield(oo_load_mh.oo_.MarginalDensity,'ModifiedHarmonicMean')
                oo_.MarginalDensity.ModifiedHarmonicMean=oo_load_mh.oo_.MarginalDensity.ModifiedHarmonicMean;
521
            end
522
523
524
525
            % field set by GetPosteriorMeanVariance
            if isfield(oo_load_mh.oo_,'posterior') && isfield(oo_load_mh.oo_.posterior,'metropolis')
                oo_.posterior.metropolis=oo_load_mh.oo_.posterior.metropolis;
            end
526
        end
527
        [error_flag,~,options_]= metropolis_draw(1,options_,estim_params_,M_);
528
529
530
531
532
533
        if ~(~isempty(options_.sub_draws) && options_.sub_draws==0)
            if options_.bayesian_irf
                if error_flag
                    error('Estimation::mcmc: I cannot compute the posterior IRFs!')
                end
                PosteriorIRF('posterior');
Stéphane Adjemian's avatar
Stéphane Adjemian committed
534
            end
535
536
537
538
539
            if options_.moments_varendo
                if error_flag
                    error('Estimation::mcmc: I cannot compute the posterior moments for the endogenous variables!')
                end
                oo_ = compute_moments_varendo('posterior',options_,M_,oo_,var_list_);
Stéphane Adjemian's avatar
Stéphane Adjemian committed
540
            end
541
542
543
544
            if options_.smoother || ~isempty(options_.filter_step_ahead) || options_.forecast
                if error_flag
                    error('Estimation::mcmc: I cannot compute the posterior statistics!')
                end
545
546
547
548
549
                if options_.order==1 && ~options_.particle.status
                    prior_posterior_statistics('posterior',dataset_,dataset_info); %get smoothed and filtered objects and forecasts
                else
                    error('Estimation::mcmc: Particle Smoothers are not yet implemented.')
                end
Stéphane Adjemian's avatar
Stéphane Adjemian committed
550
            end
551
552
        else
            fprintf('Estimation:mcmc: sub_draws was set to 0. Skipping posterior computations.')
553
        end
554
        xparam1 = get_posterior_parameters('mean',M_,estim_params_,oo_,options_);
555
        M_ = set_all_parameters(xparam1,estim_params_,M_);
556
557
    end
end
558

559
if options_.particle.status
560
561
    %reset qz_criterium
    options_.qz_criterium=qz_criterium_old;
562
563
564
    return
end

565
if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.pshape> 0) && options_.load_mh_file)) ...
566
    || ~options_.smoother ) && ~options_.partial_information  % to be fixed
567
    %% ML estimation, or posterior mode without Metropolis-Hastings or Metropolis without Bayesian smoothes variables
568
    [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_);
569
    [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);
570

571
    if ~options_.nograph
572
        [nbplt,nr,nc,lr,lc,nstar] = pltorg(M_.exo_nbr);
573
574
        if ~exist([M_.dname '/graphs'],'dir')
            mkdir(M_.dname,'graphs');
575
        end
576
        if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
577
            fidTeX = fopen([M_.dname, '/graphs/' M_.fname '_SmoothedShocks.tex'],'w');
578
            fprintf(fidTeX,'%% TeX eps-loader file generated by dynare_estimation_1.m (Dynare).\n');
579
580
            fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
            fprintf(fidTeX,' \n');
581
        end
582
        for plt = 1:nbplt
583
            fh = dyn_figure(options_.nodisplay,'Name','Smoothed shocks');
584
585
            NAMES = [];
            if options_.TeX, TeXNAMES = []; end
586
            nstar0=min(nstar,M_.exo_nbr-(plt-1)*nstar);
587
588
589
590
591
592
593
            if gend==1
                marker_string{1,1}='-ro';
                marker_string{2,1}='-ko';
            else
                marker_string{1,1}='-r';
                marker_string{2,1}='-k';
            end
594
            for i=1:nstar0
595
596
                k = (plt-1)*nstar+i;
                subplot(nr,nc,i);
597
                plot([1 gend],[0 0],marker_string{1,1},'linewidth',.5)
598
                hold on
599
                plot(1:gend,innov(k,:),marker_string{2,1},'linewidth',1)
600
                hold off
601
                name = M_.exo_names{k};
602
603
604
                if ~isempty(options_.XTick)
                    set(gca,'XTick',options_.XTick)
                    set(gca,'XTickLabel',options_.XTickLabel)
605
                end
606
607
608
                if gend>1
                    xlim([1 gend])
                end
609
                if options_.TeX
610
611
612
                    title(['$' M_.exo_names_tex{k} '$'],'Interpreter','latex')
                else
                    title(name,'Interpreter','none')
613
                end
614
            end
615
            dyn_saveas(fh,[M_.dname, '/graphs/' M_.fname '_SmoothedShocks' int2str(plt)],options_.nodisplay,options_.graph_format);
616
            if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
617
618
                fprintf(fidTeX,'\\begin{figure}[H]\n');
                fprintf(fidTeX,'\\centering \n');
619
                fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_SmoothedShocks%s}\n',options_.figures.textwidth*min(i/nc,1),[M_.dname, '/graphs/' M_.fname],int2str(plt));
620
621
622
623
                fprintf(fidTeX,'\\caption{Smoothed shocks.}');
                fprintf(fidTeX,'\\label{Fig:SmoothedShocks:%s}\n',int2str(plt));
                fprintf(fidTeX,'\\end{figure}\n');
                fprintf(fidTeX,'\n');
624
            end
625
        end
626
        if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
627
628
629
            fprintf(fidTeX,'\n');
            fprintf(fidTeX,'%% End of TeX file.\n');
            fclose(fidTeX);
630
        end
631
632
633
634
    end
    if nvn
        number_of_plots_to_draw = 0;
        index = [];
635
        for obs_iter=1:n_varobs
636
            if max(abs(measurement_error(obs_iter,:))) > options_.ME_plot_tol;
637
                number_of_plots_to_draw = number_of_plots_to_draw + 1;
638
                index = cat(1,index,obs_iter);
639
640
            end
        end
641
642
        if ~options_.nograph
            [nbplt,nr,nc,lr,lc,nstar] = pltorg(number_of_plots_to_draw);
643
            if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
644
                fidTeX = fopen([M_.dname, '/graphs/' M_.fname '_SmoothedObservationErrors.tex'],'w');
645
                fprintf(fidTeX,'%% TeX eps-loader file generated by dynare_estimation_1.m (Dynare).\n');
646
647
                fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
                fprintf(fidTeX,' \n');
648
            end
649
            for plt = 1:nbplt
650
                fh = dyn_figure(options_.nodisplay,'Name','Smoothed observation errors');
651
                nstar0=min(nstar,number_of_plots_to_draw-(plt-1)*nstar);
652
653
654
655
656
657
658
                if gend==1
                    marker_string{1,1}='-ro';
                    marker_string{2,1}='-ko';
                else
                    marker_string{1,1}='-r';
                    marker_string{2,1}='-k';
                end
659
                for i=1:nstar0
660
661
                    k = (plt-1)*nstar+i;
                    subplot(nr,nc,i);
662
                    plot([1 gend],[0 0],marker_string{1,1},'linewidth',.5)
663
                    hold on
664
                    plot(1:gend,measurement_error(index(k),:),marker_string{2,1},'linewidth',1)
665
                    hold off
666
                    name = options_.varobs{index(k)};
667
668
669
                    if gend>1
                        xlim([1 gend])
                    end
670
671
672
673
674
                    if ~isempty(options_.XTick)
                        set(gca,'XTick',options_.XTick)
                        set(gca,'XTickLabel',options_.XTickLabel)
                    end
                    if options_.TeX
675
                        idx = strmatch(options_.varobs{index(k)}, M_.endo_names, 'exact');
676
677
678
                        title(['$' M_.endo_names_tex{idx} '$'],'Interpreter','latex')
                    else
                        title(name,'Interpreter','none')
679
                    end
680
                end
681
                dyn_saveas(fh,[M_.dname, '/graphs/' M_.fname '_SmoothedObservationErrors' int2str(plt)],options_.nodisplay,options_.graph_format);
682
                if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
683
684
                    fprintf(fidTeX,'\\begin{figure}[H]\n');
                    fprintf(fidTeX,'\\centering \n');
685
                    fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_SmoothedObservationErrors%s}\n',options_.figures.textwidth*min(i/nc,1),[M_.dname, '/graphs/' M_.fname],int2str(plt));
686
687
                    fprintf(fidTeX,'\\caption{Smoothed observation errors.}');
                    fprintf(fidTeX,'\\label{Fig:SmoothedObservationErrors:%s}\n',int2str(plt));
688
689
                    fprintf(fidTeX,'\\end{figure}\n');
                    fprintf(fidTeX,'\n');
690
                end
691
            end
692
            if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
693
694
695
696
                fprintf(fidTeX,'\n');
                fprintf(fidTeX,'%% End of TeX file.\n');
                fclose(fidTeX);
            end
697
        end
698
    end
699
700
701
    %%
    %%  Historical and smoothed variabes
    %%
702
    if ~options_.nograph
703
704
        [nbplt,nr,nc,lr,lc,nstar] = pltorg(n_varobs);
        if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
705
            fidTeX = fopen([M_.dname, '/graphs/' M_.fname '_HistoricalAndSmoothedVariables.tex'],'w');
706
707
708
            fprintf(fidTeX,'%% TeX eps-loader file generated by dynare_estimation_1.m (Dynare).\n');
            fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
            fprintf(fidTeX,' \n');
709
        end
710
711
712
713
714
715
716
        for plt = 1:nbplt
            fh = dyn_figure(options_.nodisplay,'Name','Historical and smoothed variables');
            NAMES = [];
            nstar0=min(nstar,n_varobs-(plt-1)*nstar);
            if gend==1
                marker_string{1,1}='-ro';
                marker_string{2,1}='--ko';
717
            else
718
719
                marker_string{1,1}='-r';
                marker_string{2,1}='--k';
720
            end
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
            for i=1:nstar0
                k = (plt-1)*nstar+i;
                subplot(nr,nc,i);
                plot(1:gend,yf(k,:),marker_string{1,1},'linewidth',1)
                hold on
                plot(1:gend,rawdata(:,k),marker_string{2,1},'linewidth',1)
                hold off
                name = options_.varobs{k};
                if ~isempty(options_.XTick)
                    set(gca,'XTick',options_.XTick)
                    set(gca,'XTickLabel',options_.XTickLabel)
                end
                if gend>1
                    xlim([1 gend])
                end
                if options_.TeX
737
                    idx = strmatch(options_.varobs{k}, M_.endo_names,'exact');
738
739
740
                    title(['$' M_.endo_names_tex{idx} '$'],'Interpreter','latex')
                else
                    title(name,'Interpreter','none')
741
                end
742
            end
743
            dyn_saveas(fh,[M_.dname, '/graphs/' M_.fname '_HistoricalAndSmoothedVariables' int2str(plt)],options_.nodisplay,options_.graph_format);
744
745
746
            if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
                fprintf(fidTeX,'\\begin{figure}[H]\n');
                fprintf(fidTeX,'\\centering \n');
747
                fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_HistoricalAndSmoothedVariables%s}\n',options_.figures.textwidth*min(i/nc,1),[M_.dname, '/graphs/' M_.fname],int2str(plt));
748
749
750
751
                fprintf(fidTeX,'\\caption{Historical and smoothed variables.}');
                fprintf(fidTeX,'\\label{Fig:HistoricalAndSmoothedVariables:%s}\n',int2str(plt));
                fprintf(fidTeX,'\\end{figure}\n');
                fprintf(fidTeX,'\n');
752
753
            end
        end
754
        if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
755
            fprintf(fidTeX,'\n');
756
757
            fprintf(fidTeX,'%% End of TeX file.\n');
            fclose(fidTeX);
758
759
        end
    end
760
end
761

762
if options_.forecast > 0 && options_.mh_replic == 0 && ~options_.load_mh_file
763
    oo_.forecast = dyn_forecast(var_list_,M_,options_,oo_,'smoother',dataset_info);
764
765
766
767
768
769
770
end

if np > 0
    pindx = estim_params_.param_vals(:,1);
    save([M_.fname '_pindx.mat'] ,'pindx');
end

771
772
%reset qz_criterium
options_.qz_criterium=qz_criterium_old;
773
774
775
776
777

if reset_options_related_to_estimation
    options_.mode_compute = mode_compute_o;
    options_.mh_replic = mh_replic_o;
end
778
779
if first_obs_nan_indicator
    options_.first_obs=NaN;
780
end