

% Hessian
% It returns the standard errors in sevec
% [sevec,nodif,HH]=hessy(thetstar,nonev,se,namellcode,accuracy,varargin);
% thetstar=optimum
% nonev=parameter to be held constant (ie w.r. to which we s.e. are
% not calculated); type [] if want all s.e.
% nodif is a 0-1 vector that can be used to extract the estimates
% for which s.e. have been calculated 
% eg xo(nodif==1,nodif==1) will return the optima (xo) for which s.e.
% exist. This allows to compute t-ratio easily 
% t_ratio=(eg xo(nodif==1,nodif==1)./sevec
% se=1 if transformed estimates, 0 if just estimates
% G. Lombardo, Bundesbank 2002

% below, llfnsub is the likelihood evaluation algorithm

function [sevec,nodif,HH]=hessy(thetstar,nonev,se1,llfnsub,accur,varargin)
ob=varargin{1:end};
global bigpar se
se=se1; 
eee = accur;
disp(['step of numerical differentiation = ',num2str(eee)])
nodif=ones(length(thetstar),1);
nodif(nonev)=0;
scalvec=ones(size(thetstar,1),1);%
%abs(thetstar./max(abs(thetstar))).^(-1);;%%abs(thetstory./max(abs(thetstory))).^(-1);%
%scalvec=[1 10 10 10 100 10 1 1000]';
%it should be so that the increase in the parameters is of equal size
% yet then the evaluation is off the optimum...
  tstar= real(thetstar);
  scalinv = inv(diag(scalvec));
  tstars = tstar;%diag(scalvec)*tstar;% this is undone in llfn3
 eval([' fstar = ',llfnsub,'(tstars,ob)']);
e =diag(max(diag(diag(eee)),diag(diag(eee))*thetstar));
%e=eee*eye(length(scalvec));

nth=length(tstars);
Hs=zeros(nth,nth);
%for i=1:nth
%   Ht(i,i)=2*(llfn3(tstars+e(:,i))-fstar)/e(i,i)^2;
%end
% Judd p.36 two sided derivative
u=1;
h = waitbar(0,'Please wait...');
for i=1:nth
   %if nodif(i)
     eval([' Hs(i,i)=(',llfnsub, '(tstars+e(:,i),ob)+' , llfnsub ,' (tstars-e(:,i),ob)-2*fstar)/e(i,i)^2;']);
      %end
waitbar(u/(2*nth-1))
u=u+1;
%Hs(i,i)
end
H=Hs;
for i=1:nth-1
   %if nodif(i)
      for j=i+1:nth
    %     if nodif(j)
           eval([' H(i,j)=(' , llfnsub, '(tstars+e(:,i)+e(:,j),ob)-',... 
               llfnsub,'(tstars+e(:,i),ob)-',...
               llfnsub,'(tstars+e(:,j),ob)+fstar)/(e(i,i)*e(j,j));'])
            %    end
waitbar(u/(2*nth-1))

end
u=u+1;
   %end
end
close(h)
HU = H'-diag(diag(H));

  HH = H+ HU;
HH=HH(nodif==1,nodif==1);
  bighx = scalinv(nodif==1,nodif==1)*inv(HH)*scalinv(nodif==1,nodif==1)';

sevec = sqrt(diag(bighx));
se=0;
lm=eig(HH);
lneg=lm(lm<0);
if length(lneg)>0
   disp('The Hessian is not positive semidefinite')
end
 
