Added new optimization routine (particleswarm).

Only available under Matlab if the Global Optimization Toolbox is installed.

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

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/global_initialization.m

@@ -645,6 +645,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/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 ;
-     [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)

tests/Makefile.am

@@ -309,6 +309,7 @@ MODFILES = \
 	optimizers/fs2000_4_with_optim.mod \
 	optimizers/fs2000_5.mod \
 	optimizers/fs2000_7.mod \
+	optimizers/fs2000_12.mod \
 	optimizers/fs2000_101.mod \
 	optimizers/fs2000_102.mod \
 	optimizers/fs2000_w.mod \

tests/optimizers/fs2000_12.mod

+@#include "fs2000.common.inc"
+
+estimation(mode_compute=12,order=1, datafile='../fs2000/fsdat_simul', nobs=192, mh_replic=0);