Add code that conducts the posterior mode finding if desired

PF_caller.m

@@ -43,10 +43,10 @@
 x_std=NaN(T,1);
 Eff_particle=NaN(T,1);
 
-rho_sigma=draw(1,1);
-rho=draw(1,2);
-eta_sigma=draw(1,3);
-sigma_bar=draw(1,4);
+rho_sigma=draw(1);
+rho=draw(2);
+eta_sigma=draw(3);
+sigma_bar=draw(4);
 
 [rho_sigma,rho,eta_sigma,sigma_bar]=par_retransform_AR1(rho_sigma,rho,eta_sigma,sigma_bar);
 

PF_caller_minimizer.m

@@ -0,0 +1,29 @@
+function [minus_logpost,x_hat,x_std,Eff_particle]=PF_caller_minimizer(varargin)
+% [minus_logpost, flag, xhat_PF0,x_std_PF0,Eff_particle]=PF_caller_minimizer(varargin)
+%Inputs:
+%   varargin            cell array        all required input arguments
+%
+% Output
+%  logpost              scalar            The log-posterior
+%  x_hat                [T by 1] vector   the posterior state estimate
+%                                           x_t|T
+%  x_std                [T by 1] vector   standard deviation of the posterior
+%  Eff_particle     	[T by 1] vector   A measure of the effective sample size
+%
+% Copyright (C) 2013-2015 Benjamin Born + Johannes Pfeifer
+%
+% This 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.
+%
+% It 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.
+%
+% See <http://www.gnu.org/licenses/>.
+
+
+[logpost,x_hat,x_std,Eff_particle]=PF_caller(varargin{:});
+minus_logpost=-logpost;

run_filter_and_smoother_AR1.m

@@ -5,7 +5,10 @@
 % positive definite matrix will still yield draws from the ergodic distribution 
 % (see Chib/Greenberg (1995))
 % 
-% Copyright (C) 2013-2014 Benjamin Born + Johannes Pfeifer
+% To perform mode-finding using the CMA-ES algorithm, uncomment the code in lines 101 following and 
+% download the cmaes.m from the provided link.
+% 
+% Copyright (C) 2013-2015 Benjamin Born + Johannes Pfeifer
 %
 % This is free software: you can redistribute it and/or modify
 % it under the terms of the GNU General Public License as published by
@@ -72,7 +75,6 @@
 PFstreamrand = RandStream('mt19937ar','Seed',2);
 randnr             = rand(PFstreamrand,periods,num_sim_filter);
 
-
 %initialize matrices
 draws=zeros(MH_draws+burnin,npar);
 likelihood=zeros(MH_draws+burnin,1);
@@ -88,17 +90,33 @@
 scale_mh = 2; %scale factor for proposal density to get acceptance rate of 23-40 percent
 accept=0;
 
-%set startig value, starts a true values and makes mode-finding redundant
+%set startig value, starts at true values and makes mode-finding redundant
 draws(1,:)=[rho_sigma_trans,rho_trans,eta_sigma_trans,sigma_bar_trans];
-old_posterior = PF_caller(draws(1,:),observable_series,num_sim_filter,num_sim_smoother,shocks,randnr,1); %!!! reset 1 to 0
-likelihood(1,1)=old_posterior;
-
 
 % set jumping matrix for proposal density
 % uses simple identity matrix instead of Hessian at the mode
 inv_Hessian = 1e-3*eye(4); 
 Sigma_chol = cholcov(inv_Hessian)';
 
+%% Maximize posterior using the CMA-ES algorithm of Hansen et al. (for example included in Dynare and available at https://www.lri.fr/~hansen/cmaes_inmatlab.html)
+% vm=0.1*ones(npar,1);
+% opts=struct('MaxFunEvals',1000,'SaveVariables','off','DispFinal','off','WarnOnEqualFunctionValues','no','DispModulo','100','LogModulo','0','LogTime','0');
+% [x, FMIN, COUNTEVAL, STOPFLAG, OUT, BESTEVER] = cmaes('PF_caller_minimizer',draws(1,:)',vm,opts,observable_series,num_sim_filter,dimy,shocks,randnr,0);
+% draws(1,:)=BESTEVER.x;
+
+%% Compute the Hessian at the mode, using Dynare's hessian.m; Note that this is not advocated as the numerical differentiation will usually not work 
+%% due to the non-differentiability of the likelihood function
+% gstep = ones(2,1);
+% gstep(1) = 1e-2;
+% gstep(2) = 1.0;
+% hessian_mat = reshape(hessian('PF_caller',xopt,gstep,observable_series,num_sim_filter,dimy,shocks,randnr,0),length(xopt),length(xopt));
+% Sigma_chol = cholcov(inv(hessian_mat))';
+
+
+old_posterior = PF_caller(draws(1,:),observable_series,num_sim_filter,num_sim_smoother,shocks,randnr,0);
+likelihood(1,1)=old_posterior;
+
+
 
 proposal_draws = scale_mh*Sigma_chol*randn(4,MH_draws+burnin);
 for ii=2:burnin+MH_draws