function H = hessian(f,x,varargin)
% PURPOSE: Computes finite difference Hessian
% -------------------------------------------------------
% Usage:  H = hessian(func,x,varargin)
% Where: func = function name, fval = func(x,varargin)
%           x = vector of parameters (n x 1)
%    varargin = optional arguments passed to the function
% -------------------------------------------------------
% RETURNS:
%           H = finite differnce hessian
% -------------------------------------------------------
% ACKNOWLEDGEMENT: The original ARFIMA code was written by Jurgen A Doornik
% -----------------------------------------------------

% The matlab port written and changed by
% Ying Ma
% yingma@email.com

eps = 1e-5;

n = size(x,1);
fx = feval(f,x,varargin{:});
 

h = eps.^(1/3)*max(abs(x),1e-2);
xh = x+h;
h = xh-x;    
e = sparse(1:n,1:n,h,n,n);
 

g = zeros(n,1);
for i=1:n
  g(i) = feval(f,x+e(:,i),varargin{:});
end
   
H=h*h';

for i=1:n
for j=i:n
  H(i,j) = (feval(f,x+e(:,i)+e(:,j),varargin{:})-g(i)-g(j)+fx)/H(i,j);
  H(j,i) = H(i,j);
end
end
