diff --git a/doc/dynare.texi b/doc/dynare.texi
index 1ca3bc42e..0dc42dacf 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -5410,6 +5410,10 @@ Uses the simpsa algorithm, based on the combination of the non-linear simplex an
@item 11
This is not strictly speaking an optimization algorithm. The (estimated) parameters are treated as state variables and estimated jointly with the original state variables of the model using a nonlinear filter. The algorithm implemented in Dynare is described in @cite{Liu and West (2001)}.
+@item 12
+Uses @code{particleswarm} optimization routine (available under MATLAB if
+the Global Optimization Toolbox is installed; not available under Octave).
+
@item 101
Uses the SolveOpt algorithm for local nonlinear optimization problems proposed by
@cite{Kuntsevich and Kappel (1997)}.
@@ -5513,7 +5517,7 @@ A list of @var{NAME} and @var{VALUE} pairs. Can be used to set options for the o
@table @code
-@item 1, 3, 7
+@item 1, 3, 7, 12
Available options are given in the documentation of the MATLAB Optimization Toolbox or in Octave's documentation.
@item 2
@@ -6165,6 +6169,16 @@ This is the convergence criterion used in the fixed point Lyapunov solver. Its d
@anchor{lyapunov_doubling_tol}
This is the convergence criterion used in the doubling algorithm to solve the Lyapunov equation. Its default value is 1e-16.
+@item use_penalized_objective_for_hessian
+Use the penalized objective instead of the objective function to compute
+numerically the hessian matrix at the mode. The penalties decrease the value of
+the posterior density (or likelihood) when, for some perturbations, Dynare is
+not able to solve the model (issues with steady state existence, Blanchard and
+Kahn conditions, ...). In pratice, the penalized and original
+objectives will only differ if the posterior mode is found to be near a region
+where the model is ill-behaved. By default the original objective function is
+used.
+
@item analytic_derivation
Triggers estimation with analytic gradient. The final hessian is also
computed analytically. Only works for stationary models without
diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m
index 8a4a7d6b6..d911adc6e 100644
--- a/matlab/dsge_likelihood.m
+++ b/matlab/dsge_likelihood.m
@@ -141,7 +141,10 @@ trend_coeff = [];
exit_flag = 1;
info = zeros(4,1);
DLIK = [];
-Hess = [];
+Hess = [];
+
+% Ensure that xparam1 is a column vector.
+xparam1 = xparam1(:);
if DynareOptions.estimation_dll
[fval,exit_flag,SteadyState,trend_coeff,info,params,H,Q] ...
diff --git a/matlab/dsge_var_likelihood.m b/matlab/dsge_var_likelihood.m
index 646120c5f..b97beef26 100644
--- a/matlab/dsge_var_likelihood.m
+++ b/matlab/dsge_var_likelihood.m
@@ -37,7 +37,7 @@ function [fval,info,exit_flag,grad,hess,SteadyState,trend_coeff,PHI_tilde,SIGMA_
% SPECIAL REQUIREMENTS
% None.
-% Copyright (C) 2006-2016 Dynare Team
+% Copyright (C) 2006-2017 Dynare Team
%
% This file is part of Dynare.
%
@@ -68,6 +68,9 @@ iXX = [];
prior = [];
trend_coeff=[];
+% Ensure that xparam1 is a column vector.
+xparam1 = xparam1(:);
+
% Initialization of of the index for parameter dsge_prior_weight in Model.params.
if isempty(dsge_prior_weight_idx)
dsge_prior_weight_idx = strmatch('dsge_prior_weight',Model.param_names);
diff --git a/matlab/dynare_config.m b/matlab/dynare_config.m
index fc89f1539..c70e0a89f 100644
--- a/matlab/dynare_config.m
+++ b/matlab/dynare_config.m
@@ -91,6 +91,15 @@ if isoctave
p{end+1} = '/missing/ordeig';
end
+if isoctave && ~isequal(supported_octave_version(), version())
+ skipline()
+ warning(['This version of Octave is not supported. Consider installing ' ...
+ 'version %s of Octave,\notherwise m files will be used instead ' ...
+ 'of precompiled mex files and some features, like solution\n' ...
+ 'of models approximated at third order, will not be available.'], supported_octave_version())
+ skipline()
+end
+
% ilu is missing in Octave < 4.0
if isoctave && octave_ver_less_than('4.0')
p{end+1} = '/missing/ilu';
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index 6679f0cee..6077501ca 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -220,10 +220,14 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
options_.analytic_derivation = ana_deriv_old;
elseif ~isnumeric(options_.mode_compute) || ~(isequal(options_.mode_compute,5) && newratflag~=1),
- % with flag==0, we force to use the hessian from outer
- % product gradient of optimizer 5
- hh = reshape(penalty_hessian(objective_function,xparam1,fval, ...
- options_.gstep,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_),nx,nx);
+ % with flag==0, we force to use the hessian from outer product gradient of optimizer 5
+ if options_.hessian.use_penalized_objective
+ penalized_objective_function = str2func('penalty_objective_function');
+ hh = hessian(penalized_objective_function, xparam1, options_.gstep, objective_function, fval, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_);
+ else
+ hh = hessian(objective_function, xparam1, options_.gstep, dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds,oo_);
+ end
+ hh = reshape(hh, nx, nx);
elseif isnumeric(options_.mode_compute) && isequal(options_.mode_compute,5)
% other numerical hessian options available with optimizer 5
%
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index e65ca4a62..e8d1627e2 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -565,6 +565,11 @@ options_.homotopy_mode = 0;
options_.homotopy_steps = 1;
options_.homotopy_force_continue = 0;
+% numerical hessian
+hessian.use_penalized_objective = false;
+
+options_.hessian = hessian;
+
%csminwel optimization routine
csminwel.tolerance.f=1e-7;
csminwel.maxiter=1000;
@@ -645,6 +650,23 @@ options_.saopt.nt=10;
options_.saopt.step_length_c=0.1;
options_.saopt.initial_step_length=1;
+% particleswarm (global optimization toolbox needed)
+particleswarm.Display = 'iter';
+particleswarm.DisplayInterval = 1;
+particleswarm.FunctionTolerance = 1e-6;
+particleswarm.FunValCheck = 'on';
+particleswarm.HybridFcn = [];
+particleswarm.InertiaRange = [0.1, 1.1];
+particleswarm.MaxIterations = 100000;
+particleswarm.MaxStallIterations = 20;
+particleswarm.MaxStallTime = Inf;
+particleswarm.MaxTime = Inf;
+particleswarm.MinNeighborsFraction = .25;
+particleswarm.ObjectiveLimit = -Inf;
+particleswarm.UseParallel = false;
+particleswarm.UseVectorized = false;
+options_.particleswarm = particleswarm;
+
% prior analysis
options_.prior_mc = 20000;
options_.prior_analysis_endo_var_list = [];
diff --git a/matlab/hessian.m b/matlab/hessian.m
index 13ddc2f44..4faa4c500 100644
--- a/matlab/hessian.m
+++ b/matlab/hessian.m
@@ -1,4 +1,4 @@
-function hessian_mat = hessian(func,x,gstep,varargin) % --*-- Unitary tests --*--
+function hessian_mat = hessian(func,x, gstep, varargin) % --*-- Unitary tests --*--
% Computes second order partial derivatives
%
@@ -26,7 +26,7 @@ function hessian_mat = hessian(func,x,gstep,varargin) % --*-- Unitary tests --*-
% none
%
-% Copyright (C) 2001-2014 Dynare Team
+% Copyright (C) 2001-2017 Dynare Team
%
% This file is part of Dynare.
%
@@ -46,48 +46,55 @@ function hessian_mat = hessian(func,x,gstep,varargin) % --*-- Unitary tests --*-
if ~isa(func, 'function_handle')
func = str2func(func);
end
-n=size(x,1);
-h1=max(abs(x),sqrt(gstep(1))*ones(n,1))*eps^(1/6)*gstep(2);
-h_1=h1;
-xh1=x+h1;
-h1=xh1-x;
-xh1=x-h_1;
-h_1=x-xh1;
-xh1=x;
-f0=feval(func,x,varargin{:});
-f1=zeros(size(f0,1),n);
-f_1=f1;
+
+n = size(x,1);
+h1 = max(abs(x), sqrt(gstep(1))*ones(n, 1))*eps^(1/6)*gstep(2);
+h_1 = h1;
+xh1 = x+h1;
+h1 = xh1-x;
+xh1 = x-h_1;
+h_1 = x-xh1;
+xh1 = x;
+f0 = feval(func, x, varargin{:});
+f1 = zeros(size(f0, 1), n);
+f_1 = f1;
+
for i=1:n
%do step up
- xh1(i)=x(i)+h1(i);
- f1(:,i)=feval(func,xh1,varargin{:});
- %do step up
- xh1(i)=x(i)-h_1(i);
- f_1(:,i)=feval(func,xh1,varargin{:});
- xh1(i)=x(i);%reset parameter
+ xh1(i) = x(i)+h1(i);
+ f1(:,i) = feval(func, xh1, varargin{:});
+ %do step down
+ xh1(i) = x(i)-h_1(i);
+ f_1(:,i) = feval(func, xh1, varargin{:});
+ %reset parameter
+ xh1(i) = x(i);
end
-xh_1=xh1;
-hessian_mat = zeros(size(f0,1),n*n);
-temp=f1+f_1-f0*ones(1,n); %term f_(1,0)+f_(-1,0)-f_(0,0) used later
+
+xh_1 = xh1;
+temp = f1+f_1-f0*ones(1, n); %term f_(1,0)+f_(-1,0)-f_(0,0) used later
+
+hessian_mat = zeros(size(f0,1), n*n);
+
for i=1:n
- if i > 1 %fill symmetric part of Hessian based on previously computed results
- k=[i:n:n*(i-1)];
- hessian_mat(:,(i-1)*n+1:(i-1)*n+i-1)=hessian_mat(:,k);
+ if i > 1
+ %fill symmetric part of Hessian based on previously computed results
+ k = [i:n:n*(i-1)];
+ hessian_mat(:,(i-1)*n+1:(i-1)*n+i-1) = hessian_mat(:,k);
end
- hessian_mat(:,(i-1)*n+i)=(f1(:,i)+f_1(:,i)-2*f0)./(h1(i)*h_1(i)); %formula 25.3.23
+ hessian_mat(:,(i-1)*n+i) = (f1(:,i)+f_1(:,i)-2*f0)./(h1(i)*h_1(i)); %formula 25.3.23
for j=i+1:n
%step in up direction
- xh1(i)=x(i)+h1(i);
- xh1(j)=x(j)+h_1(j);
+ xh1(i) = x(i)+h1(i);
+ xh1(j) = x(j)+h_1(j);
%step in down direction
- xh_1(i)=x(i)-h1(i);
- xh_1(j)=x(j)-h_1(j);
- hessian_mat(:,(i-1)*n+j)=-(-feval(func,xh1,varargin{:})-feval(func,xh_1,varargin{:})+temp(:,i)+temp(:,j))./(2*h1(i)*h_1(j)); %formula 25.3.27
+ xh_1(i) = x(i)-h1(i);
+ xh_1(j) = x(j)-h_1(j);
+ hessian_mat(:,(i-1)*n+j) =-(-feval(func, xh1, varargin{:})-feval(func, xh_1, varargin{:})+temp(:,i)+temp(:,j))./(2*h1(i)*h_1(j)); %formula 25.3.27
%reset grid points
- xh1(i)=x(i);
- xh1(j)=x(j);
- xh_1(i)=x(i);
- xh_1(j)=x(j);
+ xh1(i) = x(i);
+ xh1(j) = x(j);
+ xh_1(i) = x(i);
+ xh_1(j) = x(j);
end
end
diff --git a/matlab/non_linear_dsge_likelihood.m b/matlab/non_linear_dsge_likelihood.m
index f8ac92e23..17d448927 100644
--- a/matlab/non_linear_dsge_likelihood.m
+++ b/matlab/non_linear_dsge_likelihood.m
@@ -107,7 +107,7 @@ function [fval,info,exit_flag,DLIK,Hess,ys,trend_coeff,Model,DynareOptions,Bayes
%! @end deftypefn
%@eod:
-% Copyright (C) 2010-2016 Dynare Team
+% Copyright (C) 2010-2017 Dynare Team
%
% This file is part of Dynare.
%
@@ -124,9 +124,6 @@ function [fval,info,exit_flag,DLIK,Hess,ys,trend_coeff,Model,DynareOptions,Bayes
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see .
-% AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT fr
-% frederic DOT karame AT univ DASH lemans DOT fr
-
% Declaration of the penalty as a persistent variable.
persistent init_flag
persistent restrict_variables_idx observed_variables_idx state_variables_idx mf0 mf1
@@ -140,6 +137,9 @@ exit_flag = 1;
DLIK = [];
Hess = [];
+% Ensure that xparam1 is a column vector.
+xparam1 = xparam1(:);
+
% Issue an error if loglinear option is used.
if DynareOptions.loglinear
error('non_linear_dsge_likelihood: It is not possible to use a non linear filter with the option loglinear!')
diff --git a/matlab/optimization/dynare_minimize_objective.m b/matlab/optimization/dynare_minimize_objective.m
index 4724fc91a..0aebc6b50 100644
--- a/matlab/optimization/dynare_minimize_objective.m
+++ b/matlab/optimization/dynare_minimize_objective.m
@@ -401,8 +401,58 @@ switch minimizer_algorithm
[LB, UB]=set_bounds_to_finite_values(bounds, options_.huge_number);
[opt_par_values, fval, exitflag] = simpsa(func2str(objective_function),start_par_value,LB,UB,simpsaOptions,varargin{:});
case 11
- options_.cova_compute = 0 ;
- [opt_par_values,stdh,lb_95,ub_95,med_param] = online_auxiliary_filter(start_par_value,varargin{:}) ;
+ options_.cova_compute = 0;
+ [opt_par_values, stdh, lb_95, ub_95, med_param] = online_auxiliary_filter(start_par_value, varargin{:});
+ case 12
+ [LB, UB] = set_bounds_to_finite_values(bounds, options_.huge_number);
+ tmp = transpose([fieldnames(options_.particleswarm), struct2cell(options_.particleswarm)]);
+ particleswarmOptions = optimoptions(@particleswarm);
+ particleswarmOptions = optimoptions(particleswarmOptions, tmp{:});
+ if ~isempty(options_.optim_opt)
+ options_list = read_key_value_string(options_.optim_opt);
+ SupportedListOfOptions = {'CreationFcn', 'Display', 'DisplayInterval', 'FunctionTolerance', ...
+ 'FunValCheck', 'HybridFcn', 'InertiaRange', 'InitialSwarmMatrix', 'InitialSwarmSpan', ...
+ 'MaxIterations', 'MaxStallIterations', 'MaxStallTime', 'MaxTime', ...
+ 'MinNeighborsFraction', 'ObjectiveLimit', 'OutputFcn', 'PlotFcn', 'SelfAdjustmentWeight', ...
+ 'SocialAdjustmentWeight', 'SwarmSize', 'UseParallel', 'UseVectorized'};
+ for i=1:rows(options_list)
+ if ismember(options_list{i,1}, SupportedListOfOptions)
+ particleswarmOptions = optimoptions(particleswarmOptions, options_list{i,1}, options_list{i,2});
+ else
+ warning(['particleswarm: Unknown option (' options_list{i,1} ')!'])
+ end
+ end
+ end
+ % Get number of instruments.
+ numberofvariables = length(start_par_value);
+ % Set objective function.
+ objfun = @(x) objective_function(x, varargin{:});
+ if ischar(particleswarmOptions.SwarmSize)
+ eval(['particleswarmOptions.SwarmSize = ' particleswarmOptions.SwarmSize ';'])
+ end
+ if isempty(particleswarmOptions.InitialSwarmMatrix)
+ particleswarmOptions.InitialSwarmMatrix = zeros(particleswarmOptions.SwarmSize, numberofvariables);
+ p = 1;
+ FVALS = zeros(particleswarmOptions.SwarmSize, 1);
+ while p<=particleswarmOptions.SwarmSize
+ candidate = rand(numberofvariables, 1).*(UB-LB)+LB;
+ [fval, info, exit_flag] = objfun(candidate);
+ if exit_flag
+ particleswarmOptions.InitialSwarmMatrix(p,:) = transpose(candidate);
+ FVALS(p) = fval;
+ p = p + 1;
+ end
+ end
+ end
+ % Set penalty to the worst value of the objective function.
+ TMP = [particleswarmOptions.InitialSwarmMatrix, FVALS];
+ TMP = sortrows(TMP, length(start_par_value)+1);
+ penalty = TMP(end,end);
+ % Define penalized objective.
+ objfun = @(x) penalty_objective_function(x, objective_function, penalty, varargin{:});
+ % Minimize the penalized objective (note that the penalty is not updated).
+ [opt_par_values, fval, exitflag, output] = particleswarm(objfun, length(start_par_value), LB, UB, particleswarmOptions);
+ opt_par_values = opt_par_values(:);
case 101
solveoptoptions = options_.solveopt;
if ~isempty(options_.optim_opt)
diff --git a/matlab/optimization/penalty_hessian.m b/matlab/optimization/penalty_hessian.m
deleted file mode 100644
index 8da7bf6a9..000000000
--- a/matlab/optimization/penalty_hessian.m
+++ /dev/null
@@ -1,94 +0,0 @@
-function hessian_mat = penalty_hessian(func,x,penalty,gstep,varargin) % --*-- Unitary tests --*--
-
-% Computes second order partial derivatives with penalty_objective_function
-%
-% INPUTS
-% func [string] name of the function
-% x [double] vector, the Hessian of "func" is evaluated at x.
-% penalty [double] penalty base used if function fails
-% gstep [double] scalar, size of epsilon.
-% varargin [void] list of additional arguments for "func".
-%
-% OUTPUTS
-% hessian_mat [double] Hessian matrix
-%
-% ALGORITHM
-% Uses Abramowitz and Stegun (1965) formulas 25.3.23
-% \[
-% \frac{\partial^2 f_{0,0}}{\partial {x^2}} = \frac{1}{h^2}\left( f_{1,0} - 2f_{0,0} + f_{ - 1,0} \right)
-% \]
-% and 25.3.27 p. 884
-%
-% \[
-% \frac{\partial ^2f_{0,0}}{\partial x\partial y} = \frac{-1}{2h^2}\left(f_{1,0} + f_{-1,0} + f_{0,1} + f_{0,-1} - 2f_{0,0} - f_{1,1} - f_{-1,-1} \right)
-% \]
-%
-% SPECIAL REQUIREMENTS
-% none
-%
-
-% Copyright (C) 2001-2014 Dynare Team
-%
-% 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 .
-
-if ~isa(func, 'function_handle')
- func = str2func(func);
-end
-n=size(x,1);
-h1=max(abs(x),sqrt(gstep(1))*ones(n,1))*eps^(1/6)*gstep(2);
-h_1=h1;
-xh1=x+h1;
-h1=xh1-x;
-xh1=x-h_1;
-h_1=x-xh1;
-xh1=x;
-f0=penalty_objective_function(x,func,penalty,varargin{:});
-f1=zeros(size(f0,1),n);
-f_1=f1;
-for i=1:n
- %do step up
- xh1(i)=x(i)+h1(i);
- f1(:,i)=penalty_objective_function(xh1,func,penalty,varargin{:});
- %do step up
- xh1(i)=x(i)-h_1(i);
- f_1(:,i)=penalty_objective_function(xh1,func,penalty,varargin{:});
- xh1(i)=x(i);%reset parameter
-end
-xh_1=xh1;
-hessian_mat = zeros(size(f0,1),n*n);
-temp=f1+f_1-f0*ones(1,n); %term f_(1,0)+f_(-1,0)-f_(0,0) used later
-for i=1:n
- if i > 1 %fill symmetric part of Hessian based on previously computed results
- k=[i:n:n*(i-1)];
- hessian_mat(:,(i-1)*n+1:(i-1)*n+i-1)=hessian_mat(:,k);
- end
- hessian_mat(:,(i-1)*n+i)=(f1(:,i)+f_1(:,i)-2*f0)./(h1(i)*h_1(i)); %formula 25.3.23
- for j=i+1:n
- %step in up direction
- xh1(i)=x(i)+h1(i);
- xh1(j)=x(j)+h_1(j);
- %step in down direction
- xh_1(i)=x(i)-h1(i);
- xh_1(j)=x(j)-h_1(j);
- hessian_mat(:,(i-1)*n+j)=-(-penalty_objective_function(xh1,func,penalty,varargin{:})-penalty_objective_function(xh_1,func,penalty,varargin{:})+temp(:,i)+temp(:,j))./(2*h1(i)*h_1(j)); %formula 25.3.27
- %reset grid points
- xh1(i)=x(i);
- xh1(j)=x(j);
- xh_1(i)=x(i);
- xh_1(j)=x(j);
- end
-end
-
diff --git a/matlab/optimization/penalty_objective_function.m b/matlab/optimization/penalty_objective_function.m
index 639a248aa..c245800c2 100644
--- a/matlab/optimization/penalty_objective_function.m
+++ b/matlab/optimization/penalty_objective_function.m
@@ -1,7 +1,19 @@
-function [fval,exit_flag,arg1,arg2] = penalty_objective_function(x0,fcn,penalty,varargin)
-%function [fval,exit_flag,arg1,arg2] = penalty_objective_function(x0,fcn,penalty,varargin)
+function [fval, exit_flag, arg1, arg2] = penalty_objective_function(x, fcn, base_penalty, varargin)
-% Copyright (C) 2016 Dynare Team
+% Encapsulates an objective function to be minimized, adding a penalty if necessary.
+%
+% INPUTS
+% - x [double] n*1 vector of instrument values.
+% - fcn [fhandle] objective function.
+% - base_penalty [double] scalar, base of the penality (typically the value of the objective at the previous iteration).
+% - varagin [cell] additional parameters for fcn.
+%
+% OUTPUTS
+% - fval [double] scalar, value of the objective function at x.
+% - exit_flag [integer] scalar, flag returned by fcn (third output).
+% - arg1, arg2 fourth and fifth output arguments of the objective function.
+
+% Copyright (C) 2016-2017 Dynare Team
%
% This file is part of Dynare.
%
@@ -18,9 +30,8 @@ function [fval,exit_flag,arg1,arg2] = penalty_objective_function(x0,fcn,penalty,
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see .
- [fval,info,exit_flag,arg1,arg2] = fcn(x0,varargin{:});
+[fval, info, exit_flag, arg1, arg2] = fcn(x, varargin{:});
- if info(1) ~= 0
- fval = penalty + info(4);
- end
+if info(1)
+ fval = base_penalty + info(4);
end
\ No newline at end of file
diff --git a/matlab/resol.m b/matlab/resol.m
index 2ad99cf2c..31ce81ede 100644
--- a/matlab/resol.m
+++ b/matlab/resol.m
@@ -109,9 +109,10 @@ if info(1)
end
if options.loglinear
+ threshold = 1e-16;
% Find variables with non positive steady state. Skip auxiliary
% variables for lagges/leaded exogenous variables
- idx = find(dr.ys(get_all_variables_but_lagged_leaded_exogenous(M)) < 1e-9);
+ idx = find(dr.ys(get_all_variables_but_lagged_leaded_exogenous(M)).
+
+v = '4.2.0';
\ No newline at end of file
diff --git a/mex/build/bytecode.am b/mex/build/bytecode.am
index 11936ea5b..3fff91434 100644
--- a/mex/build/bytecode.am
+++ b/mex/build/bytecode.am
@@ -1,6 +1,6 @@
mex_PROGRAMS = bytecode
-bytecode_CPPFLAGS = $(AM_CPPFLAGS) -I$(top_srcdir)/../../sources -I$(top_srcdir)/../../sources/bytecode -I$(top_srcdir)/../../../preprocessor
+bytecode_CPPFLAGS = -Wno-maybe-uninitialized $(AM_CPPFLAGS) -I$(top_srcdir)/../../sources -I$(top_srcdir)/../../sources/bytecode -I$(top_srcdir)/../../../preprocessor
TOPDIR = $(top_srcdir)/../../sources/bytecode
diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc
index 8240e9fc4..bc9aa02a1 100644
--- a/preprocessor/ComputingTasks.cc
+++ b/preprocessor/ComputingTasks.cc
@@ -470,7 +470,7 @@ RamseyConstraintsStatement::writeOutput(ostream &output, const string &basename,
break;
default:
cerr << "Ramsey constraints: this shouldn't happen." << endl;
- exit(1);
+ exit(EXIT_FAILURE);
}
output << "', '";
it->expression->writeOutput(output);
@@ -1094,9 +1094,7 @@ EstimatedParamsStatement::writeOutput(ostream &output, const string &basename, b
<< "estim_params_.corrn = [];" << endl
<< "estim_params_.param_vals = [];" << endl;
- vector::const_iterator it;
-
- for (it = estim_params_list.begin(); it != estim_params_list.end(); it++)
+ for (vector::const_iterator it = estim_params_list.begin(); it != estim_params_list.end(); it++)
{
int symb_id = symbol_table.getTypeSpecificID(it->name) + 1;
SymbolType symb_type = symbol_table.getType(it->name);
@@ -1214,9 +1212,7 @@ EstimatedParamsInitStatement::writeOutput(ostream &output, const string &basenam
if (use_calibration)
output << "options_.use_calibration_initialization = 1;" << endl;
- vector::const_iterator it;
-
- for (it = estim_params_list.begin(); it != estim_params_list.end(); it++)
+ for (vector::const_iterator it = estim_params_list.begin(); it != estim_params_list.end(); it++)
{
int symb_id = symbol_table.getTypeSpecificID(it->name) + 1;
SymbolType symb_type = symbol_table.getType(it->name);
@@ -1313,9 +1309,7 @@ EstimatedParamsBoundsStatement::EstimatedParamsBoundsStatement(const vector::const_iterator it;
-
- for (it = estim_params_list.begin(); it != estim_params_list.end(); it++)
+ for (vector::const_iterator it = estim_params_list.begin(); it != estim_params_list.end(); it++)
{
int symb_id = symbol_table.getTypeSpecificID(it->name) + 1;
SymbolType symb_type = symbol_table.getType(it->name);
@@ -1436,10 +1430,7 @@ void
ObservationTrendsStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
{
output << "options_.trend_coeff = {};" << endl;
-
- trend_elements_t::const_iterator it;
-
- for (it = trend_elements.begin(); it != trend_elements.end(); it++)
+ for (trend_elements_t::const_iterator it = trend_elements.begin(); it != trend_elements.end(); it++)
{
SymbolType type = symbol_table.getType(it->first);
if (type == eEndogenous)
@@ -1450,7 +1441,7 @@ ObservationTrendsStatement::writeOutput(ostream &output, const string &basename,
output << "';" << endl;
}
else
- cout << "Error : Non-variable symbol used in TREND_COEFF: " << it->first << endl;
+ cerr << "Warning : Non-variable symbol used in observation_trends: " << it->first << endl;
}
}
diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc
index a665a5779..213ead141 100644
--- a/preprocessor/DynamicModel.cc
+++ b/preprocessor/DynamicModel.cc
@@ -800,7 +800,7 @@ DynamicModel::writeModelEquationsCode(string &file_name, const string &bin_basen
code_file.open(main_name.c_str(), ios::out | ios::binary | ios::ate);
if (!code_file.is_open())
{
- cout << "Error : Can't open file \"" << main_name << "\" for writing\n";
+ cerr << "Error : Can't open file \"" << main_name << "\" for writing" << endl;
exit(EXIT_FAILURE);
}
@@ -1076,7 +1076,7 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
code_file.open(main_name.c_str(), ios::out | ios::binary | ios::ate);
if (!code_file.is_open())
{
- cout << "Error : Can't open file \"" << main_name << "\" for writing\n";
+ cerr << "Error : Can't open file \"" << main_name << "\" for writing" << endl;
exit(EXIT_FAILURE);
}
//Temporary variables declaration
@@ -1768,7 +1768,7 @@ DynamicModel::Write_Inf_To_Bin_File_Block(const string &dynamic_basename, const
SaveCode.open((bin_basename + "_dynamic.bin").c_str(), ios::out | ios::binary);
if (!SaveCode.is_open())
{
- cout << "Error : Can't open file \"" << bin_basename << "_dynamic.bin\" for writing\n";
+ cerr << "Error : Can't open file \"" << bin_basename << "_dynamic.bin\" for writing" << endl;
exit(EXIT_FAILURE);
}
u_count_int = 0;
@@ -3252,8 +3252,8 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
if (!no_tmp_terms)
computeTemporaryTermsOrdered();
int k = 0;
- equation_block = vector(equation_number());
- variable_block_lead_lag = vector< pair< int, pair< int, int> > >(equation_number());
+ equation_block = vector(equations.size());
+ variable_block_lead_lag = vector< pair< int, pair< int, int> > >(equations.size());
for (unsigned int i = 0; i < getNbBlocks(); i++)
{
for (unsigned int j = 0; j < getBlockSize(i); j++)
@@ -4210,10 +4210,10 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con
<< "% from model file (.mod)" << endl << endl
<< model_local_vars_output.str()
<< model_output.str()
- << "rp = zeros(" << equation_number() << ", "
+ << "rp = zeros(" << equations.size() << ", "
<< symbol_table.param_nbr() << ");" << endl
<< jacobian_output.str()
- << "gp = zeros(" << equation_number() << ", " << dynJacobianColsNbr << ", " << symbol_table.param_nbr() << ");" << endl
+ << "gp = zeros(" << equations.size() << ", " << dynJacobianColsNbr << ", " << symbol_table.param_nbr() << ");" << endl
<< hessian_output.str()
<< "if nargout >= 3" << endl
<< "rpp = zeros(" << residuals_params_second_derivatives.size() << ",4);" << endl
@@ -4238,10 +4238,10 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con
<< "ss_param_deriv, ss_param_2nd_deriv)" << endl
<< model_local_vars_output.str()
<< model_output.str()
- << "rp = zeros(" << equation_number() << ", "
+ << "rp = zeros(" << equations.size() << ", "
<< symbol_table.param_nbr() << ");" << endl
<< jacobian_output.str()
- << "gp = zeros(" << equation_number() << ", " << dynJacobianColsNbr << ", " << symbol_table.param_nbr() << ");" << endl
+ << "gp = zeros(" << equations.size() << ", " << dynJacobianColsNbr << ", " << symbol_table.param_nbr() << ");" << endl
<< hessian_output.str()
<< "rpp = zeros(" << residuals_params_second_derivatives.size() << ",4);" << endl
<< hessian1_output.str()
@@ -5012,7 +5012,7 @@ DynamicModel::writeFirstDerivativesC_csr(const string &basename, bool cuda) cons
break;
default:
std::cerr << "This case shouldn't happen" << std::endl;
- exit(1);
+ exit(EXIT_FAILURE);
}
derivative deriv(col_id + eq*cols_nbr,col_id,eq,it->second);
D.push_back(deriv);
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index bbae98758..7f6d28834 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -113,6 +113,7 @@ class ParsingDriver;
%token CPF_WEIGHTS AMISANOTRISTANI MURRAYJONESPARSLOW
%token FILTER_ALGORITHM PROPOSAL_APPROXIMATION CUBATURE UNSCENTED MONTECARLO DISTRIBUTION_APPROXIMATION
%token NAME
+%token USE_PENALIZED_OBJECTIVE_FOR_HESSIAN
%token NAN_CONSTANT NO_STATIC NOBS NOCONSTANT NODISPLAY NOCORR NODIAGNOSTIC NOFUNCTIONS NO_HOMOTOPY
%token NOGRAPH POSTERIOR_NOGRAPH POSTERIOR_GRAPH NOMOMENTS NOPRINT NORMAL_PDF SAVE_DRAWS
%token OBSERVATION_TRENDS OPTIM OPTIM_WEIGHTS ORDER OSR OSR_PARAMS MAX_DIM_COVA_GROUP ADVANCED OUTFILE OUTVARS OVERWRITE
@@ -1825,6 +1826,7 @@ estimation_options : o_datafile
| o_posterior_sampling_method
| o_posterior_sampler_options
| o_keep_kalman_algo_if_singularity_is_detected
+ | o_use_penalized_objective_for_hessian
;
list_optim_option : QUOTED_STRING COMMA QUOTED_STRING
@@ -3232,6 +3234,7 @@ o_mcmc_jumping_covariance : MCMC_JUMPING_COVARIANCE EQUAL HESSIAN
| MCMC_JUMPING_COVARIANCE EQUAL filename
{ driver.option_str("MCMC_jumping_covariance", $3); }
;
+o_use_penalized_objective_for_hessian : USE_PENALIZED_OBJECTIVE_FOR_HESSIAN { driver.option_num("hessian.use_penalized_objective","true"); };
o_irf_plot_threshold : IRF_PLOT_THRESHOLD EQUAL non_negative_number { driver.option_num("impulse_responses.plot_threshold", $3); };
o_dr_display_tol : DR_DISPLAY_TOL EQUAL non_negative_number { driver.option_num("dr_display_tol", $3); };
o_consider_all_endogenous : CONSIDER_ALL_ENDOGENOUS { driver.option_str("endo_vars_for_moment_computations_in_estimation", "all_endogenous_variables"); };
diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll
index 8d00b2874..e1e0aa4c5 100644
--- a/preprocessor/DynareFlex.ll
+++ b/preprocessor/DynareFlex.ll
@@ -400,6 +400,7 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2
distribution_approximation {return token::DISTRIBUTION_APPROXIMATION;}
proposal_distribution {return token::PROPOSAL_DISTRIBUTION;}
no_posterior_kernel_density {return token::NO_POSTERIOR_KERNEL_DENSITY;}
+use_penalized_objective_for_hessian {return token::USE_PENALIZED_OBJECTIVE_FOR_HESSIAN;}
alpha {
yylval->string_val = new string(yytext);
diff --git a/preprocessor/ModelTree.cc b/preprocessor/ModelTree.cc
index 02b4c93d8..d1a270208 100644
--- a/preprocessor/ModelTree.cc
+++ b/preprocessor/ModelTree.cc
@@ -36,7 +36,7 @@ using namespace MFS;
bool
ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, bool verbose)
{
- const int n = equation_number();
+ const int n = equations.size();
assert(n == symbol_table.endo_nbr());
@@ -96,7 +96,7 @@ ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, boo
#endif
// Create the resulting map, by copying the n first elements of mate_map, and substracting n to them
- endo2eq.resize(equation_number());
+ endo2eq.resize(equations.size());
transform(mate_map.begin(), mate_map.begin() + n, endo2eq.begin(), bind2nd(minus(), n));
#ifdef DEBUG
@@ -143,7 +143,7 @@ ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian
cout << "Normalizing the model..." << endl;
- int n = equation_number();
+ int n = equations.size();
// compute the maximum value of each row of the contemporaneous Jacobian matrix
//jacob_map normalized_contemporaneous_jacobian;
@@ -329,7 +329,7 @@ ModelTree::writeRevXrefs(ostream &output, const map > &xrefmap, co
void
ModelTree::computeNormalizedEquations(multimap &endo2eqs) const
{
- for (int i = 0; i < equation_number(); i++)
+ for (int i = 0; i < equations.size(); i++)
{
VariableNode *lhs = dynamic_cast(equations[i]->get_arg1());
if (lhs == NULL)
@@ -417,11 +417,11 @@ ModelTree::evaluateAndReduceJacobian(const eval_context_t &eval_context, jacob_m
void
ModelTree::computePrologueAndEpilogue(const jacob_map_t &static_jacobian_arg, vector &equation_reordered, vector &variable_reordered)
{
- vector eq2endo(equation_number(), 0);
- equation_reordered.resize(equation_number());
- variable_reordered.resize(equation_number());
+ vector eq2endo(equations.size(), 0);
+ equation_reordered.resize(equations.size());
+ variable_reordered.resize(equations.size());
bool *IM;
- int n = equation_number();
+ int n = equations.size();
IM = (bool *) calloc(n*n, sizeof(bool));
int i = 0;
for (vector::const_iterator it = endo2eq.begin(); it != endo2eq.end(); it++, i++)
@@ -1696,7 +1696,7 @@ ModelTree::Write_Inf_To_Bin_File(const string &basename,
SaveCode.open(bin_basename.c_str(), ios::out | ios::binary);
if (!SaveCode.is_open())
{
- cout << "Error : Can't open file \"" << bin_basename << "\" for writing\n";
+ cerr << "Error : Can't open file \"" << bin_basename << "\" for writing" << endl;
exit(EXIT_FAILURE);
}
u_count_int = 0;
@@ -1799,7 +1799,7 @@ ModelTree::addEquation(expr_t eq, int lineno)
void
ModelTree::addEquation(expr_t eq, int lineno, vector > &eq_tags)
{
- int n = equation_number();
+ int n = equations.size();
for (size_t i = 0; i < eq_tags.size(); i++)
equation_tags.push_back(make_pair(n, eq_tags[i]));
addEquation(eq, lineno);
@@ -1843,7 +1843,7 @@ ModelTree::addNonstationaryVariables(vector nonstationary_vars, bool log_de
void
ModelTree::initializeVariablesAndEquations()
{
- for (int j = 0; j < equation_number(); j++)
+ for (int j = 0; j < equations.size(); j++)
{
equation_reordered.push_back(j);
variable_reordered.push_back(j);
diff --git a/preprocessor/StaticModel.cc b/preprocessor/StaticModel.cc
index e38231675..35c9f516e 100644
--- a/preprocessor/StaticModel.cc
+++ b/preprocessor/StaticModel.cc
@@ -412,7 +412,7 @@ StaticModel::writeModelEquationsCode(const string file_name, const string bin_ba
code_file.open(main_name.c_str(), ios::out | ios::binary | ios::ate);
if (!code_file.is_open())
{
- cout << "Error : Can't open file \"" << main_name << "\" for writing\n";
+ cerr << "Error : Can't open file \"" << main_name << "\" for writing" << endl;
exit(EXIT_FAILURE);
}
int count_u;
@@ -596,7 +596,7 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
code_file.open(main_name.c_str(), ios::out | ios::binary | ios::ate);
if (!code_file.is_open())
{
- cout << "Error : Can't open file \"" << main_name << "\" for writing\n";
+ cerr << "Error : Can't open file \"" << main_name << "\" for writing" << endl;
exit(EXIT_FAILURE);
}
//Temporary variables declaration
@@ -990,7 +990,7 @@ StaticModel::Write_Inf_To_Bin_File_Block(const string &static_basename, const st
SaveCode.open((bin_basename + "_static.bin").c_str(), ios::out | ios::binary);
if (!SaveCode.is_open())
{
- cout << "Error : Can't open file \"" << bin_basename << "_static.bin\" for writing\n";
+ cerr << "Error : Can't open file \"" << bin_basename << "_static.bin\" for writing" << endl;
exit(EXIT_FAILURE);
}
u_count_int = 0;
@@ -2368,10 +2368,10 @@ StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) cons
<< "% from model file (.mod)" << endl << endl
<< model_local_vars_output.str()
<< model_output.str()
- << "rp = zeros(" << equation_number() << ", "
+ << "rp = zeros(" << equations.size() << ", "
<< symbol_table.param_nbr() << ");" << endl
<< jacobian_output.str()
- << "gp = zeros(" << equation_number() << ", " << symbol_table.endo_nbr() << ", "
+ << "gp = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ", "
<< symbol_table.param_nbr() << ");" << endl
<< hessian_output.str()
<< "if nargout >= 3" << endl
@@ -2396,10 +2396,10 @@ StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) cons
<< "function params_derivs(y, x, params)" << endl
<< model_local_vars_output.str()
<< model_output.str()
- << "rp = zeros(" << equation_number() << ", "
+ << "rp = zeros(" << equations.size() << ", "
<< symbol_table.param_nbr() << ");" << endl
<< jacobian_output.str()
- << "gp = zeros(" << equation_number() << ", " << symbol_table.endo_nbr() << ", "
+ << "gp = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ", "
<< symbol_table.param_nbr() << ");" << endl
<< hessian_output.str()
<< "rpp = zeros(" << residuals_params_second_derivatives.size() << ",4);" << endl
diff --git a/preprocessor/StaticModel.hh b/preprocessor/StaticModel.hh
index 75ab0b3f8..d770896ed 100644
--- a/preprocessor/StaticModel.hh
+++ b/preprocessor/StaticModel.hh
@@ -30,9 +30,6 @@ using namespace std;
class StaticModel : public ModelTree
{
private:
- //! Temporary terms for the file containing parameters dervicatives
- temporary_terms_t params_derivs_temporary_terms;
-
//! global temporary terms for block decomposed models
vector > v_temporary_terms;
diff --git a/tests/optimizers/fs2000.common.inc b/tests/optimizers/fs2000.common.inc
index 51c871bca..ee8510a3b 100644
--- a/tests/optimizers/fs2000.common.inc
+++ b/tests/optimizers/fs2000.common.inc
@@ -1,7 +1,7 @@
/*
- * This file replicates the estimation of the cash in advance model described
- * Frank Schorfheide (2000): "Loss function-based evaluation of DSGE models",
- * Journal of Applied Econometrics, 15(6), 645-670.
+ * This file replicates (except for a small difference in the priors)the estimation of
+ * the cash in advance model described Frank Schorfheide (2000): "Loss function-based
+ * evaluation of DSGE models", Journal of Applied Econometrics, 15(6), 645-670.
*
* The data are in file "fs2000/fsdat_simul.m", and have been artificially generated.
* They are therefore different from the original dataset used by Schorfheide.
@@ -17,7 +17,7 @@
*/
/*
- * Copyright (C) 2004-2010 Dynare Team
+ * Copyright (C) 2004-2017 Dynare Team
*
* This file is part of Dynare.
*
@@ -106,7 +106,16 @@ alp, beta_pdf, 0.356, 0.02;
bet, beta_pdf, 0.993, 0.002;
gam, normal_pdf, 0.0085, 0.003;
mst, normal_pdf, 1.0002, 0.007;
-rho, beta_pdf, 0.129, 0.223;
+/*
+ * The prior on the autoregressive parameter rho is different from the one considered in the original paper.
+ * We lowered the prior variance to ensure that the prior does not have an asymptote at zero, which causes
+ * troubles for some optimizers. For instance, mode_compute=12 (which relies on Mathworks' particleswarm
+ * routine) finds a posterior at zero for rho (with a much higher posterior density than the one reported
+ * by csminwell or other optimization routines). Because zero is on the boundary of the prior, this is an
+ * an issue when we have to compute the hessian matrix (which is not full rank) at the mode or if we try to
+ * use this solution as an initial guess for csminwell.
+ */
+rho, beta_pdf, 0.129, 0.100;
psi, beta_pdf, 0.65, 0.05;
del, beta_pdf, 0.01, 0.005;
stderr e_a, inv_gamma_pdf, 0.035449, inf;
diff --git a/tests/optimizers/fs2000_12.mod b/tests/optimizers/fs2000_12.mod
new file mode 100644
index 000000000..306d2aa80
--- /dev/null
+++ b/tests/optimizers/fs2000_12.mod
@@ -0,0 +1,3 @@
+@#include "fs2000.common.inc"
+
+estimation(mode_compute=12,order=1, datafile='../fs2000/fsdat_simul', nobs=192, mh_replic=0);