hessian

 function hessian_mat = hessian(func,x,gstep,varargin)
 Computes second order partial derivatives

 INPUTS
    func        [string]   name of the function
    x           [double]   vector, the Hessian of "func" is evaluated at x.
    gstep       [double]   scalar, size of epsilon.
    varargin    [void]     list of additional arguments for "func".

 OUTPUTS
    hessian_mat [double]   Hessian matrix

 ALGORITHM
    Uses Abramowitz and Stegun (1965) formulas 25.3.24 and 25.3.27 p. 884

 SPECIAL REQUIREMENTS
    none

0001 function hessian_mat = hessian(func,x,gstep,varargin)
0002 % function hessian_mat = hessian(func,x,gstep,varargin)
0003 % Computes second order partial derivatives
0004 %
0005 % INPUTS
0006 %    func        [string]   name of the function
0007 %    x           [double]   vector, the Hessian of "func" is evaluated at x.
0008 %    gstep       [double]   scalar, size of epsilon.
0009 %    varargin    [void]     list of additional arguments for "func".
0010 %
0011 % OUTPUTS
0012 %    hessian_mat [double]   Hessian matrix
0013 %
0014 % ALGORITHM
0015 %    Uses Abramowitz and Stegun (1965) formulas 25.3.24 and 25.3.27 p. 884
0016 %
0017 % SPECIAL REQUIREMENTS
0018 %    none
0019 %
0020 
0021 % Copyright (C) 2001-2009 Dynare Team
0022 %
0023 % This file is part of Dynare.
0024 %
0025 % Dynare is free software: you can redistribute it and/or modify
0026 % it under the terms of the GNU General Public License as published by
0027 % the Free Software Foundation, either version 3 of the License, or
0028 % (at your option) any later version.
0029 %
0030 % Dynare is distributed in the hope that it will be useful,
0031 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0032 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0033 % GNU General Public License for more details.
0034 %
0035 % You should have received a copy of the GNU General Public License
0036 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0037 
0038 if ~isa(func, 'function_handle') 
0039     func = str2func(func);
0040 end
0041 n=size(x,1);
0042 h1=max(abs(x),sqrt(gstep(1))*ones(n,1))*eps^(1/6)*gstep(2);
0043 h_1=h1;
0044 xh1=x+h1;
0045 h1=xh1-x;
0046 xh1=x-h_1;
0047 h_1=x-xh1;
0048 xh1=x;
0049 f0=feval(func,x,varargin{:});
0050 f1=zeros(size(f0,1),n);
0051 f_1=f1;
0052 for i=1:n    
0053     xh1(i)=x(i)+h1(i);
0054     f1(:,i)=feval(func,xh1,varargin{:});
0055     xh1(i)=x(i)-h_1(i);
0056     f_1(:,i)=feval(func,xh1,varargin{:});
0057     xh1(i)=x(i);
0058     i=i+1;
0059 end
0060 xh_1=xh1;
0061 hessian_mat = zeros(size(f0,1),n*n);
0062 for i=1:n    
0063     if i > 1        
0064         k=[i:n:n*(i-1)];
0065         hessian_mat(:,(i-1)*n+1:(i-1)*n+i-1)=hessian_mat(:,k);
0066     end     
0067     hessian_mat(:,(i-1)*n+i)=(f1(:,i)+f_1(:,i)-2*f0)./(h1(i)*h_1(i));
0068     temp=f1+f_1-f0*ones(1,n);
0069     for j=i+1:n        
0070         xh1(i)=x(i)+h1(i);
0071         xh1(j)=x(j)+h_1(j);
0072         xh_1(i)=x(i)-h1(i);
0073         xh_1(j)=x(j)-h_1(j);
0074         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));
0075         xh1(i)=x(i);
0076         xh1(j)=x(j);
0077         xh_1(i)=x(i);
0078         xh_1(j)=x(j);
0079         j=j+1;
0080     end    
0081     i=i+1;
0082 end