# Funções para obtenção matricial de estimativas
getBeta <- function(X, y) {
    return(solve(t(X) %*% X) %*% t(X) %*% y)
}

getResiduals <- function(X, y, beta) {
    return(y - X %*% beta)
}

getVarBeta <- function(X, omega){
    return(solve(t(X) %*% X) %*% t(X) %*% omega %*% 
             X %*% solve(t(X) %*% X))
}

Introdução

Na presença de heterocedasticidade:

O estimador de Mínimos Quadrados Ordinários dos coeficientes de regressão permanece não-viesado e consistente.
O estimador da matriz de covariância por MQO é viesado e inconsistente.

O modelo

O modelo considerado é:

\[y = X\beta + \epsilon\] Se os erros (\(\epsilon\)) são homocedásticos:

\[ \hat{\beta} = (X'X)^{-1}X'y\] \[ \hat{Var(\hat{\beta})} = (X'X)^{-1}X'\hat{\Omega} X (X'X)^{-1}\] onde \(\hat{\Omega} = \hat{\sigma}^2I_n\).

O estimador de Halbert White

A proposta de White é utilizar:

\[\hat{\Omega} = \begin{pmatrix} \hat{e_1}^2 & 0 & \cdots & 0 \\ 0 & \hat{e_2}^2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \hat{e_n}^2 \end{pmatrix} \]

em que \(\hat{e} = (I_n - X(X'X)^{-1}X')y\). Ou seja, encontra-se uma matriz diagonal formada a partir do vetor contendo os quadrados dos resíduos de mínimos quadrados.

Vantagens:

O estimador é consistente quando os erros são homocedásticos e quando há heterocedasticidade desconhecida.

Desvantagens:

O estimador pode ser muito viesado em amostras finitas, conduzindo a testes quase-t liberais (tende a rejeitar mais \(H_0\)).
Pontos de alta alavancagem tem grande influência sobre o desempenho dos estimadores consistentes e testes associados.

Exemplo

set.seed(14112017)
x <- seq(1, 50, length.out = 1000)
y <- 10 + 0.027 * x + 
    rnorm(1000, mean = 0, sd = seq(3, 15, length.out = 1000))

X <- model.matrix(y ~ x)
hat <- X %*% solve(t(X) %*% X) %*% t(X)
beta <- getBeta(X, y)
pred <- X %*% beta
res <- (y - pred)
sigma2 <- sum(res^2)/(length(res) - 2)

xyplot(y ~ x,
       main = "Modelo Heterocedástico",
       panel = function(x, y) {
           panel.xyplot(x, y, pch = 20)
           panel.abline(beta, col = 2)
       }, key = list(space = "bottom", columns = 1, 
        text = list(lab = c("Reta regressão")),
        lines = list(lty = 1, col = 2)))

xyplot(res ~ pred, main = "Resíduos", pch = 20,
       panel = function(x, y, ...) {
           panel.xyplot(x, y, ...)
           panel.abline(h = 0, col = "#424242", lty = 2, lwd = 2)
       }, xlab = "Preditos", ylab = "Resíduos")

# Estimação da matriz de variância
(varbeta <- getVarBeta(X, sigma2 * diag(1, 
                                       nrow = length(res))))

##             (Intercept)          x
## (Intercept)   0.3918480 -0.0117453
## x            -0.0117453  0.0004606

## Estimando a matriz com o método de White:
omega <- diag(1, nrow = length(res))
omega[omega == 1] <- res^2

varbetawhite <- getVarBeta(X, omega)
varbetawhite

##              (Intercept)             x
## (Intercept)  0.199099701 -0.0085725271
## x           -0.008572527  0.0005078929

qt(0.975, length(res) - 2)

## [1] 1.962344

beta[2]/sqrt(varbeta[2, 2])

## [1] 1.985689

beta[2]/sqrt(varbetawhite[2, 2])

## [1] 1.890981

O estimador HC3

A proposta do HC3 é utilizar:

\[\hat{\Omega} = \begin{pmatrix} \frac{\hat{e_1}^2}{(1 - h_1)^2} & 0 & \cdots & 0 \\ 0 & \frac{\hat{e_1}^2}{(1 - h_2)^2} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \frac{\hat{e_n}^2}{(1 - h_n)^2} \end{pmatrix} \]

onde \(h_i\) é o i-ésimo elemento da “matriz chapéu”, \(H = X(X'X)^{-1}X'\).

Vantagens:

Geralmente possui melhor desempenho em pequenas amostras se comparado ao estimador de White.
O estimador conduz a testes quase-t que não são marcadamente liberais e, consequêntemente, a inferências mais confiáveis.

Desvantagens:

É uma aproximação do estimador jackknife considerado por MacKinnon e White.

Exemplo

omega <- diag(1, nrow = length(res))
omega[omega == 1] <- res^2/(1 - diag(hat))^2

varbetahc3 <- getVarBeta(X, omega)
varbetahc3

##              (Intercept)             x
## (Intercept)  0.200134355 -0.0086221229
## x           -0.008622123  0.0005107905

qt(0.975, length(res) - 2)

## [1] 1.962344

beta[2]/sqrt(varbeta[2, 2])

## [1] 1.985689

beta[2]/sqrt(varbetawhite[2, 2])

## [1] 1.890981

beta[2]/sqrt(varbetahc3[2, 2])

## [1] 1.88561

O estimador de bootstrap

Se utilizarmos o bootstrap em sua forma mais simples (como proposto por Bradley Efron), não consideramos a heterocedasticidade:

Exemplo 1

set.seed(15112017)
e.new <- sample(res, length(res), replace = TRUE)
y.new <- X %*% beta + e.new
xyplot(y.new ~ x, pch = 20)

Bootstrap ponderado

A equação derivada da proposta por Wu é:

\[y_i^* = X\beta + t_i^* \frac{\hat{e_i}}{(1 - h_i)}\]

em que \[t_i^*\] é um número gerado aleatoriamente de uma distribuição com média 0 e variância 1, levando em consideração a possível não-constância das variâncias dos erros.

Exemplo 2

set.seed(16112017)
tast <- rnorm(length(res))
y.new <- X %*% beta + tast * (res/(1 - diag(hat)))
y.new.or <- X %*% beta + tast * (res/sqrt(1 - diag(hat)))

gp1 <- xyplot(y.new ~ x, pch = 20, col = 1)
gp2 <- xyplot(y.new.or ~ x, col = 3, pch = 20)

grid.arrange(grobs = list(gp1, gp2), ncol = 2)

Os valores gerados pela equação de Wu original e a equação derivada são bastante similares. A diferença entre eles pode ser vista no gráfico abaixo:

xyplot(y.new - y.new.or ~ 1:length(y.new), 
       scale = list(y = list(tick.number = 18)))

boot <- replicate(2000, {
    tast <- rnorm(length(res))
    y.new <- X %*% beta + tast * (res/(1 - diag(hat)))
    XA <- model.matrix(y.new ~ x)
    as.vector(solve(t(XA) %*% XA) %*% t(XA) %*% y.new)
})

varBoot <- apply(boot, 1, var)
names(varBoot) <- c("(intercept)", "x")
varBoot

##  (intercept)            x 
## 0.1997023512 0.0005165218

qt(0.975, length(res) - 2)

## [1] 1.962344

beta[2]/sqrt(varbeta[2, 2])

## [1] 1.985689

beta[2]/sqrt(varbetawhite[2, 2])

## [1] 1.890981

beta[2]/sqrt(varbetahc3[2, 2])

## [1] 1.88561

beta[2]/sqrt(varBoot[2])

##        x 
## 1.875119

Vantagens:

É mais precisa do que aquela obtida a partir de sua aproximação assintótica de primeira ordem.

Desvantagens:

O custo computacional é alto.

O estimador HC4

A proposta do HC4 é utilizar:

\[\hat{\Omega} = \begin{pmatrix} \frac{\hat{e_1}^2}{(1 - h_1)^{\delta_1}} & 0 & \cdots & 0 \\ 0 & \frac{\hat{e_1}^2}{(1 - h_2)^{\delta_2}} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \frac{\hat{e_n}^2}{(1 - h_n)^{\delta_n}} \end{pmatrix} \]

\(\delta_i = n \frac{h_i}{p}\), em que p, que é o posto da matriz X, dado por \(\sum_{j = 1}^{n} h_i\).

Exemplo

omega <- diag(1, nrow = length(res))
delta <- length(res) * diag(hat) / 2
omega[omega == 1] <- res^2/(1 - diag(hat))^delta

varbetahc4 <- getVarBeta(X, omega)
varbetahc4

##              (Intercept)             x
## (Intercept)  0.199853804 -0.0086106656
## x           -0.008610666  0.0005100973

qt(0.975, length(res) - 2)

## [1] 1.962344

beta[2]/sqrt(varbeta[2, 2])

## [1] 1.985689

beta[2]/sqrt(varbetawhite[2, 2])

## [1] 1.890981

beta[2]/sqrt(varbetahc3[2, 2])

## [1] 1.88561

beta[2]/sqrt(varBoot[2])

##        x 
## 1.875119

beta[2]/sqrt(varbetahc4[2, 2])

## [1] 1.886891

Aqui, o expoente que controla o grau de desconto para a observação \(i\) é dado pela razão entre o valor de \(h_i\) e a média dos \(h_i\).

Teste quase-t

boot <- replicate(1999, {
    tast <- rnorm(length(res))
    y.new <- X %*% c(beta[1], 0) + tast * res/(1 - diag(hat))
    beta.new <- getBeta(X, y.new)
    res.new <- getResiduals(X, y.new, beta.new)
    s2 <- sum(res.new^2)/(length(res.new) - 2)
    vb <- diag(getVarBeta(X, s2 * diag(1, nrow = length(res.new))))
    as.vector(beta.new/sqrt(vb))
})

boot <- cbind(beta/sqrt(diag(varbeta)), boot)

plot(density(boot[2,]), 
     main = expression("Dist. empírica de " * beta[1]))
abline(v = boot[2, 1], col = 2, lty = 2)

p.value <- sum(boot[2, ] >= boot[2, 1] | 
                   boot[2, ] <= -boot[2, 1])/dim(boot)[2]
p.value

## [1] 0.0615

Bootstrap duplo

boot <- replicate(999, {
    tast <- rnorm(length(res))
    y.new <- X %*% c(beta[1], 0) + tast * res/(1 - diag(hat))
    beta.new <- getBeta(X, y.new)
    res.new <- getResiduals(X, y.new, beta.new)
    s2 <- sum(res.new^2)/(length(res.new) - 2)
    vb <- diag(getVarBeta(X, s2 * diag(1, nrow = length(res.new))))
    replicate(20, {
        tast2 <- rnorm(length(res))
        y.new2 <- X %*% c(beta.new[1], 0) + 
            tast2 * res.new/(1 - diag(hat))
        beta.new2 <- getBeta(X, y.new2)
        res.new2 <- getResiduals(X, y.new2, beta.new2)
        s22 <- sum(res.new2^2)/(length(res.new2) - 2)
        vb2 <- diag(getVarBeta(X, s22 * diag(1,
                                             nrow = length(res.new2))))
        as.vector(beta.new2/sqrt(vb2))
    })
})

bs <- NULL
kk <- apply(boot, 3, function(x) {
    bs <<- cbind(bs, x)
})

bs <- cbind(beta/sqrt(diag(varbeta)), bs)

plot(density(bs[2,]), 
    main = expression("Dist. empírica de " * beta[1]))
abline(v = bs[2, 1], col = 2, lty = 2)

p.value <- sum(bs[2, ] >= bs[2, 1] | 
                   bs[2, ] <= -bs[2, 1])/dim(bs)[2]
p.value

## [1] 0.0622091

Vantagens:

A vantagem da utilização do estimador de bootstrap duplo reside em seu bom desempenho geral.

Desvantagens:

Alto custo computacional.

Algumas conclusões

O estimador de White conduz a testes liberais: podemos encontrar significância para variáveis que não os são;
O estimador HC3 possui melhor desempenho em pequenas amostras, semelhante ao estimador de bootstrap ponderado;
Porém, ambos conduzem a testes liberais quando há pontos de alavancagem: usar o HC4, que é substancialmente mais simples do ponto de vista computacional, em comparação ao bootstrap duplo.

Inferência em Modelos Heterocedásticos

Alcides Conte Neto & Bruna Wundervald

20 de novembro de 2017

Introdução

O modelo

O estimador de Halbert White

Exemplo

O estimador HC3

Exemplo

O estimador de bootstrap

Exemplo 1

Bootstrap ponderado

Exemplo 2

O estimador HC4

Exemplo

Teste quase-t

Bootstrap duplo

Algumas conclusões