dynare_estimation_1.m 53.6 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-2013 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_
33

34
35
% Set particle filter flag.
if options_.order > 1
36
37
38
39
    if options_.particle.status && options_.order==2
        disp(' ')
        disp('Estimation using a non linear filter!')
        disp(' ')
40
    elseif options_.particle.status && options_.order>2
41
42
        error(['Non linear filter are not implemented with order ' int2str(options_.order) ' approximation of the model!'])
    elseif ~options_.particle.status && options_.order==2
43
        error('For estimating the model with a second order approximation using a non linear filter, one should have options_.particle.status=1;')
44
45
46
    else
        error(['Cannot estimate a model with an order ' int2str(options_.order) ' approximation!'])
    end
47
48
end

49
if ~options_.dsge_var
50
    if options_.particle.status
51
52
53
54
        objective_function = str2func('non_linear_dsge_likelihood');
    else
        objective_function = str2func('dsge_likelihood');
    end
55
56
57
58
else
    objective_function = str2func('DsgeVarLikelihood');
end

59
[dataset_,xparam1, M_, options_, oo_, estim_params_,bayestopt_] = dynare_estimation_init(var_list_, dname, [], M_, options_, oo_, estim_params_, bayestopt_);
60

61
% Set sigma_e_is_diagonal flag (needed if the shocks block is not declared in the mod file).
62
M_.sigma_e_is_diagonal = 1;
63
if estim_params_.ncx || ~isequal(nnz(M_.Sigma_e),length(M_.Sigma_e))
64
65
66
    M_.sigma_e_is_diagonal = 0;
end

67
68
69
% Set the correlation matrix if necessary.
if ~isequal(estim_params_.ncx,nnz(tril(M_.Sigma_e,-1)))
    M_.Correlation_matrix = diag(1./sqrt(diag(M_.Sigma_e)))*M_.Sigma_e*diag(1./sqrt(diag(M_.Sigma_e)));
70
71
72
73
74
75
76
77
    % Remove NaNs appearing because of variances calibrated to zero.
    if any(isnan(M_.Correlation_matrix))
        zero_variance_idx = find(~diag(M_.Sigma_e));
        for i=1:length(zero_variance_idx)
            M_.Correlation_matrix(zero_variance_idx(i),:) = 0;
            M_.Correlation_matrix(:,zero_variance_idx(i)) = 0;
        end
    end
78
79
end

80
81
data = dataset_.data;
rawdata = dataset_.rawdata;
MichelJuillard's avatar
MichelJuillard committed
82
83
data_index = dataset_.missing.aindex;
missing_value = dataset_.missing.state;
84

85
% Set number of observations
86
gend = options_.nobs;
87
% Set the number of observed variables.
88
n_varobs = size(options_.varobs,1);
89
% Get the number of parameters to be estimated.
90
91
92
93
94
95
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.
96
97
98
99
100
%% Set the names of the priors.
pnames = ['     ';'beta ';'gamm ';'norm ';'invg ';'unif ';'invg2'];
%% Set parameters bounds
lb = bayestopt_.lb;
ub = bayestopt_.ub;
101

102
dr = oo_.dr;
103

104
%% load mode file is necessary
105
if ~isempty(options_.mode_file) && ~options_.mh_posterior_mode_estimation
106
    load(options_.mode_file);
107
108
109
110

    if length(xparam1) ~= nx
        error([ 'ESTIMATION: the posterior mode file ' options_.mode_file ' has been generated using another specification. Please delete it and recompute the posterior mode.'])
    end
111
112
end

113
114
115
116
117
118
119
120
121
122
123
124
125
126
%% load optimal_mh_scale parameter if previous run was with
%% mode_compute=6
mh_scale_fname = [M_.fname '_optimal_mh_scale_parameter.mat'];
if exist(mh_scale_fname)
    if options_.mode_compute == 0
        tmp = load(mh_scale_fname,'Scale');
        bayestopt_.mh_jscale = tmp.Scale;
        clear tmp;
    else
        % remove the file if mode_compute ~= 0
        delete('mh_scale_fname')
    end
end

127
128
129
if ~isempty(estim_params_)
    set_parameters(xparam1);
end
130

131
% compute sample moments if needed (bvar-dsge)
132
if options_.dsge_var
133
    if dataset_.missing.state
134
135
136
137
138
        error('I cannot estimate a DSGE-VAR model with missing observations!')
    end
    if options_.noconstant
        evalin('base',...
               ['[mYY,mXY,mYX,mXX,Ydata,Xdata] = ' ...
139
                'var_sample_moments(options_.first_obs,' ...
140
                'options_.first_obs+options_.nobs-1,options_.dsge_varlag,-1,' ...
141
142
143
144
                'options_.datafile, options_.varobs,options_.xls_sheet,options_.xls_range);'])
    else% The steady state is non zero ==> a constant in the VAR is needed!
        evalin('base',['[mYY,mXY,mYX,mXX,Ydata,Xdata] = ' ...
                       'var_sample_moments(options_.first_obs,' ...
145
                       'options_.first_obs+options_.nobs-1,options_.dsge_varlag,0,' ...
146
147
148
149
150
                       'options_.datafile, options_.varobs,options_.xls_sheet,options_.xls_range);'])
    end
end


151
oo_ = initial_estimation_checks(objective_function,xparam1,dataset_,M_,estim_params_,options_,bayestopt_,oo_);
152

153
if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && options_.mh_posterior_mode_estimation==0
154
    if options_.smoother == 1
155
        [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp] = DsgeSmoother(xparam1,gend,data,data_index,missing_value);
156
157
        oo_.Smoother.SteadyState = ys;
        oo_.Smoother.TrendCoeffs = trend_coeff;
158
        if options_.filter_covariance
159
            oo_.Smoother.Variance = P;
160
        end
161
        i_endo = bayestopt_.smoother_saved_var_list;
162
        if options_.nk ~= 0
163
164
165
166
167
168
169
170
171
172
            oo_.FilteredVariablesKStepAhead = ...
                aK(options_.filter_step_ahead,i_endo,:);
            if ~isempty(PK)
                oo_.FilteredVariablesKStepAheadVariances = ...
                    PK(options_.filter_step_ahead,i_endo,i_endo,:);
            end
            if ~isempty(decomp)
                oo_.FilteredVariablesShockDecomposition = ...
                    decomp(options_.filter_step_ahead,i_endo,:,:);
            end
173
        end
174
        for i=bayestopt_.smoother_saved_var_list'
175
            i1 = dr.order_var(bayestopt_.smoother_var_list(i));
176
177
178
179
180
181
            eval(['oo_.SmoothedVariables.' deblank(M_.endo_names(i1,:)) ...
                  ' = atT(i,:)'';']);
            if options_.nk > 0
                eval(['oo_.FilteredVariables.' deblank(M_.endo_names(i1,:)) ...
                      ' = squeeze(aK(1,i,:));']);
            end
182
            eval(['oo_.UpdatedVariables.' deblank(M_.endo_names(i1,:)) ' = updated_variables(i,:)'';']);
183
184
185
186
        end
        for i=1:M_.exo_nbr
            eval(['oo_.SmoothedShocks.' deblank(M_.exo_names(i,:)) ' = innov(i,:)'';']);
        end
187
    end
188
    return
189
190
end

191

192
%% Estimation of the posterior mode or likelihood mode
193
if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
194
195
    switch options_.mode_compute
      case 1
196
197
198
199
200
201
        if exist('OCTAVE_VERSION')
            error('Option mode_compute=1 is not available under Octave')
        elseif ~user_has_matlab_license('optimization_toolbox')
            error('Option mode_compute=1 requires the Optimization Toolbox')
        end

202
203
204
205
        optim_options = optimset('display','iter','LargeScale','off', ...
                                 'MaxFunEvals',100000,'TolFun',1e-8,'TolX',1e-6);
        if isfield(options_,'optim_opt')
            eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']);
206
207
        end
        if options_.analytic_derivation,
208
            optim_options = optimset(optim_options,'GradObj','on','TolX',1e-7);
209
210
        end
            [xparam1,fval,exitflag,output,lamdba,grad,hessian_fmincon] = ...
211
                fmincon(objective_function,xparam1,[],[],[],[],lb,ub,[],optim_options,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
212
213
214
      case 2
        error('ESTIMATION: mode_compute=2 option (Lester Ingber''s Adaptive Simulated Annealing) is no longer available')
      case 3
215
216
217
218
219
220
        if exist('OCTAVE_VERSION') && ~user_has_octave_forge_package('optim')
            error('Option mode_compute=3 requires the optim package')
        elseif ~exist('OCTAVE_VERSION') && ~user_has_matlab_license('optimization_toolbox')
            error('Option mode_compute=3 requires the Optimization Toolbox')
        end

221
222
223
224
        optim_options = optimset('display','iter','MaxFunEvals',100000,'TolFun',1e-8,'TolX',1e-6);
        if isfield(options_,'optim_opt')
            eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']);
        end
225
226
227
        if options_.analytic_derivation,
            optim_options = optimset(optim_options,'GradObj','on');
        end
228
229
230
231
232
233
234
235
        if ~exist('OCTAVE_VERSION')
            [xparam1,fval,exitflag] = fminunc(objective_function,xparam1,optim_options,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
        else
            % Under Octave, use a wrapper, since fminunc() does not have a 4th arg
            func = @(x) objective_function(x, dataset_,options_,M_,estim_params_,bayestopt_,oo_);
            [xparam1,fval,exitflag] = fminunc(func,xparam1,optim_options);
        end

236
237
238
239
240
      case 4
        H0 = 1e-4*eye(nx);
        crit = 1e-7;
        nit = 1000;
        verbose = 2;
241
242
243
244
245
246
        if options_.analytic_derivation,
            analytic_grad=1;
        else
            analytic_grad=[];
        end

247
            [fval,xparam1,grad,hessian_csminwel,itct,fcount,retcodehat] = ...
248
                csminwel1(objective_function,xparam1,H0,analytic_grad,crit,nit,options_.gradient_method,options_.gradient_epsilon,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
249
250
251
252
253
254
255
256
257
258
            disp(sprintf('Objective function at mode: %f',fval))
      case 5
        if isfield(options_,'hess')
            flag = options_.hess;
        else
            flag = 1;
        end
        if isfield(options_,'ftol')
            crit = options_.ftol;
        else
259
            crit = 1.e-5;
260
        end
261
262
263
264
265
266
267
268
269
        if options_.analytic_derivation,
            analytic_grad=1;
            ana_deriv = options_.analytic_derivation;
            options_.analytic_derivation = -1;
            crit = 1.e-7;
            flag = 0;
        else
            analytic_grad=0;
        end
270
271
272
273
274
        if isfield(options_,'nit')
            nit = options_.nit;
        else
            nit=1000;
        end
275
276
277
278
        [xparam1,hh,gg,fval,invhess] = newrat(objective_function,xparam1,analytic_grad,crit,nit,flag,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
        if options_.analytic_derivation,
            options_.analytic_derivation = ana_deriv;
        end
279
280
281
        parameter_names = bayestopt_.name;
        save([M_.fname '_mode.mat'],'xparam1','hh','gg','fval','invhess','parameter_names');
      case 6
282
        fval = feval(objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
283
        OldMode = fval;
284
        if ~exist('MeanPar','var')
285
286
            MeanPar = xparam1;
        end
287
        if exist('hh','var')
288
289
290
291
            CovJump = inv(hh);
        else% The covariance matrix is initialized with the prior
            % covariance (a diagonal matrix) %%Except for infinite variances ;-)
            varinit = 'prior';
292
            if strcmpi(varinit,'prior')
293
294
295
296
297
298
                stdev = bayestopt_.p2;
                indx = find(isinf(stdev));
                stdev(indx) = ones(length(indx),1)*sqrt(10);
                vars = stdev.^2;
                CovJump = diag(vars);
            elseif strcmpi(varinit,'eye')
299
300
                vars = ones(length(bayestopt_.p2),1)*0.1;
                CovJump = diag(vars);
301
302
303
304
305
306
307
            else
                disp('gmhmaxlik :: Error!')
                return
            end
        end
        OldPostVar = CovJump;
        Scale = options_.mh_jscale;
308
        for i=1:options_.gmhmaxlik.iterations
309
            if i == 1
310
                if options_.gmhmaxlik.iterations>1
311
312
313
314
                    flag = '';
                else
                    flag = 'LastCall';
                end
315
                [xparam1,PostVar,Scale,PostMean] = ...
316
                    gmhmaxlik(objective_function,xparam1,[lb ub],options_.gmhmaxlik,Scale,flag,MeanPar,CovJump,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
317
                fval = feval(objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
318
319
                options_.mh_jscale = Scale;
                mouvement = max(max(abs(PostVar-OldPostVar)));
320
321
322
323
324
325
                disp(' ')
                disp('========================================================== ')
                disp(['   Change in the covariance matrix = ' num2str(mouvement) '.'])
                disp(['   Mode improvement = ' num2str(abs(OldMode-fval))])
                disp(['   New value of jscale = ' num2str(Scale)])
                disp('========================================================== ')
326
327
328
                OldMode = fval;
            else
                OldPostVar = PostVar;
329
                if i<options_.gmhmaxlik.iterations
330
331
332
333
                    flag = '';
                else
                    flag = 'LastCall';
                end
334
335
                [xparam1,PostVar,Scale,PostMean] = ...
                    gmhmaxlik(objective_function,xparam1,[lb ub],...
336
                              options_.gmhmaxlik,Scale,flag,PostMean,PostVar,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
337
                fval = feval(objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
338
339
                options_.mh_jscale = Scale;
                mouvement = max(max(abs(PostVar-OldPostVar)));
340
341
342
343
344
345
346
                fval = feval(objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
                disp(' ')
                disp('========================================================== ')
                disp(['   Change in the covariance matrix = ' num2str(mouvement) '.'])
                disp(['   Mode improvement = ' num2str(abs(OldMode-fval))])
                disp(['   New value of jscale = ' num2str(Scale)])
                disp('========================================================== ')
347
348
349
350
                OldMode = fval;
            end
            hh = inv(PostVar);
            save([M_.fname '_mode.mat'],'xparam1','hh');
351
            save([M_.fname '_optimal_mh_scale_parameter.mat'],'Scale');
352
            bayestopt_.jscale = ones(length(xparam1),1)*Scale;
353
        end
354
355
356
357
358
        disp(' ')
        disp(['Optimal value of the scale parameter = ' num2str(Scale)])
        disp(' ')
        disp(['Final value of the log posterior (or likelihood): ' num2str(fval)])
        disp(' ')
359
360
        parameter_names = bayestopt_.name;
        save([M_.fname '_mode.mat'],'xparam1','hh','parameter_names');
361
      case 7
362
        % Matlab's simplex (Optimization toolbox needed).
363
364
365
366
367
368
        if exist('OCTAVE_VERSION') && ~user_has_octave_forge_package('optim')
            error('Option mode_compute=7 requires the optim package')
        elseif ~exist('OCTAVE_VERSION') && ~user_has_matlab_license('optimization_toolbox')
            error('Option mode_compute=7 requires the Optimization Toolbox')
        end

369
370
371
372
        optim_options = optimset('display','iter','MaxFunEvals',1000000,'MaxIter',6000,'TolFun',1e-8,'TolX',1e-6);
        if isfield(options_,'optim_opt')
            eval(['optim_options = optimset(optim_options,' options_.optim_opt ');']);
        end
373
        [xparam1,fval,exitflag] = fminsearch(objective_function,xparam1,optim_options,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
374
      case 8
375
        % Dynare implementation of the simplex algorithm.
376
377
        [xparam1,fval,exitflag] = simplex_optimization_routine(objective_function,xparam1,options_.simplex,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
      case 9
378
379
380
        H0 = 1e-4*ones(nx,1);
        warning('off','CMAES:NonfinitenessRange');
        warning('off','CMAES:InitialSigma');
381
        [x, fval, COUNTEVAL, STOPFLAG, OUT, BESTEVER] = cmaes(func2str(objective_function),xparam1,H0,options_.cmaes,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
382
        xparam1=BESTEVER.x;
383
384
        disp(sprintf('\n Objective function at mode: %f',fval))
      case 101
385
386
387
388
        myoptions=soptions;
        myoptions(2)=1e-6; %accuracy of argument
        myoptions(3)=1e-6; %accuracy of function (see Solvopt p.29)
        myoptions(5)= 1.0;
389
        [xparam1,fval]=solvopt(xparam1,objective_function,[],myoptions,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
390
391
392
393
394
395
396
      case 102
        %simulating annealing
        %        LB=zeros(size(xparam1))-20;
        % UB=zeros(size(xparam1))+20;
        LB = xparam1 - 1;
        UB = xparam1 + 1;
        neps=10;
397
        %  Set input parameters.
398
        maxy=0;
399
        epsilon=1.0e-9;
400
401
402
403
404
405
406
        rt_=.10;
        t=15.0;
        ns=10;
        nt=10;
        maxevl=100000000;
        idisp =1;
        npar=length(xparam1);
407
408

        disp(['size of param',num2str(length(xparam1))])
409
410
        c=.1*ones(npar,1);
        %*  Set input values of the input/output parameters.*
411

412
413
        vm=1*ones(npar,1);
        disp(['number of parameters= ' num2str(npar) 'max= '  num2str(maxy) 't=  ' num2str(t)]);
414
        disp(['rt_=  '  num2str(rt_) 'eps=  '  num2str(epsilon) 'ns=  '  num2str(ns)]);
415
416
417
418
419
420
421
422
423
        disp(['nt=  '  num2str(nt) 'neps= '   num2str(neps) 'maxevl=  '  num2str(maxevl)]);
        %      disp(['iprint=   '   num2str(iprint) 'seed=   '   num2str(seed)]);
        disp '  ';
        disp '  ';
        disp(['starting values(x)  ' num2str(xparam1')]);
        disp(['initial step length(vm)  '  num2str(vm')]);
        disp(['lower bound(lb)', 'initial conditions', 'upper bound(ub)' ]);
        disp([LB xparam1 UB]);
        disp(['c vector   ' num2str(c')]);
424

425
426
        [xparam1, fval, nacc, nfcnev, nobds, ier, t, vm] = sa(objective_function,xparam1,maxy,rt_,epsilon,ns,nt ...
                                                              ,neps,maxevl,LB,UB,c,idisp ,t,vm,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
427
428
      case 'prior'
        hh = diag(bayestopt_.p2.^2);
429
430
      otherwise
        if ischar(options_.mode_compute)
431
            [xparam1, fval, retcode ] = feval(options_.mode_compute,objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
432
        else
433
            error(['dynare_estimation:: mode_compute = ' int2str(options_.mode_compute) ' option is unknown!'])
434
435
        end
    end
436
    if ~isequal(options_.mode_compute,6) && ~isequal(options_.mode_compute,'prior')
437
        if options_.cova_compute == 1
438
439
440
441
442
443
444
445
446
447
448
            if options_.analytic_derivation && strcmp(func2str(objective_function),'dsge_likelihood'),
                ana_deriv = options_.analytic_derivation;
                options_.analytic_derivation = 2;
                [junk1, junk2, hh] = feval(objective_function,xparam1, ...
                    dataset_,options_,M_,estim_params_,bayestopt_,oo_);
                options_.analytic_derivation = ana_deriv;
                
            else
                hh = reshape(hessian(objective_function,xparam1, ...
                    options_.gstep,dataset_,options_,M_,estim_params_,bayestopt_,oo_),nx,nx);
            end
449
        end
450
451
452
453
454
455
456
    end
    parameter_names = bayestopt_.name;
    if options_.cova_compute
        save([M_.fname '_mode.mat'],'xparam1','hh','parameter_names');
    else
        save([M_.fname '_mode.mat'],'xparam1','parameter_names');
    end
457
458
459
end

if options_.cova_compute == 0
460
    hh = [];%NaN(length(xparam1),length(xparam1));
461
462
end

463
if ~options_.mh_posterior_mode_estimation && options_.cova_compute
464
465
466
467
468
469
470
471
472
473
474
    try
        chol(hh);
    catch
        disp(' ')
        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.')
        disp('You should  try  to change the initial values of the parameters using')
        disp('the estimated_params_init block, or use another optimization routine.')
        warning('The results below are most likely wrong!');
    end
475
476
end

477
if options_.mode_check == 1 && ~options_.mh_posterior_mode_estimation
478
479
    ana_deriv = options_.analytic_derivation;
    options_.analytic_derivation = 0;
Marco Ratto's avatar
Marco Ratto committed
480
    mode_check(objective_function,xparam1,hh,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
481
    options_.analytic_derivation = ana_deriv;
482
483
end

484
oo_.posterior.optimization.mode = xparam1;
485
oo_.posterior.optimization.Variance = [];
486
487
488
489
if ~options_.mh_posterior_mode_estimation
    if options_.cova_compute
        invhess = inv(hh);
        stdh = sqrt(diag(invhess));
490
        oo_.posterior.optimization.Variance = invhess;
491
    end
492
else
493
494
495
496
497
498
499
500
    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);
    xparam1 = transpose(xparam1);
501
502
end

503
504
505
if ~options_.cova_compute
    stdh = NaN(length(xparam1),1);
end
506

507
if any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
508
509
510
511
512
513
    disp(' ')
    disp('RESULTS FROM POSTERIOR MAXIMIZATION')
    tstath = zeros(nx,1);
    for i = 1:nx
        tstath(i) = abs(xparam1(i))/stdh(i);
    end
514

515
    header_width = row_header_width(M_,estim_params_,bayestopt_);
516

517
518
519
520
521
522
523
524
525
526
527
528
529
530
    tit1 = sprintf('%-*s %7s %8s %7s %6s %4s %6s\n',header_width-2,' ','prior mean', ...
                   'mode','s.d.','t-stat','prior','pstdev');
    if np
        ip = nvx+nvn+ncx+ncn+1;
        disp('parameters')
        disp(tit1)
        for i=1:np
            name = bayestopt_.name{ip};
            disp(sprintf('%-*s %7.3f %8.4f %7.4f %7.4f %4s %6.4f', ...
                         header_width,name, ...
                         bayestopt_.p1(ip),xparam1(ip),stdh(ip),tstath(ip), ...
                         pnames(bayestopt_.pshape(ip)+1,:), ...
                         bayestopt_.p2(ip)));
            eval(['oo_.posterior_mode.parameters.' name ' = xparam1(ip);']);
531
            eval(['oo_.posterior_std.parameters.' name ' = stdh(ip);']);
532
533
534
535
536
537
538
539
540
541
542
543
544
            ip = ip+1;
        end
    end
    if nvx
        ip = 1;
        disp('standard deviation of shocks')
        disp(tit1)
        for i=1:nvx
            k = estim_params_.var_exo(i,1);
            name = deblank(M_.exo_names(k,:));
            disp(sprintf('%-*s %7.3f %8.4f %7.4f %7.4f %4s %6.4f', ...
                         header_width,name,bayestopt_.p1(ip),xparam1(ip), ...
                         stdh(ip),tstath(ip),pnames(bayestopt_.pshape(ip)+1,:), ...
545
                         bayestopt_.p2(ip)));
546
547
            M_.Sigma_e(k,k) = xparam1(ip)*xparam1(ip);
            eval(['oo_.posterior_mode.shocks_std.' name ' = xparam1(ip);']);
548
            eval(['oo_.posterior_std.shocks_std.' name ' = stdh(ip);']);
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
            ip = ip+1;
        end
    end
    if nvn
        disp('standard deviation of measurement errors')
        disp(tit1)
        ip = nvx+1;
        for i=1:nvn
            name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
            disp(sprintf('%-*s %7.3f %8.4f %7.4f %7.4f %4s %6.4f', ...
                         header_width,name,bayestopt_.p1(ip), ...
                         xparam1(ip),stdh(ip),tstath(ip), ...
                         pnames(bayestopt_.pshape(ip)+1,:), ...
                         bayestopt_.p2(ip)));
            eval(['oo_.posterior_mode.measurement_errors_std.' name ' = xparam1(ip);']);
564
            eval(['oo_.posterior_std.measurement_errors_std.' name ' = stdh(ip);']);
565
566
567
568
569
570
571
572
573
574
575
576
            ip = ip+1;
        end
    end
    if ncx
        disp('correlation of shocks')
        disp(tit1)
        ip = nvx+nvn+1;
        for i=1:ncx
            k1 = estim_params_.corrx(i,1);
            k2 = estim_params_.corrx(i,2);
            name = [deblank(M_.exo_names(k1,:)) ',' deblank(M_.exo_names(k2,:))];
            NAME = [deblank(M_.exo_names(k1,:)) '_' deblank(M_.exo_names(k2,:))];
577
578
            disp(sprintf('%-*s %7.3f %8.4f %7.4f %7.4f %4s %6.4f', ...
                         header_width,name,bayestopt_.p1(ip),xparam1(ip),stdh(ip),tstath(ip),  ...
579
580
581
582
                         pnames(bayestopt_.pshape(ip)+1,:), bayestopt_.p2(ip)));
            M_.Sigma_e(k1,k2) = xparam1(ip)*sqrt(M_.Sigma_e(k1,k1)*M_.Sigma_e(k2,k2));
            M_.Sigma_e(k2,k1) = M_.Sigma_e(k1,k2);
            eval(['oo_.posterior_mode.shocks_corr.' NAME ' = xparam1(ip);']);
583
            eval(['oo_.posterior_std.shocks_corr.' NAME ' = stdh(ip);']);
584
585
586
587
588
589
590
591
592
593
594
595
            ip = ip+1;
        end
    end
    if ncn
        disp('correlation of measurement errors')
        disp(tit1)
        ip = nvx+nvn+ncx+1;
        for i=1:ncn
            k1 = estim_params_.corrn(i,1);
            k2 = estim_params_.corrn(i,2);
            name = [deblank(M_.endo_names(k1,:)) ',' deblank(M_.endo_names(k2,:))];
            NAME = [deblank(M_.endo_names(k1,:)) '_' deblank(M_.endo_names(k2,:))];
596
597
            disp(sprintf('%-*s %7.3f %8.4f %7.4f %7.4f %4s %6.4f', ...
                         header_width,name,bayestopt_.p1(ip),xparam1(ip),stdh(ip),tstath(ip), ...
598
599
                         pnames(bayestopt_.pshape(ip)+1,:), bayestopt_.p2(ip)));
            eval(['oo_.posterior_mode.measurement_errors_corr.' NAME ' = xparam1(ip);']);
600
            eval(['oo_.posterior_std.measurement_errors_corr.' NAME ' = stdh(ip);']);
601
602
603
604
            ip = ip+1;
        end
    end
    %% Laplace approximation to the marginal log density:
605
606
607
608
    if options_.cova_compute
        estim_params_nbr = size(xparam1,1);
        scale_factor = -sum(log10(diag(invhess)));
        log_det_invhess = -estim_params_nbr*log(scale_factor)+log(det(scale_factor*invhess));
609
610
        likelihood = feval(objective_function,xparam1,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
        oo_.MarginalDensity.LaplaceApproximation = .5*estim_params_nbr*log(2*pi) + .5*log_det_invhess - likelihood;
611
        disp(' ')
612
        disp(sprintf('Log data density [Laplace approximation] is %f.',oo_.MarginalDensity.LaplaceApproximation))
613
        disp(' ')
614
    end
615
elseif ~any(bayestopt_.pshape > 0) && ~options_.mh_posterior_mode_estimation
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
    disp(' ')
    disp('RESULTS FROM MAXIMUM LIKELIHOOD')
    tstath = zeros(nx,1);
    for i = 1:nx
        tstath(i) = abs(xparam1(i))/stdh(i);
    end
    header_width = row_header_width(M_,estim_params_,bayestopt_);
    tit1 = sprintf('%-*s %10s %7s %6s\n',header_width-2,' ','Estimate','s.d.','t-stat');
    if np
        ip = nvx+nvn+ncx+ncn+1;
        disp('parameters')
        disp(tit1)
        for i=1:np
            name = bayestopt_.name{ip};
            disp(sprintf('%-*s %8.4f %7.4f %7.4f', ...
                         header_width,name,xparam1(ip),stdh(ip),tstath(ip)));
            eval(['oo_.mle_mode.parameters.' name ' = xparam1(ip);']);
633
            eval(['oo_.mle_std.parameters.' name ' = stdh(ip);']);
634
635
            ip = ip+1;
        end
636
    end
637
638
639
640
641
642
643
644
645
646
    if nvx
        ip = 1;
        disp('standard deviation of shocks')
        disp(tit1)
        for i=1:nvx
            k = estim_params_.var_exo(i,1);
            name = deblank(M_.exo_names(k,:));
            disp(sprintf('%-*s %8.4f %7.4f %7.4f',header_width,name,xparam1(ip),stdh(ip),tstath(ip)));
            M_.Sigma_e(k,k) = xparam1(ip)*xparam1(ip);
            eval(['oo_.mle_mode.shocks_std.' name ' = xparam1(ip);']);
647
            eval(['oo_.mle_std.shocks_std.' name ' = stdh(ip);']);
648
649
650
651
652
653
654
655
656
657
658
            ip = ip+1;
        end
    end
    if nvn
        disp('standard deviation of measurement errors')
        disp(tit1)
        ip = nvx+1;
        for i=1:nvn
            name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
            disp(sprintf('%-*s %8.4f %7.4f %7.4f',header_width,name,xparam1(ip),stdh(ip),tstath(ip)))
            eval(['oo_.mle_mode.measurement_errors_std.' name ' = xparam1(ip);']);
659
            eval(['oo_.mle_std.measurement_errors_std.' name ' = stdh(ip);']);
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
            ip = ip+1;
        end
    end
    if ncx
        disp('correlation of shocks')
        disp(tit1)
        ip = nvx+nvn+1;
        for i=1:ncx
            k1 = estim_params_.corrx(i,1);
            k2 = estim_params_.corrx(i,2);
            name = [deblank(M_.exo_names(k1,:)) ',' deblank(M_.exo_names(k2,:))];
            NAME = [deblank(M_.exo_names(k1,:)) '_' deblank(M_.exo_names(k2,:))];
            disp(sprintf('%-*s %8.4f %7.4f %7.4f', header_width,name,xparam1(ip),stdh(ip),tstath(ip)));
            M_.Sigma_e(k1,k2) = xparam1(ip)*sqrt(M_.Sigma_e(k1,k1)*M_.Sigma_e(k2,k2));
            M_.Sigma_e(k2,k1) = M_.Sigma_e(k1,k2);
            eval(['oo_.mle_mode.shocks_corr.' NAME ' = xparam1(ip);']);
676
            eval(['oo_.mle_std.shocks_corr.' NAME ' = stdh(ip);']);
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
            ip = ip+1;
        end
    end
    if ncn
        disp('correlation of measurement errors')
        disp(tit1)
        ip = nvx+nvn+ncx+1;
        for i=1:ncn
            k1 = estim_params_.corrn(i,1);
            k2 = estim_params_.corrn(i,2);
            name = [deblank(M_.endo_names(k1,:)) ',' deblank(M_.endo_names(k2,:))];
            NAME = [deblank(M_.endo_names(k1,:)) '_' deblank(M_.endo_names(k2,:))];
            disp(sprintf('%-*s %8.4f %7.4f %7.4f',header_width,name,xparam1(ip),stdh(ip),tstath(ip)));
            eval(['oo_.mle_mode.measurement_error_corr.' NAME ' = xparam1(ip);']);
            eval(['oo_.mle_std.measurement_error_corr.' NAME ' = stdh(ip);']);
            ip = ip+1;
        end
    end
end


698
OutputDirectoryName = CheckPath('Output',M_.dname);
699

700
if any(bayestopt_.pshape > 0) && options_.TeX %% Bayesian estimation (posterior mode) Latex output
701
702
703
704
705
706
707
708
    if np
        filename = [OutputDirectoryName '/' M_.fname '_Posterior_Mode_1.TeX'];
        fidTeX = fopen(filename,'w');
        fprintf(fidTeX,'%% TeX-table generated by dynare_estimation (Dynare).\n');
        fprintf(fidTeX,'%% RESULTS FROM POSTERIOR MAXIMIZATION (parameters)\n');
        fprintf(fidTeX,['%% ' datestr(now,0)]);
        fprintf(fidTeX,' \n');
        fprintf(fidTeX,' \n');
709
710
        fprintf(fidTeX,'\\begin{center}\n');
        fprintf(fidTeX,'\\begin{longtable}{l|lcccc} \n');
711
        fprintf(fidTeX,'\\caption{Results from posterior maximization (parameters)}\n ');
712
        fprintf(fidTeX,'\\label{Table:Posterior:1}\\\\\n');
713
714
        fprintf(fidTeX,'\\hline\\hline \\\\ \n');
        fprintf(fidTeX,'  & Prior distribution & Prior mean  & Prior s.d. & Posterior mode & s.d. \\\\ \n');
715
716
717
718
719
720
721
722
        fprintf(fidTeX,'\\hline \\endfirsthead \n');
        fprintf(fidTeX,'\\caption{(continued)}\n ');
        fprintf(fidTeX,'\\label{Table:Posterior:1}\\\\\n');
        fprintf(fidTeX,'\\hline\\hline \\\\ \n');
        fprintf(fidTeX,'  & Prior distribution & Prior mean  & Prior s.d. & Posterior mode & s.d. \\\\ \n');
        fprintf(fidTeX,'\\hline \\endhead \n');
        fprintf(fidTeX,'\\hline \\multicolumn{6}{r}{(Continued on next page)} \\\\ \\hline \\endfoot \n');
        fprintf(fidTeX,'\\hline \\hline \\endlastfoot \n');
723
724
725
726
727
728
729
730
731
        ip = nvx+nvn+ncx+ncn+1;
        for i=1:np
            fprintf(fidTeX,'$%s$ & %s & %7.3f & %6.4f & %8.4f & %7.4f \\\\ \n',...
                    M_.param_names_tex(estim_params_.param_vals(i,1),:),...
                    deblank(pnames(bayestopt_.pshape(ip)+1,:)),...
                    bayestopt_.p1(ip),...
                    bayestopt_.p2(ip),...
                    xparam1(ip),...
                    stdh(ip));
732
            ip = ip + 1;
733
        end
734
735
        fprintf(fidTeX,'\\end{longtable}\n ');    
        fprintf(fidTeX,'\\end{center}\n');
736
737
738
739
740
741
742
743
744
745
746
        fprintf(fidTeX,'%% End of TeX file.\n');
        fclose(fidTeX);
    end
    if nvx
        TeXfile = [OutputDirectoryName '/' M_.fname '_Posterior_Mode_2.TeX'];
        fidTeX = fopen(TeXfile,'w');
        fprintf(fidTeX,'%% TeX-table generated by dynare_estimation (Dynare).\n');
        fprintf(fidTeX,'%% RESULTS FROM POSTERIOR MAXIMIZATION (standard deviation of structural shocks)\n');
        fprintf(fidTeX,['%% ' datestr(now,0)]);
        fprintf(fidTeX,' \n');
        fprintf(fidTeX,' \n');
747
748
        fprintf(fidTeX,'\\begin{center}\n');
        fprintf(fidTeX,'\\begin{longtable}{l|lcccc} \n');
749
        fprintf(fidTeX,'\\caption{Results from posterior maximization (standard deviation of structural shocks)}\n ');
750
        fprintf(fidTeX,'\\label{Table:Posterior:2}\\\\\n');
751
        fprintf(fidTeX,'\\hline\\hline \\\\ \n');
Marco Ratto's avatar
Marco Ratto committed
752
        fprintf(fidTeX,'  & Prior distribution & Prior mean  & Prior s.d. & Posterior mode & s.d. \\\\ \n');
753
754
755
756
757
758
759
760
        fprintf(fidTeX,'\\hline \\endfirsthead \n');
        fprintf(fidTeX,'\\caption{(continued)}\n ');
        fprintf(fidTeX,'\\label{Table:Posterior:2}\\\\\n');
        fprintf(fidTeX,'\\hline\\hline \\\\ \n');
        fprintf(fidTeX,'  & Prior distribution & Prior mean  & Prior s.d. & Posterior mode & s.d. \\\\ \n');
        fprintf(fidTeX,'\\hline \\endhead \n');
        fprintf(fidTeX,'\\hline \\multicolumn{6}{r}{(Continued on next page)} \\\\ \\hline \\endfoot \n');
        fprintf(fidTeX,'\\hline \\hline \\endlastfoot \n');
761
762
763
764
765
766
767
768
769
        ip = 1;
        for i=1:nvx
            k = estim_params_.var_exo(i,1);
            fprintf(fidTeX,[ '$%s$ & %4s & %7.3f & %6.4f & %8.4f & %7.4f \\\\ \n'],...
                    deblank(M_.exo_names_tex(k,:)),...
                    deblank(pnames(bayestopt_.pshape(ip)+1,:)),...
                    bayestopt_.p1(ip),...
                    bayestopt_.p2(ip),...
                    xparam1(ip), ...
770
                    stdh(ip));
771
772
            ip = ip+1;
        end
773
774
        fprintf(fidTeX,'\\end{longtable}\n ');    
        fprintf(fidTeX,'\\end{center}\n');
775
776
777
778
779
780
781
782
783
784
785
        fprintf(fidTeX,'%% End of TeX file.\n');
        fclose(fidTeX);
    end
    if nvn
        TeXfile = [OutputDirectoryName '/' M_.fname '_Posterior_Mode_3.TeX'];
        fidTeX  = fopen(TeXfile,'w');
        fprintf(fidTeX,'%% TeX-table generated by dynare_estimation (Dynare).\n');
        fprintf(fidTeX,'%% RESULTS FROM POSTERIOR MAXIMIZATION (standard deviation of measurement errors)\n');
        fprintf(fidTeX,['%% ' datestr(now,0)]);
        fprintf(fidTeX,' \n');
        fprintf(fidTeX,' \n');
786
787
        fprintf(fidTeX,'\\begin{center}\n');
        fprintf(fidTeX,'\\begin{longtable}{l|lcccc} \n');
788
        fprintf(fidTeX,'\\caption{Results from posterior maximization (standard deviation of measurement errors)}\n ');
789
        fprintf(fidTeX,'\\label{Table:Posterior:3}\\\\\n');
790
        fprintf(fidTeX,'\\hline\\hline \\\\ \n');
791
792
793
794
795
796
797
798
799
        fprintf(fidTeX,'  & Prior distribution & Prior mean  & Prior s.d. &  Posterior mode & s.d. \\\\ \n');
        fprintf(fidTeX,'\\hline \\endfirsthead \n');
        fprintf(fidTeX,'\\caption{(continued)}\n ');
        fprintf(fidTeX,'\\label{Table:Posterior:3}\\\\\n');
        fprintf(fidTeX,'\\hline\\hline \\\\ \n');
        fprintf(fidTeX,'  & Prior distribution & Prior mean  & Prior s.d. &  Posterior mode & s.d. \\\\ \n');
        fprintf(fidTeX,'\\hline \\endhead \n');
        fprintf(fidTeX,'\\hline \\multicolumn{6}{r}{(Continued on next page)} \\\\ \\hline \\endfoot \n');
        fprintf(fidTeX,'\\hline \\hline \\endlastfoot \n');
800
801
802
803
804
        ip = nvx+1;
        for i=1:nvn
            idx = strmatch(options_.varobs(estim_params_.var_endo(i,1),:),M_.endo_names);
            fprintf(fidTeX,'$%s$ & %4s & %7.3f & %6.4f & %8.4f & %7.4f \\\\ \n',...
                    deblank(M_.endo_names_tex(idx,:)), ...
805
                    deblank(pnames(bayestopt_.pshape(ip)+1,:)), ...
806
                    bayestopt_.p1(ip), ...
807
                    bayestopt_.p2(ip),...
808
                    xparam1(ip),...
809
                    stdh(ip));
810
811
            ip = ip+1;
        end
812
813
        fprintf(fidTeX,'\\end{longtable}\n ');    
        fprintf(fidTeX,'\\end{center}\n');
814
815
816
817
818
819
820
821
822
823
824
        fprintf(fidTeX,'%% End of TeX file.\n');
        fclose(fidTeX);
    end
    if ncx
        TeXfile = [OutputDirectoryName '/' M_.fname '_Posterior_Mode_4.TeX'];
        fidTeX = fopen(TeXfile,'w');
        fprintf(fidTeX,'%% TeX-table generated by dynare_estimation (Dynare).\n');
        fprintf(fidTeX,'%% RESULTS FROM POSTERIOR MAXIMIZATION (correlation of structural shocks)\n');
        fprintf(fidTeX,['%% ' datestr(now,0)]);
        fprintf(fidTeX,' \n');
        fprintf(fidTeX,' \n');
825
826
        fprintf(fidTeX,'\\begin{center}\n');
        fprintf(fidTeX,'\\begin{longtable}{l|lcccc} \n');
827
        fprintf(fidTeX,'\\caption{Results from posterior parameters (correlation of structural shocks)}\n ');
828
829
830
831
832
833
        fprintf(fidTeX,'\\label{Table:Posterior:4}\\\\\n');
        fprintf(fidTeX,'\\hline\\hline \\\\ \n');
        fprintf(fidTeX,'  & Prior distribution & Prior mean  & Prior s.d. &  Posterior mode & s.d. \\\\ \n');
        fprintf(fidTeX,'\\hline \\endfirsthead \n');
        fprintf(fidTeX,'\\caption{(continued)}\n ');
        fprintf(fidTeX,'\\label{Table:Posterior:4}\\\\\n');
834
        fprintf(fidTeX,'\\hline\\hline \\\\ \n');
835
836
837
838
        fprintf(fidTeX,'  & Prior distribution & Prior mean  & Prior s.d. &  Posterior mode & s.d. \\\\ \n');
        fprintf(fidTeX,'\\hline \\endhead \n');
        fprintf(fidTeX,'\\hline \\multicolumn{6}{r}{(Continued on next page)} \\\\ \\hline \\endfoot \n');
        fprintf(fidTeX,'\\hline \\hline \\endlastfoot \n');
839
840
841
842
843
844
845
846
847
848
849
850
851
        ip = nvx+nvn+1;
        for i=1:ncx
            k1 = estim_params_.corrx(i,1);
            k2 = estim_params_.corrx(i,2);
            fprintf(fidTeX,[ '$%s$ & %s & %7.3f & %6.4f & %8.4f & %7.4f \\\\ \n'],...
                    [deblank(M_.exo_names_tex(k1,:)) ',' deblank(M_.exo_names_tex(k2,:))], ...
                    deblank(pnames(bayestopt_.pshape(ip)+1,:)), ...
                    bayestopt_.p1(ip), ...
                    bayestopt_.p2(ip), ...
                    xparam1(ip), ...
                    stdh(ip));
            ip = ip+1;
        end
852
853
        fprintf(fidTeX,'\\end{longtable}\n ');    
        fprintf(fidTeX,'\\end{center}\n');
854
855
856
857
858
859
860
861
862
863
864
        fprintf(fidTeX,'%% End of TeX file.\n');
        fclose(fidTeX);
    end
    if ncn
        TeXfile = [OutputDirectoryName '/' M_.fname '_Posterior_Mode_5.TeX'];
        fidTeX = fopen(TeXfile,'w');
        fprintf(fidTeX,'%% TeX-table generated by dynare_estimation (Dynare).\n');
        fprintf(fidTeX,'%% RESULTS FROM POSTERIOR MAXIMIZATION (correlation of measurement errors)\n');
        fprintf(fidTeX,['%% ' datestr(now,0)]);
        fprintf(fidTeX,' \n');
        fprintf(fidTeX,' \n');
865
866
        fprintf(fidTeX,'\\begin{center}\n');
        fprintf(fidTeX,'\\begin{longtabe}{l|lcccc} \n');
867
        fprintf(fidTeX,'\\caption{Results from posterior parameters (correlation of measurement errors)}\n ');
868
869
870
871
872
873
        fprintf(fidTeX,'\\label{Table:Posterior:5}\\\\\n');
        fprintf(fidTeX,'\\hline\\hline \\\\ \n');
        fprintf(fidTeX,'  & Prior distribution & Prior mean  & Prior s.d. &  Posterior mode & s.d. \\\\ \n');
        fprintf(fidTeX,'\\hline \\endfirsthead \n');
        fprintf(fidTeX,'\\caption{(continued)}\n ');
        fprintf(fidTeX,'\\label{Table:Posterior:5}\\\\\n');
874
        fprintf(fidTeX,'\\hline\\hline \\\\ \n');
875
876
877
878
        fprintf(fidTeX,'  & Prior distribution & Prior mean  & Prior s.d. &  Posterior mode & s.d. \\\\ \n');
        fprintf(fidTeX,'\\hline \\endhead \n');
        fprintf(fidTeX,'\\hline \\multicolumn{6}{r}{(Continued on next page)} \\\\ \\hline \\endfoot \n');
        fprintf(fidTeX,'\\hline \\hline \\endlastfoot \n');
879
880
881
882
883
884
885
886
887
888
889
890
891
        ip = nvx+nvn+ncx+1;
        for i=1:ncn
            k1 = estim_params_.corrn(i,1);
            k2 = estim_params_.corrn(i,2);
            fprintf(fidTeX,'$%s$ & %s & %7.3f & %6.4f & %8.4f & %7.4f \\\\ \n',...
                    [deblank(M_.endo_names_tex(k1,:)) ',' deblank(M_.endo_names_tex(k2,:))], ...
                    pnames(bayestopt_.pshape(ip)+1,:), ...
                    bayestopt_.p1(ip), ...
                    bayestopt_.p2(ip), ...
                    xparam1(ip), ...
                    stdh(ip));
            ip = ip+1;
        end
892
893
        fprintf(fidTeX,'\\end{longtable}\n ');    
        fprintf(fidTeX,'\\end{center}\n');
894
895
896
897
898
899
900
901
902
903
        fprintf(fidTeX,'%% End of TeX file.\n');
        fclose(fidTeX);
    end
end

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

904
905
if (any(bayestopt_.pshape  >0 ) && options_.mh_replic) || ...
        (any(bayestopt_.pshape >0 ) && options_.load_mh_file)  %% not ML estimation
906
    bounds = prior_bounds(bayestopt_,options_);
907
908
909
910
    bounds(:,1)=max(bounds(:,1),lb);
    bounds(:,2)=min(bounds(:,2),ub);
    bayestopt_.lb = bounds(:,1);
    bayestopt_.ub = bounds(:,2);
911
    if any(xparam1 < bounds(:,1)) || any(xparam1 > bounds(:,2))
912
913
914
915
916
917
        find(xparam1 < bounds(:,1))
        find(xparam1 > bounds(:,2))
        error('Mode values are outside prior bounds. Reduce prior_trunc.')
    end
    % runs MCMC
    if options_.mh_replic
918
        if options_.load_mh_file && options_.use_mh_covariance_matrix
919
920
            invhess = compute_mh_covariance_matrix;
        end
921
922
        ana_deriv = options_.analytic_derivation;
        options_.analytic_derivation = 0;
923
        if options_.cova_compute
924
            oo_.MC_record=feval(options_.posterior_sampling_method,objective_function,options_.proposal_distribution,xparam1,invhess,bounds,dataset_,options_,M_,estim_params_,bayestopt_,oo_);
925
        else
926
            error('I Cannot start the MCMC because the hessian of the posterior kernel at the mode was not computed.')
927
        end
928
        options_.analytic_derivation = ana_deriv;
929
    end
930
931
932
    if options_.mh_posterior_mode_estimation
        CutSample(M_, options_, estim_params_);
        return
933
    else
934
        if ~options_.nodiagnostic && options_.mh_replic > 1000 && options_.mh_nblck > 1
935
            McMCDiagnostics(options_, estim_params_, M_);
936
        end
937
938
939
940
941
942
        %% Here i discard first half of the draws:
        CutSample(M_, options_, estim_params_);
        %% Estimation of the marginal density from the Mh draws:
        if options_.mh_replic
            [marginal,oo_] = marginal_density(M_, options_, estim_params_, oo_);
            oo_ = GetPosteriorParametersStatistics(estim_params_, M_, options_, bayestopt_, oo_);
943
944
945
            if ~options_.nograph
                oo_ = PlotPosteriorDistributions(estim_params_, M_, options_, bayestopt_, oo_);
            end
946
            [oo_.posterior.metropolis.mean,oo_.posterior.metropolis.Variance] ...
947
                = GetPosteriorMeanVariance(M_,options_.mh_drop);
948
949
        else
            load([M_.fname '_results'],'oo_');
950
        end
951
952
953
        metropolis_draw(1);
        if options_.bayesian_irf
            PosteriorIRF('posterior');
954
        end
955
956
        if options_.moments_varendo
            oo_ = compute_moments_varendo('posterior',options_,M_,oo_,var_list_);
957
        end
958
        if options_.smoother || ~isempty(options_.filter_step_ahead) || options_.forecast
959
            prior_posterior_statistics('posterior',dataset_);
960
961
        end
        xparam = get_posterior_parameters('mean');
962
        M_ = set_all_parameters(xparam,estim_params_,M_);
963
964
    end
end
965

966
967
968
969
if options_.particle.status
    return
end

970
971
if (~((any(bayestopt_.pshape > 0) && options_.mh_replic) || (any(bayestopt_.pshape ...
                                                      > 0) && options_.load_mh_file)) ...
MichelJuillard's avatar
MichelJuillard committed
972
    || ~options_.smoother ) && options_.partial_information == 0  % to be fixed
973
    %% ML estimation, or posterior mode without metropolis-hastings or metropolis without bayesian smooth variable
974
    [atT,innov,measurement_error,updated_variables,ys,trend_coeff,aK,T,R,P,PK,decomp] = DsgeSmoother(xparam1,dataset_.info.ntobs,dataset_.data,dataset_.missing.aindex,dataset_.missing.state);
975
976
    oo_.Smoother.SteadyState = ys;
    oo_.Smoother.TrendCoeffs = trend_coeff;
977
    oo_.Smoother.Variance = P;
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
    i_endo = bayestopt_.smoother_saved_var_list;
    if options_.nk ~= 0
        oo_.FilteredVariablesKStepAhead = aK(options_.filter_step_ahead, ...
                                             i_endo,:);
        if isfield(options_,'kalman_algo')
            if ~isempty(PK)
                oo_.FilteredVariablesKStepAheadVariances = ...
                    PK(options_.filter_step_ahead,i_endo,i_endo,:);
            end
            if ~isempty(decomp)
                oo_.FilteredVariablesShockDecomposition = ...
                    decomp(options_.filter_step_ahead,i_endo,:,:);
            end
        end
    end
    for i=bayestopt_.smoother_saved_var_list'
994
        i1 = dr.order_var(bayestopt_.smoother_var_list(i));
995
996
997
998
999
1000
        eval(['oo_.SmoothedVariables.' deblank(M_.endo_names(i1,:)) ' = ' ...
                            'atT(i,:)'';']);
        if options_.nk > 0
            eval(['oo_.FilteredVariables.' deblank(M_.endo_names(i1,:)) ...
                  ' = squeeze(aK(1,i,:));']);
        end
1001
1002
1003
        eval(['oo_.UpdatedVariables.' deblank(M_.endo_names(i1,:)) ...
              ' = updated_variables(i,:)'';']);
    end
1004
1005
1006
    for i=1:M_.exo_nbr
        eval(['oo_.SmoothedShocks.' deblank(M_.exo_names(i,:)) ' = innov(i,:)'';']);
    end
1007
1008
    if ~options_.nograph,
        [nbplt,nr,nc,lr,lc,nstar] = pltorg(M_.exo_nbr);
1009
        if options_.TeX
1010
1011
1012
1013
            fidTeX = fopen([M_.fname '_SmoothedShocks.TeX'],'w');
            fprintf(fidTeX,'%% TeX eps-loader file generated by dynare_estimation.m (Dynare).\n');
            fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
            fprintf(fidTeX,' \n');
1014
        end
1015
1016
        for plt = 1:nbplt,
            hh = dyn_figure(options_,'Name','Smoothed shocks');
1017
1018
            NAMES = [];
            if options_.TeX, TeXNAMES = []; end
1019
1020
            nstar0=min(nstar,M_.exo_nbr-(plt-1)*nstar);
            for i=1:nstar0,
1021
1022
1023
1024
1025
1026
1027
                k = (plt-1)*nstar+i;
                subplot(nr,nc,i);
                plot([1 gend],[0 0],'-r','linewidth',.5)
                hold on
                plot(1:gend,innov(k,:),'-k','linewidth',1)
                hold off
                name = deblank(M_.exo_names(k,:));
1028
1029
1030
1031
1032
                if isempty(NAMES)
                    NAMES = name;
                else
                    NAMES = char(NAMES,name);
                end
1033
1034