Registered S3 method overwritten by 'quantmod':
method from
as.zoo.data.frame zoo
library(forecast)library(prophet)
Carregando pacotes exigidos: Rcpp
Carregando pacotes exigidos: rlang
library(xts)
Carregando pacotes exigidos: zoo
Attaching package: 'zoo'
The following object is masked from 'package:tsibble':
index
The following objects are masked from 'package:base':
as.Date, as.Date.numeric
######################### Warning from 'xts' package ##########################
# #
# The dplyr lag() function breaks how base R's lag() function is supposed to #
# work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or #
# source() into this session won't work correctly. #
# #
# Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
# conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
# dplyr from breaking base R's lag() function. #
# #
# Code in packages is not affected. It's protected by R's namespace mechanism #
# Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
# #
###############################################################################
Attaching package: 'xts'
The following objects are masked from 'package:dplyr':
first, last
Para a resolução dos exercícios, será utilizada a base de dados aus_retail, que consiste em dados sobre as vendas no varejo australiano, uma série temporal mensal de diferentes categorias de varejo.
Abaixo esta listada todas as categorias presentes na base de dados:
unique(aus_retail$Industry)
[1] "Cafes, restaurants and catering services"
[2] "Cafes, restaurants and takeaway food services"
[3] "Clothing retailing"
[4] "Clothing, footwear and personal accessory retailing"
[5] "Department stores"
[6] "Electrical and electronic goods retailing"
[7] "Food retailing"
[8] "Footwear and other personal accessory retailing"
[9] "Furniture, floor coverings, houseware and textile goods retailing"
[10] "Hardware, building and garden supplies retailing"
[11] "Household goods retailing"
[12] "Liquor retailing"
[13] "Newspaper and book retailing"
[14] "Other recreational goods retailing"
[15] "Other retailing"
[16] "Other retailing n.e.c."
[17] "Other specialised food retailing"
[18] "Pharmaceutical, cosmetic and toiletry goods retailing"
[19] "Supermarket and grocery stores"
[20] "Takeaway food services"
Para realizar o trabalho iremos utilizar a categoria de “Department stores”:
Estes dados referem-se a lojas de departamentos, que são grandes estabelecimentos de varejo que oferecem uma ampla variedade de mercadorias em diferentes seções ou departamentos dentro da mesma loja.
Pre-Precessamento
Da forma que o dataframe está organizado, existem diferentes estados:
unique(department_stores_data$State)
[1] "Australian Capital Territory" "New South Wales"
[3] "Queensland" "South Australia"
[5] "Victoria" "Western Australia"
Para a resolução dos exercícios da lista, serão utilizados os dados da Austrália inteira, sendo necessária a seguinte transformação:
Com esta transformação, os dados são agrupados pela coluna “Month”, somando os turnovers de todos os estados conforme essa coluna. O resultado final é a receita total de lojas de departamentos para toda a Austrália.
Resolução:
Questão 1
Plote os gráficos da série e das ACF/PACF e discorra sobre o aspecto da série. Faça o teste de raíz unitária e explique as conclusões. Existe tendência estocástica?
O primeiro gráfico análisado será o Time Plot:
department_stores_group %>%autoplot(Turnover_total) +labs(y="Turnover", title="Australian Department Stores Turnover")
O gráfico plotado acima demonstra um aparente padrão sazonal na série, com um pico de receita que se repete ao longo dos anos. Além disso, parece haver uma tendência de aumento de receita ao longo da série.
Outro gráfico analisado será o Seasonal Plot:
department_stores_group %>%gg_season(Turnover_total, period ="year")
Observando o gráfico plotado, é possível identificar um padrão sazonal de aumento de vendas no mês de dezembro, representado pelo pico no gráfico. Além disso, após esse pico de vendas, há uma queda aparente em fevereiro, que se mostra mais acentuada nos anos mais recentes.
Para uma melhor compreensão sobre a possível tendência de aumento de receita ao longo dos anos, podemos observar o gráfico de subplots abaixo:
Como podemos observar pelo gráfico discriminado por meses acima, a sazonalidade de aumento de vendas em dezembro está bem representada. Além disso, nota-se que, de maneira geral, houve um aparente aumento acentuado de receita em todos os meses entre os anos de 1990 e o início dos anos 2000, estabilizando-se nos anos seguintes.
Para reforçar essas premissas, podemos analisar as funções de autocorrelação e autocorrelação parcial.
# Gráfico de Autocorrelaçãodepartment_stores_group %>%ACF(Turnover_total, lag_max =100) %>%autoplot()
Como esperado, a função de autocorrelação revela fortes correlações entre os lags múltiplos de 12 (12, 24, 36…) indicando a sazonalidade associada aos meses de dezembro. Além disso, é possível observar uma aparente tendência positiva, evidenciada pela significativa autocorrelação nos lags menores, que diminui gradualmente à medida que os lags aumentam.
# Gráfico de Autocorrelação Parcialdepartment_stores_group %>%PACF(Turnover_total, lag_max =100) %>%autoplot()
Analisando o gráfico gerado pela PACF, percebemos um pico positivo seguido por um pico negativo, ambos muito acima das bandas do intervalo de confiança, nos lags 12 e 13. Isso indica uma forte sazonalidade neste período do ano, conforme observado nos gráficos anteriores.
Ao que tudo indica, esta série temporal possui raiz unitária e não é estacionária. Para obtermos mais uma medida que confirme esta hipótese, podemos realizar o teste de raiz unitária de Dickey-Fuller abaixo:
# Realizar o teste ADFadf.test(department_stores_group$Turnover_total, alternative =c("stationary", "explosive"), k =trunc((length(department_stores_group$Turnover_total)-1)^(1/3)))
Warning in adf.test(department_stores_group$Turnover_total, alternative =
c("stationary", : p-value smaller than printed p-value
Augmented Dickey-Fuller Test
data: department_stores_group$Turnover_total
Dickey-Fuller = -8.9811, Lag order = 7, p-value = 0.01
alternative hypothesis: stationary
O resultado do teste realizado, com p-valor menor do que o nível de 0,05, indica que não existem evidências significativas para rejeitar a hipótese nula de que a série não é estacionária. Este resultado contraria os resultados obtidos anteriormente pelos gráficos da série e pela função ACF. Um dos motivos para isso pode ser o fato de a série possuir uma tendência fraca que não foi capturada pelo teste ADF. Nesse caso, podemos dizer que a série possui uma fraca estacionariedade e não tem uma tendência estocástica evidente. Contudo, esses resultados podem estar sendo afetados pela forte sazonalidade que a série apresenta.
Questão 2
Use um método de decomposição estudado e explique os componentes, tendo em vista o aspecto da série. Existe sazonalidade? Qual periodicidade?
Aplicando método de decomposição dos componentes:
department_stores_group %>%model(classical_decomposition(Turnover_total, type ="additive") ) |>components() |>autoplot() +labs(title ="Decomposição Série Departaments Stores")
Warning: Removed 6 rows containing missing values or values outside the scale range
(`geom_line()`).
Como podemos observar no gráfico de decomposição, existe um padrão sazonal de picos de venda no final de cada ano da série, fenômeno que já podia ser observado nos gráficos gerados anteriormente.
Ao contrário do resultado do teste ADF realizado na questão anterior, a decomposição da série temporal sugere a presença de uma tendência positiva de aumento nas vendas ao longo dos anos, o que também pode ser observado nos gráficos ACF e Seasonal plots realizados anteriormente.
Assim, podemos concluir que existe uma sazonalidade de 12 meses ao final de cada ano.
Questão 3
Identifique um modelo SARIMA apropriado. Justifique o modelo escolhido através de critérios de informação e/ou usando a etapa de identificação da abordagem Box-Jenkins. Nesta etapa é possível utilizar um método automático de seleção das ordens do modelo.
Como foi observado anteriormente, existe uma sazonalidade anual na série. Por conta disso, iremos diferenciar a série, componente d, com base nesta sazonalidade anual para tornar a série ainda mais estacionária.
# Aplicando diferenciação de 12 mesesdepartment_stores_diff <-diff(department_stores_group$Turnover_total, lag=12)# Plotar ACFacf(department_stores_diff, main="ACF da Série Diferenciada Anualmente")
# Plotando PACfpacf(department_stores_diff, main="PACF da Série Diferenciada Anualmente")
Ao aplicar a diferenciação, percebe-se que os componentes de tendência e sazonalidade observados anteriormente foram reduzidos e que a série apresentou um resultado de ACF muito mais próximo de um padrão estacionário.
Com isso, podemos rodas os modelos abaixo:
# Definindo os dadosy <-ts(department_stores_group$Turnover_total, frequency =12) # Definindo os parâmetros para testarorders <-list(c(1,1,0), c(1,1,1), c(0,1,0), c(0,1,1))seasonal_orders <-list(c(1,1,0), c(1,1,1), c(0,1,0))# Inicializando listas para armazenar os resultadosmodels <-list()aic <-numeric()aicc <-numeric()bic <-numeric()# Loop para ajustar os modelosfor (i inseq_along(orders)) {for (j inseq_along(seasonal_orders)) { model <-Arima(y, order = orders[[i]],seasonal =list(order = seasonal_orders[[j]], period =12),method ="ML") model_label <-paste("fit", i, j, sep="_") models[[model_label]] <- model aic <-c(aic, model$aic) aicc <-c(aicc, model$aicc) bic <-c(bic, model$bic) }}# Criando uma tabela com os resultadosresults <-cbind(Model =names(models), AIC = aic, AICC = aicc, BIC = bic)results <-as.data.frame(results)# Mostrando a tabela de resultadosprint(results[order(results$AICC),] )
De todos os modelos testados o que obteve melhor desempenho de acordo com a métrica AICC foi o modelo fit_2_2:
best_model <- models[["fit_2_2"]]# Mostrando um resumo do modelosummary(best_model)
Series: y
ARIMA(1,1,1)(1,1,1)[12]
Coefficients:
ar1 ma1 sar1 sma1
-0.2217 -0.8239 0.1506 -0.4830
s.e. 0.0548 0.0320 0.1098 0.0958
sigma^2 = 2514: log likelihood = -2282.44
AIC=4574.89 AICc=4575.03 BIC=4595.18
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set -0.6230457 49.16252 36.80342 -0.4086745 3.528757 0.7334425
ACF1
Training set -0.024329
Ajustando o modelo otimizado:
fit =Arima(y, order=c(1,1,1), seasonal=list(order=c(1,1,1), period=12), method="ML")
Questão 4
Faça a verificação de diagnóstico residual de seu modelo ARIMA. Os resíduos são ruído branco?
Com o modelo otimizado ajustado, podemos realizar o diagnostico de desempenho.
Faça um backtest e verifique qual modelo tem melhor poder preditivo (calcule medidas de performance que sejam comparáveis entre os dois modelos usando diferentes janelas de tempo). Nesta etapa você pode usar apenas um split treino/teste temporal ou vários (cross-validation). Justifique qual dos modelos você usaria.
Warning: There was 1 warning in `dplyr::summarise()`.
ℹ In argument: `dplyr::across(.fns = summary_fns, ...)`.
Caused by warning:
! Using `across()` without supplying `.cols` was deprecated in dplyr 1.1.0.
ℹ Please supply `.cols` instead.
Accuracy Table
.model_id
.model_desc
.type
n
mae
mape
mase
smape
rmse
rsq
1
ARIMA(1,1,1)(1,1,1)[12]
Resamples
10
61.39
4.20
0.22
4.14
73.52
0.98
2
PROPHET
Resamples
10
71.32
4.95
0.26
4.86
84.53
0.98
Apesar do cross-validation ter piorado o resultado do SARIMA, ele continua melhor que o resultado do Prophet.
