diff --git a/matlab/csminwel1.m b/matlab/csminwel1.m
index 6b0ea0166..403b9a58e 100644
--- a/matlab/csminwel1.m
+++ b/matlab/csminwel1.m
@@ -81,6 +81,10 @@ if NumGrad
[g,badg] = numgrad3(fcn, f0, x0, epsilon, varargin{:});
case 5
[g,badg] = numgrad5(fcn, f0, x0, epsilon, varargin{:});
+ case 13
+ [g,badg] = numgrad3_(fcn, f0, x0, epsilon, varargin{:});
+ case 15
+ [g,badg] = numgrad5_(fcn, f0, x0, epsilon, varargin{:});
otherwise
error('csminwel1: Unknown method for gradient evaluation!')
end
@@ -135,6 +139,10 @@ while ~done
[g1 badg1] = numgrad3(fcn, f1, x1, epsilon, varargin{:});
case 5
[g1,badg1] = numgrad5(fcn, f1, x1, epsilon, varargin{:});
+ case 13
+ [g1,badg1] = numgrad3_(fcn, f1, x1, epsilon, varargin{:});
+ case 15
+ [g1,badg1] = numgrad5_(fcn, f1, x1, epsilon, varargin{:});
otherwise
error('csminwel1: Unknown method for gradient evaluation!')
end
@@ -174,6 +182,10 @@ while ~done
[g2 badg2] = numgrad3(fcn, f2, x2, epsilon, varargin{:});
case 5
[g2,badg2] = numgrad5(fcn, f2, x2, epsilon, varargin{:});
+ case 13
+ [g2,badg2] = numgrad3_(fcn, f2, x2, epsilon, varargin{:});
+ case 15
+ [g2,badg2] = numgrad5_(fcn, f2, x2, epsilon, varargin{:});
otherwise
error('csminwel1: Unknown method for gradient evaluation!')
end
@@ -214,6 +226,10 @@ while ~done
[g3 badg3] = numgrad3(fcn, f3, x3, epsilon, varargin{:});
case 5
[g3,badg3] = numgrad5(fcn, f3, x3, epsilon, varargin{:});
+ case 13
+ [g3,badg3] = numgrad3_(fcn, f3, x3, epsilon, varargin{:});
+ case 15
+ [g3,badg3] = numgrad5_(fcn, f3, x3, epsilon, varargin{:});
otherwise
error('csminwel1: Unknown method for gradient evaluation!')
end
@@ -283,6 +299,10 @@ while ~done
[gh,badgh] = numgrad3(fcn, fh, xh, epsilon, varargin{:});
case 5
[gh,badgh] = numgrad5(fcn, fh, xh, epsilon, varargin{:});
+ case 13
+ [gh,badgh] = numgrad3_(fcn, fh, xh, epsilon, varargin{:});
+ case 15
+ [gh,badgh] = numgrad5_(fcn, fh, xh, epsilon, varargin{:});
otherwise
error('csminwel1: Unknown method for gradient evaluation!')
end
diff --git a/matlab/numgrad3_.m b/matlab/numgrad3_.m
new file mode 100644
index 000000000..a4b4e56a6
--- /dev/null
+++ b/matlab/numgrad3_.m
@@ -0,0 +1,67 @@
+function [g, badg, f0, f1, f2] = numgrad3_(fcn,f0,x,epsilon,varargin)
+% Computes the gradient of the objective function fcn using a three points
+% formula if possible.
+%
+% Adapted from Sims' numgrad routine.
+%
+% See section 25.3.4 in Abramovitz and Stegun (1972, Tenth Printing, December) Handbook of Mathematical Functions.
+% http://www.math.sfu.ca/~cbm/aands/
+
+% Original file downloaded from:
+% http://sims.princeton.edu/yftp/optimize/mfiles/numgrad.m
+
+% Copyright (C) 1993-2007 Christopher Sims
+% Copyright (C) 2008-2012 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 .
+
+delta = epsilon;
+n = length(x);
+g = zeros(n,1);
+
+badg = 0;
+
+scale = []; % one(n,1);
+
+for i=1:n
+ xiold = x(i);
+ h = step_length_correction(xiold,scale,i)*delta;
+ x(i) = xiold + h;
+ [f1,junk1,junk2,cost_flag1] = feval(fcn, x, varargin{:});
+ if ~cost_flag1
+ fprintf('Gradient w.r.t. parameter number %3d (x=%16.8f,+h=%16.8f,f0=%16.8f,f1=%16.8f,f2=%16.8f,g0=%16.8f): penalty on the right!\n',i,xiold,h,f0,f1,f2,(f1 - f2) / (2*h))
+ end
+ x(i) = xiold - h;
+ [f2,junk2,junk2,cost_flag2] = feval(fcn, x, varargin{:});
+ if ~cost_flag2
+ fprintf('Gradient w.r.t. parameter number %3d (x=%16.8f,+h=%16.8f,f0=%16.8f,f1=%16.8f,f2=%16.8f,g0=%16.8f): penalty on the left!\n',i,xiold,h,f0,f1,f2,(f1 - f2) / (2*h))
+ end
+ if f0.
+
+delta = epsilon;
+n = length(x);
+g = zeros(n,1);
+
+badg = 0;
+scale = []; % ones(n,1);
+
+for i=1:n
+ xiold = x(i);
+ h = step_length_correction(xiold,scale,i)*delta;
+ x(i) = xiold+h;
+ [f1,junk1,junk2,cost_flag1] = feval(fcn, x, varargin{:});
+ x(i) = xiold-h;
+ [f2,junk1,junk2,cost_flag2] = feval(fcn, x, varargin{:});
+ x(i) = xiold+2*h;
+ [f3,junk1,junk2,cost_flag3] = feval(fcn, x, varargin{:});
+ x(i) = xiold-2*h;
+ [f4,junk1,junk2,cost_flag4] = feval(fcn, x, varargin{:});
+ if f0