On 26/09/2024 00:02, Alecos Papadopoulos wrote:
Thanks both for the input.
The reason I asked whether we have such functionality in gretl, is not
in order to "go along" with an impaired Hessian, but for diagnostic
purposes: to give a related example, when an mle command fails to
converge, personally it helps me to inspect in detail and trace
iteration-by-iteration how the estimator lost its way. In the same
spirit, obtaining an impaired Hessian, would tell me which coefficients
under estimation are linked to the problem, and this could possibly help
in realizing the source/cause of the failure.
I understand that what you want is a way to inspect the Hessian after
convergence has taken place. To do this, you may find the "numhess"
function useful. For example, here I use mle to estimate the parameters
of a gamma density and then calculate the Hessian.
<hansl>
function series llik(const series x, scalar theta, scalar p)
series ret = -p * log(theta) - lngamma(p) + (p-1) * log(x) - x/theta
return ret
end function
function scalar total_llik(matrix param, const series x)
series l = llik(x, param[1], param[2])
return sum(l)
end function
nulldata 128
x = -log(uniform())
a = 1
p = 1
mle loglik = llik(x, a, p)
params a p
end mle
theta = {a, p}
H = numhess(theta, "total_llik(theta, x)")
</hansl>
Having said this, there are several points to note: if the Hessian turns
out to be non-negative definite, it could be because you're on a saddle
point, like John said, which is unlikely; in common problems, it's far
more likely that this may be an artifact of numerical differentiation.
In the latter case, the best thing is ALWAYS to supply analytical
derivatives, but in case it's impossible, or too complicated, there are
several strategies that may help:
* try rescaling your data: numerical differentiation may struggle with
parameters that are "too large" or "too small";
* reparametrisation is also very helpful in some cases; for example, you
don't have to estimate a variance as-is. You could estimate its log, or
its square root.
* even if you don't have a function for the Hessian, the analytical
score can be a big help: compute the score and use the fdjac() function
on that.
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------