hess_element

 function d=hess_element(func,element1,element2,args)
 returns an entry of the finite differences approximation to the hessian of func

 INPUTS
    func       [function name]    string with name of the function
    element1   [int]              the indices showing the element within the hessian that should be returned
    element2   [int]
    args       [cell array]       arguments provided to func

 OUTPUTS
    d          [double]           the (element1,element2) entry of the hessian

 SPECIAL REQUIREMENTS
    none

0001 function d=hess_element(func,element1,element2,args)
0002 % function d=hess_element(func,element1,element2,args)
0003 % returns an entry of the finite differences approximation to the hessian of func
0004 %
0005 % INPUTS
0006 %    func       [function name]    string with name of the function
0007 %    element1   [int]              the indices showing the element within the hessian that should be returned
0008 %    element2   [int]
0009 %    args       [cell array]       arguments provided to func
0010 %
0011 % OUTPUTS
0012 %    d          [double]           the (element1,element2) entry of the hessian
0013 %
0014 % SPECIAL REQUIREMENTS
0015 %    none
0016 
0017 % Copyright (C) 2010-2011 Dynare Team
0018 %
0019 % This file is part of Dynare.
0020 %
0021 % Dynare is free software: you can redistribute it and/or modify
0022 % it under the terms of the GNU General Public License as published by
0023 % the Free Software Foundation, either version 3 of the License, or
0024 % (at your option) any later version.
0025 %
0026 % Dynare is distributed in the hope that it will be useful,
0027 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0028 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0029 % GNU General Public License for more details.
0030 %
0031 % You should have received a copy of the GNU General Public License
0032 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0033 
0034 func = str2func(func);
0035 
0036 h=10e-6;
0037 p10 = args;
0038 p01 = args;
0039 m10 = args;
0040 m01 = args;
0041 p11 = args;
0042 m11 = args;
0043 for i=1:size(args,2)
0044     if i==element1
0045         p10{i} = p10{i} + h;
0046         m10{i} = m10{i} - h;
0047 
0048         p11{i} = p11{i} + h;
0049         m11{i} = m11{i} - h;
0050     end
0051 
0052     if i==element2
0053         p01{i} = p01{i} + h;
0054         m01{i} = m01{i} - h;
0055 
0056         p11{i} = p11{i} + h;
0057         m11{i} = m11{i} - h;
0058     end
0059 end
0060 
0061 % From Abramowitz and Stegun. Handbook of Mathematical Functions (1965)
0062 % formulas 25.3.24 and 25.3.27 p. 884
0063 if element1==element2
0064     d = (16*func(p10{:})...
0065          +16*func(m10{:})...
0066          -30*func(args{:})...
0067          -func(p11{:})...
0068          -func(m11{:}))/(12*h^2);
0069 else
0070     d = (func(p10{:})...
0071          +func(m10{:})...
0072          +func(p01{:})...
0073          +func(m01{:})...
0074          -2*func(args{:})...
0075          -func(p11{:})...
0076          -func(m11{:}))/(-2*h^2);
0077 end