Estimação da Média, Método Generalizado dos Momentos e Cramer Rao
Faz muito tempo que eu não escrevo aqui. O que o primeiro ano do doutorado não faz! Eu trago um post que apresenta um problema que é extremamente simples e com uma solução não óbvia - pelo menos para mim. Eu vou apresentar o problema, discutir o método generalizado dos momentos (Generalized Method of Moments) e resolver o problema. Na sequência, eu vou discutir uma extensão do problema, e usar o limite inferior de Crámer Rao.
O modelo é bem simples: nós temos duas variáveis aleatórias, \(X\) e \(Y\). \(Y\) tem média 0 e \(X\) tem uma média desconhecida, \(\mu\). Nós também sabemos que \(X\) e \(Y\) tem variâncias finitas, \(\sigma_x^2\) e \(\sigma_y^2\), respectivamente; e covariância \(\sigma_{xy}\). Nós queremos estimar a média de \(X\), \(\mu\).
O estimador mais natural é a média amostral, \(\bar{X} = \frac{1}{n} \sum_{i=1}^n X_i\), que tem variância \(\sigma_x^2/n\). A pergunta aqui é: existe um estimador melhor que esse? A intuição é simples, e assuma que a covariância entre \(X\) e \(Y\) é positiva: sabemos que a média de \(Y\) é zero. Suponha que nós observamos uma amostra e calculamos a média amostral de \(Y\). Se a média é maior que zero, nós podemos usar a correlação positiva para afirmar que a nossa média amostral de \(X\) é maior do que o valor esperado de \(X\) e obter um estimador que é mais preciso que a média amostral de \(X\) usando a informação obtida com a média amostral de \(Y\).
Vamos formalizar isso usando o método generalizados dos momentos.
O Método Generalizado dos Momentos
Vamos começar com o que a gente já (?) sabe: o método dos momentos é uma técnica de estimação que nos diz para pegar qualquer momento teórico e substituir ele por um momento amostral. Você quer estimar a média, \(E[X]\)? O método dos momentos nos diz para usar a média amostral, \(\frac{1}{n} \sum_{i=1}^n X_i\). E a variância, \(E[X^2] - E[X]^2\)? Bom, use \(\frac{1}{n} \sum_{i=1}^n X_i^2 - \bar{X}^2\) - lembre que \(\bar{X} = \frac{1}{n}\sum_{i=1}^n X_i\).
Para motivar o método generalizado dos momentos, vamos considerar que \(X\) tem uma distribuição Poisson, que tem um único parâmetro \(\lambda\). A coisa legal da Poisson é que a média e a variância são iguais a \(\lambda\). Você pode se sentir tentado a usar a média ou a variância para estimar o parâmetro. Mas usar um ou outro parece jogar informação fora e talvez uma combinação dos dois seja o ideal?
O método generalizado dos momentos faz exatamente isso. Para combinar os momentos, nós minimizamos uma forma quadrática. Me mantendo no exemplo da Poisson, nós poderíamos fazer:
\[ \min_{\lambda} \left(\frac{1}{n} \sum_{i=1}^n X_i - \lambda\right)^2 + \left(\frac{1}{n} \sum_{i=1}^n (X_i - \lambda)^2 - \lambda\right)^2 \] Eu posso ser mais geral que isso e colocar alguns pesos! Eu posso colocar até pesos que fazem o termo entre os dois (ou mais) momentos! Eu vou escrever isso em forma matricial: seja \(g(X,\theta)\) a nossa função dos momentos e \(W\) uma matriz de pesos (ela tem que ser positiva definida - isso vai garantir que a nossa função quadrática tem um mínimo e não um máximo). O que o método generalizado dos momentos faz é:
\[ \min_{\theta} g(X,\theta)^T W g(X,\theta) \] O exemplo com os momentos da Poisson usa \(W = I\) e:
\[ g(X,\theta) = \begin{bmatrix} \frac{1}{n} \sum_{i=1}^n X_i - \lambda\\ \frac{1}{n} \sum_{i=1}^n (X_i - \lambda)^2 - \lambda \end{bmatrix} \]
O método generalizado dos momentos é muito próximo de variáveis instrumentais: variáveis instrumentais impõe uma condição de momento, \(E[Zu] = 0\). Se nós temos mais variáveis instrumentais do que regressores, nós temos sobreidentificação. A solução usual é mínimos quadrados em dois estágios, que é o método generalizado dos momentos \(G = Z^T(Y - X\beta)\) e \(W = (Z^T Z)^{-1}\).
De volta ao problema
Nosso problema é achar um estimador para a média de \(X\) quando nós temos uma amostra iid \(\{X_i,Y_i\}_{i=1}^n\) tal que a média de \(Y_i\) é 0 e nós variâncias finitas e covariância entre \(X\) e \(Y\).
Nós vamos usar o método generalizado dos momentos. Os momentos que vamos usar são as médias amostrais, então nós temos:
\[ \min_{\mu} \begin{bmatrix} \bar{X} - \mu & \bar{Y} \end{bmatrix} W \begin{bmatrix} \bar{X} - \mu \\ \bar{Y} \end{bmatrix} \tag{1} \]
E os pesos? Se eu usar \(W = I\), o problema é resolvido com \(\hat{\mu} = \bar{X}\). A matriz “ótima” de pesos - significando que minimiza a variância do estimador - é a inversa da variância dos momentos. O bom desse exemplo é que a variância dos momentos é fácil de calcular: a variância da média amostral é só a variância da variável aleatória dividida pelo tamanho da amostra. Você obtém a covarância usando a mesma ideia (e o fato de que a amostra é iid):
\[ W = \begin{bmatrix} \frac{\sigma_x^2}{n} & \frac{\sigma_{xy}}{n}\\ \frac{\sigma_{xy}}{n} & \frac{\sigma_y^2}{n} \end{bmatrix}^{-1} \]
Matriz \(2 \times 2\) é inversível na mão e é fácil o suficiente que nem eu sou capaz de errar:
\[ W = \frac{n}{\sigma_x^2 \sigma_y^2 - \sigma_{xy}^2} \begin{bmatrix} \frac{\sigma_{y}^2}{n} & -\frac{\sigma_{xy}}{n} \\ - \frac{\sigma_{xy}}{n} & \frac{\sigma_x^2}{n} \end{bmatrix} \tag{2} \]
O determinante na divisão vai multiplicar a expressão toda, e como ele é positivo, não faz nenhuma diferença pro problema em mãos. Nós iremos minimizar:
\[ \min_{\mu} \sigma_y^2(\bar{X} - \mu)^2 - 2\sigma_{xy} (\bar{X} - \mu)\bar{Y} + \sigma_x^2 \bar{Y}^2 \]
Isso é obtido usando a equação 2 na equação 1 e fazendo todas as multiplicações de matriz. Agora, resta a única coisa que um economista precisa saber: derivar e igualar a zero!:
\[ -2\sigma_y^2 (\bar{X} - \mu) + 2\sigma_{xy} \bar{Y} = 0 \therefore \sigma_y^2 (\bar{X} - \mu) - \sigma_{xy} \bar{Y}\\ \hat{\mu} = \bar{X} - \frac{\sigma_{xy}}{\sigma_y^2} \bar{Y} \]
Verificando a nossa intuição lá em cima: suponha que a covariância entre \(X\) e \(Y\) é positiva. Se a média amostral de \(Y\) for maior que zero, então a média amostral de \(X\) vai ser acima do valor esperado \(\mu\). Sem muito esforço, a gente pode verificar que este estimador tem uma variância mais baixa que a média amostral:
\[ \text{Var}(\hat{\mu}) = \text{Var}\left(\bar{X} - \frac{\sigma_{xy}}{\sigma_y^2} \bar{Y}\right) = \text{Var}(\bar{X}) + \frac{\sigma_{xy}^2}{(\sigma_y^2)^2} \text{Var}(\bar{Y}) - 2\frac{\sigma_{xy}}{\sigma_y^2} \text{Cov}(\bar{X},\bar{Y}) = \\ = \frac{\sigma_{x}^2}{n} + \frac{\sigma_{xy}^2}{(\sigma_y^2)^{\cancel{2}}} \frac{\cancel{\sigma_y^2}}{n} - 2\frac{\sigma_{xy}}{\sigma_y^2}\sigma_{xy} = \frac{\sigma_x^2}{n} - \frac{\sigma_{xy}^2}{n\sigma_y^2} \]
É claro que \(\frac{\sigma_{xy}^2}{\sigma_y^2}\) é positivo, então a gente está subtraindo alguma coisa da variância da média amostral de \(X\) e o nosso estimador tem uma variância menor que a variância da média amostral.
A estrutura que eu impus no modelo é mínima. O fascinante aqui é que ter uma variável com média conhecida que é correlacionada com a variável cuja média nós desconhecemos permite obter uma estimativa melhor do que a média amostral.
Esse estimador é não viesado:
\[ E[\hat{\mu}] = E[\bar{X} - \sigma_{xy}/\sigma_y^2 \bar{Y}] = \mu - \frac{\sigma_{xy}}{\sigma_y^2} 0 = \mu \]
Limite de Crámer Rao
Se você teve um excelente curso de estatística, você já ouviu falar do limite de Crámer Rao. A ideia é bem simples: entre todos os estimadores não viesados, qual alcança a menor variância? O limite de Cramer Rao nos diz que, se \(\hat{\theta}\) é um estimador não viesado de \(\theta\) (para uma amostra i.i.d.) e \(f(X|\theta)\) é a densidade de \(X\):
\[ \text{Var}(\hat{\theta}) \geq \left(nE\left[\left(\frac{\partial \log(f(X|\theta))}{\partial \theta}\right)^2\right]\right)^{-1} \]
Por exemplo, pro caso da normal univariada com média \(\mu\) e \(\sigma^2\), nenhum estimador (não viesado) consegue obter uma variância menor que \(\frac{\sigma^2}{n}\). Isso abarca o exemplo acima, mas o univariado está fazendo todo o trabalho aqui.
Eu não impus nenhuma distribuição no exemplo acima: eu só impus uma média finita e segundo momento finito. Será que existe algum estimador melhor do que o do método generalizado dos momentos no caso bivariado? Crámer Rao pode nos dar uma resposta, mas como depende da distribuição, a gente vai ter que impor uma distribuição. A mais conveniente é a normal. A normal bivariada tem densidade:
\[ f(X,Y|\theta) = \frac{1}{\sqrt{2\pi \sigma_x^2 \sigma_y^2 (1 - \rho^2)}} \exp \left(-\frac{(X - \mu)^2}{2\sigma_x^2 (1 - \rho^2)} + \frac{\rho(X-\mu)Y}{\sigma_x \sigma_y (1 - \rho^2)} - \frac{Y^2}{2 \sigma_y^2(1-\rho^2)}\right) \]
Eu to usando \(\theta = (\mu,\sigma_{x}^2, \sigma_y^2, \rho)\) e \(\rho\) é a correlação entre \(X\) e \(Y\), logo \(\rho = \frac{\sigma_{xy}}{\sigma_x \sigma_y}\). Vamos tirar o log:
\[ \log(f(X,Y|\theta)) = -\frac{1}{2}\log(2\pi \sigma_x^2 \sigma_y^2 (1 - \rho^2)) -\frac{(X - \mu)^2}{2\sigma_x^2 (1 - \rho^2)} + \frac{\rho(X-\mu)Y}{\sigma_x \sigma_y (1 - \rho^2)} - \frac{Y^2}{2 \sigma_y^2(1-\rho^2)} \]
Vamos tirar a primeira derivada:
\[ \frac{\partial \log(f(X,Y|\theta))}{\partial \mu} = \frac{X - \mu}{\sigma_x^2 (1 - \rho^2)} - \frac{\rho Y}{\sigma_x \sigma_y (1 - \rho^2)} \]
Eu poderia calcular o quadrado disso e obter uma expressão que depende de \((X-\mu)^2\), \(Y^2\) e \((X - \mu)Y\). Honestamente, isso dá muito trabalho! Eu digo, sem fornecer nenhuma prova, que:
\[ E\left[\left(\frac{\partial \log(f(X|\theta))}{\partial \theta}\right)^2\right] = - E\left[\frac{\partial^2 \log(f(X|\theta))}{\partial \theta \partial \theta^T}\right] \]
Obtendo a segunda derivada:
\[ \frac{\partial^2 \log(f(X,Y|\theta))}{\partial \mu^2} = -\frac{1}{\sigma_x^2(1-\rho^2)} \]
Então, o limite inferior de Cramer Rao é:
\[ \left(-nE\left[\frac{\partial^2 f(X|\theta)}{\partial \mu^2}\right]\right)^{-1} = \frac{\sigma_x^2 (1 - \rho^2)}{n} \]
Agora, use a definição de \(\rho = \frac{\sigma_{xy}}{\sigma_x \sigma_y}\), para obter:
\[ \left(-nE\left[\frac{\partial^2 f(X|\theta)}{\partial \mu^2}\right]\right)^{-1} = \frac{\sigma_x^2 (1 - \rho^2)}{n} = \frac{\sigma_x^2}{n} - \frac{1}{n} \frac{\sigma_x^2 \sigma_{xy}^2}{\sigma_x^2 \sigma_y^2} = \\ = \frac{\sigma_x^2}{n} - \frac{\sigma_{xy}^2}{n\sigma_y^2} \]
Essa é exatamente a variância obtida usando o método dos momentos. O nosso estimador obtido pelo método generalizado dos momentos é o mais eficiente se a distribuição conjunta das variáveis é normal.