Statalist


[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

Re: st: condition of Hessian


From   "David M. Drukker" <ddrukker@stata.com>
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: condition of Hessian
Date   Thu, 20 Mar 2008 09:22:16 -0500 (CDT)

Teresa Dale Nelson <tdale@utdallas.edu> asked for an example of how to
compute the condition number of the Hessian from an estimator reporting an
observed information matrix (OIM) variance.

Here is an example, using a silly Poisson regression:

. set mem 200m
(204800k)

. webuse nlswork
(National Longitudinal Survey. Young Women 14-26 years of age in 1968)

. poisson grade nev_mar not_smsa union

Iteration 0: log likelihood = -46328.9 Iteration 1: log likelihood = -46328.9

Poisson regression Number of obs = 19222
LR chi2(3) = 168.52
Prob > chi2 = 0.0000
Log likelihood = -46328.9 Pseudo R2 = 0.0018

------------------------------------------------------------------------------
grade | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
nev_mar | .0162206 .0050936 3.18 0.001 .0062374 .0262038
not_smsa | -.0507212 .0045589 -11.13 0.000 -.0596564 -.0417859
union | .0209936 .0047465 4.42 0.000 .0116907 .0302966
_cons | 2.552624 .0028711 889.06 0.000 2.546997 2.558251
------------------------------------------------------------------------------

. . mata: V1 = st_matrix("e(V)")

. mata: cond(V1)
8.835376744

This condition number is close to one, as far as condition numbers go,
indicating a well-defined problem.

This condition number is based on the 2 norm,

cond = = norm(A, 2) * norm(A^(-1), 2)

where norm(A, 2) is the largest singular value of A. This calculation will
yield the same answer for A and A^(-1), ignoring precision problems.

Here are some illustrative computations

. mata: norm(V1,2)
.0000300734

. mata: svdsv(V1)
1
+---------------+
1 | .0000300734 |
2 | .0000244259 |
3 | .0000195973 |
4 | 3.40375e-06 |
+---------------+

. mata: norm(V1,2)*norm(cholinv(V1),2)
8.835376744


Now I claimed that 8.835 was close to one. Here I construct a badly
conditioned problem and calculate the condition number of the OIM variance
estimator.

. generate prob = union + .001*uniform() (9296 missing values generated)

. poisson grade nev_mar not_smsa union prob

Iteration 0: log likelihood = -46328.894 Iteration 1: log likelihood = -46328.894

Poisson regression Number of obs = 19222
LR chi2(4) = 168.53
Prob > chi2 = 0.0000
Log likelihood = -46328.894 Pseudo R2 = 0.0018

------------------------------------------------------------------------------
grade | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
nev_mar | .0162228 .0050936 3.18 0.001 .0062395 .0262061
not_smsa | -.0507209 .0045589 -11.13 0.000 -.0596561 -.0417857
union | -.7705006 6.98085 -0.11 0.912 -14.45272 12.91171
prob | .7914966 6.980869 0.11 0.910 -12.89076 14.47375
_cons | 2.552229 .0045136 565.46 0.000 2.543383 2.561076
------------------------------------------------------------------------------

. mata: V2 = st_matrix("e(V)")

. mata: cond(V2)
30757045.59

This condition number is not close to one, indicating that the problem is
not well defined.

The data must be on the same scale to use condition numbers as a metric to
measure identification.

Here I illustrate by rescaling the union indicator variable. I compute the
correlation matrix before and after the rescaling to illustrate that the
correlations are invariant to the rescaling.


. corr grade nev_mar not_smsa union (obs=19222)

| grade nev_mar not_smsa union
-------------+------------------------------------
grade | 1.0000
nev_mar | 0.0450 1.0000
not_smsa | -0.1275 -0.0759 1.0000
union | 0.0573 0.0195 -0.0684 1.0000


. replace union = 100000*union
union was byte now long
(4510 real changes made)

. corr grade nev_mar not_smsa union (obs=19222)

| grade nev_mar not_smsa union
-------------+------------------------------------
grade | 1.0000
nev_mar | 0.0450 1.0000
not_smsa | -0.1275 -0.0759 1.0000
union | 0.0573 0.0195 -0.0684 1.0000


Now I refit the model and calculate the condition number of the OIM
variance.

. . poisson grade nev_mar not_smsa union

Iteration 0: log likelihood = -46328.9 Iteration 1: log likelihood = -46328.9

Poisson regression Number of obs = 19222
LR chi2(3) = 168.52
Prob > chi2 = 0.0000
Log likelihood = -46328.9 Pseudo R2 = 0.0018

------------------------------------------------------------------------------
grade | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
nev_mar | .0162206 .0050936 3.18 0.001 .0062374 .0262038
not_smsa | -.0507212 .0045589 -11.13 0.000 -.0596564 -.0417859
union | 2.10e-07 4.75e-08 4.42 0.000 1.17e-07 3.03e-07
_cons | 2.552624 .0028711 889.06 0.000 2.546997 2.558251
------------------------------------------------------------------------------

. . mata: V2 = st_matrix("e(V)")

. mata: cond(V2)
1.70881e+10

The rescaling has made this well-defined problem appear ill defined.

Drukker and Wiggins (2004) discuss these issues using real data and a
different norm. Drukker and Wiggins (2004) place the problem in terms of
parameter identification.

I hope that this helps.

--David

References
---------

Drukker, D. M. and V. Wiggins. 2004. "Verifying the Solution from a
Nonlinear Solver: A Case Study: Comment". The American Economic Review
94(1):397-399


*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/




© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index