Merge branch 'master' into json

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

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] ...

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);

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';

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
+                % with flag==0, we force to use the hessian from outer product gradient of optimizer 5
-                % product gradient of optimizer 5
+                if options_.hessian.use_penalized_objective
-                hh = reshape(penalty_hessian(objective_function,xparam1,fval, ...
+                    penalized_objective_function = str2func('penalty_objective_function');
-                    options_.gstep,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_),nx,nx);
+                    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
                 %

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 = [];

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);
+n   = size(x,1);
-h_1=h1;
+h1  = max(abs(x), sqrt(gstep(1))*ones(n, 1))*eps^(1/6)*gstep(2);
-xh1=x+h1;
+h_1 = h1;
-h1=xh1-x;
+xh1 = x+h1;
-xh1=x-h_1;
+h1  = xh1-x;
-h_1=x-xh1;
+xh1 = x-h_1;
-xh1=x;
+h_1 = x-xh1;
-f0=feval(func,x,varargin{:});
+xh1 = x;
-f1=zeros(size(f0,1),n);
+f0  = feval(func, x, varargin{:});
-f_1=f1;
+f1  = zeros(size(f0, 1), n);
+f_1 = f1;
+
 for i=1:n
     %do step up
-    xh1(i)=x(i)+h1(i);
+    xh1(i)   = x(i)+h1(i);
-    f1(:,i)=feval(func,xh1,varargin{:});
+    f1(:,i)  = feval(func, xh1, varargin{:});
-    %do step up
+    %do step down
-    xh1(i)=x(i)-h_1(i);
+    xh1(i)   = x(i)-h_1(i);
-    f_1(:,i)=feval(func,xh1,varargin{:});
+    f_1(:,i) = feval(func, xh1, varargin{:});
-    xh1(i)=x(i);%reset parameter
+    %reset parameter
+    xh1(i)   = x(i); 
 end
-xh_1=xh1;
+
-hessian_mat = zeros(size(f0,1),n*n);
+xh_1 = xh1;
-temp=f1+f_1-f0*ones(1,n); %term f_(1,0)+f_(-1,0)-f_(0,0) used later
+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      
+    if i > 1  
-        k=[i:n:n*(i-1)];
+        %fill symmetric part of Hessian based on previously computed results      
-        hessian_mat(:,(i-1)*n+1:(i-1)*n+i-1)=hessian_mat(:,k);
+        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(i) = x(i)+h1(i);
-        xh1(j)=x(j)+h_1(j);
+        xh1(j) = x(j)+h_1(j);
         %step in down direction
-        xh_1(i)=x(i)-h1(i);
+        xh_1(i) = x(i)-h1(i);
-        xh_1(j)=x(j)-h_1(j);
+        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
+        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(i)  = x(i);
-        xh1(j)=x(j);
+        xh1(j)  = x(j);
-        xh_1(i)=x(i);
+        xh_1(i) = x(i);
-        xh_1(j)=x(j);
+        xh_1(j) = x(j);
     end    
 end
 

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 <http://www.gnu.org/licenses/>.
 
-% 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!')

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 ;
+    options_.cova_compute = 0;
-     [opt_par_values,stdh,lb_95,ub_95,med_param] = online_auxiliary_filter(start_par_value,varargin{:}) ;
+    [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)

matlab/optimization/penalty_hessian.m

@@ -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 <http://www.gnu.org/licenses/>.
-
-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
-

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(x, fcn, base_penalty, varargin)
-%function [fval,exit_flag,arg1,arg2] = penalty_objective_function(x0,fcn,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 <http://www.gnu.org/licenses/>.
 
-    [fval,info,exit_flag,arg1,arg2] = fcn(x0,varargin{:});
+[fval, info, exit_flag, arg1, arg2] = fcn(x, varargin{:});
 
-    if info(1) ~= 0
+if info(1)
-        fval = penalty + info(4);
+    fval = base_penalty + info(4);
-    end
 end

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))<threshold);
     if length(idx)
         if options.debug
             variables_with_non_positive_steady_state = M.endo_names(idx,:);
@@ -132,7 +133,7 @@ if options.loglinear
             end
         end
         info(1)=26;
-        info(2)=sum(dr.ys(dr.ys<1e-9).^2);
+        info(2)=sum(dr.ys(dr.ys<threshold).^2);
         return
     end
 end

matlab/supported_octave_version.m

@@ -0,0 +1,20 @@
+function v = supported_octave_version()
+
+% Copyright (C) 2017 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 <http://www.gnu.org/licenses/>.
+
+v = '4.2.0';

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
 

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<EstimationParams>::const_iterator it;
+  for (vector<EstimationParams>::const_iterator it = estim_params_list.begin(); it != estim_params_list.end(); it++)
-
-  for (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<EstimationParams>::const_iterator it;
+  for (vector<EstimationParams>::const_iterator it = estim_params_list.begin(); it != estim_params_list.end(); it++)
-
-  for (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<Esti
 void
 EstimatedParamsBoundsStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
 {
-  vector<EstimationParams>::const_iterator it;
+  for (vector<EstimationParams>::const_iterator it = estim_params_list.begin(); it != estim_params_list.end(); it++)
-
-  for (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;
-
+  for (trend_elements_t::const_iterator it = trend_elements.begin(); it != trend_elements.end(); it++)
-  trend_elements_t::const_iterator it;
-
-  for (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;
     }
 }
 

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<int>(equation_number());
+      equation_block = vector<int>(equations.size());
-      variable_block_lead_lag = vector< pair< int, pair< int, int> > >(equation_number());
+      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);

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 <string_val> 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"); };

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
 <DYNARE_STATEMENT>distribution_approximation {return token::DISTRIBUTION_APPROXIMATION;}
 <DYNARE_STATEMENT>proposal_distribution {return token::PROPOSAL_DISTRIBUTION;}
 <DYNARE_STATEMENT>no_posterior_kernel_density {return token::NO_POSTERIOR_KERNEL_DENSITY;}
+<DYNARE_STATEMENT>use_penalized_objective_for_hessian {return token::USE_PENALIZED_OBJECTIVE_FOR_HESSIAN;}
 
 <DYNARE_STATEMENT>alpha {
   yylval->string_val = new string(yytext);

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<int>(), 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<int, set<int> > &xrefmap, co
 void
 ModelTree::computeNormalizedEquations(multimap<int, int> &endo2eqs) const
 {
-  for (int i = 0; i < equation_number(); i++)
+  for (int i = 0; i < equations.size(); i++)
     {
       VariableNode *lhs = dynamic_cast<VariableNode *>(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<int> &equation_reordered, vector<int> &variable_reordered)
 {
-  vector<int> eq2endo(equation_number(), 0);
+  vector<int> eq2endo(equations.size(), 0);
-  equation_reordered.resize(equation_number());
+  equation_reordered.resize(equations.size());
-  variable_reordered.resize(equation_number());
+  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<int>::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<pair<string, string> > &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<int> 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);

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

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<vector<temporary_terms_t> > v_temporary_terms;
 

tests/optimizers/fs2000.common.inc

@@ -1,7 +1,7 @@
 /*
- * This file replicates the estimation of the cash in advance model described
+ * This file replicates (except for a small difference in the priors)the estimation of
- * Frank Schorfheide (2000): "Loss function-based evaluation of DSGE models",
+ * the cash in advance model described Frank Schorfheide (2000): "Loss function-based
- * Journal of Applied Econometrics, 15(6), 645-670.
+ * 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;

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);