declare Sqrt_Eps=1e-8; /* Square root of machine precission */
declare Diff_Eps2=1e-4; /* Rice's formula: ln(diff_eps)=ln(mach_prec)/4 */

proc dFiniteDiff2(x);
local t;
t=(abs(x)+sqrt_eps)*sqrt_eps;
t=t|diff_eps2;
/*
retp(max((abs(x)+Sqrt_Eps)*Sqrt_Eps,Diff_Eps2));
*/
retp(maxc(t));
endp;

proc hessian(&f,vP);
local i,j,pei,pj,hi,hj,fpp,fmm,fpm,fmp,cp,m,f0,f:proc;

    cp=rows(vP);
    m=zeros(cp,cp);

    f0=f(vP); /* get function at vP */

    i=0;do while i<cp;i=i+1;

        pei=vP[i];
        hi=dFiniteDiff2(pei);

        j=0;do while j<i;j=j+1;            /* off-diagonal elements */

            pj=vP[j];
            hj=dFiniteDiff2(pj);

            vP[i]=pei+hi;vP[j]=pj+hj;
            fpp=f(vP);                      /* +1 +1 */
            vP[i]=pei-hi;vP[j]=pj-hj;
            fmm=f(vP);                      /* -1 -1 */
            vP[i]=pei+hi;vP[j]=pj-hj;
            fpm=f(vP);                      /* +1 -1 */
            vP[i]=pei-hi;vP[j]=pj+hj;
            fmp=f(vP);                      /* -1 +1 */

            m[i,j]=(fpp+fmm)/(4*hi*hj)-(fpm+fmp)/(4*hi*hj);
            m[j,i]=m[i,j];

            vP[j]=pj;

        endo;

        vP[i]=pei+2*hi;
        fpp=f(vP);
        vP[i]=pei-2*hi;
        fmm=f(vP);

        m[i,i]=fpp/(4*hi*hi)+fmm/(4*hi*hi)-(2*f0)/(4*hi*hi);
        vP[i]=pei;

    endo;

    retp(m);

endp;

