Exercício 1

Qual dos seguintes itens são consequência da heterocedasticidade?
a) O estimador de MQO, βhat é inconsistente.
b) A estatística F não tem mais uma distribuição F.
c) O estimador de MQO não é mais BLUE

Resposta:

A) A altenativa A está incorreta pois a heterocedasticidade afeta a variância do estimador tornando-o ineficiente, entretanto o mesmo continua sendo não viesado e consistente visto que esperança dos erros continuará sendo zero. Respeitando assim as outras suposições fundamentais para os Estimadores de MQO

B) A altenativa B está correta pois na presença de heterocedasticidade, a estatística F não segue mais uma distribuição F padrão. Em vez disso, sua distribuição torna-se distorcida, o que afeta a validade dos testes de hipóteses que dependem dessa distribuição.

C) A alternativa C está correta.Quando a heterocedasticidade está presente nos dados, a variância dos erros não é constante em relação às variáveis independentes. Isso viola uma das suposições fundamentais para a aplicação da regressão por mínimos quadrados ordinários (MQO). Como resultado, os estimadores MQO não são mais os melhores estimadores lineares não tendenciosos (BLUE) para os parâmetros do modelo. Em outras palavras, a heterocedasticidade torna os estimadores MQO não eficientes e, portanto, não BLUE.

Exercício 2

WLS é preferido ao OLS quando uma variável importante foi omitida do modelo. Esta afirmativa é verdadeira ou falsa?

Resposta:

Falso. Não podemos afirmar que WLS é melhor uma vez que ele só corrigirá a variância do estimador e não seu viés. Tanto OLS quanto WLS serão viesados.

Exercício 3

A variável “smokes” é uma variável binária que é igual a um se uma pessoa fuma e zero caso contrário. Utilizando os dados em SMOKE, estimamos um modelo de probabilidade linear para “smokes”:


A variável “white” é igual a um se o respondente é branco e zero caso contrário; as outras variáveis independentes são definidas no Exemplo 8.7. Ambos os erros padrão usuais e robustos à heterocedasticidade são relatados.
a) Existem diferenças importantes entre os dois conjuntos de erros padrão?
b) Mantendo os outros fatores fixos, se a educação aumentar em quatro anos, o que acontece com a probabilidade estimada de fumar?
c) Interprete o coeficiente na variável binária “restaurn” (uma variável dummy igual a um se a pessoa mora em um estado com restrições de fumo em restaurantes). 1
d) A pessoa número 206 no conjunto de dados tem as seguintes características: “cigpric” = 67,44, “income” = 6.500, “educ” = 16, “age” = 77, “restaurn” = 0, “white” = 0 e “smokes” = 0. Calcule a probabilidade prevista de fumar para essa pessoa e comente sobre o resultado.

Resposta:

a) A principal diferença é que o erro padrão robusto é calculado para corrigir problemas de medida relacionado a problemas de heterocedasticidade na amostra, entretanto ambas as variavéis tem como objetivo medir a variabilidade padrão do modelo.

b) Se a educação aumentar em 4 anos a probabilidade de fumar cai - $ Smokes = - 0,029 * educ = - 0,029 * 4 = - 0,116 $ - usando o conceito de Ceteris Paribus.

c) Se a pessoa mora em um estado onde há restrição de fumo em restaurantes a uma menor probabilidade de que ela seja fumante.

D) $ Smokes = 0.656 - 0.069 * log(67,44) + 0,012 * log(6500)- 0,029 * 16 + .020 * 77 - 0,00026* (77)^2 -0,101 * 0 - 0,026 *0 = 0,1040196 $

Podemos interpretar esse resultado explicando a formula, a propabilidade da pessoa ser fumante aumenta se a renda dela aumentar (poder de compra) e diminui se o preço do cigarro aumentar, ou se ela for uma pessoa educada e ciente dos problemas do fumo, ou se houver proibição de fumar em restaurantes e se ela for branca (devido a problemas estruturais da sociedade). Além disso a idade atua principalmente aumentando a probabilidade pois as gerações mais velhas tendem a fumar mais entretanto a idade também tem um leve efeito negativo na propabilidade pois as pessoas conforme amadurecem compreendem que o cigarro faz mal.

Exercício 4

Instruções Acesse o site do IBGE para baixar o dicionário da base PNAD Contínua: https://www.ibge.gov.br/estatisticas/downloads-estatisticas.html?caminho=Trabalho_e_Rendimento/Pesquisa_Nacional_por_Amostra_de_Domicilios_continua/Anual/Microdados/Visita. Através do RCran é possível acessar o as funções do pacote PNADcIBGE - que deve ser instalado para que vocês consigam fazer o download da base direto no r(o script que disponibilizei na primeira aula tem as mesmas funções).
Após baixar os dados da primeira visita de 2022, resolva:
a) Utilizando o dicionário, filtre as seguintes informações do seu banco de dados: escolaridade, idade,salário (renda da atividade principal), sexo e cor/raça.
b) Verifique com comando anyNA() a existência de NAs (dados faltantes), analise a classe das variáveis e as estatísticas descritivas.resolva:
c) Gere a regressão:
\(log(salário) = β0 + β1 ∗ Escolaridade + β2 ∗ Idade + β3 ∗ Idade2 + β4 ∗ Sexo + β5 ∗ Cor\)
Interprete os resultados obtidos pelos estimadores.
d) Os estimadores do item c) são confiáveis? Existe alguma variável que pode estar sendo omitida pelo modelo que causaria viés? E quanto à sua eficiência? Mostre ao menos dois testes vistos em aula que corroborem com a presença ou não de heterocedasticidade nos dados.
e) Caso haja heterocedasticidade nos dados, gere uma nova regressão ajustada

Resposta:

Primeiramente vamos importar as bibliotecas e definir o diretório no qual a base da pnad será Baixada:

#encoding
options(encoding = "UTF-8") #codificação dos caracteres
options(scipen = 999) #desliga a notação científica
rm(list = ls()) #limpa o environment


#packages
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 errors
library(PNADcIBGE)
library(ggplot2)
library(zoo)
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(lmtest)
library(sandwich)

#directory
directory <- setwd("C:/Users/rafael.souza/Documents/Bases/Lista")

Em seguida vamos carregar a pnad, filtrar apenas os itens necessários para utilização e renomear as colunas de acordo com as nossas preferências para salvar em um arquivo .csv ajustado.

pnad22 <- get_pnadc(
  year = 2022,
  interview = 1,
  var = c("Ano", "Trimestre", "UF","V1022","V2007","V2009","V2010","VD3005","VD4016"),
  selected = FALSE,
  defyear = 2022,
  defperiod = NULL,
  labels = TRUE,
  deflator = TRUE,
  design = TRUE,
  reload = TRUE,
  savedir = directory 
)
## Paths of files downloaded in this function at the save directory provided are:
## C:/Users/rafael.souza/Documents/R_Projects/Aula 1 ISIS/deflator_PNADC_2022.xls
## C:/Users/rafael.souza/Documents/R_Projects/Aula 1 ISIS/dicionario_PNADC_microdados_2022_visita1_20231129.xls
## C:/Users/rafael.souza/Documents/R_Projects/Aula 1 ISIS/input_PNADC_2022_visita1_20231129.txt
## C:/Users/rafael.souza/Documents/R_Projects/Aula 1 ISIS/PNADC_2022_visita1.txt
## C:/Users/rafael.souza/Documents/R_Projects/Aula 1 ISIS/PNADC_2022_visita1_20231129.zip
data <- pnad22$variables %>% 
  select("V2007","V2009","V2010","VD3005","VD4016")

#tirando as informações faltantes
data <- data[!is.na(data$VD4016), ]

#mudando o nome das colunas
colnames(data)<- c("sexo","idade","cor","escolaridade","salario")

data <- data %>% 
  mutate(escolaridade = str_replace(escolaridade," anos de estudo","")) %>% 
  mutate(escolaridade = str_replace(escolaridade," ano de estudo","")) %>%
  mutate(escolaridade = str_replace(escolaridade," anos ou mais de estudo","")) %>%
  mutate(escolaridade = str_replace(escolaridade,"Sem instrução e menos de 1","0"))

data$escolaridade <- as.double(data$escolaridade)

write.csv(data,"pnad_22_ajust.csv")

pnad <- read_csv("pnad_22_ajust.csv")
## New names:
## • `` -> `...1`
## Rows: 159372 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): sexo, cor
## dbl (4): ...1, 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.
pnad <- pnad %>% 
  filter(cor != "Ignorado") %>% 
  mutate(cor = ifelse(cor %in% c("Indígena","Parda","Preta"),"Pretos ou Pardos",cor)) %>% 
  mutate(cor = ifelse(cor == "Amarela","Branca",cor))


write.csv(pnad,"pnad_22_ajust.csv")

pnad <- read_csv("pnad_22_ajust.csv")
## New names:
## Rows: 159327 Columns: 7
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (2): sexo, cor dbl (5): ...1, ...2, 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`

Ajustamos algumas variaveis e rodamos a regressão igual a do enunciado:

#ajustando categoria binária
pnad$sexo <- as.factor(pnad$sexo)
levels(pnad$sexo) #homem é nossa categoria de referência
## [1] "Homem"  "Mulher"
pnad <- pnad %>% 
  mutate(idade2 = idade^2)

mqo <- lm(log(salario) ~ escolaridade+idade+idade2+sexo+cor, data=pnad)

mqo
## 
## Call:
## lm(formula = log(salario) ~ escolaridade + idade + idade2 + sexo + 
##     cor, data = pnad)
## 
## Coefficients:
##         (Intercept)         escolaridade                idade  
##           5.0843257            0.1049163            0.0613451  
##              idade2           sexoMulher  corPretos ou Pardos  
##          -0.0005558           -0.3741120           -0.2812203

Deste modo encontramos os parametros da regressão.

Sendo então a função:
\(log(salário) = 5,08432 + 0,1049 ∗ Escolaridade + 0,0613 ∗ Idade -0,0005 ∗ Idade2 -0,3741 ∗ Sexo -0,2812 ∗ Cor\)

Agora plotarei em um gráfico o resíduos da regressão:

pnad$residuos<-mqo$residuals
pnad$yhat<-mqo$fitted

graf<-ggplot(data=pnad,
             mapping=aes(y=residuos,x=yhat))+
  geom_point(col="purple")
graf

Observando o gráfico podemos ver que os erros não se comportam de forma homocedastica,

Rodarei então testes de inferência no modelo.

Teste de White

O teste de white tem por hipótese nula que \(δ1=δ2=0\) dada a seguinte regressão:

\(u2=δ0+δ1y^+δ2y^2\)

#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: 41585505.5837144"
print(paste("Valor p do teste F de White:", pvalor_F_white))
## [1] "Valor p do teste F de White: 0"
#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] 1968.405
pchisq(LM_estat,k,lower.tail = F)
## [1] 0
bptest(mqo, studentize = F)
## 
##  Breusch-Pagan test
## 
## data:  mqo
## BP = 4150.2, df = 5, p-value < 0.00000000000000022
#estatisitica F

F_estat <- (r2_u2/k)/(1-r2_u2)*(n-k-1)
print(F_estat)
## [1] 398.5906
pvalor_F <- df(F_estat,k,n-k-1)
print(pvalor_F)
## [1] 0
coeftest(mqo, vcov = vcovHC(mqo, type="HC0"))
## 
## t test of coefficients:
## 
##                         Estimate   Std. Error t value              Pr(>|t|)    
## (Intercept)          5.084325712  0.018282725 278.094 < 0.00000000000000022 ***
## escolaridade         0.104916340  0.000562051 186.667 < 0.00000000000000022 ***
## idade                0.061345059  0.000889658  68.954 < 0.00000000000000022 ***
## idade2              -0.000555815  0.000010837 -51.287 < 0.00000000000000022 ***
## sexoMulher          -0.374111996  0.004093028 -91.402 < 0.00000000000000022 ***
## corPretos ou Pardos -0.281220291  0.004018190 -69.987 < 0.00000000000000022 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mqo
## 
## Call:
## lm(formula = log(salario) ~ escolaridade + idade + idade2 + sexo + 
##     cor, data = pnad)
## 
## Coefficients:
##         (Intercept)         escolaridade                idade  
##           5.0843257            0.1049163            0.0613451  
##              idade2           sexoMulher  corPretos ou Pardos  
##          -0.0005558           -0.3741120           -0.2812203

Em todos os casos encontramos um p-valor menordo que 0.05 e rejeitamos a hipotese de homocedasticidade.

Por isso vamos rodar uma nova regressão ajustada com o metodo de Minimos Quadrados Generalizados

r <- 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)

mqgf
## 
## Call:
## lm(formula = mqo, data = pnad, weights = w)
## 
## Coefficients:
##         (Intercept)         escolaridade                idade  
##           5.0843257            0.1049163            0.0613451  
##              idade2           sexoMulher  corPretos ou Pardos  
##          -0.0005558           -0.3741120           -0.2812203

Sendo essa nova regressão:

\(log(salário) = 5,08432 + 0,1049 ∗ Escolaridade + 0,0613 ∗ Idade -0,0005 ∗ Idade2 -0,3741 ∗ Sexo -0,2812 ∗ Cor\)