O presente trabalho apresenta uma aplicação de um modelo SARIMA na série temporal de bilheteria total mensal doméstica nos EUA. Os dados são obtidos via webscrapping do site box office mojo, disponível em: https://www.boxofficemojo.com/.
A série contém 449 observações, de janeiro de 1982 à junho de 2019.
A partir disso, ajustamos um modelo ARIMA Sazonal e realizamos algumas comparações, com um modelo de suavização exponencial e com um modelo ajustado pela função do R auto.arima.
Série
scrapping_bilheteria <- function(i){
mojo <- paste0("https://www.boxofficemojo.com/monthly/?view=releasedate&chart=bymonth&month=",
as.character(i),
"&view=releasedate")
box_office <- read_html(mojo)
html_nodes(box_office, css = "center table") -> tabela
tabela_boxoffice <- html_table(tabela, header = T)
tabela_boxoffice <- tabela_boxoffice[[1]]
tabela_boxoffice %>%
rename(bilheteria = `Total Gross`,
Ano = Year,
filmes = Movies,
media = Avg.) %>%
mutate(Mes = month.abb[i],
bilheteria = parse_number(bilheteria),
media = parse_number(media)) %>%
select(Ano, Mes, bilheteria, filmes, media) -> tabela_boxoffice
tabela_boxoffice
}meses <- lapply(1:12, scrapping_bilheteria)
meses <- reduce(meses, rbind.data.frame)
meses %>%
mutate(Mes = ymd(paste0(Ano, "-", Mes, "-01"))) %>%
select(-Ano) %>%
arrange(Mes) -> mesesAnálise exploratória dos dados
Gráfico da série
library(dygraphs)
dygraph(series) %>% dyRangeSelector()Vemos que a série apresenta tendência, pode não ficar claro que a série apresenta sazonalidade no gráfico devido ao grande número de observações.
Vejamos os gráficos abaixo para averiguar a presença de sazonalidade.
Avaliação de sazonalidade
boxplot((meses$bilheteria) ~ month(meses$Mes))
Aparentemente as distribuições são diferentes para os meses. Como a série tem tendência, é válido averiguar esse gráfico, e outros parecidos feitos com relação a série diferenciada.
seried <- diff(series)
monthplot(seried)
boxplot(diff(meses$bilheteria) ~ month(meses$Mes)[-1])
Série diferenciada
plot(seried)
Observando o gráfico da série diferenciada vemos que esta aparenta apresentar heterocedasticidade. Daí, faz-se necessária alguma transformação para estabilizar a variância. O \(\lambda\) estimado por máxima verossimilhança para a transformação de box-cox foi 0.4, como o valor foi próximo de 0.5, usaremos a transformação raiz quadrada por conta da facilidade e intuitividade dessa transformação.
Transformação para estabilizar a variância
serie_sqrt <- sqrt(series)
plot(serie_sqrt)
serie_sq_dif <- diff(serie_sqrt)
plot(serie_sq_dif)
Agora a série aparenta ter variância constante. Verificamos abaixo utilizando o teste de levene.
Teste de levene
O teste de levene utilizando 10 grupos não rejeita a hipótese nula de que a variância da série não difere entre os grupos.
library(lawstat)
grupos <- rep(1:10, each = 45)[-c(449:450)]
levene.test(serie_sq_dif, group = grupos)##
## Modified robust Brown-Forsythe Levene-type test based on the
## absolute deviations from the median
##
## data: serie_sq_dif
## Test Statistic = 1.6194, p-value = 0.1071Função de autocorrelação
acf(serie_sq_dif, lag = 60)
Vemos a partir da função de autocorrelação que a parte sazonal é não estacionária, necessitando assim, da tomada de uma diferenciação sazonal.
FAC da série após a diferenciação sazonal
serie_dif_s <- diff(serie_sq_dif, lag = 12)
acf(serie_dif_s, lag = 60)
Função de autocorrelação parcial
pacf(serie_dif_s, lag = 60)
Obeservando a função de autocorrelação e de autocorrelação parcial da série, vemos que a autocorrelação é truncada em 1, tanto para a parte sazonal quanto para a não sazonal. E além disso, a função de autocorrelação parcial decresce exponencialmente para ambas.
Sugerindo, assim, que o modelo a ser utilizado é um \(SARIMA(0,1,1)\times(0,1,1)_{12}\)
Ajuste de modelo SARIMA
modelo <- arima(serie_sqrt, order = c(0,1,1),
seasonal = list(order = c(0,1,1), period = 12))
summary(modelo)##
## Call:
## arima(x = serie_sqrt, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1),
## period = 12))
##
## Coefficients:
## ma1 sma1
## -0.9709 -0.8158
## s.e. 0.0101 0.0396
##
## sigma^2 estimated as 8.063: log likelihood = -1082.54, aic = 2171.08
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.07905327 2.798144 2.141507 -1.462824 9.713009 0.3632378
## ACF1
## Training set -0.1874064Resíduos
residuos <- residuals(modelo)
plot(residuos)
Os resíduos não apresentam indícios de heterocedasticidade ou algum padrão.
Gráfico de autocorrelação dos resíduos
acf(residuos)
A função de autocorrelação dos resíduos apresenta um ponto consideravelmente alto em \(h = 1\), por isso vamos ajustar mais um modelo para tentar tratar desse problema.
Novo ajuste
A fim de tratar melhor a série transformada, ajustamos um \(SARIMA(1,1,1)\times(0,1,1)_{12}\).
modelo1 <- arima(serie_sqrt, order = c(1,1,1), seasonal = list(order = c(0,1,1), period = 12))
summary(modelo1)##
## Call:
## arima(x = serie_sqrt, order = c(1, 1, 1), seasonal = list(order = c(0, 1, 1),
## period = 12))
##
## Coefficients:
## ar1 ma1 sma1
## -0.1932 -0.9637 -0.8101
## s.e. 0.0477 0.0117 0.0392
##
## sigma^2 estimated as 7.779: log likelihood = -1074.52, aic = 2157.05
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.07317665 2.748451 2.105144 -1.433392 9.609341 0.3570701
## ACF1
## Training set -0.02072362Resíduos
residuos <- residuals(modelo1)
plot(residuos)
Os resíduos não apresentam indícios de heterocedasticidade.
acf(residuos)
A função de autocorrelação simples apresetna um valores pequenos e fora do intervalo em \(h = 2, 9\) e \(14\). Mas como esses valores são pouco significativos, optamos por ignorar e seguir com os modelos ajustados até agora.
Previsão
Separamos as observações de Janeiro de 2017 em diante para previsão, o que nos deixa com 29 valores para comparações posteriormente.
parcela_estimacao <- window(series, start = 1982, c(2016, 12))
modelo1 <- arima(sqrt(parcela_estimacao), order = c(1,1,1), seasonal = list(order = c(0,1,1), period = 12))
prev <- forecast(modelo1)
plot(prev)
Comparacão da previsão com os valores observados
valor_prev <- prev$mean
parcela_previsao <- window(series, start = c(2017,1), c(2018, 12))
comparacao <- cbind.data.frame(sqrt(parcela_previsao), valor_prev)
comparacao## sqrt(parcela_previsao) valor_prev
## 1 20.42058 21.14991
## 2 26.89981 26.74054
## 3 36.99054 30.51579
## 4 22.57654 26.23576
## 5 29.33258 35.02100
## 6 38.34840 35.71977
## 7 31.79780 35.62754
## 8 20.60097 28.77539
## 9 26.81977 24.96184
## 10 21.12345 26.37213
## 11 37.63642 35.82046
## 12 41.92613 38.41012
## 13 18.76166 21.72717
## 14 33.50970 26.99980
## 15 28.03034 30.82654
## 16 35.01714 26.53817
## 17 28.56046 35.32476
## 18 38.79046 36.02330
## 19 32.51923 35.93111
## 20 27.33496 29.07895
## 21 24.63940 25.26541
## 22 28.23296 26.67570
## 23 37.44329 36.12402
## 24 34.51087 38.71369Ajuste - Suavização exponencial de Holt-Winters
modelo2 <- HoltWinters(sqrt(parcela_estimacao))
plot(modelo2, main = "Ajuste do Modelo de Suavização exponencial de Holt-Winters")
Previsão - Holt-Winters
previsao <- forecast(modelo2)
plot(previsao, main = "Previsão do modelo de Suavização exponencial de Holt-Winters")
Comparação da previsão com os valores observados
val_prev <- previsao$mean
val_real <- window(series, start = c(2017,1), c(2018, 12))
comp <- cbind.data.frame(sqrt(val_real), val_prev)
comp## sqrt(val_real) val_prev
## 1 20.42058 21.11347
## 2 26.89981 26.39869
## 3 36.99054 30.47178
## 4 22.57654 25.95729
## 5 29.33258 34.48891
## 6 38.34840 35.41078
## 7 31.79780 35.55309
## 8 20.60097 28.31996
## 9 26.81977 24.67453
## 10 21.12345 25.63712
## 11 37.63642 35.66869
## 12 41.92613 38.35223
## 13 18.76166 21.32098
## 14 33.50970 26.60620
## 15 28.03034 30.67929
## 16 35.01714 26.16481
## 17 28.56046 34.69642
## 18 38.79046 35.61829
## 19 32.51923 35.76060
## 20 27.33496 28.52747
## 21 24.63940 24.88205
## 22 28.23296 25.84463
## 23 37.44329 35.87620
## 24 34.51087 38.55974Residuos
residuos <- residuals(modelo2)
plot(residuos)
acf(residuos)
Vemos que os resíduos aparentam estar distribuídos aleatóriamente em torno do zero e sua função de autocorrelação apresenta dois valores significativos.
Ajuste pelo auto.arima para série transformada
Para fins de comparação, ajustamos um modelo usando a função auto.arima para a série ajustada.
modelo3 <- auto.arima(sqrt(parcela_estimacao))
summary(modelo3)## Series: sqrt(parcela_estimacao)
## ARIMA(2,0,2)(1,1,1)[12] with drift
##
## Coefficients:
## ar1 ar2 ma1 ma2 sar1 sma1 drift
## 0.8329 -0.8175 -0.9450 0.8591 0.0071 -0.7417 0.0408
## s.e. 0.0898 0.0984 0.0911 0.0761 0.0728 0.0600 0.0030
##
## sigma^2 estimated as 7.388: log likelihood=-988.51
## AIC=1993.01 AICc=1993.38 BIC=2025.1
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.009013253 2.655894 2.020801 -1.027668 9.452849 0.7387238
## ACF1
## Training set -0.02819398Comparações
#modelo sarima inicial agora só com dados separados para estimação
modelo <- arima(sqrt(parcela_estimacao), order = c(0,1,1), seasonal = list(order = c(0,1,1), period = 12))Erros
erro <- forecast(modelo)$mean - sqrt(parcela_previsao)
erro1 <- forecast(modelo1)$mean - sqrt(parcela_previsao)
erro2 <- forecast(modelo2)$mean - sqrt(parcela_previsao)
erro3 <- forecast(modelo3)$mean - sqrt(parcela_previsao)
purrr::map(list(erro, erro1, erro2, erro3),
~c(mean(.x^2), mean(abs(.x)), mean(abs(.x/sqrt(parcela_previsao))))) -> erros
err <- reduce(erros, rbind.data.frame)
colnames(err) <- c('MSE','MAE','MAPE')
rownames(err) <- c('modelo','modelo 1', "modelo 2",'modelo 3')
round(err,3) %>% kable() %>% kable_styling(full_width = F)| MSE | MAE | MAPE | |
|---|---|---|---|
| modelo | 18.685 | 3.651 | 0.127 |
| modelo 1 | 18.520 | 3.622 | 0.126 |
| modelo 2 | 17.876 | 3.576 | 0.122 |
| modelo 3 | 19.932 | 3.696 | 0.132 |
Vemos que o modelo de suavização exponencial de Holt-Winters foi o que apresentou melhores resultados levando em conta todos dos erros acima avaliados. O modelo ajustado pela função do R auto.arima foi o que apresentou os piores resultados na tabela de previsão.
Previsão
plot(forecast(modelo)$mean^2,type="l",col=1,xlim=c(2017,2019),
ylim=c(0,2000), ylab="", xlab="")
par(new=T)
plot(forecast(modelo1)$mean^2,type="l",col=2,xlim=c(2017,2019),
ylim=c(0,2000), ylab="",xlab="")
par(new=T)
plot(forecast(modelo2)$mean^2,type="l",col=3,xlim=c(2017,2019),
ylim=c(0,2000), ylab="", xlab="")
par(new=T)
plot(forecast(modelo3)$mean^2,type="l",col=4,xlim=c(2017,2019),
ylim=c(0,2000), ylab="Previsões", xlab="Tempo")
Podemos ver no gráfico acima que as previsões para os anos de 2017 e 2018 são bem próximas nos 4 modelos, sendo que o modelo ajustado pela função auto.arima (azul) é o que mais se destaca superiormente.
Conclusão
Conseguimos modelar a série utilizando alguns modelos arima sazonais e um modelo de suavização exponencial, e comparar a performance dos modelos na parcela da série reservada para previsão.
Dentre todos os modelos ajustados, o modelo de suavização exponencial foi o que apresentou menores erros e entre os modelos SARIMA, o modelo \(SARIMA(1,1,1)\times(0,1,1)_{12}\) foi o que apresentou melhores resultados.