2:53 p.m.
Hi,
for the standard Poisson model I'd like to explore different methods of
covariance estimation.
Could you please tell me how to replicate the --hessian flag with the
'Hess' function?
(it's not included in MLE-advanced.inp, I read in the manuel 'negative
inverse of the hessian', but I make something wrong since invpd(h) does
not do the job).
Have a nice weekend
Leon
<hansl>
open poisson.gdt
poisson y 0 x1 x2
series fake_y = ln(y+1)
ols fake_y 0 x1 x2 --quiet
list xList = $xlist
matrix b = $coeff
matrix mX = {xList}'
function matrix score(matrix b, series y, matrix mX)
series e = y - exp(mX'b)
return {e} .* mX'
end function
function void Hess(matrix *H, matrix b, series y, matrix mX)
#computes the negative Hessian for Poisson model,
H = -(qform(mX, (mX'b)*(mX'b)'))
end function
matrix H = {}
mle loglik = y*Xb - exp(Xb)
series Xb = mX'b
deriv b = score(b, y, mX)
hessian Hess(&H, b, y, mX)
end mle
# How to replicate --hessian with the 'Hess' function ?
mle loglik = y*Xb - exp(Xb)
series Xb = mX'b
deriv b = score(b, y, mX)
end mle --hessian
<hansl>
1:42 p.m.
On Sat, 23 Nov 2013, Pindar wrote:
Hi,
for the standard Poisson model I'd like to explore different methods of
covariance estimation.
Could you please tell me how to replicate the --hessian flag with the 'Hess'
function?
(it's not included in MLE-advanced.inp, I read in the manuel 'negative
inverse of the hessian', but I make something wrong since invpd(h) does not
do the job).
The problem is that your "Hess" function has a bug. Have a look here:
<hansl>
open poisson.gdt
poisson y 0 x1 x2
series fake_y = ln(y+1)
ols fake_y 0 x1 x2 --quiet
list xList = $xlist
matrix b = $coeff
matrix mX = {xList}
function matrix score(series y, series m, matrix mX)
return {y - m} .* mX
end function
function void Hess(matrix *H, series m, matrix mX)
#computes the negative Hessian for Poisson model
H = mX'(mX .* {m})
end function
matrix H = {}
mle loglik = y*xb - m - lngamma(y+1)
series xb = mX*b
series m = exp(xb)
deriv b = score(y, m, mX)
hessian Hess(&H, m, mX)
end mle --hessian
</hansl>
-------------------------------------------------------
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
-------------------------------------------------------
4:24 p.m.
Am 24.11.2013 14:42, schrieb Riccardo (Jack) Lucchetti:
On Sat, 23 Nov 2013, Pindar wrote:
> Hi,
>
> for the standard Poisson model I'd like to explore different methods
> of covariance estimation.
> Could you please tell me how to replicate the --hessian flag with
> the 'Hess' function?
> (it's not included in MLE-advanced.inp, I read in the manuel
> 'negative inverse of the hessian', but I make something wrong since
> invpd(h) does not do the job).
The problem is that your "Hess" function has a bug. Have a look here:
I'm still puzzled since my crude "Hess" function produced the same
results without the --hessian flag as the right one. For didactic
reasons I wanted to use qform in order to have the X'WX syntax and
thought this is qual to "qform(mX', {m}*{m}')" ). However, now with the
plain (correct) matrix notation invpd() works and if one includes "eval
diag(sqrt(invpd(mX'(mX .* {m}))))" in the "Hess" function the SE
corresponding to the --hessian flag are printed. But I'm a right, that
it is not possible to reproduce --hessian SE without giving this flag
just by tweaking the "Hess" function since the --hessian SE are somehow
not passed when using "H = inv(mX'(mX .* {m}))" ?
<hansl>
open poisson.gdt
poisson y 0 x1 x2
series fake_y = ln(y+1)
ols fake_y 0 x1 x2 --quiet
list xList = $xlist
matrix b = $coeff
matrix mX = {xList}
function matrix score(series y, series m, matrix mX)
return {y - m} .* mX
end function
function void Hess(matrix *H, series m, matrix mX)
#computes the negative Hessian for Poisson model
H = mX'(mX .* {m})
end function
matrix H = {}
mle loglik = y*xb - m - lngamma(y+1)
series xb = mX*b
series m = exp(xb)
deriv b = score(y, m, mX)
hessian Hess(&H, m, mX)
end mle --hessian
</hansl>
-------------------------------------------------------
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
-------------------------------------------------------
_______________________________________________
Gretl-users mailing list
Gretl-users(a)lists.wfu.edu
http://lists.wfu.edu/mailman/listinfo/gretl-users
5:23 p.m.
On Sun, 24 Nov 2013, Pindar wrote:
Am 24.11.2013 14:42, schrieb Riccardo (Jack) Lucchetti:
> On Sat, 23 Nov 2013, Pindar wrote:
>
>> Hi,
>>
>> for the standard Poisson model I'd like to explore different methods of
>> covariance estimation.
>> Could you please tell me how to replicate the --hessian flag with the
>> 'Hess' function?
>> (it's not included in MLE-advanced.inp, I read in the manuel 'negative
>> inverse of the hessian', but I make something wrong since invpd(h) does
>> not do the job).
>
> The problem is that your "Hess" function has a bug. Have a look here:
>
I'm still puzzled since my crude "Hess" function produced the same results
without the --hessian flag as the right one.
Let's put things in perspective: your function had two main problems. For
a start, it didn't return a pd matrix (as the negative Hessian should be,
at least at the maximum). Second, the calculations were wrong. See here:
<hansl>
set echo off
set messages off
open poisson.gdt
series fake_y = ln(y+1)
ols fake_y 0 x1 x2 --quiet
list xList = $xlist
matrix b = $coeff
matrix mX = {xList}
series xb = mX*b
series m = exp(xb)
function void Hess_leon(matrix *H, matrix b, series y, matrix mX)
#computes the negative Hessian for Poisson model,
H = -(qform(mX', (mX*b)*(mX*b)'))
end function
function void Hess_jack(matrix *H, series m, matrix mX)
#computes the negative Hessian for Poisson model
H = mX'(mX .* {m})
end function
matrix Hj Hl
Hess_jack(&Hj, m, mX)
Hess_leon(&Hl, b, y, mX)
print Hj Hl
</hansl>
For both these reasons, mle found something that was unusable for
computing $vcv, so it switched to OPG; if you take your original code and
add the --hessian flag to the first "end mle" invocation (the one where
the "Hess" function is used) gretl tells you what's going on before
displaying the estimates.
<output>
Error: Hessian non-negative definite? (err = 44); dropping back to OPG
Tolerance = 1.81899e-12
</output>
Having said this, a word of clarification maybe is needed: supplying an
analytical Hessian does not mean it gets used in any way. Usage of the
Hessian in MLE problems happens in two independent circumstances:
a) during optimisation if Newton-Raphson is used instead of BFGS (not by
default)
b) in order to compute the $vcv, but you have to say so explicitly by
giving mle either the --hessian or the --robust options.
-------------------------------------------------------
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
-------------------------------------------------------
4054
days inactive
4055
days old
3 comments
2 participants
participants (2)
-
Pindar -
Riccardo (Jack) Lucchetti