2009-12-16 18:17:34 +01:00
function dynare_estimation_1 ( var_list_,dname)
% function dynare_estimation_1(var_list_,dname)
% runs the estimation of the model
2011-09-17 15:38:49 +02:00
%
2009-12-16 18:17:34 +01:00
% INPUTS
% var_list_: selected endogenous variables vector
% dname: alternative directory name
2011-09-17 15:38:49 +02:00
%
2009-12-16 18:17:34 +01:00
% OUTPUTS
% none
%
% SPECIAL REQUIREMENTS
% none
2013-01-07 14:51:45 +01:00
% Copyright (C) 2003-2013 Dynare Team
2009-12-16 18:17:34 +01:00
%
% 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/>.
2014-06-16 17:41:59 +02:00
global M_ options_ oo_ estim_params_ bayestopt_ dataset_ dataset_info
2009-12-16 18:17:34 +01:00
2012-03-03 16:27:25 +01:00
% Set particle filter flag.
if options_ . order > 1
2012-03-05 08:46:52 +01:00
if options_ . particle . status && options_ . order == 2
2013-07-10 16:16:32 +02:00
skipline ( )
2012-03-05 08:46:52 +01:00
disp ( ' Estimation using a non linear filter!' )
2013-07-10 16:16:32 +02:00
skipline ( )
2013-07-10 16:03:49 +02:00
if ~ options_ . nointeractive && ismember ( options_ . mode_compute , [ 1 , 3 , 4 ] ) % Known gradient-based optimizers
2013-07-10 16:15:28 +02:00
disp ( ' You are using a gradient-based mode-finder. Particle filtering introduces discontinuities in the' )
disp ( ' objective function w.r.t the parameters. Thus, should use a non-gradient based optimizer.' )
2013-07-08 19:34:31 +02:00
fprintf ( ' \nPlease choose a mode-finder:\n' )
2013-07-10 16:03:49 +02:00
fprintf ( ' \t 0 - Continue using gradient-based method (it is most likely that you will no get any sensible result).\n' )
2013-07-08 19:34:31 +02:00
fprintf ( ' \t 6 - Monte Carlo based algorithm\n' )
2013-07-10 16:03:49 +02:00
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' )
2013-07-08 19:34:31 +02:00
fprintf ( ' \t 9 - CMA-ES (Covariance Matrix Adaptation Evolution Strategy) algorithm\n' )
2013-07-10 16:03:49 +02:00
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
2013-07-08 19:34:31 +02:00
end
end
2012-03-05 14:56:02 +01:00
elseif options_ . particle . status && options_ . order > 2
2012-03-05 08:46:52 +01:00
error ( [ ' Non linear filter are not implemented with order ' int2str ( options_ . order ) ' approximation of the model!' ] )
elseif ~ options_ . particle . status && options_ . order == 2
2012-06-06 17:06:28 +02:00
error ( ' For estimating the model with a second order approximation using a non linear filter, one should have options_.particle.status=1;' )
2012-03-05 08:46:52 +01:00
else
error ( [ ' Cannot estimate a model with an order ' int2str ( options_ . order ) ' approximation!' ] )
end
2012-03-03 16:27:25 +01:00
end
2011-09-19 16:38:38 +02:00
if ~ options_ . dsge_var
2012-03-05 08:46:52 +01:00
if options_ . particle . status
2012-03-03 16:27:25 +01:00
objective_function = str2func ( ' non_linear_dsge_likelihood' ) ;
2014-09-05 18:09:49 +02:00
if options_ . particle . filter_algorithm . sis == 1
options_ . particle . algorithm = ' sequential_importance_particle_filter' ;
else
if options_ . particle . filter_algorithm . apf == 1
options_ . particle . algorithm = ' auxiliary_particle_filter' ;
else
if options_ . particle . filter_algorithm . gf == 1
options_ . particle . algorithm = ' gaussian_filter' ;
else
if options_ . particle . filter_algorithm . gmf == 1
options_ . particle . algorithm = ' gaussian_mixture_filter' ;
else
if options_ . particle . filter_algorithm . cpf == 1
options_ . particle . algorithm = ' conditional_particle_filter' ;
else
error ( ' Estimation: Unknown filter!' )
end
end
end
end
end
2012-03-03 16:27:25 +01:00
else
objective_function = str2func ( ' dsge_likelihood' ) ;
end
2011-09-19 16:38:38 +02:00
else
2013-11-16 23:33:37 +01:00
objective_function = str2func ( ' dsge_var_likelihood' ) ;
2011-09-19 16:38:38 +02:00
end
2014-07-16 17:00:17 +02:00
[ dataset_ , dataset_info , xparam1 , hh , M_ , options_ , oo_ , estim_params_ , bayestopt_ ] = ...
2014-06-16 17:41:59 +02:00
dynare_estimation_init ( var_list_ , dname , [ ] , M_ , options_ , oo_ , estim_params_ , bayestopt_ ) ;
2011-09-17 15:38:49 +02:00
2013-11-16 23:26:19 +01:00
if options_ . dsge_var
check_dsge_var_model ( M_ , estim_params_ , bayestopt_ ) ;
2014-06-23 10:55:08 +02:00
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
2012-06-06 15:58:48 +02:00
end
2013-10-30 14:26:30 +01:00
%check for calibrated covariances before updating parameters
if ~ isempty ( estim_params_ )
estim_params_ = check_for_calibrated_covariances ( xparam1 , estim_params_ , M_ ) ;
2012-06-06 18:24:03 +02:00
end
2013-10-30 21:43:33 +01:00
%%read out calibration that was set in mod-file and can be used for initialization
xparam1_calib = get_all_parameters ( estim_params_ , M_ ) ; %get calibrated parameters
if ~ any ( isnan ( xparam1_calib ) ) %all estimated parameters are calibrated
full_calibration_detected = 1 ;
else
full_calibration_detected = 0 ;
end
if options_ . use_calibration_initialization %set calibration as starting values
[ xparam1 , estim_params_ ] = do_parameter_initialization ( estim_params_ , xparam1_calib , xparam1 ) ; %get explicitly initialized parameters that have precedence to calibrated values
try
check_prior_bounds ( xparam1 , [ bayestopt_ . lb bayestopt_ . ub ] , M_ , estim_params_ , options_ , bayestopt_ ) ; %check whether calibration satisfies prior bounds
2013-11-08 16:35:52 +01:00
catch
e = lasterror ( ) ;
2013-10-30 21:43:33 +01:00
fprintf ( ' Cannot use parameter values from calibration as they violate the prior bounds.' )
2013-11-08 16:35:52 +01:00
rethrow ( e ) ;
2013-03-19 21:12:49 +01:00
end
end
2012-06-06 18:24:03 +02:00
% Set sigma_e_is_diagonal flag (needed if the shocks block is not declared in the mod file).
2012-06-06 15:58:48 +02:00
M_ . sigma_e_is_diagonal = 1 ;
2013-10-30 14:26:30 +01:00
if estim_params_ . ncx || any ( nnz ( tril ( M_ . Correlation_matrix , - 1 ) ) ) || isfield ( estim_params_ , ' calibrated_covariances' )
2012-06-06 15:58:48 +02:00
M_ . sigma_e_is_diagonal = 0 ;
end
2011-09-17 15:38:49 +02:00
data = dataset_ . data ;
2014-07-09 15:28:46 +02:00
rawdata = dataset_info . rawdata ;
2014-06-16 17:41:59 +02:00
data_index = dataset_info . missing . aindex ;
missing_value = dataset_info . missing . state ;
2009-12-16 18:17:34 +01:00
2011-09-17 15:38:49 +02:00
% Set number of observations
2014-06-16 17:41:59 +02:00
gend = dataset_ . nobs ;
2011-09-17 15:38:49 +02:00
% Set the number of observed variables.
2013-05-21 16:38:17 +02:00
n_varobs = length ( options_ . varobs ) ;
2011-09-17 15:38:49 +02:00
% Get the number of parameters to be estimated.
2009-12-16 18:17:34 +01:00
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.
2011-02-21 11:20:26 +01:00
%% Set the names of the priors.
pnames = [ ' ' ; ' beta ' ; ' gamm ' ; ' norm ' ; ' invg ' ; ' unif ' ; ' invg2' ] ;
%% Set parameters bounds
lb = bayestopt_ . lb ;
ub = bayestopt_ . ub ;
2009-12-16 18:17:34 +01:00
2011-02-21 11:20:26 +01:00
dr = oo_ . dr ;
2009-12-16 18:17:34 +01:00
2014-02-25 15:41:34 +01:00
% load optimal_mh_scale parameter if previous run was with mode_compute=6
2012-10-06 16:45:19 +02:00
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
2014-02-25 15:41:34 +01:00
delete ( mh_scale_fname )
2012-10-06 16:45:19 +02:00
end
end
2009-12-16 18:17:34 +01:00
if ~ isempty ( estim_params_ )
2013-10-30 14:26:30 +01:00
M_ = set_all_parameters ( xparam1 , estim_params_ , M_ ) ;
2009-12-16 18:17:34 +01:00
end
2011-10-12 21:46:50 +02:00
2009-12-16 18:17:34 +01:00
2013-10-30 21:43:33 +01:00
%% perform initial estimation checks;
try
2014-06-16 17:41:59 +02:00
oo_ = initial_estimation_checks ( objective_function , xparam1 , dataset_ , dataset_info , M_ , estim_params_ , options_ , bayestopt_ , oo_ ) ;
2013-11-08 15:02:30 +01:00
catch % if check fails, provide info on using calibration if present
e = lasterror ( ) ;
2013-10-30 21:43:33 +01:00
if full_calibration_detected %calibrated model present and no explicit starting values
skipline ( 1 ) ;
fprintf ( ' ESTIMATION_CHECKS: There was an error in computing the likelihood for initial parameter values.\n' )
fprintf ( ' ESTIMATION_CHECKS: You should try using the calibrated version of the model as starting values. To do\n' )
2013-11-08 16:09:16 +01:00
fprintf ( ' ESTIMATION_CHECKS: this, add an empty estimated_params_init-block with use_calibration option immediately before the estimation\n' )
2013-10-30 21:43:33 +01:00
fprintf ( ' ESTIMATION_CHECKS: command (and after the estimated_params-block so that it does not get overwritten):\n' ) ;
skipline ( 2 ) ;
end
2013-11-08 15:02:30 +01:00
rethrow ( e ) ;
2013-10-30 21:43:33 +01:00
end
2009-12-16 18:17:34 +01:00
2011-01-28 14:01:57 +01:00
if isequal ( options_ . mode_compute , 0 ) && isempty ( options_ . mode_file ) && options_ . mh_posterior_mode_estimation == 0
2009-12-16 18:17:34 +01:00
if options_ . smoother == 1
2010-04-03 08:00:22 +02:00
[ atT , innov , measurement_error , updated_variables , ys , trend_coeff , aK , T , R , P , PK , decomp ] = DsgeSmoother ( xparam1 , gend , data , data_index , missing_value ) ;
2010-01-05 11:46:10 +01:00
oo_ . Smoother . SteadyState = ys ;
oo_ . Smoother . TrendCoeffs = trend_coeff ;
2010-03-23 18:17:15 +01:00
if options_ . filter_covariance
2013-03-17 22:49:28 +01:00
oo_ . Smoother . Variance = P ;
2010-03-23 18:17:15 +01:00
end
2010-03-19 20:40:49 +01:00
i_endo = bayestopt_ . smoother_saved_var_list ;
2010-01-05 11:46:10 +01:00
if options_ . nk ~= 0
2010-03-23 18:17:15 +01:00
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
2010-01-05 11:46:10 +01:00
end
2011-01-18 17:50:10 +01:00
for i = bayestopt_ . smoother_saved_var_list '
2011-06-21 13:20:51 +02:00
i1 = dr . order_var ( bayestopt_ . smoother_var_list ( i ) ) ;
2011-11-02 14:02:12 +01:00
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
2010-03-23 18:17:15 +01:00
eval ( [ ' oo_.UpdatedVariables.' deblank ( M_ . endo_names ( i1 , : ) ) ' = updated_variables(i,:)' ' ;' ] ) ;
2010-01-05 11:46:10 +01:00
end
for i = 1 : M_ . exo_nbr
eval ( [ ' oo_.SmoothedShocks.' deblank ( M_ . exo_names ( i , : ) ) ' = innov(i,:)' ' ;' ] ) ;
end
2009-12-16 18:17:34 +01:00
end
2011-09-19 18:10:29 +02:00
return
2009-12-16 18:17:34 +01:00
end
2013-10-03 11:01:11 +02:00
% Estimation of the posterior mode or likelihood mode
2011-01-28 14:01:57 +01:00
if ~ isequal ( options_ . mode_compute , 0 ) && ~ options_ . mh_posterior_mode_estimation
2009-12-16 18:17:34 +01:00
switch options_ . mode_compute
case 1
2013-11-04 10:54:45 +01:00
if isoctave
2013-01-07 14:51:45 +01:00
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
2013-10-03 11:01:11 +02:00
% Set default optimization options for fmincon.
optim_options = optimset ( ' display' , ' iter' , ' LargeScale' , ' off' , ' MaxFunEvals' , 100000 , ' TolFun' , 1e-8 , ' TolX' , 1e-6 ) ;
2009-12-16 18:17:34 +01:00
if isfield ( options_ , ' optim_opt' )
eval ( [ ' optim_options = optimset(optim_options,' options_ . optim_opt ' );' ] ) ;
2012-04-27 15:57:58 +02:00
end
if options_ . analytic_derivation ,
2012-07-05 10:22:36 +02:00
optim_options = optimset ( optim_options , ' GradObj' , ' on' , ' TolX' , 1e-7 ) ;
2009-12-16 18:17:34 +01:00
end
2013-10-03 11:01:11 +02:00
[ xparam1 , fval , exitflag , output , lamdba , grad , hessian_fmincon ] = ...
2014-06-16 17:41:59 +02:00
fmincon ( objective_function , xparam1 , [ ] , [ ] , [ ] , [ ] , lb , ub , [ ] , optim_options , dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , oo_ ) ;
2009-12-16 18:17:34 +01:00
case 2
error ( ' ESTIMATION: mode_compute=2 option (Lester Ingber' ' s Adaptive Simulated Annealing) is no longer available' )
case 3
2013-11-04 10:54:45 +01:00
if isoctave && ~ user_has_octave_forge_package ( ' optim' )
2013-01-07 14:51:45 +01:00
error ( ' Option mode_compute=3 requires the optim package' )
2013-11-04 10:54:45 +01:00
elseif ~ isoctave && ~ user_has_matlab_license ( ' optimization_toolbox' )
2013-01-07 14:51:45 +01:00
error ( ' Option mode_compute=3 requires the Optimization Toolbox' )
end
2013-10-03 11:06:07 +02:00
% Set default optimization options for fminunc.
2009-12-16 18:17:34 +01:00
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
2012-04-27 15:57:58 +02:00
if options_ . analytic_derivation ,
optim_options = optimset ( optim_options , ' GradObj' , ' on' ) ;
end
2013-11-04 10:54:45 +01:00
if ~ isoctave
2014-06-16 17:41:59 +02:00
[ xparam1 , fval , exitflag ] = fminunc ( objective_function , xparam1 , optim_options , dataset_ , dataset_info_ , options_ , M_ , estim_params_ , bayestopt_ , oo_ ) ;
2013-01-07 15:38:45 +01:00
else
% Under Octave, use a wrapper, since fminunc() does not have a 4th arg
2014-06-16 17:41:59 +02:00
func = @ ( x ) objective_function ( x , dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , oo_ ) ;
2013-01-07 15:38:45 +01:00
[ xparam1 , fval , exitflag ] = fminunc ( func , xparam1 , optim_options ) ;
end
2009-12-16 18:17:34 +01:00
case 4
2013-10-03 12:35:06 +02:00
% Set default options.
2009-12-16 18:17:34 +01:00
H0 = 1e-4 * eye ( nx ) ;
crit = 1e-7 ;
nit = 1000 ;
2013-10-03 12:35:06 +02:00
numgrad = options_ . gradient_method ;
epsilon = options_ . gradient_epsilon ;
% Change some options.
if isfield ( options_ , ' optim_opt' )
2014-02-03 12:23:14 +01:00
options_list = read_key_value_string ( options_ . optim_opt ) ;
for i = 1 : rows ( options_list )
switch options_list { i , 1 }
case ' MaxIter'
nit = options_list { i , 2 } ;
case ' InitialInverseHessian'
H0 = eval ( options_list { i , 2 } ) ;
2014-02-21 12:23:36 +01:00
case ' TolFun'
2014-02-03 12:23:14 +01:00
crit = options_list { i , 2 } ;
case ' NumgradAlgorithm'
numgrad = options_list { i , 2 } ;
case ' NumgradEpsilon'
epsilon = options_list { i , 2 } ;
2013-10-03 12:35:06 +02:00
otherwise
2014-02-03 12:23:14 +01:00
warning ( [ ' csminwel: Unknown option (' options_list { i , 1 } ' )!' ] )
2013-10-03 12:35:06 +02:00
end
end
end
% Set flag for analytical gradient.
if options_ . analytic_derivation
2012-04-27 15:57:58 +02:00
analytic_grad = 1 ;
else
analytic_grad = [ ] ;
end
2013-10-03 12:35:06 +02:00
% Call csminwell.
[ fval , xparam1 , grad , hessian_csminwel , itct , fcount , retcodehat ] = ...
2014-06-16 17:41:59 +02:00
csminwel1 ( objective_function , xparam1 , H0 , analytic_grad , crit , nit , numgrad , epsilon , dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , oo_ ) ;
2013-10-03 12:35:06 +02:00
% Disp value at the mode.
disp ( sprintf ( ' Objective function at mode: %f' , fval ) )
2009-12-16 18:17:34 +01:00
case 5
if isfield ( options_ , ' hess' )
flag = options_ . hess ;
else
flag = 1 ;
end
if isfield ( options_ , ' ftol' )
crit = options_ . ftol ;
else
2011-02-08 15:42:38 +01:00
crit = 1.e-5 ;
2009-12-16 18:17:34 +01:00
end
2012-06-08 14:23:18 +02:00
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
2009-12-16 18:17:34 +01:00
if isfield ( options_ , ' nit' )
nit = options_ . nit ;
else
nit = 1000 ;
end
2014-07-23 18:20:16 +02:00
[ xparam1 , hh , gg , fval , invhess ] = newrat ( objective_function , xparam1 , analytic_grad , crit , nit , 0 , dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , oo_ ) ;
2012-06-08 14:23:18 +02:00
if options_ . analytic_derivation ,
options_ . analytic_derivation = ana_deriv ;
2014-07-23 18:20:16 +02:00
else
if flag == 0 ,
options_ . cova_compute = 1 ;
kalman_algo0 = options_ . kalman_algo ;
if ~ ( ( options_ . kalman_algo == 2 ) || ( options_ . kalman_algo == 4 ) ) ,
options_ . kalman_algo = 2 ;
if options_ . lik_init == 3 ,
options_ . kalman_algo = 4 ;
end
end
hh = reshape ( mr_hessian ( 0 , xparam1 , objective_function , 1 , crit , dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , oo_ ) , nx , nx ) ;
options_ . kalman_algo = kalman_algo0 ;
end
2012-06-08 14:23:18 +02:00
end
2009-12-16 18:17:34 +01:00
parameter_names = bayestopt_ . name ;
2014-07-23 18:20:16 +02:00
save ( [ M_ . fname ' _mode.mat' ] , ' xparam1' , ' hh' , ' parameter_names' ) ;
2009-12-16 18:17:34 +01:00
case 6
2013-10-04 16:12:14 +02:00
% Set default options
gmhmaxlikOptions = options_ . gmhmaxlik ;
if ~ isempty ( hh ) ;
gmhmaxlikOptions . varinit = ' previous' ;
else
gmhmaxlikOptions . varinit = ' prior' ;
end
if isfield ( options_ , ' optim_opt' )
2014-02-03 12:23:14 +01:00
options_list = read_key_value_string ( options_ . optim_opt ) ;
for i = 1 : rows ( options_list )
switch options_list { i , 1 }
case ' NumberOfMh'
gmhmaxlikOptions . iterations = options_list { i , 2 } ;
case ' ncov-mh'
gmhmaxlikOptions . number = options_list { i , 2 } ;
case ' nscale'
gmhmaxlikOptions . nscale = options_list { i , 2 } ;
case ' nclimb'
gmhmaxlikOptions . nclimb = options_list { i , 2 } ;
case ' InitialCovarianceMatrix'
switch options_list { i , 2 }
2013-10-04 16:12:14 +02:00
case ' previous'
if isempty ( hh )
error ( ' gmhmaxlik: No previous estimate of the Hessian matrix available!' )
else
gmhmaxlikOptions . varinit = ' previous'
end
case { ' prior' , ' identity' }
2014-02-03 12:23:14 +01:00
gmhmaxlikOptions . varinit = options_list { i , 2 } ;
2013-10-04 16:12:14 +02:00
otherwise
error ( ' gmhmaxlik: Unknown value for option ' ' InitialCovarianceMatrix' ' !' )
end
2014-02-03 12:23:14 +01:00
case ' AcceptanceRateTarget'
gmhmaxlikOptions . target = options_list { i , 2 } ;
2013-10-04 16:12:14 +02:00
if gmhmaxlikOptions . target > 1 || gmhmaxlikOptions . target < eps
error ( ' gmhmaxlik: The value of option AcceptanceRateTarget should be a double between 0 and 1!' )
end
otherwise
2014-02-21 19:19:36 +01:00
warning ( [ ' gmhmaxlik: Unknown option (' options_list { i , 1 } ' )!' ] )
2013-10-04 16:12:14 +02:00
end
end
2014-02-03 12:23:14 +01:00
end
2013-10-04 16:12:14 +02:00
% Evaluate the objective function.
2014-06-16 17:41:59 +02:00
fval = feval ( objective_function , xparam1 , dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , oo_ ) ;
2009-12-16 18:17:34 +01:00
OldMode = fval ;
2011-03-26 20:19:23 +01:00
if ~ exist ( ' MeanPar' , ' var' )
2009-12-16 18:17:34 +01:00
MeanPar = xparam1 ;
end
2013-10-04 16:12:14 +02:00
switch gmhmaxlikOptions . varinit
case ' previous'
2009-12-16 18:17:34 +01:00
CovJump = inv ( hh ) ;
2013-10-04 16:12:14 +02:00
case ' prior'
% The covariance matrix is initialized with the prior
2009-12-16 18:17:34 +01:00
% covariance (a diagonal matrix) %%Except for infinite variances ;-)
2013-10-04 16:12:14 +02:00
stdev = bayestopt_ . p2 ;
indx = find ( isinf ( stdev ) ) ;
stdev ( indx ) = ones ( length ( indx ) , 1 ) * sqrt ( 10 ) ;
vars = stdev .^ 2 ;
CovJump = diag ( vars ) ;
case ' identity'
vars = ones ( length ( bayestopt_ . p2 ) , 1 ) * 0.1 ;
CovJump = diag ( vars ) ;
otherwise
error ( ' gmhmaxlik: This is a bug! Please contact the developers.' )
2009-12-16 18:17:34 +01:00
end
OldPostVar = CovJump ;
Scale = options_ . mh_jscale ;
2013-10-04 16:12:14 +02:00
for i = 1 : gmhmaxlikOptions . iterations
2009-12-16 18:17:34 +01:00
if i == 1
2013-10-04 16:12:14 +02:00
if gmhmaxlikOptions . iterations > 1
2009-12-16 18:17:34 +01:00
flag = ' ' ;
else
flag = ' LastCall' ;
end
2011-09-19 18:10:29 +02:00
[ xparam1 , PostVar , Scale , PostMean ] = ...
2014-06-16 17:41:59 +02:00
gmhmaxlik ( objective_function , xparam1 , [ lb ub ] , gmhmaxlikOptions , Scale , flag , MeanPar , CovJump , dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , oo_ ) ;
fval = feval ( objective_function , xparam1 , dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , oo_ ) ;
2009-12-16 18:17:34 +01:00
options_ . mh_jscale = Scale ;
mouvement = max ( max ( abs ( PostVar - OldPostVar ) ) ) ;
2013-07-10 16:16:32 +02:00
skipline ( )
2010-06-25 14:58:53 +02:00
disp ( ' ========================================================== ' )
disp ( [ ' Change in the covariance matrix = ' num2str ( mouvement ) ' .' ] )
disp ( [ ' Mode improvement = ' num2str ( abs ( OldMode - fval ) ) ] )
disp ( [ ' New value of jscale = ' num2str ( Scale ) ] )
disp ( ' ========================================================== ' )
2009-12-16 18:17:34 +01:00
OldMode = fval ;
else
OldPostVar = PostVar ;
2013-10-04 16:12:14 +02:00
if i < gmhmaxlikOptions . iterations
2009-12-16 18:17:34 +01:00
flag = ' ' ;
else
flag = ' LastCall' ;
end
2011-09-19 18:10:29 +02:00
[ xparam1 , PostVar , Scale , PostMean ] = ...
gmhmaxlik ( objective_function , xparam1 , [ lb ub ] , ...
2014-06-16 17:41:59 +02:00
gmhmaxlikOptions , Scale , flag , PostMean , PostVar , dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , oo_ ) ;
fval = feval ( objective_function , xparam1 , dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , oo_ ) ;
2009-12-16 18:17:34 +01:00
options_ . mh_jscale = Scale ;
mouvement = max ( max ( abs ( PostVar - OldPostVar ) ) ) ;
2014-06-16 17:41:59 +02:00
fval = feval ( objective_function , xparam1 , dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , oo_ ) ;
2013-07-10 16:16:32 +02:00
skipline ( )
2012-09-29 00:06:44 +02:00
disp ( ' ========================================================== ' )
disp ( [ ' Change in the covariance matrix = ' num2str ( mouvement ) ' .' ] )
disp ( [ ' Mode improvement = ' num2str ( abs ( OldMode - fval ) ) ] )
disp ( [ ' New value of jscale = ' num2str ( Scale ) ] )
disp ( ' ========================================================== ' )
2009-12-16 18:17:34 +01:00
OldMode = fval ;
end
hh = inv ( PostVar ) ;
2013-10-04 16:12:14 +02:00
parameter_names = bayestopt_ . name ;
save ( [ M_ . fname ' _mode.mat' ] , ' xparam1' , ' hh' , ' parameter_names' ) ;
2010-10-18 14:39:48 +02:00
save ( [ M_ . fname ' _optimal_mh_scale_parameter.mat' ] , ' Scale' ) ;
2011-09-19 18:10:29 +02:00
bayestopt_ . jscale = ones ( length ( xparam1 ) , 1 ) * Scale ;
2011-09-17 15:38:49 +02:00
end
2013-07-10 16:16:32 +02:00
skipline ( )
2012-09-29 00:06:44 +02:00
disp ( [ ' Optimal value of the scale parameter = ' num2str ( Scale ) ] )
2013-07-10 16:16:32 +02:00
skipline ( )
2012-09-29 00:06:44 +02:00
disp ( [ ' Final value of the log posterior (or likelihood): ' num2str ( fval ) ] )
2013-07-10 16:16:32 +02:00
skipline ( )
2012-09-27 14:55:57 +02:00
parameter_names = bayestopt_ . name ;
save ( [ M_ . fname ' _mode.mat' ] , ' xparam1' , ' hh' , ' parameter_names' ) ;
2009-12-16 18:17:34 +01:00
case 7
2011-09-19 18:10:29 +02:00
% Matlab's simplex (Optimization toolbox needed).
2013-11-04 10:54:45 +01:00
if isoctave && ~ user_has_octave_forge_package ( ' optim' )
2013-01-07 14:51:45 +01:00
error ( ' Option mode_compute=7 requires the optim package' )
2013-11-04 10:54:45 +01:00
elseif ~ isoctave && ~ user_has_matlab_license ( ' optimization_toolbox' )
2013-01-07 14:51:45 +01:00
error ( ' Option mode_compute=7 requires the Optimization Toolbox' )
end
2009-12-16 18:17:34 +01:00
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
2014-06-16 17:41:59 +02:00
[ xparam1 , fval , exitflag ] = fminsearch ( objective_function , xparam1 , optim_options , dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , oo_ ) ;
2011-05-10 13:51:11 +02:00
case 8
2011-09-19 18:10:29 +02:00
% Dynare implementation of the simplex algorithm.
2013-10-08 12:55:11 +02:00
simplexOptions = options_ . simplex ;
if isfield ( options_ , ' optim_opt' )
2014-02-03 12:23:14 +01:00
options_list = read_key_value_string ( options_ . optim_opt ) ;
for i = 1 : rows ( options_list )
switch options_list { i , 1 }
case ' MaxIter'
simplexOptions . maxiter = options_list { i , 2 } ;
case ' TolFun'
simplexOptions . tolerance . f = options_list { i , 2 } ;
case ' TolX'
simplexOptions . tolerance . x = options_list { i , 2 } ;
case ' MaxFunEvals'
simplexOptions . maxfcall = options_list { i , 2 } ;
case ' MaxFunEvalFactor'
simplexOptions . maxfcallfactor = options_list { i , 2 } ;
case ' InitialSimplexSize'
simplexOptions . delta_factor = options_list { i , 2 } ;
2013-10-08 12:55:11 +02:00
otherwise
2014-02-03 12:23:14 +01:00
warning ( [ ' simplex: Unknown option (' options_list { i , 1 } ' )!' ] )
2013-10-08 12:55:11 +02:00
end
end
end
2014-06-16 17:41:59 +02:00
[ xparam1 , fval , exitflag ] = simplex_optimization_routine ( objective_function , xparam1 , simplexOptions , dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , oo_ ) ;
2012-03-03 16:26:47 +01:00
case 9
2013-10-08 15:53:55 +02:00
% Set defaults
2012-03-03 15:30:55 +01:00
H0 = 1e-4 * ones ( nx , 1 ) ;
2013-10-08 15:53:55 +02:00
cmaesOptions = options_ . cmaes ;
% Modify defaults
if isfield ( options_ , ' optim_opt' )
2014-02-03 12:23:14 +01:00
options_list = read_key_value_string ( options_ . optim_opt ) ;
for i = 1 : rows ( options_list )
switch options_list { i , 1 }
case ' MaxIter'
cmaesOptions . MaxIter = options_list { i , 2 } ;
case ' TolFun'
cmaesOptions . TolFun = options_list { i , 2 } ;
case ' TolX'
cmaesOptions . TolX = options_list { i , 2 } ;
case ' MaxFunEvals'
cmaesOptions . MaxFunEvals = options_list { i , 2 } ;
2013-10-08 15:53:55 +02:00
otherwise
2014-02-03 12:23:14 +01:00
warning ( [ ' cmaes: Unknown option (' options_list { i , 1 } ' )!' ] )
2013-10-08 15:53:55 +02:00
end
end
end
2012-03-03 15:30:55 +01:00
warning ( ' off' , ' CMAES:NonfinitenessRange' ) ;
warning ( ' off' , ' CMAES:InitialSigma' ) ;
2014-06-16 17:41:59 +02:00
[ x , fval , COUNTEVAL , STOPFLAG , OUT , BESTEVER ] = cmaes ( func2str ( objective_function ) , xparam1 , H0 , cmaesOptions , dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , oo_ ) ;
2012-03-03 15:30:55 +01:00
xparam1 = BESTEVER . x ;
2012-03-03 16:26:47 +01:00
disp ( sprintf ( ' \n Objective function at mode: %f' , fval ) )
2013-10-02 15:30:55 +02:00
case 10
2013-10-08 15:18:14 +02:00
simpsaOptions = options_ . simpsa ;
if isfield ( options_ , ' optim_opt' )
2014-02-03 12:23:14 +01:00
options_list = read_key_value_string ( options_ . optim_opt ) ;
for i = 1 : rows ( options_list )
switch options_list { i , 1 }
case ' MaxIter'
simpsaOptions . MAX_ITER_TOTAL = options_list { i , 2 } ;
case ' TolFun'
simpsaOptions . TOLFUN = options_list { i , 2 } ;
case ' TolX'
tolx = options_list { i , 2 } ;
2013-10-08 16:18:54 +02:00
if tolx < 0
2014-02-03 12:23:14 +01:00
simpsaOptions = rmfield ( simpsaOptions , ' TOLX' ) ; % Let simpsa choose the default.
2013-10-08 16:18:54 +02:00
else
simpsaOptions . TOLX = tolx ;
end
2014-02-03 12:23:14 +01:00
case ' EndTemparature'
simpsaOptions . TEMP_END = options_list { i , 2 } ;
case ' MaxFunEvals'
simpsaOptions . MAX_FUN_EVALS = options_list { i , 2 } ;
2013-10-08 15:18:14 +02:00
otherwise
2014-02-03 12:23:14 +01:00
warning ( [ ' simpsa: Unknown option (' options_list { i , 1 } ' )!' ] )
2013-10-08 15:18:14 +02:00
end
end
end
simpsaOptionsList = options2cell ( simpsaOptions ) ;
simpsaOptions = simpsaset ( simpsaOptionsList { : } ) ;
2014-06-16 17:41:59 +02:00
[ xparam1 , fval , exitflag ] = simpsa ( func2str ( objective_function ) , xparam1 , lb , ub , simpsaOptions , dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , oo_ ) ;
2014-02-03 12:23:14 +01:00
case 11
2013-06-18 16:18:55 +02:00
options_ . cova_compute = 0 ;
2014-06-16 17:41:59 +02:00
[ xparam1 , stdh , lb_95 , ub_95 , med_param ] = online_auxiliary_filter ( xparam1 , dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , oo_ ) ;
2012-03-03 16:26:47 +01:00
case 101
2009-12-16 18:17:34 +01:00
myoptions = soptions ;
myoptions ( 2 ) = 1e-6 ; %accuracy of argument
myoptions ( 3 ) = 1e-6 ; %accuracy of function (see Solvopt p.29)
myoptions ( 5 ) = 1.0 ;
2014-06-16 17:41:59 +02:00
[ xparam1 , fval ] = solvopt ( xparam1 , objective_function , [ ] , myoptions , dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , oo_ ) ;
2009-12-16 18:17:34 +01:00
case 102
%simulating annealing
% LB=zeros(size(xparam1))-20;
% UB=zeros(size(xparam1))+20;
LB = xparam1 - 1 ;
UB = xparam1 + 1 ;
neps = 10 ;
2011-09-17 15:38:49 +02:00
% Set input parameters.
2009-12-16 18:17:34 +01:00
maxy = 0 ;
2010-12-23 08:34:13 +01:00
epsilon = 1.0e-9 ;
2009-12-16 18:17:34 +01:00
rt_ = . 10 ;
t = 15.0 ;
ns = 10 ;
nt = 10 ;
maxevl = 100000000 ;
idisp = 1 ;
npar = length ( xparam1 ) ;
2011-09-17 15:38:49 +02:00
disp ( [ ' size of param' , num2str ( length ( xparam1 ) ) ] )
2009-12-16 18:17:34 +01:00
c = . 1 * ones ( npar , 1 ) ;
%* Set input values of the input/output parameters.*
2011-09-17 15:38:49 +02:00
2009-12-16 18:17:34 +01:00
vm = 1 * ones ( npar , 1 ) ;
disp ( [ ' number of parameters= ' num2str ( npar ) ' max= ' num2str ( maxy ) ' t= ' num2str ( t ) ] ) ;
2010-12-23 08:34:13 +01:00
disp ( [ ' rt_= ' num2str ( rt_ ) ' eps= ' num2str ( epsilon ) ' ns= ' num2str ( ns ) ] ) ;
2009-12-16 18:17:34 +01:00
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 ' ) ] ) ;
2011-09-17 15:38:49 +02:00
2011-09-19 18:10:29 +02:00
[ xparam1 , fval , nacc , nfcnev , nobds , ier , t , vm ] = sa ( objective_function , xparam1 , maxy , rt_ , epsilon , ns , nt ...
2014-06-16 17:41:59 +02:00
, neps , maxevl , LB , UB , c , idisp , t , vm , dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , oo_ ) ;
2009-12-16 18:17:34 +01:00
otherwise
if ischar ( options_ . mode_compute )
2014-06-16 17:41:59 +02:00
[ xparam1 , fval , retcode ] = feval ( options_ . mode_compute , objective_function , xparam1 , dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , oo_ ) ;
2009-12-16 18:17:34 +01:00
else
2011-09-19 18:10:29 +02:00
error ( [ ' dynare_estimation:: mode_compute = ' int2str ( options_ . mode_compute ) ' option is unknown!' ] )
2009-12-16 18:17:34 +01:00
end
end
2013-10-27 10:59:04 +01:00
if ~ isequal ( options_ . mode_compute , 6 ) %always already computes covariance matrix
if options_ . cova_compute == 1 %user did not request covariance not to be computed
2012-04-27 15:57:58 +02:00
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 , ...
2014-06-16 17:41:59 +02:00
dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , oo_ ) ;
2012-04-27 15:57:58 +02:00
options_ . analytic_derivation = ana_deriv ;
2014-07-23 18:20:16 +02:00
elseif ~ ( isequal ( options_ . mode_compute , 5 ) && flag == 0 ) ,
% with flag==0, we force to use the hessian from outer
% product gradient of optimizer 5
2012-04-27 15:57:58 +02:00
hh = reshape ( hessian ( objective_function , xparam1 , ...
2014-06-16 17:41:59 +02:00
options_ . gstep , dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , oo_ ) , nx , nx ) ;
2012-04-27 15:57:58 +02:00
end
2009-12-16 18:17:34 +01:00
end
2011-02-04 17:17:48 +01:00
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
2009-12-16 18:17:34 +01:00
end
if options_ . cova_compute == 0
2011-09-22 11:17:31 +02:00
hh = [ ] ; %NaN(length(xparam1),length(xparam1));
2009-12-16 18:17:34 +01:00
end
2013-10-27 10:59:04 +01:00
switch options_ . MCMC_jumping_covariance
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.' )
2014-02-03 12:23:14 +01:00
else
2013-10-27 10:59:04 +01:00
hh = diag ( 1. / ( bayestopt_ . p2 .^ 2 ) ) ;
end
case ' identity_matrix' %Use identity
2014-02-03 12:23:14 +01:00
hh = eye ( nx ) ;
2013-10-27 10:59:04 +01:00
otherwise %user specified matrix in file
2014-02-03 12:23:14 +01:00
try
2013-10-27 10:59:04 +01:00
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 ) ;
catch
error ( [ ' Specified jumping_covariance is not positive definite' ] )
end
end
2011-05-13 12:32:23 +02:00
if ~ options_ . mh_posterior_mode_estimation && options_ . cova_compute
2010-09-01 22:15:47 +02:00
try
chol ( hh ) ;
catch
2013-07-10 16:16:32 +02:00
skipline ( )
2010-09-01 22:15:47 +02:00
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.' )
2013-10-31 08:07:55 +01:00
disp ( ' You should try to change the initial values of the parameters using' )
2010-09-01 22:15:47 +02:00
disp ( ' the estimated_params_init block, or use another optimization routine.' )
2013-07-06 00:23:09 +02:00
params_at_bound = find ( xparam1 == ub | xparam1 == lb ) ;
if ~ isempty ( params_at_bound )
for ii = 1 : length ( params_at_bound )
2013-09-30 17:10:39 +02:00
params_at_bound_name { ii , 1 } = get_the_name ( params_at_bound ( ii ) , 0 , M_ , estim_params_ , options_ ) ;
2013-07-06 00:23:09 +02:00
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' )
fprintf ( ' - Use a different mode_compute like 6 or 9.\n' )
fprintf ( ' - Check whether the parameters estimated are identified.\n' )
fprintf ( ' - Increase the informativeness of the prior.\n' )
end
2010-09-01 22:15:47 +02:00
warning ( ' The results below are most likely wrong!' ) ;
end
2009-12-16 18:17:34 +01:00
end
2013-06-20 17:19:44 +02:00
if options_ . mode_check . status == 1 && ~ options_ . mh_posterior_mode_estimation
2012-07-05 10:22:36 +02:00
ana_deriv = options_ . analytic_derivation ;
options_ . analytic_derivation = 0 ;
2014-06-16 17:41:59 +02:00
mode_check ( objective_function , xparam1 , hh , dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , oo_ ) ;
2012-07-05 10:22:36 +02:00
options_ . analytic_derivation = ana_deriv ;
2009-12-16 18:17:34 +01:00
end
2012-04-21 21:28:03 +02:00
oo_ . posterior . optimization . mode = xparam1 ;
2013-03-17 22:49:28 +01:00
oo_ . posterior . optimization . Variance = [ ] ;
2011-05-13 14:34:53 +02:00
if ~ options_ . mh_posterior_mode_estimation
if options_ . cova_compute
invhess = inv ( hh ) ;
stdh = sqrt ( diag ( invhess ) ) ;
2013-03-17 22:49:28 +01:00
oo_ . posterior . optimization . Variance = invhess ;
2011-05-13 14:34:53 +02:00
end
2009-12-16 18:17:34 +01:00
else
2010-09-01 22:15:47 +02:00
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 ) ;
2009-12-16 18:17:34 +01:00
end
2011-05-13 12:32:23 +02:00
if ~ options_ . cova_compute
stdh = NaN ( length ( xparam1 ) , 1 ) ;
end
2009-12-16 18:17:34 +01:00
2011-02-10 15:54:23 +01:00
if any ( bayestopt_ . pshape > 0 ) && ~ options_ . mh_posterior_mode_estimation
2013-05-21 16:38:17 +02:00
% display results table and store parameter estimates and standard errors in results
2013-05-10 18:58:08 +02:00
oo_ = display_estimation_results_table ( xparam1 , stdh , M_ , options_ , estim_params_ , bayestopt_ , oo_ , pnames , ' Posterior' , ' posterior' ) ;
2013-05-21 16:38:17 +02:00
% Laplace approximation to the marginal log density:
2011-05-13 14:55:10 +02:00
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 ) ) ;
2014-06-16 17:41:59 +02:00
likelihood = feval ( objective_function , xparam1 , dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , oo_ ) ;
2011-09-22 11:17:31 +02:00
oo_ . MarginalDensity . LaplaceApproximation = . 5 * estim_params_nbr * log ( 2 * pi ) + . 5 * log_det_invhess - likelihood ;
2013-07-10 16:16:32 +02:00
skipline ( )
2011-09-22 11:17:31 +02:00
disp ( sprintf ( ' Log data density [Laplace approximation] is %f.' , oo_ . MarginalDensity . LaplaceApproximation ) )
2013-07-10 16:16:32 +02:00
skipline ( )
2009-12-16 18:17:34 +01:00
end
2012-03-15 15:32:31 +01:00
elseif ~ any ( bayestopt_ . pshape > 0 ) && ~ options_ . mh_posterior_mode_estimation
2013-05-10 18:58:08 +02:00
oo_ = display_estimation_results_table ( xparam1 , stdh , M_ , options_ , estim_params_ , bayestopt_ , oo_ , pnames , ' Maximum Likelihood' , ' mle' ) ;
2009-12-16 18:17:34 +01:00
end
2011-12-15 11:56:23 +01:00
OutputDirectoryName = CheckPath ( ' Output' , M_ . dname ) ;
2009-12-16 18:17:34 +01:00
if np > 0
pindx = estim_params_ . param_vals ( : , 1 ) ;
save ( [ M_ . fname ' _params.mat' ] , ' pindx' ) ;
end
2011-02-10 15:54:23 +01:00
if ( any ( bayestopt_ . pshape > 0 ) && options_ . mh_replic ) || ...
( any ( bayestopt_ . pshape > 0 ) && options_ . load_mh_file ) %% not ML estimation
2011-09-14 23:40:23 +02:00
bounds = prior_bounds ( bayestopt_ , options_ ) ;
2009-12-16 18:17:34 +01:00
bounds ( : , 1 ) = max ( bounds ( : , 1 ) , lb ) ;
bounds ( : , 2 ) = min ( bounds ( : , 2 ) , ub ) ;
bayestopt_ . lb = bounds ( : , 1 ) ;
bayestopt_ . ub = bounds ( : , 2 ) ;
2013-07-06 00:20:07 +02:00
outside_bound_pars = find ( xparam1 < bounds ( : , 1 ) | xparam1 > bounds ( : , 2 ) ) ;
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 , : } ] ;
2013-05-12 13:27:55 +02:00
end
2013-11-27 12:27:51 +01:00
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
2009-12-16 18:17:34 +01:00
end
% runs MCMC
if options_ . mh_replic
2011-02-10 15:54:23 +01:00
if options_ . load_mh_file && options_ . use_mh_covariance_matrix
2009-12-16 18:17:34 +01:00
invhess = compute_mh_covariance_matrix ;
end
2012-04-27 15:57:58 +02:00
ana_deriv = options_ . analytic_derivation ;
options_ . analytic_derivation = 0 ;
2011-05-13 14:55:10 +02:00
if options_ . cova_compute
2014-06-16 17:41:59 +02:00
feval ( options_ . posterior_sampling_method , objective_function , options_ . proposal_distribution , xparam1 , invhess , bounds , dataset_ , dataset_info , options_ , M_ , estim_params_ , bayestopt_ , oo_ ) ;
2009-12-16 18:17:34 +01:00
else
2013-06-13 14:24:17 +02:00
error ( ' I Cannot start the MCMC because the Hessian of the posterior kernel at the mode was not computed.' )
2009-12-16 18:17:34 +01:00
end
2012-04-27 15:57:58 +02:00
options_ . analytic_derivation = ana_deriv ;
2009-12-16 18:17:34 +01:00
end
2010-09-01 22:15:47 +02:00
if options_ . mh_posterior_mode_estimation
CutSample ( M_ , options_ , estim_params_ ) ;
return
2010-02-11 14:26:30 +01:00
else
2014-02-25 17:24:34 +01:00
if ~ options_ . nodiagnostic && options_ . mh_replic > 0
2013-09-16 17:19:56 +02:00
oo_ = McMCDiagnostics ( options_ , estim_params_ , M_ , oo_ ) ;
2009-12-16 18:17:34 +01:00
end
2010-09-01 22:15:47 +02:00
%% 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
2013-11-03 11:40:36 +01:00
[ marginal , oo_ ] = marginal_density ( M_ , options_ , estim_params_ , ...
oo_ ) ;
% Store posterior statistics by parameter name
2010-09-01 22:15:47 +02:00
oo_ = GetPosteriorParametersStatistics ( estim_params_ , M_ , options_ , bayestopt_ , oo_ ) ;
2012-11-29 14:26:39 +01:00
if ~ options_ . nograph
oo_ = PlotPosteriorDistributions ( estim_params_ , M_ , options_ , bayestopt_ , oo_ ) ;
end
2013-11-03 11:40:36 +01:00
% Store posterior mean in a vector and posterior variance in
% a matrix
2013-03-17 22:49:28 +01:00
[ oo_ . posterior . metropolis . mean , oo_ . posterior . metropolis . Variance ] ...
2012-04-21 21:28:03 +02:00
= GetPosteriorMeanVariance ( M_ , options_ . mh_drop ) ;
2010-09-01 22:15:47 +02:00
else
load ( [ M_ . fname ' _results' ] , ' oo_' ) ;
2009-12-16 18:17:34 +01:00
end
2013-12-18 16:39:41 +01:00
error_flag = metropolis_draw ( 1 ) ;
2010-09-01 22:15:47 +02:00
if options_ . bayesian_irf
2013-12-18 16:39:41 +01:00
if error_flag
error ( ' Estimation::mcmc: I cannot compute the posterior IRFs!' )
end
2010-09-01 22:15:47 +02:00
PosteriorIRF ( ' posterior' ) ;
2009-12-16 18:17:34 +01:00
end
2010-09-01 22:15:47 +02:00
if options_ . moments_varendo
2013-12-18 16:39:41 +01:00
if error_flag
error ( ' Estimation::mcmc: I cannot compute the posterior moments for the endogenous variables!' )
end
2010-09-01 22:15:47 +02:00
oo_ = compute_moments_varendo ( ' posterior' , options_ , M_ , oo_ , var_list_ ) ;
2009-12-16 18:17:34 +01:00
end
2011-02-10 15:54:23 +01:00
if options_ . smoother || ~ isempty ( options_ . filter_step_ahead ) || options_ . forecast
2013-12-18 16:39:41 +01:00
if error_flag
error ( ' Estimation::mcmc: I cannot compute the posterior statistics!' )
end
2014-06-23 12:24:36 +02:00
prior_posterior_statistics ( ' posterior' , dataset_ , dataset_info ) ;
2010-09-01 22:15:47 +02:00
end
xparam = get_posterior_parameters ( ' mean' ) ;
2012-06-07 14:15:59 +02:00
M_ = set_all_parameters ( xparam , estim_params_ , M_ ) ;
2010-09-13 21:34:15 +02:00
end
end
2010-09-01 22:15:47 +02:00
2012-06-06 17:04:27 +02:00
if options_ . particle . status
return
end
2011-02-10 15:54:23 +01:00
if ( ~ ( ( any ( bayestopt_ . pshape > 0 ) && options_ . mh_replic ) || ( any ( bayestopt_ . pshape ...
> 0 ) && options_ . load_mh_file ) ) ...
2012-04-18 21:15:18 +02:00
|| ~ options_ . smoother ) && options_ . partial_information == 0 % to be fixed
2010-09-13 21:34:15 +02:00
%% ML estimation, or posterior mode without metropolis-hastings or metropolis without bayesian smooth variable
2014-06-16 17:41:59 +02:00
[ atT , innov , measurement_error , updated_variables , ys , trend_coeff , aK , T , R , P , PK , decomp ] = DsgeSmoother ( xparam1 , dataset_ . nobs , transpose ( dataset_ . data ) , dataset_info . missing . aindex , dataset_info . missing . state ) ;
2010-09-13 21:34:15 +02:00
oo_ . Smoother . SteadyState = ys ;
oo_ . Smoother . TrendCoeffs = trend_coeff ;
2013-03-17 22:49:28 +01:00
oo_ . Smoother . Variance = P ;
2010-09-13 21:34:15 +02:00
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 '
2011-06-21 13:20:51 +02:00
i1 = dr . order_var ( bayestopt_ . smoother_var_list ( i ) ) ;
2011-11-02 14:02:12 +01:00
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
2010-09-13 21:34:15 +02:00
eval ( [ ' oo_.UpdatedVariables.' deblank ( M_ . endo_names ( i1 , : ) ) ...
' = updated_variables(i,:)' ' ;' ] ) ;
end
2012-08-27 17:38:22 +02:00
for i = 1 : M_ . exo_nbr
eval ( [ ' oo_.SmoothedShocks.' deblank ( M_ . exo_names ( i , : ) ) ' = innov(i,:)' ' ;' ] ) ;
end
2012-02-13 11:58:28 +01:00
if ~ options_ . nograph ,
[ nbplt , nr , nc , lr , lc , nstar ] = pltorg ( M_ . exo_nbr ) ;
2010-09-13 21:34:15 +02:00
if options_ . TeX
2012-02-13 11:58:28 +01:00
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' ) ;
2010-09-13 21:34:15 +02:00
end
2012-02-13 11:58:28 +01:00
for plt = 1 : nbplt ,
2013-09-10 17:40:56 +02:00
fh = dyn_figure ( options_ , ' Name' , ' Smoothed shocks' ) ;
2010-09-13 21:34:15 +02:00
NAMES = [ ] ;
if options_ . TeX , TeXNAMES = [ ] ; end
2012-02-13 11:58:28 +01:00
nstar0 = min ( nstar , M_ . exo_nbr - ( plt - 1 ) * nstar ) ;
2013-08-15 16:35:01 +02:00
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
2012-02-13 11:58:28 +01:00
for i = 1 : nstar0 ,
2010-09-13 21:34:15 +02:00
k = ( plt - 1 ) * nstar + i ;
subplot ( nr , nc , i ) ;
2013-08-15 16:35:01 +02:00
plot ( [ 1 gend ] , [ 0 0 ] , marker_string { 1 , 1 } , ' linewidth' , . 5 )
2010-09-13 21:34:15 +02:00
hold on
2013-08-15 16:35:01 +02:00
plot ( 1 : gend , innov ( k , : ) , marker_string { 2 , 1 } , ' linewidth' , 1 )
2010-09-13 21:34:15 +02:00
hold off
name = deblank ( M_ . exo_names ( k , : ) ) ;
2010-09-26 13:42:41 +02:00
if isempty ( NAMES )
NAMES = name ;
else
NAMES = char ( NAMES , name ) ;
end
2010-09-13 21:34:15 +02:00
if ~ isempty ( options_ . XTick )
set ( gca , ' XTick' , options_ . XTick )
set ( gca , ' XTickLabel' , options_ . XTickLabel )
2009-12-16 18:17:34 +01:00
end
2013-08-15 16:35:01 +02:00
if gend > 1
xlim ( [ 1 gend ] )
end
2010-09-13 21:34:15 +02:00
if options_ . TeX
texname = M_ . exo_names_tex ( k , : ) ;
2010-09-26 13:42:41 +02:00
if isempty ( TeXNAMES )
TeXNAMES = [ ' $ ' deblank ( texname ) ' $' ] ;
else
TeXNAMES = char ( TeXNAMES , [ ' $ ' deblank ( texname ) ' $' ] ) ;
end
2011-09-17 15:38:49 +02:00
end
2010-09-13 21:34:15 +02:00
title ( name , ' Interpreter' , ' none' )
2009-12-16 18:17:34 +01:00
end
2013-09-10 17:40:56 +02:00
dyn_saveas ( fh , [ M_ . fname ' _SmoothedShocks' int2str ( plt ) ] , options_ ) ;
2009-12-16 18:17:34 +01:00
if options_ . TeX
2010-09-13 21:34:15 +02:00
fprintf ( fidTeX , ' \\begin{figure}[H]\n' ) ;
2012-02-13 11:58:28 +01:00
for jj = 1 : nstar0
2010-09-13 21:34:15 +02:00
fprintf ( fidTeX , ' \\psfrag{%s}[1][][0.5][0]{%s}\n' , deblank ( NAMES ( jj , : ) ) , deblank ( TeXNAMES ( jj , : ) ) ) ;
2011-09-17 15:38:49 +02:00
end
2010-09-13 21:34:15 +02:00
fprintf ( fidTeX , ' \\centering \n' ) ;
fprintf ( fidTeX , ' \\includegraphics[scale=0.5]{%s_SmoothedShocks%s}\n' , M_ . fname , int2str ( plt ) ) ;
fprintf ( fidTeX , ' \\caption{Smoothed shocks.}' ) ;
fprintf ( fidTeX , ' \\label{Fig:SmoothedShocks:%s}\n' , int2str ( plt ) ) ;
fprintf ( fidTeX , ' \\end{figure}\n' ) ;
fprintf ( fidTeX , ' \n' ) ;
2011-09-17 15:38:49 +02:00
end
2010-09-13 21:34:15 +02:00
end
if options_ . TeX
fprintf ( fidTeX , ' \n' ) ;
fprintf ( fidTeX , ' %% End of TeX file.\n' ) ;
fclose ( fidTeX ) ;
2011-09-17 15:38:49 +02:00
end
2010-09-13 21:34:15 +02:00
end
%%
%% Smooth observational errors...
%%
2012-05-11 19:58:16 +02:00
if options_ . prefilter == 1
2014-06-16 17:41:59 +02:00
yf = atT ( bayestopt_ . mf , : ) + repmat ( dataset_info . descriptive . mean ' , 1 , gend ) ;
2012-05-11 19:58:16 +02:00
elseif options_ . loglinear == 1
2014-06-16 17:41:59 +02:00
gend
bayestopt_ . mfys
ys
2012-05-11 19:58:16 +02:00
yf = atT ( bayestopt_ . mf , : ) + repmat ( log ( ys ( bayestopt_ . mfys ) ) , 1 , gend ) + ...
trend_coeff * [ 1 : gend ] ;
2010-09-13 21:34:15 +02:00
else
2012-05-11 19:58:16 +02:00
yf = atT ( bayestopt_ . mf , : ) + repmat ( ys ( bayestopt_ . mfys ) , 1 , gend ) + ...
trend_coeff * [ 1 : gend ] ;
2010-09-13 21:34:15 +02:00
end
if nvn
number_of_plots_to_draw = 0 ;
index = [ ] ;
for i = 1 : n_varobs
2011-11-02 15:49:02 +01:00
if max ( abs ( measurement_error ) ) > 0.000000001
2010-09-13 21:34:15 +02:00
number_of_plots_to_draw = number_of_plots_to_draw + 1 ;
index = cat ( 1 , index , i ) ;
end
2013-05-21 16:38:17 +02:00
eval ( [ ' oo_.SmoothedMeasurementErrors.' options_ . varobs { i } ' = measurement_error(i,:)' ' ;' ] ) ;
2010-09-13 21:34:15 +02:00
end
2012-02-13 11:58:28 +01:00
if ~ options_ . nograph
[ nbplt , nr , nc , lr , lc , nstar ] = pltorg ( number_of_plots_to_draw ) ;
2010-09-13 21:34:15 +02:00
if options_ . TeX
2012-02-13 11:58:28 +01:00
fidTeX = fopen ( [ M_ . fname ' _SmoothedObservationErrors.TeX' ] , ' w' ) ;
fprintf ( fidTeX , ' %% TeX eps-loader file generated by dynare_estimation.m (Dynare).\n' ) ;
fprintf ( fidTeX , [ ' %% ' datestr ( now , 0 ) ' \n' ] ) ;
fprintf ( fidTeX , ' \n' ) ;
2010-09-13 21:34:15 +02:00
end
2011-11-02 15:49:02 +01:00
for plt = 1 : nbplt
2013-09-10 17:40:56 +02:00
fh = dyn_figure ( options_ , ' Name' , ' Smoothed observation errors' ) ;
2009-12-16 18:17:34 +01:00
NAMES = [ ] ;
if options_ . TeX , TeXNAMES = [ ] ; end
2013-03-19 21:12:49 +01:00
nstar0 = min ( nstar , number_of_plots_to_draw - ( nbplt - 1 ) * nstar ) ;
2013-08-15 16:35:01 +02:00
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
2013-03-19 21:12:49 +01:00
for i = 1 : nstar0
2010-09-13 21:34:15 +02:00
k = ( plt - 1 ) * nstar + i ;
subplot ( nr , nc , i ) ;
2013-08-15 16:35:01 +02:00
plot ( [ 1 gend ] , [ 0 0 ] , marker_string { 1 , 1 } , ' linewidth' , . 5 )
2009-12-16 18:17:34 +01:00
hold on
2013-08-15 16:35:01 +02:00
plot ( 1 : gend , measurement_error ( index ( k ) , : ) , marker_string { 2 , 1 } , ' linewidth' , 1 )
2009-12-16 18:17:34 +01:00
hold off
2013-05-21 16:38:17 +02:00
name = options_ . varobs { index ( k ) } ;
2013-08-15 16:35:01 +02:00
if gend > 1
xlim ( [ 1 gend ] )
end
2010-09-26 13:42:41 +02:00
if isempty ( NAMES )
NAMES = name ;
else
NAMES = char ( NAMES , name ) ;
end
2009-12-16 18:17:34 +01:00
if ~ isempty ( options_ . XTick )
set ( gca , ' XTick' , options_ . XTick )
set ( gca , ' XTickLabel' , options_ . XTickLabel )
end
if options_ . TeX
2013-05-21 16:38:17 +02:00
idx = strmatch ( options_ . varobs { index ( k ) } , M_ . endo_names , ' exact' ) ;
2010-09-13 21:34:15 +02:00
texname = M_ . endo_names_tex ( idx , : ) ;
2010-09-26 13:42:41 +02:00
if isempty ( TeXNAMES )
TeXNAMES = [ ' $ ' deblank ( texname ) ' $' ] ;
else
TeXNAMES = char ( TeXNAMES , [ ' $ ' deblank ( texname ) ' $' ] ) ;
end
end
2009-12-16 18:17:34 +01:00
title ( name , ' Interpreter' , ' none' )
end
2013-09-10 17:40:56 +02:00
dyn_saveas ( fh , [ M_ . fname ' _SmoothedObservationErrors' int2str ( plt ) ] , options_ ) ;
2009-12-16 18:17:34 +01:00
if options_ . TeX
fprintf ( fidTeX , ' \\begin{figure}[H]\n' ) ;
2013-03-19 21:12:49 +01:00
for jj = 1 : nstar0
2009-12-16 18:17:34 +01:00
fprintf ( fidTeX , ' \\psfrag{%s}[1][][0.5][0]{%s}\n' , deblank ( NAMES ( jj , : ) ) , deblank ( TeXNAMES ( jj , : ) ) ) ;
2011-09-17 15:38:49 +02:00
end
2009-12-16 18:17:34 +01:00
fprintf ( fidTeX , ' \\centering \n' ) ;
2010-09-13 21:34:15 +02:00
fprintf ( fidTeX , ' \\includegraphics[scale=0.5]{%s_SmoothedObservationErrors%s}\n' , M_ . fname , int2str ( plt ) ) ;
fprintf ( fidTeX , ' \\caption{Smoothed observation errors.}' ) ;
fprintf ( fidTeX , ' \\label{Fig:SmoothedObservationErrors:%s}\n' , int2str ( plt ) ) ;
2009-12-16 18:17:34 +01:00
fprintf ( fidTeX , ' \\end{figure}\n' ) ;
fprintf ( fidTeX , ' \n' ) ;
2011-09-17 15:38:49 +02:00
end
2009-12-16 18:17:34 +01:00
end
2012-02-13 11:58:28 +01:00
if options_ . TeX
fprintf ( fidTeX , ' \n' ) ;
fprintf ( fidTeX , ' %% End of TeX file.\n' ) ;
fclose ( fidTeX ) ;
end
2010-09-13 21:34:15 +02:00
end
2011-09-17 15:38:49 +02:00
end
2010-09-13 21:34:15 +02:00
%%
%% Historical and smoothed variabes
%%
2012-02-13 11:58:28 +01:00
if ~ options_ . nograph
2010-09-13 21:34:15 +02:00
[ nbplt , nr , nc , lr , lc , nstar ] = pltorg ( n_varobs ) ;
if options_ . TeX
fidTeX = fopen ( [ M_ . fname ' _HistoricalAndSmoothedVariables.TeX' ] , ' w' ) ;
fprintf ( fidTeX , ' %% TeX eps-loader file generated by dynare_estimation.m (Dynare).\n' ) ;
fprintf ( fidTeX , [ ' %% ' datestr ( now , 0 ) ' \n' ] ) ;
fprintf ( fidTeX , ' \n' ) ;
2011-09-17 15:38:49 +02:00
end
2012-02-13 11:58:28 +01:00
for plt = 1 : nbplt ,
2013-09-10 17:40:56 +02:00
fh = dyn_figure ( options_ , ' Name' , ' Historical and smoothed variables' ) ;
2010-09-13 21:34:15 +02:00
NAMES = [ ] ;
if options_ . TeX , TeXNAMES = [ ] ; end
2012-02-13 11:58:28 +01:00
nstar0 = min ( nstar , n_varobs - ( plt - 1 ) * nstar ) ;
2013-08-15 16:35:01 +02:00
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
2012-02-13 11:58:28 +01:00
for i = 1 : nstar0 ,
k = ( plt - 1 ) * nstar + i ;
2010-09-13 21:34:15 +02:00
subplot ( nr , nc , i ) ;
2013-08-15 16:35:01 +02:00
plot ( 1 : gend , yf ( k , : ) , marker_string { 1 , 1 } , ' linewidth' , 1 )
2010-09-13 21:34:15 +02:00
hold on
2013-08-15 16:35:01 +02:00
plot ( 1 : gend , rawdata ( : , k ) , marker_string { 2 , 1 } , ' linewidth' , 1 )
2010-09-13 21:34:15 +02:00
hold off
2013-05-21 16:38:17 +02:00
name = options_ . varobs { k } ;
2010-09-26 13:42:41 +02:00
if isempty ( NAMES )
NAMES = name ;
else
NAMES = char ( NAMES , name ) ;
end
2010-09-13 21:34:15 +02:00
if ~ isempty ( options_ . XTick )
set ( gca , ' XTick' , options_ . XTick )
set ( gca , ' XTickLabel' , options_ . XTickLabel )
2009-12-16 18:17:34 +01:00
end
2013-08-15 16:35:01 +02:00
if gend > 1
xlim ( [ 1 gend ] )
end
2010-09-13 21:34:15 +02:00
if options_ . TeX
2013-05-21 16:38:17 +02:00
idx = strmatch ( options_ . varobs { k } , M_ . endo_names , ' exact' ) ;
2010-09-13 21:34:15 +02:00
texname = M_ . endo_names_tex ( idx , : ) ;
2010-09-26 13:42:41 +02:00
if isempty ( TeXNAMES )
TeXNAMES = [ ' $ ' deblank ( texname ) ' $' ] ;
else
TeXNAMES = char ( TeXNAMES , [ ' $ ' deblank ( texname ) ' $' ] ) ;
end
2010-09-13 21:34:15 +02:00
end
2012-02-13 11:58:28 +01:00
title ( name , ' Interpreter' , ' none' )
2010-09-13 21:34:15 +02:00
end
2013-09-10 17:40:56 +02:00
dyn_saveas ( fh , [ M_ . fname ' _HistoricalAndSmoothedVariables' int2str ( plt ) ] , options_ ) ;
2010-09-13 21:34:15 +02:00
if options_ . TeX
fprintf ( fidTeX , ' \\begin{figure}[H]\n' ) ;
2012-02-13 11:58:28 +01:00
for jj = 1 : nstar0 ,
2010-09-13 21:34:15 +02:00
fprintf ( fidTeX , ' \\psfrag{%s}[1][][0.5][0]{%s}\n' , deblank ( NAMES ( jj , : ) ) , deblank ( TeXNAMES ( jj , : ) ) ) ;
2011-09-17 15:38:49 +02:00
end
2010-09-13 21:34:15 +02:00
fprintf ( fidTeX , ' \\centering \n' ) ;
2012-02-13 11:58:28 +01:00
fprintf ( fidTeX , ' \\includegraphics[scale=0.5]{%s_HistoricalAndSmoothedVariables%s}\n' , M_ . fname , int2str ( plt ) ) ;
2010-09-13 21:34:15 +02:00
fprintf ( fidTeX , ' \\caption{Historical and smoothed variables.}' ) ;
2012-02-13 11:58:28 +01:00
fprintf ( fidTeX , ' \\label{Fig:HistoricalAndSmoothedVariables:%s}\n' , int2str ( plt ) ) ;
2010-09-13 21:34:15 +02:00
fprintf ( fidTeX , ' \\end{figure}\n' ) ;
fprintf ( fidTeX , ' \n' ) ;
2009-12-16 18:17:34 +01:00
end
end
2012-02-13 11:58:28 +01:00
if options_ . TeX
fprintf ( fidTeX , ' \n' ) ;
fprintf ( fidTeX , ' %% End of TeX file.\n' ) ;
fclose ( fidTeX ) ;
end
end
2010-09-01 22:15:47 +02:00
end
2010-09-13 21:34:15 +02:00
2011-09-17 15:38:49 +02:00
if options_ . forecast > 0 && options_ . mh_replic == 0 && ~ options_ . load_mh_file
2011-12-11 11:15:38 +01:00
dyn_forecast ( var_list_ , ' smoother' ) ;
2009-12-16 18:17:34 +01:00
end
if np > 0
pindx = estim_params_ . param_vals ( : , 1 ) ;
save ( [ M_ . fname ' _pindx.mat' ] , ' pindx' ) ;
end