mr_hessian

  [hessian_mat, gg, htol1, ihh, hh_mat0, hh1] = mr_hessian(init,x,func,hflag,htol0,varargin)

  numerical gradient and Hessian, with 'automatic' check of numerical
  error

 adapted from Michel Juillard original rutine hessian.m

  func =  function handle. The function must give two outputs:
    - the log-likelihood AND the single contributions at times t=1,...,T
    of the log-likelihood to compute outer product gradient
  x = parameter values
  hflag = 0, Hessian computed with outer product gradient, one point
  increments for partial derivatives in  gradients
  hflag = 1, 'mixed' Hessian: diagonal elements computed with numerical second order derivatives
             with correlation structure as from outer product gradient;
             two point evaluation of derivatives for partial derivatives
             in gradients
  hflag = 2, full numerical Hessian, computes second order partial derivatives
          uses Abramowitz and Stegun (1965) formulas 25.3.24 and 25.3.27
          p. 884.
  htol0 = 'precision' of increment of function values for numerical
  derivatives

  varargin: other parameters of func

0001 function [hessian_mat, gg, htol1, ihh, hh_mat0, hh1] = mr_hessian(init,x,func,hflag,htol0,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults)
0002 %  [hessian_mat, gg, htol1, ihh, hh_mat0, hh1] = mr_hessian(init,x,func,hflag,htol0,varargin)
0003 %
0004 %  numerical gradient and Hessian, with 'automatic' check of numerical
0005 %  error
0006 %
0007 % adapted from Michel Juillard original rutine hessian.m
0008 %
0009 %  func =  function handle. The function must give two outputs:
0010 %    - the log-likelihood AND the single contributions at times t=1,...,T
0011 %    of the log-likelihood to compute outer product gradient
0012 %  x = parameter values
0013 %  hflag = 0, Hessian computed with outer product gradient, one point
0014 %  increments for partial derivatives in  gradients
0015 %  hflag = 1, 'mixed' Hessian: diagonal elements computed with numerical second order derivatives
0016 %             with correlation structure as from outer product gradient;
0017 %             two point evaluation of derivatives for partial derivatives
0018 %             in gradients
0019 %  hflag = 2, full numerical Hessian, computes second order partial derivatives
0020 %          uses Abramowitz and Stegun (1965) formulas 25.3.24 and 25.3.27
0021 %          p. 884.
0022 %  htol0 = 'precision' of increment of function values for numerical
0023 %  derivatives
0024 %
0025 %  varargin: other parameters of func
0026 
0027 % Copyright (C) 2004-2011 Dynare Team
0028 %
0029 % This file is part of Dynare.
0030 %
0031 % Dynare is free software: you can redistribute it and/or modify
0032 % it under the terms of the GNU General Public License as published by
0033 % the Free Software Foundation, either version 3 of the License, or
0034 % (at your option) any later version.
0035 %
0036 % Dynare is distributed in the hope that it will be useful,
0037 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0038 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0039 % GNU General Public License for more details.
0040 %
0041 % You should have received a copy of the GNU General Public License
0042 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0043 
0044 persistent h1 htol
0045 
0046 n=size(x,1);
0047 if init
0048     gstep_=DynareOptions.gstep;
0049     htol = 1.e-4;
0050     h1=DynareOptions.gradient_epsilon*ones(n,1);
0051     return
0052 end
0053 
0054 [f0, ff0]=feval(func,x,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
0055 h2=BayesInfo.ub-BayesInfo.lb;
0056 hmax=BayesInfo.ub-x;
0057 hmax=min(hmax,x-BayesInfo.lb);
0058 
0059 h1 = min(h1,0.5.*hmax);
0060 
0061 if htol0<htol
0062     htol=htol0;
0063 end
0064 xh1=x;
0065 f1=zeros(size(f0,1),n);
0066 f_1=f1;
0067 ff1=zeros(size(ff0));
0068 ff_1=ff1;
0069 ggh=zeros(size(ff0,1),n);
0070 
0071 i=0;
0072 while i<n
0073     i=i+1;
0074     h10=h1(i);
0075     hcheck=0;
0076     xh1(i)=x(i)+h1(i);
0077     try
0078         [fx, ffx]=feval(func,xh1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
0079     catch
0080         fx=1.e8;
0081     end
0082     it=1;
0083     dx=(fx-f0);
0084     ic=0;
0085     icount = 0;
0086     h0=h1(i);
0087     while (abs(dx(it))<0.5*htol || abs(dx(it))>(3*htol)) && icount<10 && ic==0
0088         icount=icount+1;
0089         if abs(dx(it))<0.5*htol
0090             if abs(dx(it)) ~= 0,
0091                 h1(i)=min(max(1.e-10,0.3*abs(x(i))), 0.9*htol/abs(dx(it))*h1(i));
0092             else
0093                 h1(i)=2.1*h1(i);
0094             end
0095             h1(i) = min(h1(i),0.5*hmax(i));
0096             xh1(i)=x(i)+h1(i);
0097             try
0098                 [fx, ffx]=feval(func,xh1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
0099             catch
0100                 fx=1.e8;
0101             end
0102         end
0103         if abs(dx(it))>(3*htol)
0104             h1(i)= htol/abs(dx(it))*h1(i);
0105             xh1(i)=x(i)+h1(i);
0106             try
0107                 [fx, ffx]=feval(func,xh1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
0108             catch
0109                 fx=1.e8;
0110             end
0111             while (fx-f0)==0
0112                 h1(i)= h1(i)*2;
0113                 xh1(i)=x(i)+h1(i);
0114                 [fx, ffx]=feval(func,xh1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
0115                 ic=1;
0116             end
0117         end
0118         it=it+1;
0119         dx(it)=(fx-f0);
0120         h0(it)=h1(i);
0121         if (h1(i)<1.e-12*min(1,h2(i)) && h1(i)<0.5*hmax(i))% || (icount==10 &&  abs(dx(it))>(3*htol)),
0122             ic=1;
0123             hcheck=1;
0124         end
0125     end
0126     f1(:,i)=fx;
0127     if any(isnan(ffx))
0128         ff1=ones(size(ff0)).*fx/length(ff0);
0129     else
0130         ff1=ffx;
0131     end
0132     xh1(i)=x(i)-h1(i);
0133     [fx, ffx]=feval(func,xh1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
0134     f_1(:,i)=fx;
0135     if any(isnan(ffx))
0136         ff_1=ones(size(ff0)).*fx/length(ff0);
0137     else
0138         ff_1=ffx;
0139     end
0140     ggh(:,i)=(ff1-ff_1)./(2.*h1(i));
0141     xh1(i)=x(i);
0142     if hcheck && htol<1
0143         htol=min(1,max(min(abs(dx))*2,htol*10));
0144         h1(i)=h10;
0145         i=0;
0146     end
0147 end
0148 
0149 h_1=h1;
0150 xh1=x;
0151 xh_1=xh1;
0152 
0153 gg=(f1'-f_1')./(2.*h1);
0154 
0155 if hflag==2
0156     gg=(f1'-f_1')./(2.*h1);
0157     hessian_mat = zeros(size(f0,1),n*n);
0158     for i=1:n
0159         if i > 1
0160             k=[i:n:n*(i-1)];
0161             hessian_mat(:,(i-1)*n+1:(i-1)*n+i-1)=hessian_mat(:,k);
0162         end
0163         hessian_mat(:,(i-1)*n+i)=(f1(:,i)+f_1(:,i)-2*f0)./(h1(i)*h_1(i));
0164         temp=f1+f_1-f0*ones(1,n);
0165         for j=i+1:n
0166             xh1(i)=x(i)+h1(i);
0167             xh1(j)=x(j)+h_1(j);
0168             xh_1(i)=x(i)-h1(i);
0169             xh_1(j)=x(j)-h_1(j);
0170             temp1 = feval(func,xh1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
0171             temp2 = feval(func,xh_1,DynareDataset,DynareOptions,Model,EstimatedParameters,BayesInfo,DynareResults);
0172             hessian_mat(:,(i-1)*n+j)=-(-temp1 -temp2+temp(:,i)+temp(:,j))./(2*h1(i)*h_1(j));
0173             xh1(i)=x(i);
0174             xh1(j)=x(j);
0175             xh_1(i)=x(i);
0176             xh_1(j)=x(j);
0177             j=j+1;
0178         end
0179         i=i+1;
0180     end
0181 elseif hflag==1
0182     hessian_mat = zeros(size(f0,1),n*n);
0183     for i=1:n
0184         dum = (f1(:,i)+f_1(:,i)-2*f0)./(h1(i)*h_1(i));
0185         if dum>eps
0186             hessian_mat(:,(i-1)*n+i)=dum;
0187         else
0188             hessian_mat(:,(i-1)*n+i)=max(eps, gg(i)^2);
0189         end
0190     end
0191 end
0192 
0193 gga=ggh.*kron(ones(size(ff1)),2.*h1');  % re-scaled gradient
0194 hh_mat=gga'*gga;  % rescaled outer product hessian
0195 hh_mat0=ggh'*ggh;  % outer product hessian
0196 A=diag(2.*h1);  % rescaling matrix
0197 % igg=inv(hh_mat);  % inverted rescaled outer product hessian
0198 ihh=A'*(hh_mat\A);  % inverted outer product hessian
0199 if hflag>0 && min(eig(reshape(hessian_mat,n,n)))>0
0200     hh0 = A*reshape(hessian_mat,n,n)*A';  %rescaled second order derivatives
0201     hh = reshape(hessian_mat,n,n);  %rescaled second order derivatives
0202     sd0=sqrt(diag(hh0));   %rescaled 'standard errors' using second order derivatives
0203     sd=sqrt(diag(hh_mat));  %rescaled 'standard errors' using outer product
0204     hh_mat=hh_mat./(sd*sd').*(sd0*sd0');  %rescaled inverse outer product with 'true' std's
0205     igg=inv(hh_mat);   % rescaled outer product hessian with 'true' std's
0206     ihh=A'*(hh_mat\A);  % inverted outer product hessian
0207     hh_mat0=inv(A)'*hh_mat*inv(A);  % outer product hessian with 'true' std's
0208     sd=sqrt(diag(ihh));   %standard errors
0209     sdh=sqrt(1./diag(hh));   %diagonal standard errors
0210     for j=1:length(sd)
0211         sd0(j,1)=min(BayesInfo.p2(j), sd(j));  %prior std
0212         sd0(j,1)=10^(0.5*(log10(sd0(j,1))+log10(sdh(j,1))));
0213     end
0214     ihh=ihh./(sd*sd').*(sd0*sd0');  %inverse outer product with modified std's
0215     igg=inv(A)'*ihh*inv(A);  % inverted rescaled outer product hessian with modified std's
0216     hh_mat=inv(igg);   % outer product rescaled hessian with modified std's
0217     hh_mat0=inv(A)'*hh_mat*inv(A);  % outer product hessian with modified std's
0218                                     %     sd0=sqrt(1./diag(hh0));   %rescaled 'standard errors' using second order derivatives
0219                                     %     sd=sqrt(diag(igg));  %rescaled 'standard errors' using outer product
0220                                     %     igg=igg./(sd*sd').*(sd0*sd0');  %rescaled inverse outer product with 'true' std's
0221                                     %     hh_mat=inv(igg);   % rescaled outer product hessian with 'true' std's
0222                                     %     ihh=A'*igg*A;  % inverted outer product hessian
0223                                     %     hh_mat0=inv(A)'*hh_mat*inv(A);  % outer product hessian with 'true' std's
0224 end
0225 if hflag<2
0226     hessian_mat=hh_mat0(:);
0227 end
0228 
0229 if any(isnan(hessian_mat))
0230     hh_mat0=eye(length(hh_mat0));
0231     ihh=hh_mat0;
0232     hessian_mat=hh_mat0(:);
0233 end
0234 hh1=h1;
0235 htol1=htol;
0236 save hess.mat hessian_mat