Séries Temporais: Rotina para GARCH com o pacote garchmodels

Abstract

This is an undergrad student level instruction for class use.

Licença

This work is licensed under the Creative Commons Attribution-ShareAlike 4.0 International License. To view a copy of this license, visit http://creativecommons.org/licenses/by-sa/4.0/ or send a letter to Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.

License: CC BY-SA 4.0

Citação

Sugestão de citação: FIGUEIREDO, Adriano Marcos Rodrigues. Séries Temporais: Rotina para GARCH com o pacote garchmodels. Campo Grande-MS,Brasil: RStudio/Rpubs, 2021. Disponível em <http://rpubs.com/amrofi/garchmodels_exemplo>.

1 Introdução

Este exemplo tem por base a vinheta explicativa do pacote garchmodels, adaptado de Almuiña (2021a,b). A sugestão é utilizar o arcabouço do timetk e o garchmodels invés das alternativas via fable,forecast,rugarch.

2 Pacotes

library(garchmodels)
library(timetk)
library(tidyverse)
library(tidymodels)

O exemplo original usa um dataset previamente disponibilizado no pacote garchmodels, a saber, o de retornos da IBM: rIBM. São retornos diários da empresa IBM, 3523 linhas e 2 variáveis em formato tibble. É possível ver que apresentam valores positivos e negativos, e uma coluna date, de 2007-01-03 a 2020-12-29, com gaps implícitos para dias não negociados.

3 Organização dos Dados

rIBM<-rIBM
rIBM

tail(rIBM)

Modificarei a dataset para o formato longo, e dividiremos em amostra treino e estenderemos este dataset ao longo do horizonte de previsão. O dataset rIBM_extended tem 3526 observações, ou seja 3 dias a mais, os quais estão em rIBM_future e o treino com 3523 observações em rIBM_train. O leitor pode verificar o chunk em que cria 3 observacoes por meio de timetk::future_frame, depois, o código que divide e cria a série rIBM_train, e o que filtra e cria o dataset rIBM_future.

set.seed(123456)
rIBM_extended <- rIBM %>%
    future_frame(.length_out = 3, .bind_data = TRUE) 

rIBM_train  <- rIBM_extended %>% drop_na()
rIBM_future <- rIBM_extended %>% filter(is.na(daily_returns))

3.1 Teste ARCH

# teste para efeitos ARCH no {feasts}
feasts::stat_arch_lm(rIBM_train$daily_returns, lags = 30)

stat_arch_lm 
   0.1600181 

4 Especificação do modelo GARCH

O próximo passo é especificar o modelo garch. Marcaremos quais parâmetros vamos ajustar. Também é necessário usar o argumento tune_by em que devemos especificar sigmaFor ou seriesFor para nos referirmos às previsões de série (series forecast) ou previsões da variância (sigma forecast). Com base nessas previsões obtidas no processo de ajuste, é que as métricas serão calculadas. Especificamos o modo de regressão com ajuste da ordem da parte do ARCH (nível da série) e da parte do GARCH (variância). O mecanismo de estimação usará o pacote `rugarch` em background.

# Model Spec
model_spec <-garchmodels::garch_reg(mode = "regression",
                                    arch_order = tune(),
                                    garch_order = tune(),
                                    tune_by = "sigmaFor") %>%
    set_engine("rugarch") 

A próxima etapa é criar uma receita (recipe), uma regra na qual especificamos a fórmula que usaremos. Especificamos os retornos diários com base na data, para o dataset treino:

# Recipe spec
recipe_spec <- recipe(daily_returns ~ date, data = rIBM_train) 

Reunimos tudo em um fluxo de trabalho:

# Workflow
wflw <- workflow() %>%
    add_recipe(recipe_spec) %>%
    add_model(model_spec)

5 Modelo univariado

A seguir, criamos o modelo GARCH para a variância:

model_fit_garch <-garchmodels::garch_reg(mode = "regression",
                                         arch_order = 1,
                                         garch_order = 1,
                                         ma_order = 0,
                                         ar_order = 0) %>%
                  set_engine("rugarch", mean.model = list(include.mean = FALSE)) %>%
                  fit(daily_returns ~ date, data = rIBM_train)
print(model_fit_garch)

parsnip model object

Fit time:  891ms 

*---------------------------------*
*          GARCH Model Fit        *
*---------------------------------*

Conditional Variance Dynamics   
-----------------------------------
GARCH Model : sGARCH(1,1)
Mean Model  : ARFIMA(0,0,0)
Distribution    : norm 

Optimal Parameters
------------------------------------
        Estimate  Std. Error  t value Pr(>|t|)
omega   0.000014    0.000001   21.001        0
alpha1  0.115846    0.006528   17.745        0
beta1   0.818369    0.009052   90.408        0

Robust Standard Errors:
        Estimate  Std. Error  t value Pr(>|t|)
omega   0.000014    0.000002   7.8280        0
alpha1  0.115846    0.015520   7.4641        0
beta1   0.818369    0.025162  32.5246        0

LogLikelihood : 10191.5 

Information Criteria
------------------------------------
                    
Akaike       -5.7840
Bayes        -5.7787
Shibata      -5.7840
Hannan-Quinn -5.7821

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                        statistic p-value
Lag[1]                      1.427  0.2323
Lag[2*(p+q)+(p+q)-1][2]     1.640  0.3301
Lag[4*(p+q)+(p+q)-1][5]     2.623  0.4802
d.o.f=0
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.3874  0.5337
Lag[2*(p+q)+(p+q)-1][5]    1.5028  0.7392
Lag[4*(p+q)+(p+q)-1][9]    2.4485  0.8453
d.o.f=2

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale P-Value
ARCH Lag[3]    0.9991 0.500 2.000  0.3175
ARCH Lag[5]    1.5830 1.440 1.667  0.5708
ARCH Lag[7]    2.0297 2.315 1.543  0.7110

Nyblom stability test
------------------------------------
Joint Statistic:  23.7215
Individual Statistics:             
omega  1.6697
alpha1 0.2575
beta1  0.1887

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:         0.846 1.01 1.35
Individual Statistic:    0.35 0.47 0.75

Sign Bias Test
------------------------------------
                   t-value   prob sig
Sign Bias           0.6816 0.4956    
Negative Sign Bias  0.8871 0.3751    
Positive Sign Bias  0.1652 0.8688    
Joint Effect        3.3802 0.3366    


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20     165.7    1.955e-25
2    30     182.5    3.461e-24
3    40     188.2    1.946e-21
4    50     224.8    2.351e-24


Elapsed time : 0.2044561 

5.1 Plot da série estimada

Representaremos graficamente a série com limites de VaR (Value-at-risk) de 1% e o SD condicional (vs | retorns |):

plot(model_fit_garch$fit$models$model_1, which = 2)

please wait...calculating quantiles...

plot(model_fit_garch$fit$models$model_1, which = 3)

Agora faremos a previsão (predict):

predict(model_fit_garch, rIBM_future)

Vamos agora criar um modelo combinado usando um modelo GARCH para a variância e um modelo ARMA para a média:

model_fit_garch2 <-garchmodels::garch_reg(mode = "regression",
                                         arch_order = 1,
                                         garch_order = 1,
                                         ma_order = 2,
                                         ar_order = 2) %>%
                  set_engine("rugarch") %>%
                  fit(daily_returns ~ date, data = rIBM_train)
print(model_fit_garch2)

parsnip model object

Fit time:  1.2s 

*---------------------------------*
*          GARCH Model Fit        *
*---------------------------------*

Conditional Variance Dynamics   
-----------------------------------
GARCH Model : sGARCH(1,1)
Mean Model  : ARFIMA(2,0,2)
Distribution    : norm 

Optimal Parameters
------------------------------------
        Estimate  Std. Error    t value Pr(>|t|)
mu      0.000362    0.000205     1.7618 0.078099
ar1    -0.271525    0.006396   -42.4506 0.000000
ar2    -0.987467    0.001171  -843.1570 0.000000
ma1     0.267126    0.004876    54.7809 0.000000
ma2     0.992600    0.000076 13021.0202 0.000000
omega   0.000014    0.000001    17.1715 0.000000
alpha1  0.120510    0.006581    18.3114 0.000000
beta1   0.813198    0.009300    87.4383 0.000000

Robust Standard Errors:
        Estimate  Std. Error   t value Pr(>|t|)
mu      0.000362    0.000213    1.7035 0.088466
ar1    -0.271525    0.007100  -38.2423 0.000000
ar2    -0.987467    0.002143 -460.8064 0.000000
ma1     0.267126    0.005348   49.9471 0.000000
ma2     0.992600    0.000102 9702.0988 0.000000
omega   0.000014    0.000002    5.9675 0.000000
alpha1  0.120510    0.017233    6.9930 0.000000
beta1   0.813198    0.027910   29.1368 0.000000

LogLikelihood : 10197.46 

Information Criteria
------------------------------------
                    
Akaike       -5.7845
Bayes        -5.7705
Shibata      -5.7845
Hannan-Quinn -5.7795

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                         statistic p-value
Lag[1]                      0.9806  0.3220
Lag[2*(p+q)+(p+q)-1][11]    5.5126  0.7874
Lag[4*(p+q)+(p+q)-1][19]    8.6984  0.6883
d.o.f=4
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.3362  0.5620
Lag[2*(p+q)+(p+q)-1][5]    1.5925  0.7171
Lag[4*(p+q)+(p+q)-1][9]    2.5995  0.8228
d.o.f=2

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale P-Value
ARCH Lag[3]     1.116 0.500 2.000  0.2907
ARCH Lag[5]     1.699 1.440 1.667  0.5417
ARCH Lag[7]     2.182 2.315 1.543  0.6787

Nyblom stability test
------------------------------------
Joint Statistic:  26.7383
Individual Statistics:              
mu     0.26152
ar1    0.21834
ar2    0.05954
ma1    0.17730
ma2    0.09443
omega  1.53005
alpha1 0.25073
beta1  0.18857

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:         1.89 2.11 2.59
Individual Statistic:    0.35 0.47 0.75

Sign Bias Test
------------------------------------
                   t-value   prob sig
Sign Bias           0.3681 0.7128    
Negative Sign Bias  1.0545 0.2917    
Positive Sign Bias  0.2667 0.7897    
Joint Effect        2.9874 0.3936    


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20     170.9    1.853e-26
2    30     184.7    1.324e-24
3    40     200.9    1.144e-23
4    50     215.1    1.105e-22


Elapsed time : 0.8607299 

Vamos representar a série:

plot(model_fit_garch2$fit$models$model_1, which = 1)
plot(model_fit_garch2$fit$models$model_1, which = "all")

please wait...calculating quantiles...

6 Criando subamostras para validação cruzada

A rotina de estimação considera resamples (reamostragens). O procedimento por timetk::time_series_cv faz a cross-validation da série temporal, produzindo um plano amostral iniciado com a observação mais recente e “rolando”para trás.

O argumento initial indica o número de amostras usadas para análise / modelagem na reamostragem inicial. O argumento assess indica o número de amostras usadas para cada reamostragem de avaliação. O argumento skip indica o número inteiro que indica quantas (se houver) reamostragens adicionais pular para diminuir a quantidade total de pontos de dados na reamostragem de análise.

Por se tratar de um problema de série temporal, a ordenação por data é importante entre observações, portanto, usaremos a função timetk e visualizaremos graficamente estas reamostras criadas.

resamples <- timetk::time_series_cv(rIBM_train, 
                                    date_var = date, 
                                    initial = "6 years", 
                                    assess = "24 months",
                                    skip = "24 months",
                                    cumulative = TRUE,
                                    slice_limit = 3)

timetk::plot_time_series_cv_plan(resamples, .date_var = date, .value = daily_returns)

7 Ajuste fino e estimação

Por fim, usamos a função tune_grid() para aplicar o fluxo de trabalho nas reamostragens de maneira paralela.

tune_results <- tune_grid(
    object     = wflw,
    resamples  = resamples,
    param_info = parameters(wflw),
    grid       = 5,
    control    = control_grid(verbose = TRUE, allow_par = TRUE, parallel_over = "everything")
)

ugarchfilter-->error: parameters names do not match specification
Expected Parameters are: mu ar1 ma1 omega

ugarchfilter-->error: parameters names do not match specification
Expected Parameters are: mu ar1 ma1 omega alpha1 beta1

ugarchfilter-->error: parameters names do not match specification
Expected Parameters are: mu ar1 ma1 omega

ugarchfilter-->error: parameters names do not match specification
Expected Parameters are: mu ar1 ma1 omega

#> 
#> ugarchfilter-->error: parameters names do not match specification
#> Expected Parameters are: mu ar1 ma1 omega
#> 
#> ugarchfilter-->error: parameters names do not match specification
#> Expected Parameters are: mu ar1 ma1 omega alpha1 beta1
#> 
#> ugarchfilter-->error: parameters names do not match specification
#> Expected Parameters are: mu ar1 ma1 omega
#> 
#> ugarchfilter-->error: parameters names do not match specification
#> Expected Parameters are: mu ar1 ma1 omega

Por fim, podemos ver quais são os melhores resultados ordenados pelo RMSE:

tune_results %>% show_best(metric = "rmse")

Esta rotina de ajuste fino acaba sendo interessante para avaliar possibilidades distintas da ordem do ARCH-GARCH.

8 GJR-GARCH

Agora farei a opção com o gjrGARCH:

model_garch_fit3 <-garchmodels::garch_reg(mode = "regression",
                                        arch_order = 2,
                                        garch_order = 2) %>%
    set_engine("rugarch", variance.model = list(model='gjrGARCH', 
                                                garchOrder=c(1,1)),
               mean.model     = list(armaOrder=c(0,0))) %>%
    fit(daily_returns ~ date, data = rIBM_train)

predict(model_garch_fit3, rIBM_future)

plot(model_garch_fit3$fit$models$model_1, which = "all")

please wait...calculating quantiles...

print(model_garch_fit3)

parsnip model object

Fit time:  1.8s 

*---------------------------------*
*          GARCH Model Fit        *
*---------------------------------*

Conditional Variance Dynamics   
-----------------------------------
GARCH Model : gjrGARCH(2,2)
Mean Model  : ARFIMA(1,0,1)
Distribution    : norm 

Optimal Parameters
------------------------------------
        Estimate  Std. Error   t value Pr(>|t|)
mu      0.000243    0.000201  1.207981 0.227055
ar1    -0.354502    0.455085 -0.778978 0.435992
ma1     0.328176    0.459617  0.714019 0.475215
omega   0.000015    0.000001 16.168282 0.000000
alpha1  0.065235    0.026952  2.420441 0.015502
alpha2  0.004847    0.028877  0.167841 0.866708
beta1   0.816588    0.112206  7.277607 0.000000
beta2   0.000003    0.104563  0.000027 0.999978
gamma1  0.102956    0.043392  2.372715 0.017658
gamma2 -0.019374    0.046895 -0.413146 0.679500

Robust Standard Errors:
        Estimate  Std. Error   t value Pr(>|t|)
mu      0.000243    0.000252  0.963378 0.335358
ar1    -0.354502    0.397052 -0.892834 0.371946
ma1     0.328176    0.397529  0.825539 0.409066
omega   0.000015    0.000003  4.921118 0.000001
alpha1  0.065235    0.060895  1.071266 0.284050
alpha2  0.004847    0.053615  0.090399 0.927970
beta1   0.816588    0.179501  4.549222 0.000005
beta2   0.000003    0.169577  0.000017 0.999987
gamma1  0.102956    0.102708  1.002417 0.316142
gamma2 -0.019374    0.098255 -0.197186 0.843682

LogLikelihood : 10200.59 

Information Criteria
------------------------------------
                    
Akaike       -5.7852
Bayes        -5.7677
Shibata      -5.7852
Hannan-Quinn -5.7789

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                        statistic p-value
Lag[1]                   0.002558  0.9597
Lag[2*(p+q)+(p+q)-1][5]  1.075541  1.0000
Lag[4*(p+q)+(p+q)-1][9]  3.805704  0.7357
d.o.f=2
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                         statistic p-value
Lag[1]                      0.5454  0.4602
Lag[2*(p+q)+(p+q)-1][11]    3.3964  0.8187
Lag[4*(p+q)+(p+q)-1][19]    6.2561  0.8527
d.o.f=4

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale P-Value
ARCH Lag[5]    0.6716 0.500 2.000  0.4125
ARCH Lag[7]    1.3861 1.473 1.746  0.6522
ARCH Lag[9]    2.3588 2.402 1.619  0.6883

Nyblom stability test
------------------------------------
Joint Statistic:  299.8073
Individual Statistics:              
mu     0.18874
ar1    0.04978
ma1    0.04765
omega  1.44117
alpha1 0.28401
alpha2 0.29277
beta1  0.21439
beta2  0.21733
gamma1 0.35347
gamma2 0.39375

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:         2.29 2.54 3.05
Individual Statistic:    0.35 0.47 0.75

Sign Bias Test
------------------------------------
                   t-value   prob sig
Sign Bias           0.6748 0.4999    
Negative Sign Bias  0.2144 0.8302    
Positive Sign Bias  0.2464 0.8054    
Joint Effect        0.8974 0.8261    


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20     172.5    9.234e-27
2    30     195.5    1.296e-26
3    40     206.7    1.054e-24
4    50     223.9    3.446e-24


Elapsed time : 1.092078 

Referências

ALMUIÑA, Alberto (2021a). Tuning Univariate Garch Models. Disponível em: https://albertoalmuinha.github.io/garchmodels/articles/tuning_univariate_algorithms.html#example-1. Accessed on 07 Jun 2021.

ALMUIÑA, Alberto (2021b). garchmodels: The ‘Tidymodels’ Extension for GARCH Models. R package version 0.1.1. https://CRAN.R-project.org/package=garchmodels. Accessed on 07 Jun 2021.

ALMUIÑA, Alberto (2021c). Getting Started with garchmodels. Disponível em: https://albertoalmuinha.github.io/garchmodels/articles/getting-started.html. Accessed on 07 Jun 2021.

HYNDMAN, Rob J. (2018). fpp2: Data for “Forecasting: Principles and Practice” (2nd Edition). R package version 2.3. Disponível em: https://CRAN.R-project.org/package=fpp2. Accessed on 20 May 2021.

HYNDMAN, Rob J. (2019). fpp3: Data for “Forecasting: Principles and Practice” (3rd Edition). R package. Disponível em: https://github.com/robjhyndman/fpp3-package, https://OTexts.org/fpp3/. Accessed on 20 May 2021.

HYNDMAN, R.J.; ATHANASOPOULOS, G. (2020) Forecasting: principles and practice, 3rd edition, OTexts: Melbourne, Australia. Disponível em: https://otexts.com/fpp3/. Accessed on 20 May 2021.

O’HARA-WILD, Mitchell; HYNDMAN, Rob J.; WANG, Earo. (2021). feasts: Feature Extraction and Statistics for Time Series. R package version 0.2.1. Disponível em: https://CRAN.R-project.org/package=feasts. Accessed on 20 May 2021.