diff --git a/doc/dynare.texi b/doc/dynare.texi
index 1ca3bc42e460743d3771584944d96faa7291d009..bf63caf7097bb3feae9a49ab4f827c099ead15fd 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
diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m
index 8a4a7d6b61bf6b8e09dc72bf4ac93a87746e58db..d911adc6e0031556eb62fe1fdaeee6c73a8e1c9a 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 646120c5f6b44ec6be72f4d1ebab038d1641616a..b97beef262e40ad6ccd799c7572ca06fc82beb69 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/global_initialization.m b/matlab/global_initialization.m
index e65ca4a62a53ff320a98718a10f0a4dca76da1ae..aa1531bb42b39899331548812d1647665afb7c51 100644
--- a/matlab/global_initialization.m
+++ b/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 = [];
diff --git a/matlab/non_linear_dsge_likelihood.m b/matlab/non_linear_dsge_likelihood.m
index f8ac92e230f6602c133ce7c7fcdd045c3b153ed6..17d448927ad82bd7ee77def16b93904239b1624c 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 <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!')
diff --git a/matlab/optimization/dynare_minimize_objective.m b/matlab/optimization/dynare_minimize_objective.m
index 4724fc91a47ffe81de22e6c23415e8b29f2dd278..0aebc6b50621e19577bd7879fb6af594f1c298d6 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/tests/Makefile.am b/tests/Makefile.am
index c2b935217b4ff535bdf56430572df9e888be4919..a85b41085eaef28c9acc65910819c34657f9f068 100644
--- a/tests/Makefile.am
+++ b/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 \
diff --git a/tests/optimizers/fs2000_12.mod b/tests/optimizers/fs2000_12.mod
new file mode 100644
index 0000000000000000000000000000000000000000..306d2aa8011ee3f01c92ab40a0886bae54a23b75
--- /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);