1function d=hess_element(func,element1,element2,args) 2% function d=hess_element(func,element1,element2,args) 3% returns an entry of the finite differences approximation to the hessian of func 4% 5% INPUTS 6% func [function name] string with name of the function 7% element1 [int] the indices showing the element within the hessian that should be returned 8% element2 [int] 9% args [cell array] arguments provided to func 10% 11% OUTPUTS 12% d [double] the (element1,element2) entry of the hessian 13% 14% SPECIAL REQUIREMENTS 15% none 16 17% Copyright (C) 2010-2011 Dynare Team 18% 19% This file is part of Dynare. 20% 21% Dynare is free software: you can redistribute it and/or modify 22% it under the terms of the GNU General Public License as published by 23% the Free Software Foundation, either version 3 of the License, or 24% (at your option) any later version. 25% 26% Dynare is distributed in the hope that it will be useful, 27% but WITHOUT ANY WARRANTY; without even the implied warranty of 28% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 29% GNU General Public License for more details. 30% 31% You should have received a copy of the GNU General Public License 32% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 33 34func = str2func(func); 35 36h=10e-6; 37p10 = args; 38p01 = args; 39m10 = args; 40m01 = args; 41p11 = args; 42m11 = args; 43for i=1:size(args,2) 44 if i==element1 45 p10{i} = p10{i} + h; 46 m10{i} = m10{i} - h; 47 48 p11{i} = p11{i} + h; 49 m11{i} = m11{i} - h; 50 end 51 52 if i==element2 53 p01{i} = p01{i} + h; 54 m01{i} = m01{i} - h; 55 56 p11{i} = p11{i} + h; 57 m11{i} = m11{i} - h; 58 end 59end 60 61% From Abramowitz and Stegun. Handbook of Mathematical Functions (1965) 62% formulas 25.3.24 and 25.3.27 p. 884 63if element1==element2 64 d = (16*func(p10{:})... 65 +16*func(m10{:})... 66 -30*func(args{:})... 67 -func(p11{:})... 68 -func(m11{:}))/(12*h^2); 69else 70 d = (func(p10{:})... 71 +func(m10{:})... 72 +func(p01{:})... 73 +func(m01{:})... 74 -2*func(args{:})... 75 -func(p11{:})... 76 -func(m11{:}))/(-2*h^2); 77end 78