Aula 1 — Função de Expectativa Condicional e Projeção Linear
A Aula 1 tem uma função muito importante dentro do curso: separar três ideias que muitas vezes aparecem misturadas na prática aplicada: predição, projeção e causalidade. A Função de Expectativa Condicional, ou CEF, descreve a média de \(Y\) para cada valor de \(X\). O Melhor Preditor Linear, ou BLP, é a melhor aproximação linear de \(Y\) usando \(X\), mesmo quando a verdadeira CEF não é linear. Já a interpretação causal exige hipóteses adicionais, como a CIA, e não decorre automaticamente de uma regressão.
O fio condutor desta aula é: primeiro entendemos o que significa prever \(Y\) a partir de \(X\); depois vemos o que acontece quando restringimos a previsão a uma forma linear; por fim, discutimos por que uma relação preditiva ou associativa não é, por si só, uma relação causal.
1.1 Intuição: o que é a CEF?
A Função de Expectativa Condicional é definida por
\[m(x) \equiv \mathbb{E}[Y \mid X=x].\]
Ela responde à pergunta: qual é o valor médio esperado de \(Y\) quando observamos \(X=x\)? Portanto, a CEF não descreve cada observação individual, mas sim o comportamento médio de \(Y\) dentro de cada valor, ou faixa de valores, de \(X\).
A CEF é central porque é o melhor preditor de \(Y\) dado \(X\) quando usamos erro quadrático médio como critério. Isto é, entre todos os preditores possíveis da forma \(g(X)\), a função que minimiza
\[\mathbb{E}[(Y-g(X))^2]\]
é justamente
\[g(X)=\mathbb{E}[Y\mid X].\]
Assim, a CEF é o objeto populacional que concentra toda a informação preditiva de \(X\) sobre a média de \(Y\). Se definirmos o erro da CEF como
\[e=Y-m(X),\]
então vale
\[\mathbb{E}[e\mid X]=0.\]
Essa propriedade significa que, depois de usar a melhor previsão condicional possível, não sobra erro sistemático previsível por \(X\). O erro pode continuar existindo, mas sua média condicional é zero. O slide da Aula 1 destaca justamente essa propriedade: o erro da CEF é ortogonal, em média, a qualquer função integrável de \(X\).
Vamos verificar essa ideia empiricamente. Simulamos dados em que a CEF verdadeira é conhecida:
set.seed(42)n <-1000X <-runif(n, -3, 3)Y <- X^2+rnorm(n)dados <-data.frame(X = X, Y = Y)dados$faixa <-cut(X, breaks =20)cef_empirica <- dados |>group_by(faixa) |>summarise(X_medio =mean(X),Y_medio =mean(Y),.groups ="drop" )ggplot(dados, aes(x = X, y = Y)) +geom_point(alpha =0.2, color ="steelblue") +geom_line(data = cef_empirica,aes(x = X_medio, y = Y_medio),color ="red", linewidth =1.2) +stat_function(fun =function(x) x^2,color ="black", linetype ="dashed", linewidth =1) +labs(title ="CEF empírica (vermelho) vs. CEF verdadeira (preto tracejado)",x ="X", y ="Y") +theme_minimal()
CEF verdadeira vs. médias condicionais empíricas
A linha vermelha representa médias empíricas de \(Y\) dentro de faixas de \(X\). Ela é uma aproximação amostral da CEF. A linha preta tracejada é a CEF verdadeira, \(m(x)=x^2\). Como a amostra é grande, as médias por faixa se aproximam bem da curva verdadeira.
1.2 CEF vs BLP: o que acontece quando a CEF não é linear?
A CEF pode ter qualquer formato: linear, quadrático, descontínuo, com interações, com efeitos heterogêneos etc. O Melhor Preditor Linear impõe uma restrição adicional: ele procura a melhor aproximação de \(Y\) entre funções lineares de \(X\).
A diferença central é: a CEF é o melhor preditor geral, enquanto o BLP é o melhor preditor dentro da classe linear. Se a CEF for linear, o BLP recupera a CEF. Mas se a CEF for não linear, o BLP apenas fornece a melhor reta no sentido de erro quadrático médio.
blp <-lm(Y ~ X, data = dados)ggplot(dados, aes(x = X, y = Y)) +geom_point(alpha =0.15, color ="steelblue") +stat_function(fun =function(x) x^2,color ="black", linetype ="dashed", linewidth =1) +geom_abline(intercept =coef(blp)[1], slope =coef(blp)[2],color ="blue", linewidth =1.2) +labs(title ="BLP linear (azul) falha em capturar a CEF quadrática (tracejado)",x ="X", y ="Y") +theme_minimal()
BLP (azul) vs CEF verdadeira (preto tracejado)
cat("Coeficientes do BLP:\n")
Coeficientes do BLP:
cat("Intercepto:", round(coef(blp)[1], 4), "\n")
Intercepto: 3.027
cat("Inclinação:", round(coef(blp)[2], 4), "\n")
Inclinação: -0.0128
cat("\nNota: inclinação próxima de zero — X~U(-3,3) é simétrico,","então E[XY] = E[X·X²] = E[X³] = 0 por simetria.\n")
Nota: inclinação próxima de zero — X~U(-3,3) é simétrico, então E[XY] = E[X·X²] = E[X³] = 0 por simetria.
Neste exemplo, a inclinação estimada fica próxima de zero. Isso não significa que \(X\) não tenha relação com \(Y\). Na verdade, a relação é fortíssima: \(Y\) cresce quando \(|X|\) cresce. O problema é que essa relação é quadrática e simétrica, enquanto o BLP procura uma reta.
Como \(X\sim U(-3,3)\) é simétrico, temos aproximadamente
Por isso a melhor reta é praticamente horizontal. Esse exemplo é fundamental: uma inclinação linear perto de zero não implica ausência de relação; pode significar apenas que a forma linear é inadequada.
Observação importante:\(\mathbb{E}[Xe]=0\) é uma propriedade da projeção linear. Ela vale por construção para o erro do BLP. Isso não deve ser confundido com exogeneidade causal.
Aqui, a regressão linear tem uma interpretação forte: ela descreve a própria média condicional de \(Y\). Já quando a CEF não é linear, ainda podemos calcular o BLP, mas ele deixa de ser a CEF e passa a ser apenas a melhor aproximação linear.
Essa distinção é uma das mais importantes da primeira aula:
\[\text{CEF linear} \Rightarrow \text{BLP coincide com a CEF},\]
mas
\[\text{BLP existe} \nRightarrow \text{CEF é linear}.\]
1.3 Viés de Variável Omitida (OVB)
O viés de variável omitida aparece quando estimamos uma regressão curta, omitindo um regressor relevante que está correlacionado com os regressores incluídos.
cat("\n→ Os dois coincidem: a fórmula funciona exatamente.\n")
→ Os dois coincidem: a fórmula funciona exatamente.
1.4 CIA e identificação causal
A Aula 1 termina com uma mensagem essencial: regressão, CEF e BLP são objetos de predição e associação. Para transformar diferenças condicionais em efeitos causais, precisamos de hipóteses adicionais.
Considere um tratamento binário \(T\in\{0,1\}\). O resultado observado é
cat("Diferença de médias ingênua: ", round(diff_ingenua, 4), "\n")
Diferença de médias ingênua: 1.5798
cat("Viés de seleção: ", round(diff_ingenua - ace_verdadeiro, 4), "\n")
Viés de seleção: 1.0798
cat("\n→ A diferença de médias superestima o efeito causal.\n")
→ A diferença de médias superestima o efeito causal.
cat(" Quem faz mestrado já teria salário maior (U alto) mesmo sem mestrado.\n")
Quem faz mestrado já teria salário maior (U alto) mesmo sem mestrado.
A simulação mostra um caso de seleção positiva. Pessoas com maior \(U\), interpretado como habilidade não observada, têm maior probabilidade de fazer mestrado. Como \(U\) também aumenta o salário, a diferença de médias entre quem fez e quem não fez mestrado mistura o efeito causal verdadeiro com a diferença prévia de habilidade entre os grupos.
Aula 2 — Álgebra do MQO: Matrizes P, M, FWL e Alavancagem
A Aula 2 muda o foco da população para a amostra. Na Aula 1, o objeto central era o BLP populacional,
cat("\n→ Idênticos. O lm() usa decomposição QR (mais estável),","mas o resultado é o mesmo.\n")
→ Idênticos. O lm() usa decomposição QR (mais estável), mas o resultado é o mesmo.
2.3 Matrizes \(P\) e \(M\)
A matriz de projeção \(P=X(X^\top X)^{-1}X^\top\) transforma \(Y\) em sua projeção no espaço coluna de \(X\): \(\hat Y=PY\). A matriz residualizadora \(M=I_n-P\) remove essa parte projetada: \(\hat e=MY\).
Como \(PX=X\), temos \(MX=(I-P)X=0\), o que implica a ortogonalidade dos resíduos: \(X^\top\hat e=0\).
P <- X %*%solve(XtX) %*%t(X)M <-diag(n) - PY_hat <- P %*% Ye_hat <- M %*% Ycat("=== Propriedades de P ===\n")
No modelo \(Y=X_1\beta_1+X_2\beta_2+e\), o coeficiente de \(X_1\) mede a associação entre \(Y\) e a parte de \(X_1\) que não é explicada por \(X_2\). O procedimento é:
remover de \(Y\) a parte explicada por \(X_2\);
remover de \(X_1\) a parte explicada por \(X_2\);
regressar os resíduos de \(Y\) nos resíduos de \(X_1\).
O coeficiente obtido é exatamente o mesmo da regressão múltipla original.
cat("max(h_ii):", round(max(h_ii), 4), " | deve ser <= 1\n")
max(h_ii): 0.1787 | deve ser <= 1
cat("sum(h_ii):", round(sum(h_ii), 4), " | deve ser k =", k, "\n\n")
sum(h_ii): 3 | deve ser k = 3
data.frame(h = h_ii, e =as.numeric(e_hat)) |>ggplot(aes(x = h, y = e)) +geom_point(aes(size =abs(h_ii * e_loo), color =abs(h_ii * e_loo)),alpha =0.7) +scale_color_gradient(low ="steelblue", high ="red") +geom_hline(yintercept =0, linetype ="dashed") +geom_vline(xintercept =2*k/n, linetype ="dashed", color ="gray50") +labs(title ="Alavancagem × Resíduo",subtitle ="Tamanho/cor = influência |h_ii × ẽ_i| | linha = limiar 2k/n",x ="Alavancagem (h_ii)", y ="Resíduo (ê_i)",color ="Influência", size ="Influência") +theme_minimal()
Alavancagem vs resíduo — identificando observações influentes
No gráfico, pontos mais à direita têm maior alavancagem. Pontos com bolhas maiores e mais avermelhadas são mais influentes segundo a medida \(|h_{ii}\tilde e_i|\). A influência nasce da combinação de alta alavancagem com erro leave-one-out grande.
Fechamento da Aula 2
A Aula 2 mostra que o MQO é uma projeção ortogonal em \(\mathbb{R}^n\). A matriz \(P\) projeta \(Y\) no espaço coluna de \(X\), gerando \(\hat Y\), enquanto \(M=I-P\) gera os resíduos \(\hat e\). As equações normais \(X^\top\hat e=0\) dizem que os resíduos são ortogonais a todos os regressores. O Teorema FWL mostra que coeficientes de regressão múltipla podem ser entendidos como associações parciais entre resíduos. Por fim, alavancagem e leave-one-out permitem estudar observações potencialmente influentes e antecipam temas posteriores como HC3, jackknife e validação cruzada.
Aula 3 — Propriedades em Amostra Finita e Erros-Padrão Robustos
A Aula 3 estuda o comportamento do MQO em amostras finitas, isto é, sem depender ainda de argumentos assintóticos. A pergunta central é: dado um desenho amostral \(X\), o que conseguimos dizer sobre \(\hat\beta\)?
A aula parte de três hipóteses fundamentais:
amostragem aleatória/i.i.d.;
exogeneidade em média, \(\mathbb{E}[e\mid X]=0\);
homocedasticidade, quando aplicável, \(\mathbb{E}[e^2\mid X]=\sigma^2\).
Sob essas condições, o MQO é não viesado condicionalmente, sua variância pode ser escrita explicitamente e, sob homocedasticidade, ele é BLUE pelo Teorema de Gauss-Markov. Quando há heterocedasticidade, o MQO continua não viesado sob exogeneidade, mas a matriz de variância clássica deixa de ser válida; por isso entram os estimadores robustos HC0–HC3.
3.1 Setup
Vamos simular um modelo linear em que a verdade é conhecida:
\[Y=X\beta+e, \qquad \beta=(1,2,-1)^\top.\]
Neste primeiro desenho, os erros são homocedásticos: \(e_i\sim N(0,1)\).
library(ggplot2)set.seed(123)n <-200k <-3X1 <-rnorm(n)X2 <-0.5* X1 +rnorm(n)X <-cbind(1, X1, X2)beta_true <-c(1, 2, -1)e_true <-rnorm(n)Y <- X %*% beta_true + e_trueXtX <-t(X) %*% Xbeta_hat <-solve(XtX) %*%t(X) %*% Ye_hat <- Y - X %*% beta_hatP <- X %*%solve(XtX) %*%t(X)h_ii <-diag(P)
3.2 Não viés do MQO: verificação por simulação
O Teorema C1 afirma que, sob amostragem i.i.d. e exogeneidade em média,
\[\mathbb{E}[\hat\beta\mid X]=\beta.\]
Como
\[\hat\beta = (X^\top X)^{-1}X^\top Y = \beta+(X^\top X)^{-1}X^\top e,\]
Na simulação abaixo, mantemos \(X\) fixo e geramos muitos vetores de erro \(e\). Assim, estamos verificando a ideia de média condicional: ao repetir o experimento muitas vezes com o mesmo desenho \(X\), a média de \(\hat\beta\) se aproxima de \(\beta\).
B <-2000beta1_sim <-numeric(B)for (b in1:B) { e_b <-rnorm(n) Y_b <- X %*% beta_true + e_b b_hat <-solve(t(X) %*% X) %*%t(X) %*% Y_b beta1_sim[b] <- b_hat[2]}cat("=== Verificação do não-viés ===\n\n")
=== Verificação do não-viés ===
cat("beta1 verdadeiro: ", beta_true[2], "\n")
beta1 verdadeiro: 2
cat("Média de β̂₁ (B=2000): ", round(mean(beta1_sim), 4), "\n")
ggplot(data.frame(b1 = beta1_sim), aes(x = b1)) +geom_histogram(aes(y =after_stat(density)),bins =50, fill ="steelblue", alpha =0.7) +geom_vline(xintercept = beta_true[2],color ="red", linewidth =1.2, linetype ="dashed") +geom_vline(xintercept =mean(beta1_sim),color ="black", linewidth =1) +labs(title ="Distribuição de β̂₁ em 2000 simulações",subtitle ="Vermelho = valor verdadeiro | Preto = média simulada",x ="β̂₁", y ="Densidade") +theme_minimal()
Distribuição de β̂₁ em 2000 amostras — centragem no valor verdadeiro
A média simulada de \(\hat\beta_1\) fica próxima de \(2\), o valor verdadeiro. Isso não significa que cada amostra individual produz exatamente \(2\). Significa que, repetindo o processo muitas vezes sob as hipóteses do modelo, o MQO fica centrado no parâmetro verdadeiro.
3.3 Variância do MQO e Teorema de Gauss-Markov
A variância condicional geral do MQO é
\[\operatorname{Var}(\hat\beta\mid X) = (X^\top X)^{-1}X^\top D X(X^\top X)^{-1},\]
onde \(D=\operatorname{Var}(e\mid X)\). No caso homocedástico, \(D=\sigma^2 I_n\), então
Essa fórmula mostra que a precisão do MQO depende de duas coisas: da variância do erro \(\sigma^2\) e da geometria dos regressores, capturada por \((X^\top X)^{-1}\).
O Teorema de Gauss-Markov diz que, sob homocedasticidade, o MQO é o melhor estimador linear não viesado: nenhum outro estimador linear e não viesado tem matriz de variância menor. Nos slides, isso aparece como o resultado BLUE.
O problema é que os resíduos MQO são menores, em média, que os erros verdadeiros, porque a regressão escolhe \(\hat\beta\) justamente para minimizar a soma dos quadrados. Sob homocedasticidade,
Logo, \(\hat\sigma^2\) é viesado para baixo. A correção clássica é
\[s^2=\frac{1}{n-k}\sum_{i=1}^n \hat e_i^2,\]
que satisfaz \(\mathbb{E}[s^2\mid X]=\sigma^2\). Essa é a origem dos graus de liberdade residuais \(n-k\): a regressão estimou \(k\) parâmetros antes de calcular a soma de quadrados residual, consumindo \(k\) graus de liberdade.
cat("Estimador com divisor n: ", round(sigma2_hat_n, 4), "\n")
Estimador com divisor n: 0.9228
cat("Estimador com divisor n-k (s²): ", round(sigma2_hat_df, 4), "\n")
Estimador com divisor n-k (s²): 0.9369
cat("\nFator esperado de viés com divisor n: (n-k)/n =",round((n-k)/n, 4), "\n")
Fator esperado de viés com divisor n: (n-k)/n = 0.985
3.5 Estimadores HC0 a HC3 — construídos à mão
Quando há heterocedasticidade, a matriz correta de variância condicional do MQO é
\[V_{\hat\beta} = (X^\top X)^{-1}(X^\top D X)(X^\top X)^{-1},\]
onde \(D=\operatorname{diag}(\sigma_1^2,\ldots,\sigma_n^2)\). Como \(D\) é desconhecida, os estimadores HC substituem os erros verdadeiros \(e_i^2\) por funções dos resíduos \(\hat e_i^2\):
for (i in1:4) { tipo <-c("HC0","HC1","HC2","HC3")[i] man <-c(se_HC0, se_HC1, se_HC2, se_HC3)[i] pkg <- se_pkg[i]cat(sprintf("%-6s %8.6f %8.6f %8.2e\n", tipo, man, pkg, man - pkg))}
e portanto, para uma mesma observação, \(\text{HC0} \le \text{HC2} \le \text{HC3}\). Essa desigualdade é especialmente relevante quando \(h_{ii}\) é grande. O HC3 aumenta bastante o peso de observações com alta alavancagem, por isso costuma ser mais conservador em amostras finitas.
pesos <-data.frame(obs =1:n,HC0 =as.numeric(e_hat)^2,HC2 = w_HC2,HC3 = w_HC3,h_ii = h_ii)pesos_long <- tidyr::pivot_longer(pesos, cols =c(HC0, HC2, HC3),names_to ="tipo", values_to ="peso")ggplot(pesos_long, aes(x = h_ii, y = peso, color = tipo)) +geom_point(alpha =0.5) +facet_wrap(~tipo, scales ="free_y") +labs(title ="Pesos dos estimadores HC em função da alavancagem",subtitle ="HC2 e HC3 inflam os pesos para obs. com alta alavancagem",x ="Alavancagem (h_ii)", y ="Peso") +theme_minimal() +theme(legend.position ="none")
Pesos dos estimadores HC em função da alavancagem
3.7 Simulação com heterocedasticidade: erro-padrão clássico vs. robusto
Agora geramos um desenho em que a variância do erro depende de \(X_1\):
de modo que \(\operatorname{Var}(e_i\mid X_i)=(1+|X_{1i}|)^2\). O MQO continua centrado em \(\beta\), porque mantemos \(\mathbb{E}[e_i\mid X_i]=0\). Porém, a fórmula homocedástica da variância deixa de ser apropriada.
set.seed(321)eps_het <-rnorm(n)e_het <- (1+abs(X1)) * eps_hetY_het <- X %*% beta_true + e_hetmod_het <-lm(as.numeric(Y_het) ~ X1 + X2)se_classico <-sqrt(vcov(mod_het)[2, 2])se_hc1 <-sqrt(vcovHC(mod_het, type ="HC1")[2, 2])se_hc2 <-sqrt(vcovHC(mod_het, type ="HC2")[2, 2])se_hc3 <-sqrt(vcovHC(mod_het, type ="HC3")[2, 2])cat("=== Erros-padrão sob heterocedasticidade ===\n\n")
Quando há heterocedasticidade, o erro-padrão clássico pode ficar subestimado ou superestimado, dependendo do desenho. Os erros-padrão robustos estimam diretamente
\[V_{\hat\beta} = (X^\top X)^{-1}(X^\top D X)(X^\top X)^{-1}\]
sem impor \(D=\sigma^2 I_n\).
3.8 Multicolinearidade: quando \(X^\top X\) quase não é invertível
Na multicolinearidade estrita, existe \(\alpha\neq 0\) tal que \(X\alpha=0\). Nesse caso, \(X^\top X\) é singular e \((X^\top X)^{-1}\) não existe. O MQO não é definido de forma única.
Na quase multicolinearidade, \(X^\top X\) ainda é invertível, mas é mal condicionado. Sob homocedasticidade,
Quando \((X^\top X)^{-1}\) tem entradas muito grandes, os erros-padrão também crescem. A quase multicolinearidade afeta a precisão do MQO exclusivamente por meio de \((X^\top X)^{-1}\).
Mesmo que os coeficientes verdadeiros estejam bem definidos, a amostra tem dificuldade em separar o efeito individual de \(X_1\) e \(X_2\) quando eles carregam quase a mesma informação. O resultado típico é: coeficientes instáveis, erros-padrão grandes e baixa precisão individual.
3.9 Amostragem por clusters e dependência intracluster
Em muitos bancos de dados, as observações vêm agrupadas em clusters: alunos em escolas, trabalhadores em firmas, indivíduos em municípios. A hipótese apropriada é independência entre clusters, mas dependência dentro do cluster.
Para o cluster \(g\), escrevemos \(Y_g=X_g\beta+e_g\) e permitimos \(\Sigma_g=\mathbb{E}[e_ge_g^\top\mid X_g]\). O MQO continua sendo \((X^\top X)^{-1}X^\top Y\), mas a matriz de variância muda para
A mensagem central: com clusters, o tamanho efetivo da amostra é mais próximo de \(G\), o número de clusters, do que de \(n\), o número total de observações.
cat("EP clássico de beta1: ", round(se_classico_cl, 4), "\n")
EP clássico de beta1: 0.0755
cat("EP cluster-robusto de beta1: ", round(se_cluster_cl, 4), "\n")
EP cluster-robusto de beta1: 0.0603
Como os erros compartilham um componente comum dentro de cada cluster, tratar todas as observações como independentes tende a subestimar a incerteza. O erro-padrão cluster-robusto corrige a matriz de variância para permitir dependência arbitrária dentro dos clusters.
Fechamento da Aula 3
A Aula 3 mostra que, em amostra finita, o MQO tem propriedades muito fortes sob exogeneidade em média. Condicionalmente em \(X\), ele é não viesado:
\[\mathbb{E}[\hat\beta\mid X]=\beta.\]
Sua variância geral é \((X^\top X)^{-1}X^\top D X(X^\top X)^{-1}\) e, sob homocedasticidade, simplifica para \(\sigma^2(X^\top X)^{-1}\). O Teorema de Gauss-Markov garante que o MQO é BLUE sob homocedasticidade. Quando ela falha, entram os estimadores robustos HC0–HC3. A aula também mostra que \(\sigma^2\) precisa da correção \(n-k\), que multicolinearidade prejudica precisão via \((X^\top X)^{-1}\), e que em dados agrupados a inferência precisa respeitar a dependência intracluster.
Aula 4 — Regressões Normais: Inferência Exata
A Aula 4 introduz o modelo de regressão normal. Até aqui, nas aulas anteriores, conseguimos obter propriedades importantes do MQO, como não viés, variância condicional e Teorema de Gauss-Markov. Mas esses resultados ainda não davam, em geral, a distribuição exata das estatísticas de interesse.
A hipótese adicional desta aula é:
\[e\mid X\sim \mathcal{N}(0,\sigma^2 I_n).\]
Essa hipótese permite sair de resultados apenas sobre média e variância e obter distribuições exatas em amostra finita. Assim, conseguimos justificar formalmente:
A grande mensagem da Aula 4 é que, sob normalidade dos erros, a inferência clássica da regressão linear deixa de ser apenas aproximada e passa a ser exata em amostras finitas.
4.0 As distribuições da inferência exata
A normal padrão é \(Z\sim \mathcal{N}(0,1)\). Toda normal escalar pode ser escrita como transformação afim da normal padrão:
\[X=\mu+\sigma Z \quad\Rightarrow\quad X\sim \mathcal{N}(\mu,\sigma^2).\]
Na versão multivariada, se \(Z\sim \mathcal{N}(0,I_k)\), então seus componentes são independentes e normais padrão. Além disso, transformações afins preservam normalidade:
Essas propriedades são fundamentais porque o MQO é uma transformação linear de \(Y\). Logo, se \(Y\mid X\) é normal, então \(\hat\beta\mid X\) também será normal.
As distribuições \(\chi^2\), \(t\) e \(F\) aparecem a partir da normal:
Se \(Z\sim \mathcal{N}(0,1)\), \(Q\sim \chi^2_k\) e \(Z\perp Q\), então
\[\frac{Z}{\sqrt{Q/k}}\sim t_k.\]
E se \(U\sim\chi^2_q\), \(V\sim\chi^2_m\) e \(U\perp V\), então
\[\frac{U/q}{V/m}\sim F_{q,m}.\]
Essas três relações são a base matemática de toda a inferência exata da regressão normal.
4.1 Setup: modelo normal
Vamos trabalhar com o modelo \(Y=X\beta+e\), com \(e\mid X\sim\mathcal{N}(0,\sigma^2I_n)\).
A normalidade dos erros implica que, condicionalmente em \(X\), o vetor \(Y\) também é normal:
\[Y\mid X\sim \mathcal{N}(X\beta,\sigma^2I_n).\]
Essa hipótese é mais forte do que as hipóteses da Aula 3. Na Aula 3, bastava ter \(\mathbb{E}[e\mid X]=0\) e, para a fórmula clássica de variância, homocedasticidade. Agora, além disso, assumimos uma forma distributiva completa para os erros. O ganho é que podemos obter distribuições exatas em amostras finitas.
Para \(\sigma^2\) fixo, maximizar \(\ell_n(\beta,\sigma^2)\) em relação a \(\beta\) é o mesmo que minimizar \(\sum_{i=1}^n(Y_i-X_i^\top\beta)^2\). Portanto, no modelo normal, o MLE de \(\beta\) coincide com o MQO:
\[\hat\beta_{MLE}=\hat\beta_{MQO}.\]
Essa é uma conexão importante: o MQO pode ser visto tanto como uma projeção geométrica quanto como um estimador de máxima verossimilhança sob normalidade.
sigma2_mle <-as.numeric(t(e_hat) %*% e_hat) / ncat("=== MLE no modelo normal ===\n\n")
=== MLE no modelo normal ===
cat("beta_hat MQO/MLE:\n")
beta_hat MQO/MLE:
print(round(beta_hat, 4))
[,1]
0.6995
X1 2.1548
X2 -0.9570
cat("\nMLE de sigma² usa divisor n: ", round(sigma2_mle, 4), "\n")
MLE de sigma² usa divisor n: 3.3346
cat("Estimador não viesado s² usa n-k: ", round(s2, 4), "\n")
Estimador não viesado s² usa n-k: 3.5474
Note que o MLE de \(\sigma^2\) usa divisor \(n\), enquanto o estimador clássico não viesado usa divisor \(n-k\). Assim, no modelo normal, o MLE de \(\beta\) coincide com o MQO, mas o MLE de \(\sigma^2\) não coincide com \(s^2\).
4.2 Distribuição exata de \(\hat\beta\): verificação por simulação
Como \(\hat\beta = \beta+(X^\top X)^{-1}X^\top e\) e \(e\mid X\sim\mathcal{N}(0,\sigma^2I_n)\), segue que \(\hat\beta\mid X\) é uma transformação afim de um vetor normal. Logo,
\[\hat\beta\mid X \sim \mathcal{N}(\beta,\sigma^2(X^\top X)^{-1}).\]
Esse é o Teorema D.4 dos slides. A diferença em relação à Aula 3 é que lá sabíamos apenas a média e a variância de \(\hat\beta\). Agora sabemos sua distribuição completa.
x_grid <-seq(min(beta1_sim), max(beta1_sim), length.out =200)normal_teo <-dnorm(x_grid, mean = beta_true[2], sd = sd_beta1_teo)ggplot(data.frame(b1 = beta1_sim), aes(x = b1)) +geom_histogram(aes(y =after_stat(density)),bins =60, fill ="steelblue", alpha =0.6) +geom_line(data =data.frame(x = x_grid, y = normal_teo),aes(x = x, y = y), color ="red", linewidth =1.2) +labs(title ="β̂₁ ~ N(β₁, σ²[(X'X)⁻¹]₁₁) — Teorema D.4",subtitle ="Histograma simulado (azul) vs Normal teórica (vermelho)",x ="β̂₁", y ="Densidade") +theme_minimal()
Distribuição simulada de β̂₁ vs Normal teórica
O histograma simulado de \(\hat\beta_1\) coincide muito bem com a densidade normal teórica. Como estamos mantendo \(X\) fixo e simulando apenas novos erros normais, essa é uma ilustração direta da distribuição condicional \(\hat\beta\mid X\).
4.2.1 Distribuição dos resíduos e independência entre \(\hat\beta\) e \(\hat{e}\)
Da Aula 2, sabemos que \(\hat{e}=MY\) com \(M=I-X(X^\top X)^{-1}X^\top\). Como \(MX=0\), temos \(\hat{e}=Me\). Como \(e\mid X\sim\mathcal{N}(0,\sigma^2I_n)\), os resíduos também são normais:
\[\hat{e}\mid X\sim\mathcal{N}(0,\sigma^2M).\]
Além disso, no modelo normal, \(\hat\beta\perp \hat{e}\). Essa independência é crucial: a estatística \(t\) nasce justamente da combinação de uma parte que depende de \(\hat\beta\) com outra que depende de \(s^2\), que é função dos resíduos. A normalidade garante que essas partes sejam independentes.
cat("Correlação simulada entre β̂₁ e s²:", round(corr_beta1_s2, 4), "\n")
Correlação simulada entre β̂₁ e s²: 0.0147
cat("Sob normalidade, a independência teórica implica correlação zero.\n")
Sob normalidade, a independência teórica implica correlação zero.
4.3 Distribuição exata de \(s^2\): qui-quadrado
O estimador não viesado da variância do erro é \(s^2=\frac{1}{n-k}\hat{e}^\top\hat{e}\). Como \(\hat{e}=Me\) e \(M\) é simétrica, idempotente e tem posto \(n-k\), temos sob normalidade:
\[\frac{(n-k)s^2}{\sigma^2}\sim\chi^2_{n-k}.\]
Esse resultado é o Teorema D.6 e é a base para construir tanto a estatística \(t\) quanto intervalos de confiança para \(\sigma^2\).
stat_sim <- (n - k) * s2_sim / sigma2_truex_grid2 <-seq(0, max(stat_sim), length.out =300)chi2_teo <-dchisq(x_grid2, df = n - k)ggplot(data.frame(stat = stat_sim), aes(x = stat)) +geom_histogram(aes(y =after_stat(density)),bins =60, fill ="steelblue", alpha =0.6) +geom_line(data =data.frame(x = x_grid2, y = chi2_teo),aes(x = x, y = y), color ="red", linewidth =1.2) +labs(title ="(n-k)s²/σ² ~ χ²(n-k) — Teorema D.6",subtitle =paste("n-k =", n-k, "graus de liberdade"),x ="(n-k)s²/σ²", y ="Densidade") +theme_minimal()
O teste-\(F\) serve para testar hipóteses conjuntas, como \(H_0:\beta_1=\beta_2=0\). Comparamos dois modelos: o restrito (apenas intercepto) e o irrestrito (intercepto, \(X_1\) e \(X_2\)). A estatística é
\[F = \frac{(SSE_R-SSE_U)/q}{SSE_U/(n-k)},\]
onde \(SSE_R\) é a soma dos quadrados residual do modelo restrito, \(SSE_U\) é a do modelo irrestrito, \(q\) é o número de restrições e \(n-k\) são os graus de liberdade residuais do modelo irrestrito.
A intuição é simples: o numerador mede quanto o ajuste piora ao impor \(H_0\), por restrição. O denominador mede o ruído residual médio do modelo irrestrito. No modelo normal, sob \(H_0\), \(F\sim F_{q,n-k}\).
X_rest <-matrix(1, n, 1)b_rest <-solve(t(X_rest) %*% X_rest) %*%t(X_rest) %*% Ye_rest <- Y - X_rest %*% b_restSSE_rest <-as.numeric(t(e_rest) %*% e_rest)SSE_irr <-as.numeric(t(e_hat) %*% e_hat)q <-2F_obs <- ((SSE_rest - SSE_irr) / q) / (SSE_irr / (n - k))pval_F <-pf(F_obs, df1 = q, df2 = n - k, lower.tail =FALSE)cat("=== Teste F manual ===\n\n")
cat("F via anova(): ", round(anova_res$F[2], 4), "\n")
F via anova(): 36.7867
cat("p via anova(): ", round(anova_res$`Pr(>F)`[2], 6), "\n")
p via anova(): 0
cat("\n→ Idênticos. F construído à mão = F do lm().\n")
→ Idênticos. F construído à mão = F do lm().
O resultado manual coincide com anova(). Como \(\beta_1=2\) e \(\beta_2=-1\), a hipótese \(H_0:\beta_1=\beta_2=0\) é falsa; por isso o \(F\) observado é grande e o p-valor é praticamente zero.
4.6 Intervalo de confiança para \(\beta_1\)
O intervalo de confiança exato para \(\beta_j\) usa o quantil da distribuição \(t_{n-k}\):
Um intervalo de confiança de 95% não significa que há 95% de probabilidade de o parâmetro fixo estar naquele intervalo específico depois de observado. A interpretação frequentista é: se repetíssemos o experimento muitas vezes, aproximadamente 95% dos intervalos construídos conteriam o parâmetro verdadeiro. No modelo normal, usando a distribuição \(t_{n-k}\), essa cobertura é exata em amostra finita.
Como \(\operatorname{Var}(\hat\beta\mid X)=\sigma^2(X^\top X)^{-1}\), o MQO atinge o bound de Cramér–Rao para \(\beta\). Portanto, no modelo normal, \(\hat\beta\) é eficiente no sentido de Cramér–Rao: nenhum estimador não viesado de \(\beta\) pode ter matriz de variância menor.
V_beta_cr <- sigma2_true * XtX_invcat("=== Cramér-Rao para beta ===\n\n")
=== Cramér-Rao para beta ===
cat("Variância exata de beta_hat = bloco de beta em I^{-1}:\n")
Variância exata de beta_hat = bloco de beta em I^{-1}:
A Aula 4 mostra que a normalidade dos erros transforma a inferência do MQO em uma teoria exata de amostra finita. Sob \(e\mid X\sim \mathcal{N}(0,\sigma^2I_n)\), temos
Esses resultados justificam os testes \(t\), os intervalos de confiança, os p-valores e os testes \(F\) clássicos. A aula também mostra que, no modelo normal, o MQO coincide com o MLE de \(\beta\) e atinge o bound de Cramér–Rao. Portanto, nessa estrutura, o MQO é não apenas BLUE, mas também eficiente em um sentido mais forte.
Aula 5 — Normalidade Assintótica do MQO
A Aula 5 muda o regime da inferência. Na Aula 4, a distribuição exata de \(\hat\beta\) vinha da hipótese forte de normalidade dos erros:
\[e\mid X\sim \mathcal{N}(0,\sigma^2 I_n).\]
Agora, essa hipótese é removida. Em vez de exigir normalidade exata dos erros, passamos a usar resultados assintóticos: Lei dos Grandes Números, Teorema Central do Limite, Teorema do Mapeamento Contínuo, Slutsky e Método Delta.
Ou seja, o MQO é consistente e, depois de reescalado por \(\sqrt{n}\), tem distribuição aproximadamente normal para amostras grandes. A Aula 5 mostra que essa normalidade não vem mais da normalidade dos erros, mas do Teorema Central do Limite aplicado ao termo
\[\frac{1}{\sqrt{n}}\sum_{i=1}^n X_i e_i.\]
5.1 A ideia central: de exato para assintótico
Na Aula 4, a normalidade de \(\hat\beta\) veio da hipótese \(e \mid X \sim \mathcal{N}(0,\sigma^2 I_n)\). Aqui removemos essa hipótese. Os erros podem ter qualquer distribuição com média zero e variância finita. A normalidade de \(\hat\beta\) emerge assintoticamente pelo TCL — mas só vale aproximadamente para \(n\) finito.
library(ggplot2)set.seed(42)k <-3gera_dados <-function(n, dist_erro ="normal") { X1 <-rnorm(n) X2 <-0.5* X1 +rnorm(n) X <-cbind(1, X1, X2) e <-switch(dist_erro,"normal"=rnorm(n),"chi2"= (rchisq(n, df =2) -2) /2,"uniforme"=runif(n, -sqrt(3), sqrt(3)) ) Y <- X %*%c(1, 2, -1) + elist(X = X, Y = Y)}
5.2 Consistência: \(\hat\beta \to_p \beta\) quando \(n \to \infty\)
A consistência do MQO vem da ideia plug-in. O estimador pode ser escrito como
B <-3000distrs <-c("normal", "chi2", "uniforme")ns <-c(30, 200, 1000)resultados <-list()for (dist in distrs) {for (n in ns) { t_stats <-replicate(B, { d <-gera_dados(n, dist) b_hat <-solve(t(d$X) %*% d$X) %*%t(d$X) %*% d$Ysqrt(n) * (b_hat[2] - beta1_true) }) resultados[[paste(dist, n)]] <-data.frame(stat = t_stats,dist = dist,n = n ) }}df_plot <-do.call(rbind, resultados)df_plot$label <-factor(paste0("n=", df_plot$n), levels =paste0("n=", ns))ggplot(df_plot, aes(x = stat)) +geom_histogram(aes(y =after_stat(density)),bins =50, fill ="steelblue", alpha =0.6) +facet_grid(dist ~ label, scales ="free") +labs(title ="√n(β̂₁ - β₁): aproximação normal assintótica",subtitle ="A forma se aproxima de uma normal conforme n cresce",x ="√n(β̂₁ - β₁)", y ="Densidade") +theme_minimal(base_size =9)
√n(β̂₁ - β₁) para diferentes distribuições de erro
A normalidade assintótica não diz que \(\hat\beta\) é exatamente normal em amostras finitas. Ela diz que, para \(n\) grande, \(\sqrt{n}(\hat\beta-\beta)\) passa a ser aproximadamente normal. A qualidade dessa aproximação depende da distribuição dos erros, da presença de caudas pesadas, da heterocedasticidade e do tamanho amostral. Por isso, não existe um \(n\) universalmente “grande o suficiente” que garanta boa aproximação para todas as distribuições.
5.4 Método Delta: inferência sobre funções de \(\beta\)
O Método Delta permite obter a distribuição assintótica de funções suaves dos parâmetros. Se
Na prática, se usamos diretamente a matriz vcov() do modelo, ela já estima a variância de \(\hat\beta\), não a variância de \(\sqrt{n}(\hat\beta-\beta)\). Por isso, nesse caso, não devemos dividir novamente por \(n\).
Exemplo: razão entre os coeficientes dos dois regressores:
onde \(\hat{V}_\theta\) estima a variância assintótica de \(\sqrt{n}(\hat\theta-\theta)\). O erro comum é misturar as duas escalas: se \(\widehat{\operatorname{Var}}(\hat\theta)\) já está na escala finita-amostral, não se multiplica por \(n\).
onde \(Q_{XX}=\mathbb{E}[XX^\top]\) e \(\Omega=\mathbb{E}[XX^\top e^2]\). Sob homocedasticidade, \(\Omega=\sigma^2 Q_{XX}\), então \(V_\beta^0=\sigma^2 Q_{XX}^{-1}\).
Na prática, substituímos esses objetos por estimadores plug-in. Sob heterocedasticidade, usamos
A diferença em relação à Aula 4 é importante: aqui a distribuição normal padrão é uma aproximação assintótica, não uma distribuição exata finita-amostral. Por isso, os intervalos e testes baseados em \(1.96\) são aproximados.
Fechamento da Aula 5
A Aula 5 mostra que a inferência do MQO não depende necessariamente da normalidade exata dos erros. Sob condições de momentos, amostragem i.i.d. e posto adequado de \(Q_{XX}\), temos consistência \(\hat\beta\to_p\beta\) e normalidade assintótica:
A matriz \(V_\beta\) depende de \(\Omega=\mathbb{E}[XX^\top e^2]\), permitindo heterocedasticidade. Por isso, erros-padrão robustos HC são fundamentais para inferência assintótica válida quando a variância dos erros não é constante. A aula também mostra como usar o Método Delta para funções de parâmetros e como construir estatísticas de Wald para hipóteses conjuntas.
Aula 6 — Testes de Hipóteses e Monte Carlo
A Aula 6 organiza a inferência econométrica em forma de testes de hipóteses. Nas aulas anteriores, construímos estimadores, erros-padrão e distribuições exatas ou assintóticas. Agora, usamos esses elementos para tomar decisões estatísticas.
A lógica de um teste não é provar que \(H_0\) é verdadeira ou falsa. A lógica é perguntar:
se \(H_0\) fosse verdadeira, o resultado observado seria plausível?
Se a estatística observada fica muito distante do que esperaríamos sob \(H_0\), rejeitamos \(H_0\). Se não fica distante o suficiente, não rejeitamos \(H_0\). O fio condutor da aula é:
\[\text{estimador} + \text{erro-padrão} + \text{distribuição nula} \Rightarrow \text{regra de decisão}.\]
A aula também mostra que testes podem ser avaliados por simulação de Monte Carlo, especialmente quando queremos estudar seu comportamento em amostras finitas.
6.0 Conceitos básicos de testes de hipóteses
O parâmetro de interesse pode ser um coeficiente ou uma função dos coeficientes: \(\theta=r(\beta)\). Exemplos:
Uma hipótese pontual impõe que \(H_0:\theta=\theta_0\). Quando \(\theta=r(\beta)\), escrevemos \(H_0:r(\beta)=\theta_0\). A hipótese alternativa pode ser bilateral, \(H_1:\theta\neq\theta_0\), ou unilateral, \(H_1:\theta>\theta_0\) ou \(H_1:\theta<\theta_0\).
Um teste é uma regra de decisão que divide as possíveis amostras em duas regiões:
\[S_0=\{\text{amostras em que não rejeitamos }H_0\}, \qquad S_1=\{\text{amostras em que rejeitamos }H_0\}.\]
Na prática, resumimos a amostra por uma estatística \(T\) e rejeitamos quando ela cai na região crítica: rejeitar \(H_0\) se \(T>c\), onde \(c\) é o valor crítico.
6.1 Teste-\(t\): tamanho e poder
Para um parâmetro escalar \(\theta\), com estimador \(\hat\theta\) e erro-padrão \(s(\hat\theta)\), a estatística-\(t\) é
Para o teste bilateral \(H_0:\theta=\theta_0\) contra \(H_1:\theta\neq\theta_0\), usamos \(|T(\theta_0)|\). Sob aproximação normal assintótica, \(T(\theta_0)\Rightarrow N(0,1)\). Por isso, para um teste bilateral a 5%, a regra usual é rejeitar \(H_0\) se \(|T|>1.96\), pois \(P(|Z|>1.96)=0.05\) com \(Z\sim N(0,1)\).
Nesta seção vamos simular dois conceitos importantes:
tamanho do teste: probabilidade de rejeitar \(H_0\) quando \(H_0\) é verdadeira;
poder do teste: probabilidade de rejeitar \(H_0\) quando a alternativa é verdadeira.
library(ggplot2)set.seed(42)gera_mqo <-function(n, beta_true, sigma2 =1, dist ="normal") { X1 <-rnorm(n); X2 <-0.5*X1 +rnorm(n) X <-cbind(1, X1, X2) e <-switch(dist,"normal"=rnorm(n, sd =sqrt(sigma2)),"chi2"= (rchisq(n, 2) -2) /sqrt(2),"hetero"= (1+abs(X1)) *rnorm(n, sd =sqrt(sigma2)) ) Y <- X %*% beta_true + e XtX_inv <-solve(t(X) %*% X) b_hat <- XtX_inv %*%t(X) %*% Y e_hat <- Y - X %*% b_hat s2 <-sum(e_hat^2) / (n -ncol(X)) se_hom <-sqrt(s2 * XtX_inv[2, 2]) n_ <- n; k_ <-ncol(X) Om <-t(X) %*%diag(as.numeric(e_hat)^2) %*% X V_hc1 <- (n_/(n_-k_)) * XtX_inv %*% Om %*% XtX_inv se_rob <-sqrt(V_hc1[2, 2])list(b1 = b_hat[2], se_hom = se_hom, se_rob = se_rob)}
6.1.1 Testes bilaterais e unilaterais
No teste bilateral \(H_0:\theta=\theta_0\) contra \(H_1:\theta\neq\theta_0\), rejeitamos para valores grandes de \(|T|\). A regra assintótica a 5% é \(|T|>1.96\).
No teste unilateral à direita \(H_0:\theta=\theta_0\) contra \(H_1:\theta>\theta_0\), rejeitamos para valores grandes e positivos de \(T\). A regra a 5% é \(T>1.645\).
No teste unilateral à esquerda \(H_0:\theta=\theta_0\) contra \(H_1:\theta<\theta_0\), rejeitamos para valores muito negativos: \(T<-1.645\).
A escolha entre teste unilateral e bilateral deve ser feita antes de olhar os dados. Caso contrário, o pesquisador pode escolher a alternativa que favorece o resultado observado.
6.1.2 \(p\)-valores
O \(p\)-valor mede quão extremo foi o resultado observado sob \(H_0\). Para testes que rejeitam para valores grandes de \(T\), com distribuição nula assintótica \(G\), o p-valor é \(p=1-G(T)\).
No teste-\(t\) bilateral, usando a normal padrão:
\[p=2(1-\Phi(|t|)).\]
Se usamos a distribuição \(t\) de Student com \(\nu\) graus de liberdade: \(p=2(1-F_{t_\nu}(|t|))\).
A regra via \(p\)-valor é equivalente à regra por valor crítico: \(p<\alpha \Longleftrightarrow\) rejeitar \(H_0\) ao nível \(\alpha\). Importante: um \(p\)-valor pequeno não mede o tamanho do efeito, nem sua relevância econômica. Ele mede evidência estatística contra \(H_0\), dada a escala do erro-padrão.
O tamanho do teste é \(P(\text{rejeitar }H_0\mid H_0\text{ verdadeira})\). Se o teste é nominalmente de 5%, esperamos que, sob \(H_0\), ele rejeite aproximadamente 5% das vezes.
Nesta simulação, impomos \(\beta_1=0\), então \(H_0:\beta_1=0\) é verdadeira. Geramos erros heterocedásticos e comparamos o teste-\(t\) homocedástico convencional com o teste-\(t\) robusto à heterocedasticidade.
B <-3000ns <-c(50, 100, 300, 500)beta_nulo <-c(1, 0, -1)res_tam <-lapply(ns, function(n) { rej_hom <- rej_rob <-numeric(B)for (b in1:B) { r <-gera_mqo(n, beta_nulo, dist ="hetero") t_hom <-abs(r$b1 / r$se_hom) t_rob <-abs(r$b1 / r$se_rob) rej_hom[b] <- t_hom >qt(0.975, df = n -3) rej_rob[b] <- t_rob >1.96 }data.frame(n = n, Homocedástico =mean(rej_hom),Robusto =mean(rej_rob))})res_tam_df <-do.call(rbind, res_tam)cat("=== Taxa de rejeição sob H₀ (deveria ser 5%) ===\n")
tidyr::pivot_longer(res_tam_df, -n,names_to ="Teste",values_to ="Taxa") |>ggplot(aes(x = n, y = Taxa, color = Teste)) +geom_line(linewidth =1.2) +geom_point(size =3) +geom_hline(yintercept =0.05, linetype ="dashed", color ="gray40") +scale_color_manual(values =c("Homocedástico"="red","Robusto"="steelblue")) +labs(title ="Distorção de tamanho sob heterocedasticidade",subtitle ="Linha tracejada = nível nominal 5%",x ="n", y ="Taxa de rejeição") +theme_minimal()
Taxa de rejeição sob H₀ — distorção de tamanho
O resultado mostra que o teste homocedástico rejeita muito mais do que 5% sob heterocedasticidade — comete erro tipo I em excesso. O teste robusto fica mais próximo do tamanho nominal, especialmente conforme \(n\) aumenta. Essa é exatamente a motivação prática para usar erros-padrão robustos.
6.3 Poder do teste: detectar desvios de \(H_0\)
O poder é \(\pi(\theta)=P(\text{rejeitar }H_0\mid H_1\text{ verdadeira})\). Variamos o valor verdadeiro de \(\beta_1\): quando \(\beta_1=0\) estamos sob \(H_0\); quando \(\beta_1\neq0\), estamos sob alternativas.
n <-200B <-3000desvios <-seq(-1, 1, by =0.1)poder_hom <- poder_rob <-numeric(length(desvios))for (j inseq_along(desvios)) { beta_alt <-c(1, desvios[j], -1) rh <- rr <-numeric(B)for (b in1:B) { r <-gera_mqo(n, beta_alt, dist ="normal") rh[b] <-abs(r$b1 / r$se_hom) >qt(0.975, df = n -3) rr[b] <-abs(r$b1 / r$se_rob) >1.96 } poder_hom[j] <-mean(rh) poder_rob[j] <-mean(rr)}df_poder <-data.frame(delta = desvios, Homocedástico = poder_hom,Robusto = poder_rob)tidyr::pivot_longer(df_poder, -delta,names_to ="Teste",values_to ="Poder") |>ggplot(aes(x = delta, y = Poder, color = Teste)) +geom_line(linewidth =1.2) +geom_hline(yintercept =0.05, linetype ="dashed", color ="gray50") +geom_vline(xintercept =0, linetype ="dotted") +scale_color_manual(values =c("Homocedástico"="red","Robusto"="steelblue")) +labs(title ="Curva de poder do teste-t (n = 200)",subtitle ="Em delta = 0: taxa deve ser 5% (tamanho)",x ="β₁ verdadeiro", y ="P(rejeitar H₀)") +theme_minimal()
Curva de poder: probabilidade de rejeitar H₀ quando H₁ é verdadeira
A curva de poder é baixa perto de \(\beta_1=0\), porque ali a hipótese nula é verdadeira ou quase verdadeira. À medida que \(\beta_1\) se afasta de zero, o teste passa a rejeitar com maior probabilidade.
6.4 Teste de Wald — construção manual
O teste de Wald mede a distância entre a estimativa irrestrita e o valor imposto pela hipótese nula. Para uma hipótese vetorial \(H_0:\theta=\theta_0\) com \(\theta=r(\beta)\in\mathbb{R}^q\), a estatística é
Sob \(H_0\), \(W\Rightarrow\chi_q^2\). No caso de restrições lineares \(\theta=R\beta\), temos \(\hat\theta=R\hat\beta\) e \(\widehat{\operatorname{Var}}(\hat\theta)=R\widehat{\operatorname{Var}}(\hat\beta)R^\top\). Quando \(q=1\), o Wald é simplesmente o quadrado da estatística-\(t\): \(W=T^2\).
6.4.1 Teste-\(F\): comparação entre modelo restrito e irrestrito
O teste-\(F\) é especialmente útil para hipóteses lineares conjuntas. A estatística é
\[F = \frac{(SSE_R-SSE_U)/q}{SSE_U/(n-k)}.\]
O numerador mede quanto o ajuste piora ao impor as restrições. O denominador mede o ruído residual médio do modelo irrestrito. Para hipóteses lineares sob homocedasticidade, \(F=W^0/q\). Em amostra finita normal, usa-se \(F_{q,n-k}\). Assintoticamente, \(F\Rightarrow\chi_q^2/q\).
k <-3mod_u <-lm(as.numeric(Y) ~ X1 + X2)SSE_u <-sum(resid(mod_u)^2)Y_rest <-as.numeric(Y) -2* X1 +1* X2mod_r <-lm(Y_rest ~1)SSE_r <-sum(resid(mod_r)^2)q <-2F_sse <- ((SSE_r - SSE_u) / q) / (SSE_u / (n - k))pval_F_sse <-pf(F_sse, df1 = q, df2 = n - k, lower.tail =FALSE)cat("=== Teste F via SSE restrito e irrestrito ===\n\n")
=== Teste F via SSE restrito e irrestrito ===
cat("SSE restrito: ", round(SSE_r, 4), "\n")
SSE restrito: 319.3385
cat("SSE irrestrito: ", round(SSE_u, 4), "\n")
SSE irrestrito: 318.1061
cat("F via SSE: ", round(F_sse, 4), "\n")
F via SSE: 0.5753
cat("F via Wald/q: ", round(F_stat, 4), "\n")
F via Wald/q: 0.5753
cat("p-valor: ", round(pval_F_sse, 4), "\n")
p-valor: 0.5632
6.5 Teste de Hausman
O teste de Hausman compara dois estimadores: (1) um eficiente sob \(H_0\) mas inconsistente sob \(H_1\); (2) outro consistente sob \(H_0\) e \(H_1\), mas menos eficiente sob \(H_0\). A lógica é: se \(H_0\) é verdadeira, os dois estimadores devem ser parecidos. Se a diferença é grande demais, rejeitamos \(H_0\).
Exemplo clássico em painel: sob \(H_0\), efeitos aleatórios é consistente e eficiente; efeitos fixos também é consistente, mas menos eficiente; sob \(H_1\), efeitos aleatórios é inconsistente, enquanto efeitos fixos continua consistente. Então rejeitar \(H_0\) sugere preferir efeitos fixos.
Outro exemplo é o teste de endogeneidade Durbin-Wu-Hausman: sob exogeneidade, MQO é consistente e eficiente; sob endogeneidade, MQO é inconsistente; IV pode continuar consistente se os instrumentos forem válidos. Se MQO e IV diferem muito, há evidência contra exogeneidade.
6.6 Teste de score
O teste de score é baseado na inclinação da log-verossimilhança no ponto restrito. A ideia é: se \(H_0\) é verdadeira, o estimador restrito deve estar perto do máximo irrestrito; portanto, o score avaliado no estimador restrito deve ser pequeno. Se \(H_0\) é falsa, o ponto restrito fica longe do máximo, e o score será grande.
Comparação intuitiva entre os testes: Wald olha a distância entre o estimador irrestrito e a hipótese; Score olha a inclinação no ponto restrito; LR compara valores da log-verossimilhança; \(F\) compara perda de ajuste em soma de quadrados.
6.7 Sensibilidade do Wald a reparametrizações não lineares
A mesma hipótese escrita de formas diferentes pode dar estatísticas de Wald distintas em amostras finitas. Por exemplo, \(H_0:\beta_1/\beta_2=1\) é equivalente a \(H_0:\beta_1-\beta_2=0\) quando \(\beta_2\neq0\). Em população, descrevem o mesmo conjunto nulo. Mas em amostras finitas o Wald pode produzir estatísticas diferentes.
cat("\n→ Mesma H₀, mas distorções diferentes em n=100.\n")
→ Mesma H₀, mas distorções diferentes em n=100.
cat(" A forma linear é mais estável (evita divisão por β̂ próximo de zero).\n")
A forma linear é mais estável (evita divisão por β̂ próximo de zero).
A forma linear costuma ser mais estável, pois evita divisão por \(\hat\beta_2\). Essa é uma das mensagens centrais da Aula 6: para hipóteses não lineares, o teste de Wald pode ser sensível à parametrização.
6.8 Simulação de Monte Carlo
A simulação de Monte Carlo serve para estudar a distribuição finita-amostral de uma estatística quando essa distribuição é desconhecida. O pesquisador escolhe a distribuição dos dados \(F\), o tamanho amostral \(n\), o valor verdadeiro do parâmetro \(\theta\) e o número de repetições \(B\). Em cada réplica \(b\): gera uma amostra artificial de tamanho \(n\), calcula a estatística de interesse \(T_b\) e guarda o resultado.
Após \(B\) repetições, a coleção \(T_1,\ldots,T_B\) aproxima numericamente a distribuição finita-amostral de \(T\). Monte Carlo é útil porque a teoria assintótica pode dizer que \(T\Rightarrow N(0,1)\), mas em amostras pequenas ou médias a distribuição real de \(T\) pode ser bem diferente.
6.8.1 Características estimadas por Monte Carlo
Se queremos estudar o erro de um estimador, definimos em cada réplica \(T_b=\hat\theta_b-\theta\). O viés simulado é
O erro quadrático médio simulado é \(\widehat{MSE} = \frac{1}{B}\sum_{b=1}^B(\hat\theta_b-\theta)^2\), e a variância simulada é \(\widehat{Var} = \widehat{MSE} - \widehat{bias}^2\).
Monte Carlo também tem erro amostral. Quando calculamos uma taxa simulada de rejeição \(\hat{p}=\frac{1}{B}\sum_{b=1}^B \mathbf{1}\{\text{rejeitou na réplica }b\}\), seu erro-padrão aproximado é
Logo, maior \(B\) reduz o erro de simulação mas aumenta custo computacional. Por isso, resultados de Monte Carlo devem ser reportados com cuidado: uma diferença pequena entre duas taxas simuladas pode ser apenas ruído de simulação.
Como resultados de Monte Carlo são aleatórios, duas simulações com o mesmo desenho podem gerar resultados diferentes, especialmente quando \(B\) é pequeno. Por exemplo, se um artigo reporta, com \(B=100\), \(\hat{p}_1=0.07\), e uma replicação encontra \(\hat{p}_2=0.03\), a diferença é \(0.04\). Isso parece grande, mas com \(B=100\) o erro de simulação pode ser considerável. Se \(p\approx0.05\),
Então uma diferença de \(0.04\) não prova automaticamente falha de replicação. É preciso comparar a diferença com o erro-padrão de Monte Carlo.
Fechamento da Aula 6
A Aula 6 transforma a inferência assintótica em regras práticas de decisão. A partir de estimadores, erros-padrão e distribuições nulas, construímos testes-\(t\), testes de Wald, testes-\(F\), Hausman e score.
A estatística-\(t\) é usada para hipóteses escalares. O Wald generaliza essa lógica para hipóteses vetoriais \(H_0:\theta=\theta_0\in\mathbb{R}^q\), com \(W\Rightarrow\chi_q^2\). O teste-\(F\) compara modelos restritos e irrestritos em hipóteses lineares sob homocedasticidade. O Hausman compara dois estimadores com propriedades diferentes sob \(H_0\) e \(H_1\). O score avalia a inclinação da log-verossimilhança no ponto restrito.
A aula também mostra que simulações de Monte Carlo são uma ferramenta fundamental para estudar o comportamento finito-amostral de estimadores, testes e intervalos. Elas permitem avaliar tamanho, poder, viés, MSE, quantis e sensibilidade a hipóteses como heterocedasticidade ou não linearidade.
Aula 7 — Inferência por Reamostragem: Jackknife e Bootstrap
A Aula 7 apresenta métodos de inferência baseados em reamostragem. A ideia central é usar a própria amostra observada para aproximar a distribuição amostral de um estimador.
Se a amostra observada é \(Z_1,\ldots,Z_n\), a distribuição empírica \(\hat{F}_n\) atribui massa \(1/n\) a cada observação. Como a distribuição populacional verdadeira \(F\) é desconhecida, a reamostragem usa \(\hat{F}_n\) como aproximação de \(F\).
O objetivo é aproximar a variabilidade de \(\hat\theta\) recalculando o estimador em amostras modificadas. Há duas estratégias principais:
jackknife: remove uma observação, ou cluster, por vez;
bootstrap: sorteia novas amostras com reposição da amostra original.
A ideia é simples: se retirar uma observação muda muito o estimador, essa observação é influente e o estimador é mais instável. Se retirar qualquer observação quase não muda o resultado, o estimador é mais estável.
onde \(\bar\theta_{jack} = \frac{1}{n}\sum_{i=1}^n\hat\theta_{(-i)}\). O fator \((n-1)/n\) aparece porque os estimadores leave-one-out são calculados com amostras muito parecidas entre si; eles não são \(n\) estimadores independentes.
cat("\n→ Jackknife ≈ HC3 — confirma a equivalência do Teorema da Aula 7.\n")
→ Jackknife ≈ HC3 — confirma a equivalência do Teorema da Aula 7.
O resultado mostra que o jackknife para MQO fica numericamente muito próximo do HC3. Isso acontece porque ambos corrigem a influência de observações com alta alavancagem por meio do termo \((1-h_{ii})^{-2}\).
A intuição é que \(\tilde{e}_i = \hat{e}_i/(1-h_{ii})\) é o erro de previsão leave-one-out. Quando \(h_{ii}\) é alto, o resíduo usual \(\hat{e}_i\) pode subestimar a verdadeira influência da observação. O HC3 e o jackknife corrigem justamente esse problema.
7.2.1 Delete-cluster jackknife
Quando os dados têm clusters, remover uma observação individual pode ser inadequado. Se alunos estão agrupados em escolas, empresas em municípios ou trabalhadores em firmas, observações dentro do mesmo cluster podem compartilhar choques comuns.
Nesse caso, a unidade independente não é necessariamente a observação individual, mas o cluster. Por isso, o jackknife apropriado remove um cluster inteiro por vez.
Se há \(G\) clusters, o estimador delete-cluster é \(\hat\beta_{(-g)}\), obtido estimando o modelo sem o cluster \(g\). A matriz de variância jackknife cluster-robusta é
com \(\bar\beta=\frac{1}{G}\sum_{g=1}^G\hat\beta_{(-g)}\). A validade desse procedimento depende de \(G\) ser suficientemente grande, não apenas de \(n\). Muitos indivíduos em poucos clusters ainda podem gerar inferência frágil.
for (j in1:length(b_full)) {cat(sprintf("%-12s %10.5f\n", names(b_full)[j], se_jack_cl[j]))}
(Intercept) 0.15503
X1_cl 0.07765
X2_cl 0.07372
7.3 Bootstrap não paramétrico
O bootstrap não paramétrico sorteia amostras de tamanho \(n\), com reposição, da amostra original. Em regressão, cada unidade reamostrada deve ser o par completo \(W_i=(Y_i,X_i)\).
A amostra bootstrap é \(W_1^*,\ldots,W_n^*\sim i.i.d.\ \hat{F}_n\). Em cada amostra bootstrap, recalculamos o estimador \(\hat\theta_b^*=\hat\theta(W_1^*,\ldots,W_n^*)\). Depois de \(B\) repetições, usamos a variabilidade de \(\hat\theta_1^*,\ldots,\hat\theta_B^*\) para estimar erro-padrão, quantis e intervalos.
O sorteio precisa ser com reposição. Se sorteássemos sem reposição uma amostra de tamanho \(n\), apenas embaralharíamos os dados originais, e o estimador não mudaria.
x_grid <-seq(min(beta1_boot), max(beta1_boot), length.out =200)ggplot(data.frame(b1 = beta1_boot), aes(x = b1)) +geom_histogram(aes(y =after_stat(density)),bins =60, fill ="steelblue", alpha =0.6) +geom_line(data =data.frame(x = x_grid,y =dnorm(x_grid, mean = b_hat[2], sd = se_boot)),aes(x = x, y = y), color ="red", linewidth =1.2) +geom_vline(xintercept = b_hat[2], color ="black", linetype ="dashed") +labs(title ="Distribuição bootstrap de β̂₁ (B = 3000)",subtitle ="Vermelho = Normal com EP bootstrap | Tracejado = β̂₁ observado",x ="β̂₁*", y ="Densidade") +theme_minimal()
Distribuição bootstrap de β̂₁ vs Normal assintótica
A distribuição bootstrap de \(\hat\beta_1^*\) aproxima a distribuição amostral de \(\hat\beta_1\) usando a amostra observada como se fosse a população. No universo bootstrap, a “verdade” é \(\hat{F}_n\), e o análogo do parâmetro verdadeiro é \(\hat\beta_1\). Por isso, a distribuição relevante é frequentemente \(\hat\beta_1^*-\hat\beta_1\).
7.3.1 Por que cerca de 63,2% das observações aparecem?
Em uma amostra bootstrap de tamanho \(n\), a probabilidade de uma observação específica não ser sorteada em uma retirada é \(1-\frac{1}{n}\). A probabilidade de ela não aparecer em nenhuma das \(n\) retiradas é \(\left(1-\frac{1}{n}\right)^n\). Logo, a probabilidade de aparecer pelo menos uma vez é
7.3.2 Exemplo: salário esperado em nível a partir de regressão em log
Considere o modelo \(\log(sal_i)=\beta_1 educ_i+\beta_2+e_i\). O interesse é \(\mu=\mathbb{E}[sal\mid educ=16]\).
Como \(sal=\exp(16\beta_1+\beta_2+e)\) e, se \(e\sim N(0,\sigma^2)\), temos \(\mathbb{E}[\exp(e)]=\exp(\sigma^2/2)\), então
\[\mu=\exp(16\beta_1+\beta_2+\sigma^2/2).\]
O estimador plug-in é \(\hat\mu=\exp(16\hat\beta_1+\hat\beta_2+\hat\sigma^2/2)\). Esse parâmetro é uma função não linear de \((\hat\beta_1,\hat\beta_2,\hat\sigma^2)\). Poderíamos usar método delta, mas o bootstrap é uma alternativa prática: reamostramos os dados, reestimamos o modelo e recalculamos \(\hat\mu^*\).
7.4 Inferência para parâmetro não linear: ponto de máximo
O exemplo central da Aula 7 é a regressão quadrática em experiência:
O ponto de máximo da parábola é \(\theta = -\beta_2 / 2\beta_3\), e o estimador plug-in é \(\hat\theta=-\hat\beta_2/2\hat\beta_3\).
Esse exemplo é delicado porque \(\hat\theta\) é uma razão. Se \(\hat\beta_3\) estiver próximo de zero, o estimador pode ficar muito instável. Por isso, para testar \(H_0:\theta=\theta_0\), muitas vezes é mais estável usar a forma linear equivalente
ggplot(data.frame(theta = theta_boot), aes(x = theta)) +geom_histogram(aes(y =after_stat(density)),bins =60, fill ="steelblue", alpha =0.6) +geom_vline(xintercept = theta_true, color ="red",linewidth =1.2, linetype ="dashed") +geom_vline(xintercept = theta_hat, color ="black", linewidth =1) +geom_vline(xintercept = IC_boot, color ="darkgreen",linetype ="dotted", linewidth =1) +labs(title ="Distribuição bootstrap de θ̂ = -β̂₂/2β̂₃",subtitle ="Vermelho = θ verdadeiro | Verde = IC 95% bootstrap",x ="θ̂* (anos de experiência no máximo salarial)",y ="Densidade") +theme_minimal()
Distribuição bootstrap de θ̂ (ponto de máximo)
7.4.1 Teste da hipótese sobre o ponto de máximo
Para testar \(H_0:\theta=\theta_0\) com \(\theta=-\beta_2/2\beta_3\), podemos reescrever a hipótese como \(H_0:\beta_2+2\theta_0\beta_3=0\). Essa forma é linear em \(\beta_2\) e \(\beta_3\), e tende a ser mais estável do que trabalhar diretamente com a razão.
theta0 <- theta_truer_theta <-c(0, 0, 1, 2* theta0)R_hat <-as.numeric(r_theta %*%coef(mod2))se_R <-sqrt(as.numeric(t(r_theta) %*%vcov(mod2) %*% r_theta))T_R <- R_hat / se_Rp_R <-2* (1-pnorm(abs(T_R)))cat("=== Teste linear para H0: theta =", theta0, "===\n\n")
cat("Erro-padrão de R(theta0): ", round(se_R, 5), "\n")
Erro-padrão de R(theta0): 0.00281
cat("T: ", round(T_R, 4), "\n")
T: 0.0687
cat("p-valor: ", round(p_R, 4), "\n")
p-valor: 0.9453
7.4.2 \(p\)-valor bootstrap para a restrição linear
Também podemos usar o bootstrap para aproximar a distribuição da restrição \(R(\theta_0)=\hat\beta_2+2\theta_0\hat\beta_3\). Em cada réplica bootstrap, calculamos
Esta seção é computacionalmente mais pesada e está desativada por padrão (eval: false). Para estudar cobertura de forma adequada, o ideal é gerar novos bancos de dados a partir do processo gerador verdadeiro e, dentro de cada banco, construir os intervalos por método delta e bootstrap.
A Aula 7 mostra como usar a própria amostra para aproximar a distribuição amostral de estimadores e estatísticas. O jackknife mede sensibilidade removendo uma observação por vez; no MQO, ele se conecta diretamente aos resíduos leave-one-out, à alavancagem e ao HC3. Em dados agrupados, a remoção deve ser feita por cluster, pois a dependência intracluster torna inadequado remover observações isoladas.
O bootstrap não paramétrico sorteia amostras com reposição da distribuição empírica \(\hat{F}_n\). A distribuição dos estimadores bootstrap aproxima a distribuição amostral do estimador original. No universo bootstrap, a “verdade” é \(\hat{F}_n\), e por isso a distribuição relevante é frequentemente \(\hat\theta^*-\hat\theta\).
A aula também mostra que reamostragem é especialmente útil para parâmetros não lineares, como salários esperados em nível ou pontos de máximo de regressões quadráticas. No entanto, bootstrap não é mágica: sua validade depende de regularidade, momentos adequados e estabilidade da função estimada. Razões, divisões por coeficientes próximos de zero e estimadores na fronteira exigem cuidado.
