1function hessian_mat = hessian(func,x, gstep, varargin) % --*-- Unitary tests --*--
2
3% Computes second order partial derivatives
4%
5% INPUTS
6%    func        [string]   name of the function
7%    x           [double]   vector, the Hessian of "func" is evaluated at x.
8%    gstep       [double]   scalar, size of epsilon.
9%    varargin    [void]     list of additional arguments for "func".
10%
11% OUTPUTS
12%    hessian_mat [double]   Hessian matrix
13%
14% ALGORITHM
15%    Uses Abramowitz and Stegun (1965) formulas 25.3.23
16% \[
17%     \frac{\partial^2 f_{0,0}}{\partial {x^2}} = \frac{1}{h^2}\left( f_{1,0} - 2f_{0,0} + f_{ - 1,0} \right)
18% \]
19% and 25.3.27 p. 884
20%
21% \[
22%     \frac{\partial ^2f_{0,0}}{\partial x\partial y} = \frac{-1}{2h^2}\left(f_{1,0} + f_{-1,0} + f_{0,1} + f_{0,-1} - 2f_{0,0} - f_{1,1} - f_{-1,-1} \right)
23% \]
24%
25% SPECIAL REQUIREMENTS
26%    none
27%
28
29% Copyright (C) 2001-2017 Dynare Team
30%
31% This file is part of Dynare.
32%
33% Dynare is free software: you can redistribute it and/or modify
34% it under the terms of the GNU General Public License as published by
35% the Free Software Foundation, either version 3 of the License, or
36% (at your option) any later version.
37%
38% Dynare is distributed in the hope that it will be useful,
39% but WITHOUT ANY WARRANTY; without even the implied warranty of
40% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
41% GNU General Public License for more details.
42%
43% You should have received a copy of the GNU General Public License
44% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
45
46if ~isa(func, 'function_handle')
47    func = str2func(func);
48end
49
50n   = size(x,1);
51h1  = max(abs(x), sqrt(gstep(1))*ones(n, 1))*eps^(1/6)*gstep(2);
52h_1 = h1;
53xh1 = x+h1;
54h1  = xh1-x;
55xh1 = x-h_1;
56h_1 = x-xh1;
57xh1 = x;
58f0  = feval(func, x, varargin{:});
59f1  = zeros(size(f0, 1), n);
60f_1 = f1;
61
62for i=1:n
63    %do step up
64    xh1(i)   = x(i)+h1(i);
65    f1(:,i)  = feval(func, xh1, varargin{:});
66    %do step down
67    xh1(i)   = x(i)-h_1(i);
68    f_1(:,i) = feval(func, xh1, varargin{:});
69    %reset parameter
70    xh1(i)   = x(i);
71end
72
73xh_1 = xh1;
74temp = f1+f_1-f0*ones(1, n); %term f_(1,0)+f_(-1,0)-f_(0,0) used later
75
76hessian_mat = zeros(size(f0,1), n*n);
77
78for i=1:n
79    if i > 1
80        %fill symmetric part of Hessian based on previously computed results
81        k = [i:n:n*(i-1)];
82        hessian_mat(:,(i-1)*n+1:(i-1)*n+i-1) = hessian_mat(:,k);
83    end
84    hessian_mat(:,(i-1)*n+i) = (f1(:,i)+f_1(:,i)-2*f0)./(h1(i)*h_1(i)); %formula 25.3.23
85    for j=i+1:n
86        %step in up direction
87        xh1(i) = x(i)+h1(i);
88        xh1(j) = x(j)+h_1(j);
89        %step in down direction
90        xh_1(i) = x(i)-h1(i);
91        xh_1(j) = x(j)-h_1(j);
92        hessian_mat(:,(i-1)*n+j) =-(-feval(func, xh1, varargin{:})-feval(func, xh_1, varargin{:})+temp(:,i)+temp(:,j))./(2*h1(i)*h_1(j)); %formula 25.3.27
93                                                                                                                                          %reset grid points
94        xh1(i)  = x(i);
95        xh1(j)  = x(j);
96        xh_1(i) = x(i);
97        xh_1(j) = x(j);
98    end
99end
100
101
102%@test:1
103%$ % Create a function.
104%$ fid = fopen('exfun.m','w+');
105%$ fprintf(fid,'function [f,g,H] = exfun(xvar)\\n');
106%$ fprintf(fid,'x = xvar(1);\\n');
107%$ fprintf(fid,'y = xvar(2);\\n');
108%$ fprintf(fid,'f = x^2* log(y);\\n');
109%$ fprintf(fid,'if nargout>1\\n');
110%$ fprintf(fid,'    g = zeros(2,1);\\n');
111%$ fprintf(fid,'    g(1) = 2*x*log(y);\\n');
112%$ fprintf(fid,'    g(2) = x*x/y;\\n');
113%$ fprintf(fid,'end\\n');
114%$ fprintf(fid,'if nargout>2\\n');
115%$ fprintf(fid,'    H = zeros(2,2);\\n');
116%$ fprintf(fid,'    H(1,1) = 2*log(y);\\n');
117%$ fprintf(fid,'    H(1,2) = 2*x/y;\\n');
118%$ fprintf(fid,'    H(2,1) = H(1,2);\\n');
119%$ fprintf(fid,'    H(2,2) = -x*x/(y*y);\\n');
120%$ fprintf(fid,'    H = H(:);\\n');
121%$ fprintf(fid,'end\\n');
122%$ fclose(fid);
123%$
124%$ rehash;
125%$
126%$ t = zeros(5,1);
127%$
128%$ % Evaluate the Hessian at (1,e)
129%$ try
130%$    H = hessian('exfun',[1; exp(1)],[1e-2; 1]);
131%$    t(1) = 1;
132%$ catch
133%$    t(1) = 0;
134%$ end
135%$
136%$ % Compute the true Hessian matrix
137%$ [f, g, Htrue] = exfun([1 exp(1)]);
138%$
139%$ % Delete exfun routine from disk.
140%$ delete('exfun.m');
141%$
142%$ % Compare the values in H and Htrue
143%$ if t(1)
144%$    t(2) = dassert(abs(H(1)-Htrue(1))<1e-6,true);
145%$    t(3) = dassert(abs(H(2)-Htrue(2))<1e-6,true);
146%$    t(4) = dassert(abs(H(3)-Htrue(3))<1e-6,true);
147%$    t(5) = dassert(abs(H(4)-Htrue(4))<1e-6,true);
148%$ end
149%$ T = all(t);
150%@eof:1
151