Exemplo 1: Heterocedasticidade

Gastos per capita em escolas públicas e renda per capita por estado em 1979

data("PublicSchools")
(data <- PublicSchools) ; str(data)

               Expenditure Income
Alabama                275   6247
Alaska                 821  10851
Arizona                339   7374
Arkansas               275   6183
California             387   8850
Colorado               452   8001
Connecticut            531   8914
Delaware               424   8604
Florida                316   7505
Georgia                265   6700
Hawaii                 403   8380
Idaho                  304   6813
Illinois               437   8745
Indiana                345   7696
Iowa                   431   7873
Kansas                 355   8001
Kentucky               260   6615
Louisiana              316   6640
Maine                  327   6333
Maryland               427   8306
Massachusetts          427   8063
Michigan               466   8442
Minnesota              477   7847
Mississippi            259   5736
Missouri               274   7342
Montana                433   7051
Nebraska               294   7391
Nevada                 359   9032
New Hampshire          279   7277
New Jersey             423   8818
New Mexico             388   6505
New York               447   8267
North Carolina         335   6607
North Dakota           311   7478
Ohio                   322   7812
Oklahoma               320   6951
Oregon                 397   7839
Pennsylvania           412   7733
Rhode Island           342   7526
South Carolina         315   6242
South Dakota           321   6841
Tennessee              268   6489
Texas                  315   7697
Utah                   417   6622
Vermont                353   6541
Virginia               356   7624
Washington             415   8450
Washington DC          428  10022
West Virginia          320   6456
Wisconsin               NA   7597
Wyoming                500   9096

'data.frame':   51 obs. of  2 variables:
 $ Expenditure: int  275 821 339 275 387 452 531 424 316 265 ...
 $ Income     : int  6247 10851 7374 6183 8850 8001 8914 8604 7505 6700 ...

Não temos o gasto per capita para o estado de Wisconsin

data <- na.omit(data) ; str(data)

'data.frame':   50 obs. of  2 variables:
 $ Expenditure: int  275 821 339 275 387 452 531 424 316 265 ...
 $ Income     : int  6247 10851 7374 6183 8850 8001 8914 8604 7505 6700 ...
 - attr(*, "na.action")=Class 'omit'  Named int 50
  .. ..- attr(*, "names")= chr "Wisconsin"

names(data) <- c("gastos", "renda")

Deixando a renda na escala de 10mil obamas

data$renda <- data$renda * .0001

xyplot(gastos ~ renda, data
       , pch = 16
       , jitter.x = TRUE
       , ylab = list(rot = 0)
       , panel = panel.smoothScatter) +
  as.layer(xyplot(data$gastos ~ data$renda, type = "g"))

Iniciamos com um modelo de regressão linear de termo quadrático

model <- lm(gastos ~ renda + I(renda^2), data)

Será que o termo quadrático é realmente necessário?

Testes parciais quase-t ou quase-z para todos os coeficientes do modelo podem ser computados com a função coeftest

coeftest(model, df = Inf # df = infinito, i.e., aproximação normal é usada para computar os p-valores  
         , vcov = vcovHC(model, type = "HC0")) # HC0 = estimador de White


z test of coefficients:

            Estimate Std. Error z value Pr(>|z|)  
(Intercept)   832.91     460.89  1.8072  0.07073 .
renda       -1834.20    1243.04 -1.4756  0.14006  
I(renda^2)   1587.04     829.99  1.9121  0.05586 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

O p-valor é um pouquinho maior que 5% - a regra é clara, arnaldo - a 5% de significância o termo quadrático não é significativo

E se usarmos outro estimador, o HC4 proposto pelo Cribari-Neto, por exemplo

coeftest(model, df = Inf
         , vcov = vcovHC(model, type = "HC4"))


z test of coefficients:

            Estimate Std. Error z value Pr(>|z|)
(Intercept)   832.91    3008.01  0.2769   0.7819
renda       -1834.20    8183.19 -0.2241   0.8226
I(renda^2)   1587.04    5488.93  0.2891   0.7725

É, o termo quadrático é claramente não significativo / desnecessário

Mas por que os p-valores diferem tanto de um estimador para o outro?

pred <- predict(model
                , newdata = data.frame(renda = seq(min(data$renda), max(data$renda)
                                                   , length.out = 3 * nrow(data)))
                , interval = "confidence")
pred2 <- predict(lm(gastos ~ renda, data)
                 , newdata = data.frame(renda = seq(min(data$renda), max(data$renda)
                                                    , length.out = 3 * nrow(data)))
                 , interval = "confidence")
xyplot(gastos ~ renda, data
       , pch = 16
       , jitter.x = TRUE
       , type = c("p", "g")
       , ylab = list(rot = 0)
       , key = list(corner = c(.05, .95), lines = list(col = c(2, 3), lwd = 2)
                    , text = list(c("Modelo com termo quadrático", "Modelo sem termo quadrático")))) +
  as.layer(xyplot(pred[, 1] ~ seq(min(data$renda), max(data$renda), length.out = 3 * nrow(data))
                  , type = "l"
                  , prepanel = prepanel.cbH
                  , cty = "bands"
                  , ly = pred[, 2]
                  , uy = pred[, 3]
                  , panel = panel.cbH
                  , lwd = 2
                  , col = 2)) +
  as.layer(xyplot(pred2[, 1] ~ seq(min(data$renda), max(data$renda), length.out = 3 * nrow(data))
                  , type = "l"
                  , prepanel = prepanel.cbH
                  , cty = "bands"
                  , ly = pred2[, 2]
                  , uy = pred2[, 3]
                  , panel = panel.cbH
                  , lwd = 2
                  , col = 3))

Podemos observar no gráfico um outlier, ele corresponde ao estado do Alaska

A correção para pontos de alta alavancagem presente no estimador HC4 é a responsável por essa diferença de p-valores

Exemplo 2: Autocorrelação

Uma série temporal anual de 1963 até 1982 do produto nacional bruto nominal (GNP - gross national product) - investimento interno privado bruto nominal, um índice de preços e taxa de juros que é usada para formular um modelo que explica o investimento real pelo GNP real e juros real

data("Investment")
(data2 <- Investment)

Time Series:
Start = 1963 
End = 1982 
Frequency = 1 
        GNP Investment  Price Interest   RealGNP  RealInv     RealInt
1963  596.7       90.9 0.7167     3.23  832.5659 126.8313          NA
1964  637.7       97.4 0.7277     3.55  876.3227 133.8464  2.01518767
1965  691.1      113.5 0.7436     4.04  929.3975 152.6358  1.85503367
1966  756.0      125.7 0.7676     4.50  984.8880 163.7572  1.27245831
1967  799.6      122.8 0.7906     4.19 1011.3838 155.3251  1.19364773
1968  873.4      133.3 0.8254     5.16 1058.1536 161.4975  0.75827979
1969  944.0      149.3 0.8679     5.87 1087.6829 172.0244  0.72098134
1970  992.7      144.2 0.9145     5.95 1085.5112 157.6818  0.58071782
1971 1077.6      166.4 0.9601     4.88 1122.3831 173.3153 -0.10633133
1972 1185.9      195.0 1.0000     4.50 1185.9000 195.0000  0.34418290
1973 1326.4      229.8 1.0575     6.44 1254.2790 217.3050  0.69000000
1974 1434.2      228.7 1.1508     7.83 1246.2635 198.7313 -0.99269504
1975 1549.2      206.1 1.2579     6.25 1231.5764 163.8445 -3.05656934
1976 1718.0      257.9 1.3234     5.50 1298.1714 194.8768  0.29290882
1977 1918.3      324.1 1.4005     5.46 1369.7251 231.4174 -0.36590298
1978 2163.9      386.6 1.5042     7.46 1438.5720 257.0137  0.05550161
1979 2417.8      423.0 1.6342    10.28 1479.5007 258.8422  1.63753224
1980 2631.7      401.9 1.7842    11.77 1475.0028 225.2550  2.59119692
1981 2954.1      474.9 1.9514    13.42 1513.8362 243.3637  4.04885327
1982 3073.0      414.5 2.0688    11.02 1485.4022 200.3577  5.00380650

model2 <- lm(RealInv ~ RealGNP + RealInt, data2)

xyplot(data2[, "RealInv"], type = c("b", "g")
       , pch = 16
       , lwd = 2
       , xlab = NULL
       , scales = list(y = list(rot = 0))
       , panel = function(...){
         panel.xyplot(...)
         panel.lines(ts(fitted(model2), start = 1964), col = 3, lwd = 2)})

sc3 <- scatterplot3d(data2[, c(5, 7, 6)], type = "b"
                     , angle = 65
                     , scale.y = 1
                     , pch = 16
                     , lwd = 2) ; sc3$plane3d(model2, lty.box = "solid", col = "#0080ff", lwd = 2)

coeftest(model2, df = Inf
         , vcov = NeweyWest(
           model2 # estimador de Newey e West que propõem o uso de pesos de linear decaímento
           , lag = 4 # específica o lag máximo
                     # se lag = NULL (default)
                     # um procedimento de seleção de largura de banda não paramétrico é usado
           , prewhite = FALSE) # se TRUE usa um modelo VAR(1) (defalut)
)                              # vetor autoregressivo de ordem 1 - nas funções de estimação


z test of coefficients:

              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -12.533601  18.958298 -0.6611   0.5085    
RealGNP       0.169136   0.016751 10.0972   <2e-16 ***
RealInt      -1.001438   3.342375 -0.2996   0.7645    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

coeftest(model2, df = Inf
         , vcov = NeweyWest) # tudo no default dela


z test of coefficients:

              Estimate Std. Error z value  Pr(>|z|)    
(Intercept) -12.533601  24.374177 -0.5142    0.6071    
RealGNP       0.169136   0.023586  7.1709 7.449e-13 ***
RealInt      -1.001438   3.639935 -0.2751    0.7832    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Caso você precise usar um determinado estimador HAC várias vezes em um script você pode criar uma função para agilizar seu trabalho

parzenHAC <- function(x, ...)
  kernHAC(x, kernel = "Parzen" # kernel Parzen
          , prewhite = 2 # VAR(2) - vetor autoregressivo de ordem 2
          , adjust = FALSE # sem ajuste para amostras de tamanho finito
          , bw = bwNeweyWest, ...) # estimador de Newey e West

Os três estimadores geram erros padrão diferentes, contudo as tomadas de decisão em relação as covariáveis é sempre a mesma

coeftest(model2, df = Inf
         , vcov = parzenHAC)


z test of coefficients:

              Estimate Std. Error z value  Pr(>|z|)    
(Intercept) -12.533601  24.663944 -0.5082    0.6113    
RealGNP       0.169136   0.020835  8.1181 4.737e-16 ***
RealInt      -1.001438   3.947469 -0.2537    0.7997    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Exemplo 3:
Heterocedasticidade e autocorrelação

Estimadores de covariância podem ser utilizados em outras situações, não apenas em testes quase-t (ou quase-z)

Série temporal trimestral da taxa de juros real ex-post dos Estados Unidos de 1961 até 1986

data("RealInt")

plot(RealInt
     , las = 1
     , pch = 16
     , lwd = 2
     , col = "#0080ff"
     , xlab = NULL)
abline(v = seq(1961, 1986, length = 26), h = seq(-5, 10, length = 7), col = "gray90")

gefp: computa o processo de flutuação M empírico a partir dos scores de um modelo ajustado

(model3 <- gefp(RealInt ~ 1, fit = lm # ajuste de uma regressão na média
                , vcov = kernHAC # variância estimada com a função kernHAC
                                 # estimador HAC kernel espectral quadrático
                                 # com VAR(1) nas funções de ligação
                                 # e seleção de largura de banda automática
               ))                # baseada em uma aproximação AR(1)


Generalized Empirical M-Fluctuation Process

Call: gefp(RealInt ~ 1, fit = lm, vcov = kernHAC)


Fitted model: 
Call:
fit(formula = ..1, data = data)

Coefficients:
(Intercept)  
      1.375

plot(model3, aggregate = FALSE # valores críticos à 5% em vermelho
                               # se o processo ajustado ultrapassa a linha horizontal
                               # existe uma mudança significativa na média
     , xlab = "", main = NULL) # neste caso existe pelo menos uma forte quebra em torno de 1980

sctest(model3) # teste formal para tal verificação


    M-fluctuation test

data:  model3
f(efp) = 1.6586, p-value = 0.008159

bp <- breakpoints(RealInt ~ 1) # por default os pontos são selecionados pela partição BIC mínima

confint(bp
        , vcov = kernHAC)


     Confidence intervals for breakpoints
     of optimal 3-segment partition: 

Call:
confint.breakpointsfull(object = bp, vcov. = kernHAC)

Breakpoints at observation number:
  2.5 % breakpoints 97.5 %
1    37          47     48
2    77          79     81

Corresponding to breakdates:
  2.5 %   breakpoints 97.5 % 
1 1970(1) 1972(3)     1972(4)
2 1980(1) 1980(3)     1981(1)

plot(RealInt
     , las = 1
     , pch = 16
     , lwd = 2
     , col = "#0080ff"
     , xlab = NULL)
lines(ts(fitted(bp), start = start(RealInt), frequency = 4), col = 3, lwd = 2)
lines(confint(bp, vcov = kernHAC), col = 2, lwd = 2)
abline(v = seq(1961, 1986, length = 26), h = seq(-5, 10, length = 7), col = "gray90")

REGRESSÃO HETEROCEDÁSTICA

Henrique Aparecido Laureano [Lattes, GitLab, LEG GitLab]

Dezembro de 2015

Exemplo 1: Heterocedasticidade

Exemplo 2: Autocorrelação

Exemplo 3:
Heterocedasticidade e autocorrelação

REGRESSÃO HETEROCEDÁSTICA

Henrique Aparecido Laureano [Lattes, GitLab, LEG GitLab]

Dezembro de 2015

Exemplo 1: Heterocedasticidade

Exemplo 2: Autocorrelação

Exemplo 3: Heterocedasticidade e autocorrelação

Exemplo 3:
Heterocedasticidade e autocorrelação