O Hessiano empírico de um estimador M pode ser indefinido?

Jeffrey Wooldridge em sua Análise Econométrica de Dados de Painel e Seção Transversal (página 357) diz que o Hessian empírico "não é garantido que seja definido positivo, ou mesmo semidefinido positivo, para a amostra em particular com a qual estamos trabalhando".

Isso me parece errado, pois (problemas numéricos à parte) o Hessian deve ser positivo semidefinido como resultado da definição do estimador M como o valor do parâmetro que minimiza a função objetivo da amostra em questão e o fato bem conhecido de que no mínimo (local), o hessiano é semidefinido positivo.

Meu argumento está certo?

[EDIT: A declaração foi removida no 2º ed. do livro. Veja o comentário.]

ANTECEDENTES Suponha que seja um estimador obtido minimizando que denota a ésima observação.$\widehat \theta_N$wi

Vamos denotar o Hessian de por$q$$H$ ,

A covariância assintótica de θ n envolve E [ H ( q , q 0 ) ] onde θ 0 é o verdadeiro valor do parâmetro. Uma maneira de estimar isso é usar o método Hesssiano empírico$\widehat \theta_n$$E[H(q,\theta_0)]$$\theta_0$

É o definiteness de H que está em questão.$\widehat H$

Eu acho que você está certo. Vamos destilar seu argumento em sua essência:

1. minimiza a funçãoQdefinido comoQ(θ)=$\widehat \theta_N$$Q$$Q(\theta) = {1 \over N}\sum_{i=1}^N q(w_i,\theta).$

2. Seja o hessiano de Q , de onde H ( θ ) = ∂ 2$H$$Q$ por definição e isso, por sua vez, por linearidade de diferenciação, é igual a$H(\theta) = \frac{\partial^2 Q}{\partial \theta_i \partial \theta_j}$.$\frac{1}{N}\sum_{i=1}^N H(w_i, \theta_n)$

3. Assumindo q N reside no interior do domínio de Q , em seguida, H ( θ N ) deve ser positivo semi-definida.$\widehat \theta_N$$Q$$H(\widehat \theta_N)$

Essa é apenas uma afirmação sobre a função : como ela é definida é meramente uma distração, exceto na medida em que a diferenciabilidade assumida de q de segunda ordem em relação ao seu segundo argumento ( θ ) assegura a diferenciabilidade de Q de segunda ordem .$Q$$q$$\theta$$Q$

Encontrar estimadores-M pode ser complicado. Considere estes dados fornecidos pelo @mpiktas:

{1.168042, 0.3998378}, {1.807516, 0.5939584}, {1.384942, 3.6700205}, {1.327734, -3.3390724}, {1.602101, 4.1317608}, {1.604394, -1.9045958}, {1.124633, -3.0865249}, {1.294601, -1.8331763},{1.577610, 1.0865977}, { 1.630979, 0.7869717}

O procedimento R para encontrar o estimador M com produziu a solução ( c 1 , c 2 ) = ( - 114,91316 , - 32,54386 ) . O valor da função objetivo (a média dos q 's) nesse ponto é igual a 62,3542. Aqui está um gráfico do ajuste:$q((x,y),\theta)=(y-c_1x^{c_2})^4$$(c_1, c_2)$$(-114.91316, -32.54386)$$q$

Aqui está um gráfico da função objetivo (log) em uma vizinhança desse ajuste:

Algo é suspeito aqui: os parâmetros do ajuste estão muito distantes dos parâmetros usados ​​para simular os dados (próximo ) e não parecemos estar no mínimo: estamos em um vale extremamente raso, inclinado para valores maiores de ambos os parâmetros:$(0.3, 0.2)$

O determinante negativo do Hessian neste momento confirma que este não é um mínimo local! No entanto, quando você olha para os rótulos do eixo z, é possível ver que essa função é plana com precisão de cinco dígitos em toda a região, porque é igual a uma constante 4.1329 (o logaritmo de 62.354). Isso provavelmente levou o minimizador da função R (com suas tolerâncias padrão) a concluir que estava próximo do mínimo.

De fato, a solução está longe deste ponto. Para ter certeza de encontrá-lo, empreguei o método " Eixo Principal " computacionalmente caro, mas altamente eficaz , no Mathematica , usando precisão de 50 dígitos (base 10) para evitar possíveis problemas numéricos. Ele encontra um mínimo próximo onde a função objetivo tem o valor 58.292655: cerca de 6% menor que o "mínimo" encontrado por R. Esse mínimo ocorre em uma seção extremamente plana , mas posso fazer com que pareça (apenas por pouco) um mínimo verdadeiro, com contornos elípticos, exagerando a c $(c_1, c_2) = (0.02506, 7.55973)$$c_2$ direção na trama:

Os contornos variam de 58,29266 no meio até 58,29284 nos cantos (!). Aqui está a visualização 3D (novamente do objetivo do log):

Aqui, o hessiano é positivo-definido: seus valores próprios são 55062.02 e 0.430978. Portanto, este ponto é um mínimo local (e provavelmente um mínimo global). Aqui está o ajuste a que corresponde:

Eu acho que é melhor que o outro. Os valores dos parâmetros são certamente mais realistas e está claro que não conseguiremos fazer muito melhor com essa família de curvas.

Existem lições úteis que podemos extrair deste exemplo:

1. A otimização numérica pode ser difícil, especialmente com as funções de ajuste não linear e de perda não quadrática. Portanto:
2. Verifique os resultados da maneira mais possível, incluindo:
3. Faça um gráfico da função objetivo sempre que puder.
4. Quando os resultados numéricos parecerem violar os teoremas matemáticos, seja extremamente suspeito.
5. Quando os resultados estatísticos forem surpreendentes - como os surpreendentes valores de parâmetros retornados pelo código R -, seja extremamente suspeito.

$\hat{\theta}_N$

$\hat{\theta}_N$$\Theta$$\hat{H}$

$N^{-1}\sum_{i=1}^Nq(w_i,\theta)$$\theta_0$

$N^{-1}\sum_{i=1}^Nq(w_i,\theta)$$\Theta$

Mais adiante, em seu livro, Wooldridge dá exemplos de estimativas de Hessian que são garantidas como definidas numericamente positivas. Na prática, a definição não positiva de Hessian deve indicar que a solução está no ponto limite ou o algoritmo não conseguiu encontrar a solução. O que geralmente é mais uma indicação de que o modelo ajustado pode ser inadequado para um dado dado.

Aqui está o exemplo numérico. Gero um problema de mínimos quadrados não lineares:

$X$$[1,2]$$\varepsilon$$\sigma^2$set.seed(3)$x_i$$y_i$

Escolhi a função objetiva quadrado da função objetiva de mínimos quadrados não linear usual:

Here is the code in R for optimising function, its gradient and hessian.

##First set-up the epxressions for optimising function, its gradient and hessian.
##I use symbolic derivation of R to guard against human error    
mt <- expression((y-c1*x^c2)^4)

gradmt <- c(D(mt,"c1"),D(mt,"c2"))

hessmt <- lapply(gradmt,function(l)c(D(l,"c1"),D(l,"c2")))

##Evaluate the expressions on data to get the empirical values. 
##Note there was a bug in previous version of the answer res should not be squared.
optf <- function(p) {
    res <- eval(mt,list(y=y,x=x,c1=p[1],c2=p[2]))
    mean(res)
}

gf <- function(p) {
    evl <- list(y=y,x=x,c1=p[1],c2=p[2]) 
    res <- sapply(gradmt,function(l)eval(l,evl))
    apply(res,2,mean)
}

hesf <- function(p) {
    evl <- list(y=y,x=x,c1=p[1],c2=p[2]) 
    res1 <- lapply(hessmt,function(l)sapply(l,function(ll)eval(ll,evl)))
    res <- sapply(res1,function(l)apply(l,2,mean))
    res
}

First test that gradient and hessian works as advertised.

set.seed(3)
x <- runif(10,1,2)
y <- 0.3*x^0.2

> optf(c(0.3,0.2))
[1] 0
> gf(c(0.3,0.2))
[1] 0 0
> hesf(c(0.3,0.2))
     [,1] [,2]
[1,]    0    0
[2,]    0    0
> eigen(hesf(c(0.3,0.2)))$values
[1] 0 0

The hessian is zero, so it is positive semi-definite. Now for the values of $x$ and $y$ given in the link we get

> df <- read.csv("badhessian.csv")
> df
          x          y
1  1.168042  0.3998378
2  1.807516  0.5939584
3  1.384942  3.6700205
4  1.327734 -3.3390724
5  1.602101  4.1317608
6  1.604394 -1.9045958
7  1.124633 -3.0865249
8  1.294601 -1.8331763
9  1.577610  1.0865977
10 1.630979  0.7869717
> x <- df$x
> y <- df$y
> opt <- optim(c(1,1),optf,gr=gf,method="BFGS")  
> opt$par
[1] -114.91316  -32.54386
> gf(opt$par)
[1] -0.0005795979 -0.0002399711
> hesf(opt$par)
              [,1]         [,2]
[1,]  0.0002514806 -0.003670634
[2,] -0.0036706345  0.050998404
> eigen(hesf(opt$par))$values
[1]  5.126253e-02 -1.264959e-05

Gradient is zero, but the hessian is non positive.

Note: This is my third attempt to give an answer. I hope I finally managed to give precise mathematical statements, which eluded me in the previous versions.

O hessian é indefinido em um ponto de sela. É possível que este seja o único ponto estacionário no interior do espaço de parâmetros.

Atualização: Deixe-me elaborar. Primeiro, vamos supor que o Hessian empírico exista em toda parte.

E se $\hat{\theta}_n$ é um mínimo local (ou mesmo global) de $\sum_i q(w_i, \cdot)$ e no interior do espaço de parâmetros (assumido como um conjunto aberto), então necessariamente o Hessian $(1/N) \sum_i H(w_i, \hat{\theta}_n)$é semidefinido positivo. Se não então$\hat{\theta}_n$não é um mínimo local. Isso decorre das condições de otimização de segunda ordem - localmente$\sum_i q(w_i, \cdot)$ não deve diminuir em nenhuma direção longe do $\hat{\theta}_n$.

Uma fonte da confusão pode ser a definição "funcional" de um estimador-M. Embora, em princípio, um estimador M deva ser definido como$\arg\min_\theta \sum_i q(w_i, \theta)$, também pode ser definido como uma solução para a equação

Na prática, mesmo um hessiano definido positivo, quase singular ou mal condicionado, sugeriria que o estimador é ruim e você precisa se preocupar mais do que estimar sua variação.

Tem havido muita discussão nesse tópico sobre se o Hessian deve ser positivo (semi) definido no mínimo local. Então, farei uma declaração clara sobre isso.

Presumindo que a função objetivo e todas as funções restritivas sejam duas vezes continuamente diferenciáveis, então, em qualquer mínimo local, o hessiano do lagrangiano projetado no espaço nulo do jacobiano de restrições ativas deve ser semidefinido positivo. Ou seja, se$Z$ é uma base para o espaço nulo do jacobiano de restrições ativas, então $Z^T*(\text{Hessian of Lagrangian})*Z$deve ser positivo semidefinido. Isso deve ser definido positivamente para um mínimo local estrito.

Portanto, o hessiano da função objetivo em um problema restrito com restrição (s) ativa (s) não precisa ser positivo semidefinido se houver restrições ativas.

Notas:

1) As restrições ativas consistem em todas as restrições de igualdade, mais as restrições de desigualdade que são satisfeitas com a igualdade.

2) Veja a definição do Lagrangiano em https://www.encyclopediaofmath.org/index.php/Karush-Kuhn-Tucker_conditions .

3) Se todas as restrições são lineares, então o Hessiano da Lagrangiana = Hessiano da função objetivo, porque as segundas derivadas das funções lineares são zero. Mas você ainda precisa fazer o jazz de projeção se alguma dessas restrições estiver ativa. Observe que restrições de limite inferior ou superior são casos particulares de restrições de desigualdade linear. Se as únicas restrições ativas são restrições vinculadas, a projeção do Hessian no espaço nulo do Jacobiano de restrições ativas equivale a eliminar as linhas e colunas do Hessian correspondentes a esses componentes em seus limites.

4) Como os multiplicadores de Lagrange de restrições inativas são zero, se não houver restrições ativas, o Hessian do Lagrangiano = o Hessian da função objetivo e a matriz Identity são a base do espaço nulo do Jacobiano de restrições ativas, que resulta na simplificação do critério, sendo a condição familiar de que o hessiano da função objetivo seja semidefinido positivo em um mínimo local (positivo positivo se for um mínimo local estrito).

As respostas positivas acima são verdadeiras, mas deixam de fora a suposição crucial de identificação - se o seu modelo não for identificado (ou apenas definido), você poderá, como Wooldridge indicou corretamente, encontrar-se com um Hessian empírico não-PSD. Basta executar um modelo psicométrico / econométrico que não seja um brinquedo e ver por si mesmo.