Abstract
This is an undergrad student level instruction for class use. The reader is encouraged to enhance the model changing the train data.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
Sugestão de citação: FIGUEIREDO, Adriano Marcos Rodrigues. Séries Temporais: Estimação do método theta com série diária de ação LAME4. Campo Grande-MS,Brasil: RStudio/Rpubs, 2021. Disponível em <http://rpubs.com/amrofi/acoes_lame4_theta>.
Neste exercício, usarei a série de preços da ação das Lojas Americanas (LAME4).
Vou chamar os dados por uma função derivada do quantmod que está no tidyquant, tq_get
. Veja que ele vem em tibble. Então geramos o tsibble
. Depois preenchemos os gaps implicitos com o fill_gaps.
library(tidyquant)
lame.tbl<-tq_get("LAME4.SA")
#XTS PARA TSIBBLE
library(tsibble)
library(tsibbledata)
library(magrittr)
library(dplyr)
LAME_tbl_ts<-as_tsibble(lame.tbl)
head(LAME_tbl_ts)
fabletools::autoplot(LAME_tbl_ts, close)
# fill gaps com ultimo valor ----------
library(tidyr)
dados.full2 <- LAME_tbl_ts[-1] %>%
tsibble::fill_gaps() %>%
tidyr::fill(c(open,high,low,close,volume,adjusted), .direction = "down")
print(dados.full2)
# A tsibble: 3,806 x 7 [1D]
date open high low close volume adjusted
<date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 2011-01-03 7.90 7.94 7.81 7.85 2735138 7.33
2 2011-01-04 7.90 7.91 7.69 7.80 5385557 7.28
3 2011-01-05 7.80 8.06 7.74 8.06 8986715 7.52
4 2011-01-06 8.06 8.07 7.92 7.93 5828372 7.40
5 2011-01-07 7.94 8.03 7.91 7.99 3071662 7.45
6 2011-01-08 7.94 8.03 7.91 7.99 3071662 7.45
7 2011-01-09 7.94 8.03 7.91 7.99 3071662 7.45
8 2011-01-10 7.73 7.93 7.70 7.92 6060762 7.39
9 2011-01-11 7.93 7.98 7.80 7.98 5422818 7.45
10 2011-01-12 7.95 8.00 7.86 7.93 4718197 7.41
# ... with 3,796 more rows
O gráfico para o preço de fechamento (close) será:
fabletools::autoplot(dados.full2, close)
O leitor pode testar agora as ferramentas de forecast do fable e fabletools.
library(fpp3)
fit <- dados.full2 %>%
model(
arima = ARIMA(close),
ets = ETS(close),
theta = THETA(close))
fit_fc <- fit %>% forecast(h = 365)
fit_fc %>% autoplot(dados.full2, level = NULL) +
labs(y = "valor de X", title = "Exemplo de forecast após fill_gaps")
prophet
Agora farei um ajuste com o Facebook Prophet:
# Prophet model ------------------
# seguindo https://rpubs.com/amrofi/prophet_bitcoin
library(prophet)
library(tidyverse)
dados.ph<-dados.full2[,c(1,5)] # peguei apenas a date e a close
colnames(dados.ph)<-c("ds","y") # o prophet requer esses nomes de colunas
Model1 <- prophet(dados.ph,daily.seasonality=TRUE)
class(Model1) # objeto do prophet
[1] "prophet" "list"
## Criar dataframe com mais 365 obs a frente (fora da amostra) ---------
Future1 <- make_future_dataframe(Model1, periods = 365)
tail(Future1)
## Previsao pelo Model1 para as datas de Future1 ------------
Forecast1 <- predict(Model1, Future1)
class(Forecast1)
[1] "data.frame"
## gráficos do forecast -------------
dyplot.prophet(Model1, Forecast1)
Fica a dica para o leitor avaliar a presença de efeitos ARCH.
fable.prophet
Neste caso, conforme Hyndman e Athanasopoulos (2020), capítulo 12.2, o modelo Facebook (TAYLOR e LETHAM, 2018) é útil para estimar séries diárias, e principalmente aquelas com sazonalidades semanal e anual, além de efeitos de feriados. Ele pode ser considerado um modelo de regressão não linear da forma:
\[ y_t = g(t) + s(t) + h(t) + \varepsilon_t, \]
em que \(g(t)\) descreve a tendência linear (ou “termo de crescimento”), \(s(t)\) descreve os vários padrões sazonais, \(h(t)\) captura os efeitos do feriado, e \(\varepsilon_t\) é um termo de erro de ruído branco. A tendência linear é selecionada automaticamente se não for especificada explicitamente. O componente sazonal consiste em termos de Fourier dos períodos relevantes. Por padrão, a ordem 10 é usada para sazonalidade anual e a ordem 3 é usada para sazonalidade semanal, conforme Hyndman e Athanasopoulos (2020, tradução livre). Os efeitos de feriado são adicionados como variáveis dummies simples. O modelo é estimado usando uma abordagem bayesiana para permitir a seleção automática dos pontos de mudança e outras características do modelo.
Estimaremos juntamente aos modelos ARIMA
, ETS
, THETA
anteriormente estimados, dentro do ambiente fable
e fable.prophet
. Farei a regra 80-20 para amostras treino (3806*0.80 = 3044 dias = “2019-05-04”) e teste (3806-3044 = 762 dias = após “2019-05-05”), obtido a partir do dataset completo para os gaps implicitos.
library(fable.prophet);library(fable)
train <- dados.full2 %>% filter_index(~"2019-05-04") # 3806*0.80 = 3044
h=length(dados.full2$close)-length(train$close)
h
[1] 762
fit.all <- train %>%
model(
arima = ARIMA(close,stepwise=FALSE,approximation = FALSE),
ets = ETS(close),
theta = THETA(close),
prophet = prophet(close ~
growth("linear") +
season(period = "day", order = 10) +
season(period = "week", order = 5) +
season(period = "year", order = 3))
)%>%
mutate(combination = (ets + arima+theta+prophet) / 4)
fitall_fc <- fit.all %>% forecast(h = h)
fitall_fc %>% autoplot(dados.full2, level = NULL) +
labs(y = "valor da ação", title = "Forecast com ARIMA, ETS, THETA e FABLE.PROPHET, para o Varejo de MS")
fc_accuracy <-fitall_fc %>% accuracy(dados.full2)
# acurácia ordenada pelo RMSE
fc_accuracy %>% group_by(.model) %>% summarise(RMSE = mean(RMSE), MAE = mean(MAE),
MASE = mean(MASE)) %>% arrange(RMSE)
Forecasts:
library(fable);library(fable.prophet)
fitall_fc2 <- fit.all %>% forecast(h = h+60)
fitall_fc2 %>% autoplot(dados.full2, level = c(95)) +
labs(y = "valor da ação", title = "Forecast com ARIMA, ETS, THETA e FABLE.PROPHET, para o Varejo de MS")
Representarei apenas os forecasts do período de “2021-06-01” a “2021-06-15” para não ficar muito extenso.
fitall_fc2 %>%
filter_index("2021-06-01" ~ "2021-06-15") %>% knitr::kable()
.model | date | close | .mean |
---|---|---|---|
arima | 2021-06-01 | N(16, 41) | 15.59990 |
arima | 2021-06-02 | N(16, 41) | 15.59990 |
arima | 2021-06-03 | N(16, 41) | 15.59990 |
arima | 2021-06-04 | N(16, 41) | 15.59990 |
arima | 2021-06-05 | N(16, 41) | 15.59990 |
arima | 2021-06-06 | N(16, 41) | 15.59990 |
arima | 2021-06-07 | N(16, 41) | 15.59990 |
arima | 2021-06-08 | N(16, 41) | 15.59990 |
arima | 2021-06-09 | N(16, 41) | 15.59990 |
arima | 2021-06-10 | N(16, 41) | 15.59990 |
arima | 2021-06-11 | N(16, 41) | 15.59990 |
arima | 2021-06-12 | N(16, 42) | 15.59990 |
arima | 2021-06-13 | N(16, 42) | 15.59990 |
arima | 2021-06-14 | N(16, 42) | 15.59991 |
arima | 2021-06-15 | N(16, 42) | 15.59991 |
ets | 2021-06-01 | N(16, 63) | 15.64643 |
ets | 2021-06-02 | N(16, 63) | 15.64643 |
ets | 2021-06-03 | N(16, 63) | 15.64643 |
ets | 2021-06-04 | N(16, 63) | 15.64643 |
ets | 2021-06-05 | N(16, 63) | 15.64643 |
ets | 2021-06-06 | N(16, 63) | 15.64643 |
ets | 2021-06-07 | N(16, 63) | 15.64643 |
ets | 2021-06-08 | N(16, 63) | 15.64643 |
ets | 2021-06-09 | N(16, 63) | 15.64643 |
ets | 2021-06-10 | N(16, 64) | 15.64643 |
ets | 2021-06-11 | N(16, 64) | 15.64643 |
ets | 2021-06-12 | N(16, 64) | 15.64643 |
ets | 2021-06-13 | N(16, 64) | 15.64643 |
ets | 2021-06-14 | N(16, 64) | 15.64643 |
ets | 2021-06-15 | N(16, 64) | 15.64643 |
theta | 2021-06-01 | N(17, 45) | 17.15431 |
theta | 2021-06-02 | N(17, 45) | 17.16208 |
theta | 2021-06-03 | N(17, 45) | 17.16243 |
theta | 2021-06-04 | N(17, 45) | 17.16500 |
theta | 2021-06-05 | N(17, 45) | 17.16290 |
theta | 2021-06-06 | N(17, 45) | 17.16144 |
theta | 2021-06-07 | N(17, 45) | 17.16945 |
theta | 2021-06-08 | N(17, 45) | 17.16822 |
theta | 2021-06-09 | N(17, 45) | 17.17600 |
theta | 2021-06-10 | N(17, 45) | 17.17635 |
theta | 2021-06-11 | N(17, 45) | 17.17892 |
theta | 2021-06-12 | N(17, 45) | 17.17681 |
theta | 2021-06-13 | N(17, 46) | 17.17535 |
theta | 2021-06-14 | N(17, 46) | 17.18337 |
theta | 2021-06-15 | N(17, 46) | 17.18213 |
prophet | 2021-06-01 | sample[5000] | 20.27055 |
prophet | 2021-06-02 | sample[5000] | 20.28071 |
prophet | 2021-06-03 | sample[5000] | 20.26173 |
prophet | 2021-06-04 | sample[5000] | 20.26027 |
prophet | 2021-06-05 | sample[5000] | 20.21970 |
prophet | 2021-06-06 | sample[5000] | 20.23131 |
prophet | 2021-06-07 | sample[5000] | 20.19852 |
prophet | 2021-06-08 | sample[5000] | 20.21045 |
prophet | 2021-06-09 | sample[5000] | 20.25915 |
prophet | 2021-06-10 | sample[5000] | 20.18836 |
prophet | 2021-06-11 | sample[5000] | 20.20295 |
prophet | 2021-06-12 | sample[5000] | 20.18492 |
prophet | 2021-06-13 | sample[5000] | 20.19415 |
prophet | 2021-06-14 | sample[5000] | 20.20745 |
prophet | 2021-06-15 | sample[5000] | 20.17323 |
combination | 2021-06-01 | 17.21997 | 17.21997 |
combination | 2021-06-02 | 17.22646 | 17.22646 |
combination | 2021-06-03 | 17.21925 | 17.21925 |
combination | 2021-06-04 | 17.21928 | 17.21928 |
combination | 2021-06-05 | 17.21497 | 17.21497 |
combination | 2021-06-06 | 17.21357 | 17.21357 |
combination | 2021-06-07 | 17.21329 | 17.21329 |
combination | 2021-06-08 | 17.21882 | 17.21882 |
combination | 2021-06-09 | 17.2123 | 17.21230 |
combination | 2021-06-10 | 17.21061 | 17.21061 |
combination | 2021-06-11 | 17.21338 | 17.21338 |
combination | 2021-06-12 | 17.21147 | 17.21147 |
combination | 2021-06-13 | 17.20688 | 17.20688 |
combination | 2021-06-14 | 17.21332 | 17.21332 |
combination | 2021-06-15 | 17.20762 | 17.20762 |
Como podemos ver, as previsões ainda podem ser melhores, e podemos preocupar com os valores atípicos ao final da série.
Neste caso, supondo a série completa, estimaremos tudo junto (ETS, ARIMA, THETA, PROPHET, Combination) sem fazer o split de treino e teste.
fit.semsplit <- dados.full2 %>%
model(
arima = ARIMA(close,stepwise=FALSE,approximation = FALSE),
ets = ETS(close),
theta = THETA(close),
prophet = prophet(close ~
growth("linear") +
season(period = "day", order = 10) +
season(period = "week", order = 5) +
season(period = "year", order = 3))
)%>%
mutate(combination = (ets + arima+theta+prophet) / 4)
Colocarei um forecast de 150 dias a frente (para melhorar a visualização gráfica).
fitsem_fc <- fit.semsplit %>% forecast(h = 150)
fitsem_fc %>% autoplot(dados.full2, level = c(95)) +
labs(y = "valor da ação", title = "Forecast com ARIMA, ETS, THETA e FABLE.PROPHET, para o Varejo de MS")
Interessante observar que a previsão com o dataset completo levou o prophet a prever muito acima do nível da última observação. No caso, considerando a amostra de treino, o forecast estava mais perto do nível recente, embora pouco abaixo. Nessas situações é que a combinação aparece como interessante, mesclando as estimativas dos diferentes modelos.
Deixo ao leitor para investigar a estabilidade do modelo e propriedades dos resíduos.
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.
TAYLOR, S. J.; LETHAM, B. (2018). Forecasting at scale. The American Statistician, 72(1), 37–45.