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