The objective of this project is to make predictions of future changes in the retail sales volume index in the Brazilian state of Minas Gerais. To do this, we will use the index provided by IBGE as a time series and apply different statistical models to obtain forecasts. In particular, we will use the ARIMA (Autoregressive Integrated Moving Average) model with conditional volatility adjustment (GARCH) and the ETS model.

The ARIMA model is a commonly used approach for time series forecasting. It takes into account past patterns and trends to make future predictions. Additionally, we will apply the conditional volatility adjustment GARCH, which considers the volatility of the time series when making forecasts.

In addition to the ARIMA model with GARCH adjustment, we will also use the ETS (Error, Trend, Seasonality) model. This model takes into account the error, trend, and seasonality present in the time series to generate forecasts.

After estimating both models, we will compare and analyze the forecasts generated by each of them. Based on this analysis, we can evaluate the effectiveness of each model in predicting future changes in the retail sales volume index in Minas Gerais.

The time series used in this work has a fixed base index, referring to the monthly average of sales volumes in the year 2014. This means that the observed values are compared to this average, allowing us to assess variations over time.

In the end, we will have models that could assist companies, financial institutions, and governments in anticipating trends and preparing adequately for changes in retail sales volume.

The database/time series used:
varejoMG (Please right-click and select “Open in a new tab/window.”)

Installation and Loading of All Packages

pacotes <- c("readr","readxl","plotly","trend","tidyverse","gridExtra","forecast","TTR",
             "smooth", "tsibble", "fable","tsibbledata", "fpp3","lubridate",
             "urca", "dygraphs", "quantmod","BETS","tseries","FinTS","feasts",
             "gridExtra", "scales", "caret","xtable", "tsutils","GetBCBData", 
             "quantmod","dgof","seasonal","devtools","transformr","gganimate","rugarch")

if(sum(as.numeric(!pacotes %in% installed.packages())) != 0){
  instalador <- pacotes[!pacotes %in% installed.packages()]
  for(i in 1:length(instalador)) {
    install.packages(instalador, dependencies = T)
    break()}
  sapply(pacotes, require, character = T)
} else {
  sapply(pacotes, require, character = T)
}

Saving and analyzing the time series in R

We will use the gbcbd_get_series command to retrieve the time series we will be using directly from the Brazilian Central Bank’s website.

varejoMG=gbcbd_get_series(1472,first.date='2000-01-01')
## 
## Fetching id = 1472 [1472] from BCB-SGS with cache 
##   Found 276 observations
TSvarejoMG=ts(varejoMG[2], start = c(2000,1), end = c(2022,12), frequency = 12)
dygraph(TSvarejoMG)

_https://www3.bcb.gov.br/sgspub/localizarseries/localizarSeries.do?method=prepararTelaLocalizarSeries_

Decomposing and visualizing the components of the time series

decomposicao <- decompose(TSvarejoMG)
plot(decomposicao)

We can clearly see that there is a trend and seasonality in the series.

Window splitting

Let’s split our series into two parts. One for training the algorithm, and the other for evaluating its efficiency.

TSvarejotreinoMG=window(TSvarejoMG, start=c(2000,1), end=c(2020,12))
TSvarejotesteMG=window(TSvarejoMG,start=c(2021,1), end=c(2022,12))
length(TSvarejotesteMG)
## [1] 24

Plotting the two series together for verification.

autoplot(TSvarejoMG) +
  autolayer(TSvarejotreinoMG, series="Treino") +
  autolayer(TSvarejotesteMG, series="Teste") +
  scale_color_viridis_d() +
  theme_bw()

Analysis of the series

ggtsdisplay(TSvarejotreinoMG)

As we observe the presence of both trend and seasonality, we are likely dealing with a non-stationary series. Let’s confirm this.

Stationarity Test

ur.df(TSvarejotreinoMG)
## 
## ############################################################### 
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test # 
## ############################################################### 
## 
## The value of the test statistic is: 0.1342

The series is not stationary, as 0.1342 is greater than the critical values (for significance levels of 1%, 5%, and 10%). Therefore, it needs to be differenced or transformed before applying an ARIMA model to it. The differencing process is a commonly used technique to transform non-stationary time series into stationary time series. The basic idea behind differencing is to calculate the difference between consecutive observations in the series.

Let’s see how many differentiations will be required using the ndiffs command:

ndiffs(TSvarejotreinoMG)
## [1] 1

1 differentiation is needed to make the series stationary.

ARIMA Modeling

Now that we’ve conducted a brief analysis of the series, we can proceed with modeling. Generally, ARIMA models are suitable for series with trends and seasonality. They can capture these patterns through autoregressive (AR), moving average (MA), and differencing (I) components.

In addition to ARIMA, ETS models can also handle this type of series. They are based on exponential smoothing techniques and are useful when series patterns are additive (A) or multiplicative (M). ETS models can include additional components to capture seasonality, such as seasonal exponential smoothing (S).

Therefore, since we are dealing with a series with such characteristics, we will use R to estimate an ARIMA model for our series. Then, we will estimate an ETS model and compare the predictive capabilities of the two.

Estimating an ARIMA model

arimavarejoMG=auto.arima(TSvarejotreinoMG, trace=T)
summary(arimavarejoMG)
## Series: TSvarejotreinoMG 
## ARIMA(4,1,1)(0,1,2)[12] 
## 
## Coefficients:
##          ar1     ar2     ar3      ar4      ma1     sma1     sma2
##       0.4114  0.0958  0.2099  -0.1980  -0.8623  -0.2939  -0.1795
## s.e.  0.0933  0.0780  0.0741   0.0691   0.0759   0.0791   0.0742
## 
## sigma^2 = 7.45:  log likelihood = -577.15
## AIC=1170.31   AICc=1170.93   BIC=1198.12
## 
## Training set error measures:
##                      ME     RMSE      MAE        MPE     MAPE     MASE
## Training set 0.07765287 2.618903 1.827275 0.09886015 2.263547 0.464168
##                      ACF1
## Training set -0.001970372

We can see that R has correctly estimated a seasonal ARIMA (sARIMA) model with a seasonal period of 12. Now, the question remains whether the estimated parameters, which include 4 autoregressive terms (AR), 1 non-seasonal differencing (I), 1 non-seasonal moving average term (MA), as well as 1 seasonal differencing and 2 seasonal moving average terms, will be able to make accurate forecasts. Notice that R performed the differencing we discussed earlier.

Validation and Diagnosis

Ljung-Box Test

We will use the Ljung-Box test through the checkresiduals function to test for autocorrelation in the residuals. It is desirable for the residuals of predictive time series models not to exhibit autocorrelation. When residuals show autocorrelation, it means that subsequent observations are related to each other, and the temporal structure of the data has not been adequately captured by the model. This can lead to inaccurate and inefficient forecasts, as well as incorrect interpretations of results.

checkresiduals(arimavarejoMG)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(4,1,1)(0,1,2)[12]
## Q* = 17.121, df = 17, p-value = 0.4462
## 
## Model df: 7.   Total lags used: 24

Basically, the closer the p-values are to 1 in the Ljung-Box test, the lower the chance of autocorrelation among the residuals of such a model. In this regard, we can say that our model appears to be reasonable. However, this does not necessarily indicate that the model is bad. It’s important to remember that assessing the quality of a time series is not solely based on the absence of correlation in the residuals. Other criteria, such as the context of the time series, data stationarity, forecast accuracy, and the normality of residuals, should also be considered. Let’s perform more tests to further evaluate the model.

Kolmogorov-Smirnov

We will use the Kolmogorov-Smirnov test through the ks.test function to check whether the residuals of our models follow a normal distribution. The Kolmogorov-Smirnov test is a goodness-of-fit test used to compare a theoretical distribution with a sample of data, and in our case, it will be used to assess the normality of the residuals. We expect the residuals to follow a normal distribution because a model with residuals adhering to normality tends to indicate that the model is capable of more accurate forecasts.

ks.test(arimavarejoMG$residuals, "pnorm", mean(arimavarejoMG$residuals),
        sd(arimavarejoMG$residuals))
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  arimavarejoMG$residuals
## D = 0.090908, p-value = 0.03105
## alternative hypothesis: two-sided

We can see that at a significance level of 0.5, the residuals of our model do not appear to follow a normal (Gaussian) distribution.

Lagrange Multiplier

Finally, let’s test for the presence of ARCH effects. ARCH effects are a specific form of heteroscedasticity in which the variance of errors is modeled as an autoregressive function of the errors themselves. Such effects can lead to parameter estimation inefficiency, violation of the homoscedasticity assumption, difficulty in interpreting results, and potential residual autocorrelation.

ArchTest(arimavarejoMG$residuals)
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  arimavarejoMG$residuals
## Chi-squared = 50.372, df = 12, p-value = 1.202e-06

Based on these results, we can interpret that there is significant evidence to reject the null hypothesis, suggesting that we are dealing with a model with ARCH effects in the residuals. Despite the negative indicators, let’s test the predictive capability of the model using our test series “TSvarejotesteMG.”

Testing the ARIMA model

prevvarejo=forecast::forecast(arimavarejoMG, h=24)
autoplot(prevvarejo) +
  theme_bw()

ggplotly(
  autoplot(TSvarejotreinoMG)+
    autolayer(TSvarejotesteMG,serie="Valores Reais")+
    autolayer(prevvarejo$mean, serie="Forecast")+
    scale_colour_viridis_d()+
    scale_y_continuous(labels=scales::comma)+
    theme_bw()
)
forecast::accuracy(prevvarejo, TSvarejotesteMG)
##                       ME     RMSE      MAE         MPE     MAPE     MASE
## Training set  0.07765287 2.618903 1.827275  0.09886015 2.263547 0.464168
## Test set     -0.12693049 5.112050 4.439883 -0.03313566 4.060955 1.127828
##                      ACF1 Theil's U
## Training set -0.001970372        NA
## Test set      0.617369420 0.6267072

Interestingly, we obtained good forecasts with a MAPE of 4.06. However, since the ArchTest indicated the possible presence of conditional heteroscedasticity in the residuals, let’s try implementing a GARCH model on the residuals of our model to obtain their conditional volatility. Then, we will use this volatility to adjust the forecast values.

Adjusting the forecasts for conditional volatility

Initially, we need to store our residuals in an object:

residuos_arimavarejoMG <- residuals(arimavarejoMG)

Let’s now create a time series using the residuals:

TSresiduos_arimavarejoMG <- ts(residuos_arimavarejoMG, frequency = 12, start=c(2000,1), end=c(2020,12))

In some cases, we would need to difference the time series of residuals because GARCH models typically work with return series, not absolute value series. However, since our original time series has already been differenced, we can proceed with modeling the residuals. Let’s specify a basic GARCH(2,2) model.

GARCH <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(2, 2)),
                    mean.model = list(armaOrder = c(0, 0)))

several values were tested for (p, q), and (2,2) were the ones that yielded the best results.

Now, let’s create a model with the chosen specifications using our time series of residuals from the ARIMA model:

GARCHresiduos_arimavarejoMG <- ugarchfit(GARCH, data = TSresiduos_arimavarejoMG)

Now we can save the conditional volatility forecasts in another object:

volatilidade <- sigma(GARCHresiduos_arimavarejoMG)

Finally, we can apply the volatility to the forecasts. Let’s create a time series with the volatility for this purpose:

TSvolatilidade <- ts(volatilidade, start = start(prevvarejo$mean), frequency = frequency(prevvarejo$mean))

prevvarejo$mean_garch <- prevvarejo$mean * sqrt(TSvolatilidade)

media1 <- mean(prevvarejo$mean)
media2 <- mean(prevvarejo$mean_garch)
media3 <- media2/media1

prevvarejo$mean_garch <- (prevvarejo$mean * sqrt(TSvolatilidade)) / media3

Above, we have a commonly used formula to adjust forecasts to conditional volatilities. The subsequent process using the means is necessary to adjust and normalize the impact of volatilities.

Testing the forecasts using GARCH volatility as an adjustment.

autoplot(prevvarejo) +
  autolayer(prevvarejo$mean_garch,serie="ajuste GARCH")

ggplotly(
  autoplot(TSvarejotreinoMG)+
    autolayer(TSvarejotesteMG,serie="Valores Reais")+
    autolayer(prevvarejo$mean_garch, serie="Forecast")+
    scale_colour_viridis_d()+
    scale_y_continuous(labels=scales::comma)+
    theme_bw()
)
forecast::accuracy(prevvarejo$mean_garch, TSvarejotesteMG)
##                  ME    RMSE      MAE        MPE     MAPE      ACF1 Theil's U
## Test set -0.1269305 5.49163 4.295751 -0.2318458 3.980282 0.6783739 0.6836237

With a MAPE of 3.98, the adjusted model appears to be slightly better than the original ARIMA. Now, let’s try an ETS model and compare all the models in the end.

ETS Model

ETSvarejoMG = ets(TSvarejotreinoMG)

summary(ETSvarejoMG)
## ETS(M,Ad,M) 
## 
## Call:
##  ets(y = TSvarejotreinoMG) 
## 
##   Smoothing parameters:
##     alpha = 0.4339 
##     beta  = 0.0163 
##     gamma = 0.2898 
##     phi   = 0.9799 
## 
##   Initial states:
##     l = 49.0505 
##     b = -0.0063 
##     s = 1.3137 0.9942 1.0031 0.9641 0.9973 1.005
##            0.9605 0.9986 0.9598 0.9763 0.8956 0.9318
## 
##   sigma:  0.0332
## 
##      AIC     AICc      BIC 
## 1869.763 1872.699 1933.293 
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.1935402 2.783216 1.966726 0.2560567 2.438439 0.4995916
##                    ACF1
## Training set 0.08862719
ETSprev=forecast::forecast(ETSvarejoMG, h=24)

The “M” indicates the presence of multiplicative seasonality, while the “Ad” indicates the presence of additive trend. Therefore, R has estimated a model that is suitable for time series with a trend added to a multiplicative seasonal component.

autoplot(ETSprev) +
   theme_bw()

ggplotly(
  autoplot(TSvarejotreinoMG)+
    autolayer(TSvarejotesteMG,serie="Valores Reais")+
    autolayer(ETSprev$mean, serie="Previstos por ETS")+
    scale_colour_viridis_d()+
    scale_y_continuous(labels=scales::comma)+
    theme_bw()
)

Checking the predictive accuracy of the ETS model

forecast::accuracy(ETSprev$mean, TSvarejotesteMG)
##                 ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
## Test set -3.440922 4.565962 3.883694 -3.152128 3.558263 0.2637249  0.561279

We can see that the ETS model has had the best accuracy so far…

Ljung-Box

checkresiduals(ETSvarejoMG)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(M,Ad,M)
## Q* = 38.856, df = 7, p-value = 2.082e-06
## 
## Model df: 17.   Total lags used: 24

Despite the high accuracy, and unlike the ARIMA model, the Ljung-Box test strongly indicates the presence of autocorrelation in the residuals. This is not necessarily a bad thing, but it suggests that the model may not be capturing the entire dependency structure in the time series.

Kolmogorov-Smirnov

ks.test(ETSvarejoMG$residuals, "pnorm", mean(ETSvarejoMG$residuals),
        sd(ETSvarejoMG$residuals))
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  ETSvarejoMG$residuals
## D = 0.041501, p-value = 0.7782
## alternative hypothesis: two-sided

Unlike the ARIMA model, for a significance level of 0.05, the residuals of the ETS model follow a normal distribution.

Lagrange Multiplier

ArchTest(ETSvarejoMG$residuals)
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  ETSvarejoMG$residuals
## Chi-squared = 56.902, df = 12, p-value = 8.232e-08

In the presence of indications of ARCH effects in the residuals of an ETS model, it is possible to explore the possibility of adjusting the model to conditional volatilities, similar to what was done in the case of the ARIMA model. In this context, a similar procedure was conducted where an sGARCH(0,1) model was fitted to the conditional volatilities. However, the results showed that these adjustments did not improve the forecasts of the model. This finding suggests that even with the presence of ARCH effects, our ETS model performs better without adjusting for conditional volatilities.

Taking a look at the accuracy of the “adjusted” ETS model

forecast::accuracy(ETSprev$mean_garch, TSvarejotesteMG)
##                 ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
## Test set -3.440922 4.585305 3.921992 -3.143813 3.584261 0.2582928 0.5616741

Values of (p, q) different from (0, 1) resulted in even worse outcomes.

As we can see, the MAPE is slightly higher than the initial ETS model, and it indicates that, unlike our ARIMA model, adjusting forecasts for conditional volatilities in this case did not improve the model. For this reason, let’s consider the initial ETS model as our final ETS model..

Comparisons

Overall, we obtained three models with good predictive capability, but the ETS model had the lowest MAPE, indicating the best predictive ability. It’s worth noting that despite ARIMA showing a better fit to the data (with lower AIC, AICc, and BIC), the ETS model provided better forecasts.

Accuracies

ARIMA model

forecast::accuracy(prevvarejo, TSvarejotesteMG)
##                       ME     RMSE      MAE         MPE     MAPE     MASE
## Training set  0.07765287 2.618903 1.827275  0.09886015 2.263547 0.464168
## Test set     -0.12693049 5.112050 4.439883 -0.03313566 4.060955 1.127828
##                      ACF1 Theil's U
## Training set -0.001970372        NA
## Test set      0.617369420 0.6267072

ARIMA model adjusted with conditional volatilities

forecast::accuracy(prevvarejo$mean_garch, TSvarejotesteMG)
##                  ME    RMSE      MAE        MPE     MAPE      ACF1 Theil's U
## Test set -0.1269305 5.49163 4.295751 -0.2318458 3.980282 0.6783739 0.6836237

ETS model

forecast::accuracy(ETSprev, TSvarejotesteMG)
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set  0.1935402 2.783216 1.966726  0.2560567 2.438439 0.4995916
## Test set     -3.4409221 4.565962 3.883694 -3.1521276 3.558263 0.9865437
##                    ACF1 Theil's U
## Training set 0.08862719        NA
## Test set     0.26372494  0.561279

Conclusion

Despite ARIMA being a suitable choice for time series with trends and seasonality, it doesn’t guarantee that it will always produce the best forecasts compared to other models like ETS. In our case, the ETS model seems to have had better forecasts due to its ability to capture the trend and seasonality present in the time series. The ARIMA model, even with conditional volatility adjustment, may not have been able to capture these patterns as efficiently.

Therefore, while ARIMA is a common choice for time series with trends and seasonality, it’s important to consider other approaches like ETS to determine which model offers the best predictive performance for your specific time series.

To have confidence in forecasts, it’s always a good practice to test different models and compare them to obtain the best possible model for a given series. One justification for the model with higher AIC, AICc, and BIC values having a lower MAPE is that the time series in question may be relatively simple, with patterns that are easily captured by a less complex model.