1function hessian_mat = hessian_sparse(func,x,gstep,varargin) 2% function hessian_mat = hessian_sparse(func,x,gstep,varargin) 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, sparse] Hessian matrix 13% 14% ALGORITHM 15% Uses Abramowitz and Stegun (1965) formulas 25.3.24 and 25.3.27 p. 884 16% 17% SPECIAL REQUIREMENTS 18% none 19% 20 21% Copyright (C) 2001-2017 Dynare Team 22% 23% This file is part of Dynare. 24% 25% Dynare is free software: you can redistribute it and/or modify 26% it under the terms of the GNU General Public License as published by 27% the Free Software Foundation, either version 3 of the License, or 28% (at your option) any later version. 29% 30% Dynare is distributed in the hope that it will be useful, 31% but WITHOUT ANY WARRANTY; without even the implied warranty of 32% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 33% GNU General Public License for more details. 34% 35% You should have received a copy of the GNU General Public License 36% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 37 38if ~isa(func, 'function_handle') 39 func = str2func(func); 40end 41n=size(x,1); 42h1=max(abs(x),sqrt(gstep(1))*ones(n,1))*eps^(1/6)*gstep(2); 43h_1=h1; 44xh1=x+h1; 45h1=xh1-x; 46xh1=x-h_1; 47h_1=x-xh1; 48xh1=x; 49f0=feval(func,x,varargin{:}); 50f1=zeros(size(f0,1),n); 51f_1=f1; 52for i=1:n 53 xh1(i)=x(i)+h1(i); 54 f1(:,i)=feval(func,xh1,varargin{:}); 55 xh1(i)=x(i)-h_1(i); 56 f_1(:,i)=feval(func,xh1,varargin{:}); 57 xh1(i)=x(i); 58end 59xh_1=xh1; 60hessian_mat = sparse(size(f0,1),n*n); 61 62for i=1:n 63 % if i > 1 64 % k=[i:n:n*(i-1)]; 65 % hessian_mat(:,(i-1)*n+1:(i-1)*n+i-1)=hessian_mat(:,k); 66 % hessian_mat(:,k)=0; 67 % end 68 hessian_mat(:,(i-1)*n+i)=(f1(:,i)+f_1(:,i)-2*f0)./(h1(i)*h_1(i)); 69 temp=f1+f_1-f0*ones(1,n); 70 for j=1:i-1 71 xh1(i)=x(i)+h1(i); 72 xh1(j)=x(j)+h_1(j); 73 xh_1(i)=x(i)-h1(i); 74 xh_1(j)=x(j)-h_1(j); 75 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)); 76 xh1(i)=x(i); 77 xh1(j)=x(j); 78 xh_1(i)=x(i); 79 xh_1(j)=x(j); 80 end 81end