Aula 1 - Heterocedasticidade
Ísis Michelan
Econometria II
Saindo do mundo perfeito das hipóteses do modelo linear geral, o nosso querido estimador de mínimos quadrados ordinários (MQO) perde sua qualidade BLUE - best linear unbiased estimator.
Vamos relembrar rapidamente as hipóteses do Modelo Linear Geral(MLG)?
Separabilidade Aditiva;
Linearidade em parâmetros (linear em \(\beta\));
Exogeneidade (\(E(u | X) = 0\));
Ortogonalidade dos erros (\(u_i \neq u_j\))
Homocedasticidade (\(var(u | X) = \sigma^2\)).
Quando uma dessas hipóteses não é válida, precisamos desconfiar dos resultados que a nossa regressão linear nos apresenta. No caso da heterocedasticidade devemos nos preocupar com a variância/segundo momento, isso significa que o nosso termo de erro não possui variância constante em sua distribuição (condicional ao nosso X) , i.e, $ var(u|x) ^2 $. Muitas vezes a heterocedasticidade pode ser causada pela natureza da variável, por exemplo, quando a mesma tem características binárias.
média de \(\hat\beta_{mqo}\) sob heterocedasticidade
\[E(\hat\beta_{mqo}|X) = E[(X'X)^{-1}X'(X\beta+u)|X] \\ E(\hat\beta_{mqo}|X) = E[(X'X)^{-1}X'X\beta+(X'X)^{-1}X'u|X] \\ E(\hat\beta_{mqo}|X)= E(\beta)+ E[(X'X)^{-1}X'u|X] \\ E(\hat\beta_{mqo}|X)= \beta + (X'X)^{-1}X'E(u|X) \]
Sob a hipótese de ortogonalidade, \(E(u|X)=0\), então \(E(\hat\beta_{mqo}|X)= \beta\). Pela Lei das Espectativas Iteradas temos:
\[E[E(\hat\beta_{mqo}|X)]=E(\hat\beta_{mqo}) \Rightarrow \\ E(\beta)= E(\hat\beta_{mqo}) = \beta\]
Observe que o problema de heterocedasticidade não causa viés, i.e, a esperança do parâmetro amostral converge em média para o parâmetro populacional, de forma que os parâmetros estimados via MQO podem ser considerados como corretos; o problema recai apenas sobre a variância e desvio padrão produzidos pelo estimador, que como já vimos, é ineficiente.
É importante, portanto, que o pesquisador se atente aos testes que validam a existência ou não de heterocedasticidade.Vamos começar estimando via MQO uma das regressões mais famosas do mundo da economia, a equação de Mincer(1974):
\[ log(salário) = \beta_0 + \beta_1*Escolaridade + \beta_2*Idade + \beta_3*Idade^2+ \beta_4*Sexo + \beta_5*Cor \]
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.0 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorspnad <- read_csv("pnad_22_ajust.csv")## New names:
## Rows: 155539 Columns: 12
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (4): uf, zona, sexo, cor dbl (8): ...1, ...2, ...3, ano, trim, idade,
## escolaridade, salario
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`
## • `...1` -> `...2`
## • `...2` -> `...3`glimpse(pnad)## Rows: 155,539
## Columns: 12
## $ ...1 <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17…
## $ ...2 <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17…
## $ ...3 <dbl> 1, 2, 9, 11, 12, 17, 19, 23, 24, 32, 35, 36, 39, 41, 44, …
## $ ano <dbl> 2022, 2022, 2022, 2022, 2022, 2022, 2022, 2022, 2022, 202…
## $ trim <dbl> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, …
## $ uf <chr> "Rondônia", "Rondônia", "Rondônia", "Rondônia", "Rondônia…
## $ zona <chr> "Urbana", "Urbana", "Urbana", "Urbana", "Urbana", "Urbana…
## $ sexo <chr> "Mulher", "Homem", "Homem", "Homem", "Homem", "Mulher", "…
## $ idade <dbl> 38, 42, 20, 26, 19, 47, 60, 46, 46, 37, 19, 39, 57, 20, 4…
## $ cor <chr> "Pretos ou Pardos", "Pretos ou Pardos", "Pretos ou Pardos…
## $ escolaridade <dbl> 12, 6, 11, 10, 9, 0, 1, 9, 16, 16, 12, 16, 12, 11, 11, 12…
## $ salario <dbl> 1212, 1212, 1212, 1212, 1212, 1400, 1500, 2500, 3200, 220…summary(pnad)## ...1 ...2 ...3 ano
## Min. : 1 Min. : 1 Min. : 1 Min. :2022
## 1st Qu.: 38886 1st Qu.: 38886 1st Qu.:109185 1st Qu.:2022
## Median : 77770 Median : 77770 Median :213115 Median :2022
## Mean : 77770 Mean : 77770 Mean :203361 Mean :2022
## 3rd Qu.:116655 3rd Qu.:116655 3rd Qu.:298927 3rd Qu.:2022
## Max. :155539 Max. :155539 Max. :380928 Max. :2022
## trim uf zona sexo
## Min. :1.000 Length:155539 Length:155539 Length:155539
## 1st Qu.:2.000 Class :character Class :character Class :character
## Median :3.000 Mode :character Mode :character Mode :character
## Mean :2.541
## 3rd Qu.:4.000
## Max. :4.000
## idade cor escolaridade salario
## Min. :14.00 Length:155539 Min. : 0.00 Min. : 4
## 1st Qu.:30.00 Class :character 1st Qu.: 8.00 1st Qu.: 1200
## Median :40.00 Mode :character Median :12.00 Median : 1500
## Mean :40.73 Mean :10.67 Mean : 2389
## 3rd Qu.:51.00 3rd Qu.:13.00 3rd Qu.: 2500
## Max. :93.00 Max. :16.00 Max. :300000#ajustando categoria binária
pnad$sexo <- as.factor(pnad$sexo)
levels(pnad$sexo) #homem é nossa categoria de referência## [1] "Homem" "Mulher"#Criando a variável de idade^2
pnad <- pnad %>%
mutate(idade2 = idade^2)
mqo <- lm(log(salario) ~ escolaridade+idade+idade2+sexo+cor, data=pnad)
summary(mqo)##
## Call:
## lm(formula = log(salario) ~ escolaridade + idade + idade2 + sexo +
## cor, data = pnad)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.0039 -0.3898 0.0361 0.4454 4.8853
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.072e+00 1.728e-02 293.55 <2e-16 ***
## escolaridade 1.062e-01 4.952e-04 214.52 <2e-16 ***
## idade 6.015e-02 7.963e-04 75.54 <2e-16 ***
## idade2 -5.397e-04 9.279e-06 -58.17 <2e-16 ***
## sexoMulher -3.810e-01 4.032e-03 -94.50 <2e-16 ***
## corPretos ou Pardos -2.786e-01 4.014e-03 -69.39 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7655 on 155533 degrees of freedom
## Multiple R-squared: 0.3045, Adjusted R-squared: 0.3044
## F-statistic: 1.362e+04 on 5 and 155533 DF, p-value: < 2.2e-16Agora, precisamos testar a validade das estatísticas apresentadas. Podemos inicialmente realizar uma inspeção visual dos resíduos. No gráfico abaixo podemos ver que a distribuição não parece exatamente o que esperamos de um resíduo homocedástico:
#gráfico para heterocedasticidade
pnad$residuos<-mqo$residuals
pnad$yhat<-mqo$fitted
graf<-ggplot(data=pnad,
mapping=aes(y=residuos,x=yhat))+
geom_point(col="purple")
graf
Podemos então partir para os testes de hipótese. Iremos realizar o teste de Breusch-Pagan e White.
Teste de Breuch-Pagan(BP) e Teste de White
Neste teste, nosso objetivo é compreender se a hipótese 5 é atendida pelo estimador, assim, queremos observar se \(u^2\) é relacionado com alguma das variáveis explicativas. A hipótese nula é dada por:
\[H_0: Var(u|x_1,x_2,...,x_k)= \sigma^2 \]
Então se o p-valor do teste estiver acima do valor que usamos costumeiramente (0,05 e 0,01 de significância), isto é, \(p>0.05\), não rejeitamos \(H_0\) e portanto temos um erro homocedástico.
Teste de White
O teste de white tem por hipótese nula que \(\delta_1=\delta_2=0\) dada a seguinte regressão:
\[u^2=\delta_0+\delta_1\hat{y}+\delta_2\hat{y}^2\]
library(lmtest)## Carregando pacotes exigidos: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric#Teste de White - LM
white_test <- bptest(mqo, ~fitted(mqo)+I(fitted(mqo)^2))
#Teste F de White
F_white <- (white_test$statistic * (nrow(pnad) - length(coef(mqo)) - 1)) / length(coef(mqo))
# Graus de liberdade
graus_liberdade <- length(coef(mqo))
# Valor p do teste F
pvalor_F_white <- 1 - pf(F_white, df1 = graus_liberdade, df2 = length(residuals(mqo)))
print(paste("Estatística F de White:", F_white))## [1] "Estatística F de White: 38377395.931501"print(paste("Valor p do teste F de White:", pvalor_F_white))## [1] "Valor p do teste F de White: 0"A estatística do Multiplicador de Lagrange(LM) de Breuch-Pagan é calculada por:
\[LM=nR^2_{u^2}\] Em que \(R^2_{u^2}\) representa o \(R^2\) da regressão dos resíduos.
#manualmente
resid2 <- (mqo$residuals)^2
aux <- lm(resid2 ~ escolaridade+idade+idade2+sexo+cor, data=pnad)
r2_u2 <- summary(aux)$r.squared
k=5 #n de regressores
n = as.numeric(length(aux$fitted.values))
#Estatística LM
LM_estat <- (n*r2_u2)
print(LM_estat)## [1] 1801.847pchisq(LM_estat,k,lower.tail = F)## [1] 0#pacote
library(lmtest)
bptest(mqo, studentize = F)##
## Breusch-Pagan test
##
## data: mqo
## BP = 3833.9, df = 5, p-value < 2.2e-16Teste F
Já a estatísitica F do teste é dada por:
\[F= \frac{R^2_{u^2}/k}{(1-R^2_{u^2})/(n-k-1)}\]
#estatisitica F
F_estat <- (r2_u2/k)/(1-r2_u2)*(n-k-1)
print(F_estat)## [1] 364.579pvalor_F <- df(F_estat,k,n-k-1)
print(pvalor_F)## [1] 0Erro-Padrão Robusto
Uma maneira mais fácil de corrigir a variância do estimador de MQO sobre heterocedasticidade é utilizar a matriz de White e produzir um erro padrão robusto para sua estimação:
library(sandwich)
coeftest(mqo, vcov = vcovHC(mqo, type="HC0"))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.0723e+00 1.8080e-02 280.549 < 2.2e-16 ***
## escolaridade 1.0623e-01 5.6130e-04 189.256 < 2.2e-16 ***
## idade 6.0150e-02 8.7711e-04 68.577 < 2.2e-16 ***
## idade2 -5.3971e-04 1.0699e-05 -50.443 < 2.2e-16 ***
## sexoMulher -3.8100e-01 4.0779e-03 -93.432 < 2.2e-16 ***
## corPretos ou Pardos -2.7855e-01 3.9704e-03 -70.156 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Mínimos Quadrados Generalizados (MQGF)
O estimador de Mínimos Quadrados Generalizados Factíveis (MQGF) tem por sua vez a função de realizar o ajuste da variância através de uma matriz de ponderação sobre a heterocedasticidade:
#MQGF
mqo##
## Call:
## lm(formula = log(salario) ~ escolaridade + idade + idade2 + sexo +
## cor, data = pnad)
##
## Coefficients:
## (Intercept) escolaridade idade
## 5.0722530 0.1062293 0.0601496
## idade2 sexoMulher corPretos ou Pardos
## -0.0005397 -0.3810030 -0.2785506r <- log(resid2)
varreg <- lm(r ~ escolaridade+idade+idade2+sexo+cor, data=pnad)
#ponderação
w <- 1/exp(fitted(varreg))
#MQP
mqgf <- lm(mqo, weights = w, data = pnad)Em geral o MQGF será viesado, mas consistente e assintoticamente (para grande amostras) mais eficiente que MQO.
(Fonte: JEFFREY, M. Wooldridge, introductory econometrics—A modern approach. 2018.)