12.13 Modelos ARCH
12.13.1 ARCH
(R. F. Engle 1982) propõe o modelo ARCH (Autoregressive Conditional Heteroscedasticity), que permite modelar a variância condicional de uma série em função dos quadrados dos retornos passados. A variância inconstante ao longo do tempo é chamada heterocedasticidade.
O modelo ARCH(\(s\)) pode ser expresso por
\[\begin{equation} y_t = a_t \tag{12.34} \end{equation}\]
\[\begin{equation} a_t = \varepsilon_t \sqrt{h_t} \tag{12.35} \end{equation}\]
\[\begin{equation} h_t = \alpha_0 + \sum_{i=1}^s \alpha_i y_{t-i}^2 \tag{12.36} \end{equation}\]
\[\begin{equation} \varepsilon_t \sim WN(0,1) \tag{12.37} \end{equation}\]
Demonstra-se que \(a_t \sim WN(0,\sigma^2_a)\), i.e., \(a_t\) tem variância constante tal que
\[\begin{equation} \sigma^2_a = \frac{\alpha_0}{1-\sum_{i=1}^s \alpha_i} \tag{12.38} \end{equation}\]
Portanto \(a_t\) possui estacionariedade fraca ou de segunda ordem. Para mais detalhes recomenda-se (R. F. Engle 1982), (Shephard 1996) e (G. E. P. Box, Jenkins, and Reinsel 2008, 413).
A função tseries::garch()
(Trapletti and Hornik 2024) ajusta um modelo de série temporal GARCH\((r,s)\) via máxima verossimilhança.
Exemplo 12.22 Adaptado da documentação de tseries::garch()
, simulando e ajustando um modelo ARCH\((2)\) ou GARCH\((0,2)\).
library(tseries)
n <- 1100 # número de observações simuladas, 100 primeiras descartadas
a <- c(0.1, 0.5, 0.2) # coeficientes ARCH(2): alfa_0, alfa_1 e alfa_2
set.seed(42); e <- rnorm(n) # \varepsilon_t, Eq. 12.40
y <- double(n) # vetor para conter a série temporal y simulada (double-precision vector)
set.seed(314); y[1:2] <- rnorm(2, sd = sqrt(a[1]/(1-a[2]-a[3]))) # Eq. 12.41
for(i in 3:n){ # gera processo ARCH(2)
ht <- a[1] + a[2]*y[i-1]^2 + a[3]*y[i-2]^2 # Eq. 12.39
y[i] <- e[i]*sqrt(ht) # Eqs. 12.37 e 12.38
}
y <- ts(y[101:n]) # descartando as 100 primeiras observações
y_arch <- tseries::garch(y, order = c(0,2), trace = FALSE) # ajusta ARCH(2)
summary(y_arch) # testes de diagnóstico
##
## Call:
## tseries::garch(x = y, order = c(0, 2), trace = FALSE)
##
## Model:
## GARCH(0,2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.34285 -0.70636 -0.02587 0.63041 3.41272
##
## Coefficient(s):
## Estimate Std. Error t value Pr(>|t|)
## a0 0.095299 0.009129 10.439 <2e-16 ***
## a1 0.507829 0.060802 8.352 <2e-16 ***
## a2 0.225846 0.045248 4.991 6e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Diagnostic Tests:
## Jarque Bera Test
##
## data: Residuals
## X-squared = 1.5664, df = 2, p-value = 0.457
##
##
## Box-Ljung test
##
## data: Squared.Residuals
## X-squared = 0.15934, df = 1, p-value = 0.6898
12.13.2 GARCH
(Bollerslev 1986) estendeu o trabalho original de Engle, propondo o modelo GARCH (Generalized Autoregressive Conditional Heteroscedasticity), que permite que a variância condicional não dependa somente dos quadrados dos retornos passados, mas também da própria variância passada. O modelo GARCH\((r,s)\) é definido pelas Eq. (12.34), (12.35) e (12.37), com \(h_t\) dado pela Eq. (12.39).
\[\begin{equation} h_t = \alpha_0 + \sum_{i=1}^s \alpha_i y_{t-i}^2 + \sum_{i=1}^r \beta_i h_{t-i} \tag{12.39} \end{equation}\]
Para esta formulação \(a_t\) também possui variância constante, definida pela Eq. (12.40).
\[\begin{equation} \sigma^2_a = \frac{\alpha_0}{1-\sum_{i=1}^s \alpha_i - \sum_{i=1}^r \beta_i} \tag{12.40} \end{equation}\]
Exemplo 12.23 Adaptado da documentação de tseries::garch()
, simulando e ajustando um modelo GARCH(1,1). Lembre-se que “[t]he [a]utocovariances [d]o [n]ot [a]lways [d]ecrease” (Francq and Zakoian 2010, 51).
library(tseries)
n <- 1100 # número de observações simuladas, 100 primeiras descartadas
a <- c(0.1, 0.5) # coeficientes ARCH(1): alfa_0 e alfa_1
b <- c(0.3) # coeficiente GARCH(s,1): beta_1
set.seed(42); e <- rnorm(n) # \varepsilon_t, Eq. 12.40
y <- double(n) # vetor para conter a série temporal y simulada (double-precision vector)
h <- double(n) # vetor para conter a série temporal de ht
set.seed(314); y[1] <- rnorm(1, sd = sqrt(a[1]/(1-a[2]-b[1]))) # Eq. 12.43
for(i in 2:n){ # gera processo GARCH(1,1)
h[i] <- a[1] + a[2]*y[i-1]^2 + b[1]*h[i-1] # Eq. 12.42
y[i] <- e[i]*sqrt(h[i]) # Eqs. 12.37 e 12.38
}
y <- ts(y[101:n]) # descartando as 100 primeiras observações
fit_garch_tseries <- tseries::garch(y, order = c(1,1), trace = FALSE) # ajusta GARCH(1,1)
summary(fit_garch_tseries) # testes de diagnóstico
##
## Call:
## tseries::garch(x = y, order = c(1, 1), trace = FALSE)
##
## Model:
## GARCH(1,1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.01545 -0.70722 -0.02598 0.63435 3.22918
##
## Coefficient(s):
## Estimate Std. Error t value Pr(>|t|)
## a0 0.07349 0.01238 5.935 2.94e-09 ***
## a1 0.47695 0.05420 8.801 < 2e-16 ***
## b1 0.39325 0.04371 8.997 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Diagnostic Tests:
## Jarque Bera Test
##
## data: Residuals
## X-squared = 0.6907, df = 2, p-value = 0.708
##
##
## Box-Ljung test
##
## data: Squared.Residuals
## X-squared = 0.77497, df = 1, p-value = 0.3787
A função bayesforecast::garch()
(Alonzo and Cruz 2020) implementa um modelo GARCH\((s,k,h)\) baseado em (Ardia and Hoogerheide 2010) e (Fonseca et al. 2019). As três componentes \((s,k,h) \equiv (s,r,h)\) são a ordem do ARCH (\(s\)), a ordem do GARCH (\(k \equiv r\)) e a ordem do MGARCH (\(h\)). Para mais detalhes veja (Engle Robert and Kevin 2001) e (R. Engle 2002).
library(tseries)
library(bayesforecast)
n <- 1100 # número de observações simuladas, 100 primeiras descartadas
a <- c(0.1, 0.5) # coeficientes ARCH(1): alfa_0 e alfa_1
b <- c(0.3) # coeficiente GARCH(s,1): beta_1
set.seed(42); e <- rnorm(n) # \varepsilon_t, Eq. 12.40
y <- double(n) # vetor para conter a série temporal y simulada (double-precision vector)
h <- double(n) # vetor para conter a série temporal de ht
set.seed(314); y[1] <- rnorm(1, sd = sqrt(a[1]/(1-a[2]-b[1]))) # Eq. 12.43
for(i in 2:n){ # gera processo GARCH(1,1)
h[i] <- a[1] + a[2]*y[i-1]^2 + b[1]*h[i-1] # Eq. 12.42
y[i] <- e[i]*sqrt(h[i]) # Eqs. 12.37 e 12.38
}
y <- ts(y[101:n]) # descartando as 100 primeiras observações
fit_garch_bayes <- bayesforecast::garch(y, order = c(1,1,0)) # ajusta GARCH(1,1)
fit_garch_varstan <- bayesforecast::varstan(fit_garch_bayes, iter = 500, chains = 1)
##
## SAMPLING FOR MODEL 'tgarch' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.000912 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 9.12 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 500 [ 0%] (Warmup)
## Chain 1: Iteration: 50 / 500 [ 10%] (Warmup)
## Chain 1: Iteration: 100 / 500 [ 20%] (Warmup)
## Chain 1: Iteration: 150 / 500 [ 30%] (Warmup)
## Chain 1: Iteration: 200 / 500 [ 40%] (Warmup)
## Chain 1: Iteration: 250 / 500 [ 50%] (Warmup)
## Chain 1: Iteration: 251 / 500 [ 50%] (Sampling)
## Chain 1: Iteration: 300 / 500 [ 60%] (Sampling)
## Chain 1: Iteration: 350 / 500 [ 70%] (Sampling)
## Chain 1: Iteration: 400 / 500 [ 80%] (Sampling)
## Chain 1: Iteration: 450 / 500 [ 90%] (Sampling)
## Chain 1: Iteration: 500 / 500 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 1.303 seconds (Warm-up)
## Chain 1: 0.735 seconds (Sampling)
## Chain 1: 2.038 seconds (Total)
## Chain 1:
Exercício 12.29 Considere o modelo GARCH\((r,s)\) e GARCH\((s,r,h)\).
- No Exemplo 12.22 verifique a equivalência entre os modelos ARCH\((2)\), GARCH\((0,2)\) de
tseries
e GARCH\((2,0,0)\) debayesforecast
. - Veja https://asael697.r-universe.dev/bayesforecast.
- Veja https://cran.r-universe.dev/bayesforecast/doc/manual.html.
- Veja https://cran.r-project.org/web/packages/bayesforecast/readme/README.html.
12.13.3 ARMA-GARCH
Exemplo 12.24 Considere a série temporal ‘r unlist(tsdl::meta_tsdl[20,2][[1]][[1]])
’ apresentada por r tsdl::meta_tsdl[20,1]
, disponível em tsdl::tsdl[[20]]
23.
library(forecast)
library(bayesforecast)
# ajustando um modelo ARIMA automático
y <- tsdl::tsdl[[20]]
forecast::auto.arima(y) # ARIMA(2,1,2)
## Series: y
## ARIMA(2,1,2)
##
## Coefficients:
## ar1 ar2 ma1 ma2
## 1.3467 -0.3963 -1.7710 0.8103
## s.e. 0.0303 0.0287 0.0205 0.0194
##
## sigma^2 = 243.8: log likelihood = -11745.5
## AIC=23500.99 AICc=23501.01 BIC=23530.71
# ajustando ARMA-GARCH para a série com uma diferença via bayesforecast
fit_garch_bayes <- bayesforecast::garch(y, order = c(1,1,0), arma = c(2,2)) # ajusta ARMA(2,2)-GARCH(1,1)
fit_garch_varstan <- bayesforecast::varstan(fit_garch_bayes, iter = 500, chains = 1)
##
## SAMPLING FOR MODEL 'tgarch' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.001142 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 11.42 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 500 [ 0%] (Warmup)
## Chain 1: Iteration: 50 / 500 [ 10%] (Warmup)
## Chain 1: Iteration: 100 / 500 [ 20%] (Warmup)
## Chain 1: Iteration: 150 / 500 [ 30%] (Warmup)
## Chain 1: Iteration: 200 / 500 [ 40%] (Warmup)
## Chain 1: Iteration: 250 / 500 [ 50%] (Warmup)
## Chain 1: Iteration: 251 / 500 [ 50%] (Sampling)
## Chain 1: Iteration: 300 / 500 [ 60%] (Sampling)
## Chain 1: Iteration: 350 / 500 [ 70%] (Sampling)
## Chain 1: Iteration: 400 / 500 [ 80%] (Sampling)
## Chain 1: Iteration: 450 / 500 [ 90%] (Sampling)
## Chain 1: Iteration: 500 / 500 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 41.931 seconds (Warm-up)
## Chain 1: 65.134 seconds (Sampling)
## Chain 1: 107.065 seconds (Total)
## Chain 1:
# # tomando uma diferença da série
# ydiff <- diff(y)
# fit_garch_bayes_diff <- bayesforecast::garch(ydiff, order = c(1,1,0), arma = c(2,2)) # ajusta ARMA(2,2)-GARCH(1,1)
# fit_garch_varstan_diff <- bayesforecast::varstan(fit_garch_bayes_diff, iter = 500, chains = 1)
# check_residuals(fit_garch_varstan_diff) # análise de resíduos
# bf <- bayesforecast::forecast(fit_garch_varstan_diff, h = 48)
# # desfazendo as diferenças
# bf$x <- diffinv(bf$x, xi = y[1])
# bf$mean <- diffinv(bf$mean, xi = y[length(y)])
# autoplot(bf) # projetando 48 passos
A função rugarch::ugarchspec()
(Galanos 2023) permite ajustar um modelo ARFIMA-GARCH, com variantes do modelo GARCH e variáveis externas (veja documentação).
# libs
library(forecast)
library(rugarch)
# ajustando um modelo ARIMA automático
y <- tsdl::tsdl[[37]]
(fit_arima <- forecast::auto.arima(y)) # ARIMA(1,1,1)
## Series: y
## ARIMA(1,1,1)
##
## Coefficients:
## ar1 ma1
## -0.5322 0.9476
## s.e. 0.1389 0.0657
##
## sigma^2 = 95.43: log likelihood = -199.28
## AIC=404.57 AICc=405.05 BIC=410.54
# avaliando a série ao quadrado e os resíduos ao quadrado do ARIMA(1,1,1)
par(mfrow=c(2,2))
acf(y^2)
pacf(y^2)
acf(fit_arima$residuals^2)
pacf(fit_arima$residuals^2)
# tomando uma diferença da série
ydiff <- diff(y)
# ajustando ARMA-GARCH para a série com uma diferença via rugarch
spec <- ugarchspec(
mean.model = list(armaOrder = c(2,2), arfima = FALSE),
variance.model = list(garchOrder = c(0,1))
)
fit_garch <- ugarchfit(spec = spec, data = y, solver.control = list(trace=0))
fit_garch
##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(0,1)
## Mean Model : ARFIMA(2,0,2)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 9.766902 0.729497 13.38854 0.000000
## ar1 0.973896 0.396518 2.45612 0.014045
## ar2 0.315427 0.525500 0.60024 0.548345
## ma1 -0.033163 0.048435 -0.68469 0.493542
## ma2 -1.091586 0.082842 -13.17667 0.000000
## omega 0.744145 5.858488 0.12702 0.898925
## beta1 0.999000 0.082041 12.17687 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 9.766902 1.68928 5.781703 0.000000
## ar1 0.973896 1.37412 0.708740 0.478486
## ar2 0.315427 1.83694 0.171713 0.863663
## ma1 -0.033163 0.15564 -0.213077 0.831267
## ma2 -1.091586 0.29991 -3.639651 0.000273
## omega 0.744145 23.03245 0.032309 0.974226
## beta1 0.999000 0.32684 3.056503 0.002239
##
## LogLikelihood : -190.7872
##
## Information Criteria
## ------------------------------------
##
## Akaike 7.1923
## Bayes 7.4477
## Shibata 7.1645
## Hannan-Quinn 7.2911
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.5043 0.4776
## Lag[2*(p+q)+(p+q)-1][11] 4.4613 0.9972
## Lag[4*(p+q)+(p+q)-1][19] 8.3803 0.7387
## d.o.f=4
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 6.099 0.01352
## Lag[2*(p+q)+(p+q)-1][2] 7.205 0.01052
## Lag[4*(p+q)+(p+q)-1][5] 9.123 0.01541
## d.o.f=1
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[2] 2.055 0.500 2.000 0.1517
## ARCH Lag[4] 3.418 1.397 1.611 0.2099
## ARCH Lag[6] 3.773 2.222 1.500 0.3403
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 13.3382
## Individual Statistics:
## mu 0.03699
## ar1 0.03213
## ar2 0.03141
## ma1 0.03472
## ma2 0.06053
## omega 0.58638
## beta1 0.46955
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.69 1.9 2.35
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.06206 0.950761
## Negative Sign Bias 3.02336 0.003937 ***
## Positive Sign Bias 1.64254 0.106754
## Joint Effect 14.34155 0.002475 ***
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 17.36 0.5652
## 2 30 25.18 0.6688
## 3 40 30.09 0.8463
## 4 50 44.09 0.6720
##
##
## Elapsed time : 0.139035
# previsão com o modelo
fcast <- ugarchforecast(fit_garch)
# plot(fcast, ask = FALSE)
# desfazendo a diferença
# fcast2 <- fcast
# fcast2@model$modeldata$data <- as.matrix(diffinv(fcast2@model$modeldata$data, xi = y[1]))
# fcast2@model$modeldata$data <- y[-1]
# fcast2@forecast$seriesFor <- as.matrix(diffinv(fcast2@forecast$seriesFor)[-1])
# plot(fcast2)