Séries Temporais com R: Modelo híbrido para BVSP Bovespa - B3

Abstract

This is an undergrad student level instruction for class use. We try to ilustrate time series methods with hybrid models using forecastHybrid package.

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 com R: Modelo híbrido para BVSP Bovespa - B3. Campo Grande-MS,Brasil: RStudio/Rpubs, 2019. Disponível em http://rpubs.com/amrofi/Time_Series_Hybrid.

1 Apresentação

Este artigo procura ilustrar os procedimentos para trabalhar com métodos de séries temporais para modelos híbridos utilizando o pacote forecastHybrid.

1.1 Pacotes

O primeiro passo para tal é baixar os dados da Bolsa B3, antiga Bovespa (BVSP). O pacote quantmod nos auxiliará a obter os dados eletronicamente. O pesquisador deve conhecer os tickers do papel desejado. Alguns pacotes foram retirados do CRAN e precisará pegar e instalar do arquivo (GetHFData e stockPortfolio).

if (!require(quantmod)) {
    install.packages("quantmod")
    library(quantmod)
}
if (!require(fpp2)) {
    install.packages("fpp2")
    library(fpp2)
}
if (!require(GetHFData)) {
    install.packages("GetHFData")
    library(GetHFData)
}
if (!require(forecastHybrid)) {
    install.packages("forecastHybrid")
    library(forecastHybrid)
}

2 Dados Bovespa

quantmod::getSymbols("^BVSP", src = "yahoo")  # pelo {quantmod}

[1] "^BVSP"

chartSeries(BVSP)

2.0.1 Calculando o retorno

# require(stockPortfolio)

# install.packages('quantmod') require(quantmod)

### Bovespa getSymbols('^BVSP',src='yahoo') # pelo {quantmod}
BVSP.Close <- BVSP$BVSP.Close
#### limpeza dos dados
BVSP_clean <- is.na(BVSP.Close) <- BVSP.Close == 0
BVSP_clean <- na.locf(BVSP.Close, fromLast = FALSE, coredata = NULL)
max(BVSP_clean)

[1] 130776

min(BVSP_clean)

[1] 29435

plot(BVSP_clean)

# obtendo o retorno diário
r_BVSP <- dailyReturn(BVSP_clean, type = "log")
plot(r_BVSP)

hist(r_BVSP, 100)

# janela para após 2010
r_BVSP_pos2010 <- window(r_BVSP, start = "2010-01-01")
plot(r_BVSP_pos2010)

hist(r_BVSP_pos2010, 100)

max(r_BVSP_pos2010)

[1] 0.1302228

min(r_BVSP_pos2010)

[1] -0.1599303

plot(window(BVSP_clean, start = "2010-01-01"))

3 Modelo Híbrido forecastHybrid de Shaub e Ellis (2018)

Vou usar apenas a série “limpa” e com janela para após 2010. É um objeto diário com gaps implicitos para os finais de semana (non-traded days).

x <- as.ts(r_BVSP_pos2010)
library(forecastHybrid)
modelorapido <- hybridModel(x)
fcst <- forecast(modelorapido, h = 30)
fcst

     Point Forecast       Lo 80      Hi 80       Lo 95      Hi 95
2848   0.0011424536 -0.02010487 0.02259185 -0.03086524 0.03329822
2849  -0.0001370635 -0.02112048 0.02115043 -0.03187043 0.03191081
2850   0.0003097723 -0.02019024 0.02115067 -0.03095086 0.03191105
2851   0.0005659007 -0.02036992 0.02115092 -0.03113079 0.03191129
2852   0.0002069339 -0.02031911 0.02115116 -0.03108001 0.03191153
2853   0.0004170993 -0.02033015 0.02115140 -0.03109105 0.03191178
2854   0.0003835426 -0.02032731 0.02115164 -0.03108821 0.03191202
2855   0.0001520740 -0.02032796 0.02115189 -0.03108887 0.03191226
2856   0.0002530307 -0.02032780 0.02115213 -0.03108870 0.03191251
2857   0.0004605503 -0.02032784 0.02115237 -0.03108874 0.03191275
2858   0.0002028914 -0.02032783 0.02115261 -0.03108873 0.03191299
2859   0.0002670639 -0.02032783 0.02115286 -0.03108873 0.03191323
2860   0.0002931478 -0.02032783 0.02115310 -0.03108873 0.03191348
2861   0.0002417562 -0.02032783 0.02115334 -0.03108873 0.03191372
2862   0.0002885208 -0.02032783 0.02115358 -0.03108873 0.03191396
2863   0.0002896009 -0.02032783 0.02115383 -0.03108873 0.03191420
2864   0.0002233696 -0.02032783 0.02115407 -0.03108873 0.03191445
2865   0.0002935928 -0.02032783 0.02115431 -0.03108873 0.03191469
2866   0.0002943184 -0.02032783 0.02115455 -0.03108873 0.03191493
2867   0.0002529228 -0.02032783 0.02115480 -0.03108873 0.03191517
 [ reached 'max' / getOption("max.print") -- omitted 10 rows ]

plot(forecast(modelorapido))

A função hybridModel utiliza a opção models para indicar quais os modelos a serem utilizados: - a = auto.arima;
- e = ets;
- n = nnetar;
- s = stlm;
- t = tbats;
- z = snaive.

Assim, ao especificar models = "aet", weights = "equal", estamos definindo que serão utilizados os modelos ARIMA, ETS e TBATS com pesos iguais para cada na realização da combinação dos forecasts.

3.1 Modelo 1

# opção aet indica o modelo hibrido para
hm1 <- hybridModel(y = x, models = "aet", weights = "equal")

# resumo arima
hm1$auto.arima

Series: y 
ARIMA(2,0,0) with zero mean 

Coefficients:
          ar1     ar2
      -0.0850  0.0375
s.e.   0.0187  0.0187

sigma^2 estimated as 0.0002493:  log likelihood=7771.93
AIC=-15537.86   AICc=-15537.85   BIC=-15520

# o pesquisador pode retirar dos demais modelos caso deseje
plot(forecast(hm1$auto.arima))

print(hm1)

Hybrid forecast model comprised of the following models: auto.arima, ets, tbats
############
auto.arima with weight 0.333 
############
ets with weight 0.333 
############
tbats with weight 0.333 

plot(hm1, type = "fit")

plot(hm1, type = "models")

Veja que nesse caso, o modelo “automático” levou a um ARIMA (0,0,0). A função hybridModel com opção weights = "equal" levou aos pesos 0.33 para cada modelo por serem 3 modelos (ARIMA, ETS e TBATS).

3.2 Modelo 2

hm2 <- hybridModel(y = x, models = "aet", weights = "equal")

hm2$auto.arima

Series: y 
ARIMA(2,0,0) with zero mean 

Coefficients:
          ar1     ar2
      -0.0850  0.0375
s.e.   0.0187  0.0187

sigma^2 estimated as 0.0002493:  log likelihood=7771.93
AIC=-15537.86   AICc=-15537.85   BIC=-15520

plot(forecast(hm2$auto.arima))

print(hm2)

Hybrid forecast model comprised of the following models: auto.arima, ets, tbats
############
auto.arima with weight 0.333 
############
ets with weight 0.333 
############
tbats with weight 0.333 

plot(hm2, type = "fit")

plot(hm2, type = "models")

# hm2$ets

etspetr <- ets(x)
summary(etspetr)

ETS(A,N,N) 

Call:
 ets(y = x) 

  Smoothing parameters:
    alpha = 1e-04 

  Initial states:
    l = 2e-04 

  sigma:  0.0159

      AIC      AICc       BIC 
-946.1829 -946.1745 -928.3209 

Training set error measures:
                       ME       RMSE        MAE  MPE MAPE      MASE        ACF1
Training set 3.217528e-05 0.01585554 0.01121434 -Inf  Inf 0.6826243 -0.08846717

arimabvsp <- auto.arima(x, stepwise = FALSE, approximation = FALSE)
summary(arimabvsp)

Series: x 
ARIMA(2,0,2) with non-zero mean 

Coefficients:
          ar1      ar2     ma1     ma2   mean
      -1.4211  -0.6426  1.3303  0.5598  5e-04
s.e.   0.0255   0.0762  0.0271  0.0804  3e-04

sigma^2 estimated as 0.0002473:  log likelihood=7784.79
AIC=-15557.57   AICc=-15557.54   BIC=-15521.85

Training set error measures:
                        ME       RMSE        MAE MPE MAPE      MASE       ACF1
Training set -0.0003308969 0.01571207 0.01119238 NaN  Inf 0.6812878 0.01398369

plot(hm2, type = "fit", ggplot = TRUE)

3.3 Modelo 3 - pesos do modelo

Veja que agora o modelo obteve pesos automaticamente e por acaso foram 0.33 cada.

hm3 <- hybridModel(y = x, models = "aet", weights = "insample.errors", errorMethod = "MASE")
hm3

Hybrid forecast model comprised of the following models: auto.arima, ets, tbats
############
auto.arima with weight 0.333 
############
ets with weight 0.333 
############
tbats with weight 0.334 

hm3$auto.arima

Series: y 
ARIMA(2,0,0) with zero mean 

Coefficients:
          ar1     ar2
      -0.0850  0.0375
s.e.   0.0187  0.0187

sigma^2 estimated as 0.0002493:  log likelihood=7771.93
AIC=-15537.86   AICc=-15537.85   BIC=-15520

# hm3$ets # suprimi saida apenas para reduzir espaço hm3$tbats # suprimi saida
# apenas para reduzir espaço
tbats(x)

BATS(1, {2,0}, -, -)

Call: tbats(y = x)

Parameters
  Alpha: -0.0002107319
  AR coefficients: -0.085003 0.037504

Seed States:
             [,1]
[1,] 0.0003585277
[2,] 0.0000000000
[3,] 0.0000000000

Sigma: 0.0157815
AIC: -966.835

3.4 Modelo 4 - pesos do usuário

Recorde os pesos obtidos no modelo 2:

hm2$weights

auto.arima        ets      tbats 
 0.3333333  0.3333333  0.3333333 

Vou alterar os pesos aleatoriamente treinando com outros pesos fornecidos pelo usuario.

novospesos <- c(0.4, 0.4, 0.2)
names(novospesos) <- c("auto.arima", "ets", "tbats")
hm2$weights <- novospesos
hm2

Hybrid forecast model comprised of the following models: auto.arima, ets, tbats
############
auto.arima with weight 0.4 
############
ets with weight 0.4 
############
tbats with weight 0.2 

accuracy(hm2)

                    ME       RMSE        MAE MPE MAPE        ACF1 Theil's U
Test set -1.370438e-05 0.01578899 0.01119625 NaN  Inf -0.02915247       NaN

accuracy(hm2, individual = TRUE)

$auto.arima
                       ME       RMSE        MAE MPE MAPE      MASE
Training set 0.0002278112 0.01578324 0.01120936 NaN  Inf 0.6823217
                      ACF1
Training set -0.0003974951

$ets
                       ME       RMSE        MAE  MPE MAPE      MASE        ACF1
Training set 3.217528e-05 0.01585554 0.01121434 -Inf  Inf 0.6826243 -0.08846717

$tbats
                        ME      RMSE        MAE MPE MAPE      MASE
Training set -0.0003010996 0.0157815 0.01120292 NaN  Inf 0.6819296
                      ACF1
Training set -0.0005892115

hForecast <- forecast(hm2, h = 180)
plot(hForecast)

hForecast

     Point Forecast       Lo 80      Hi 80       Lo 95      Hi 95
2848   0.0012865785 -0.02010487 0.02259185 -0.03086524 0.03329822
2849  -0.0002974781 -0.02112048 0.02054880 -0.03187043 0.03130918
2850   0.0002777778 -0.02019024 0.02098895 -0.03095086 0.03174469
2851   0.0001694804 -0.02036992 0.02080635 -0.03113079 0.03156237
2852   0.0002002562 -0.02031911 0.02085854 -0.03108001 0.03161458
2853   0.0001935793 -0.02033015 0.02084725 -0.03109105 0.03160329
2854   0.0001953008 -0.02032731 0.02085017 -0.03108821 0.03160621
2855   0.0001949041 -0.02032796 0.02084949 -0.03108887 0.03160554
2856   0.0001950024 -0.02032780 0.02084966 -0.03108870 0.03160570
2857   0.0001949792 -0.02032784 0.02084962 -0.03108874 0.03160566
2858   0.0001949848 -0.02032783 0.02084963 -0.03108873 0.03160567
2859   0.0001949835 -0.02032783 0.02084963 -0.03108873 0.03160567
2860   0.0001949838 -0.02032783 0.02084963 -0.03108873 0.03160567
2861   0.0001949837 -0.02032783 0.02084963 -0.03108873 0.03160567
2862   0.0001949837 -0.02032783 0.02084963 -0.03108873 0.03160567
2863   0.0001949837 -0.02032783 0.02084963 -0.03108873 0.03160568
2864   0.0001949837 -0.02032783 0.02084963 -0.03108873 0.03160568
2865   0.0001949837 -0.02032783 0.02084963 -0.03108873 0.03160568
2866   0.0001949837 -0.02032783 0.02084963 -0.03108873 0.03160568
2867   0.0001949837 -0.02032783 0.02084963 -0.03108873 0.03160568
 [ reached 'max' / getOption("max.print") -- omitted 160 rows ]

3.5 Modelo 4 - passando argumentos ao arima

# aet = auto.arima, ets, tbats
hm4 <- hybridModel(y = x, models = "aet", a.args = list(stepwise = FALSE, approximation = FALSE),
    weights = "insample.errors", errorMethod = "MASE")
accuracy(hm4, individual = TRUE)

$auto.arima
                        ME       RMSE        MAE MPE MAPE      MASE       ACF1
Training set -0.0003308969 0.01571207 0.01119238 NaN  Inf 0.6812878 0.01398369

$ets
                       ME       RMSE        MAE  MPE MAPE      MASE        ACF1
Training set 3.217528e-05 0.01585554 0.01121434 -Inf  Inf 0.6826243 -0.08846717

$tbats
                        ME      RMSE        MAE MPE MAPE      MASE
Training set -0.0003010996 0.0157815 0.01120292 NaN  Inf 0.6819296
                      ACF1
Training set -0.0005892115

accuracy(hm4)

                    ME      RMSE        MAE MPE MAPE        ACF1 Theil's U
Test set -0.0002000603 0.0157584 0.01118714 NaN  Inf -0.02371715       NaN

hm4$weights

auto.arima        ets      tbats 
 0.3336558  0.3330025  0.3333417 

Referências

HYNDMAN, Rob J. R packages for forecast combinations. 2016. Disponível em:https://robjhyndman.com/hyndsight/forecast-combinations/

HYNDMAN, Rob J.; ATHANASOPOULOS, George . Forecasting: principles and practice. 2nd ed. Otexts, 2018. Disponível em: https://otexts.org/fpp2/.

SHAUB, David; ELLIS, Peter. forecastHybrid: Convenient Functions for Ensemble Time Series Forecasts. R package version 3.0.14. 2018. Disponível em: https://CRAN.R-project.org/package=forecastHybrid.

SHAUB, David; ELLIS, Peter.forecastHybrid: vignette. 2016. Disponível em: https://cran.r-project.org/web/packages/forecastHybrid/vignettes/forecastHybrid.html.