A Comparative Analysis of U.S. Retail Sales Forecasting Models: Implications for Investment Decisions

Md Mahmudul Hasan

Assignment:

Context:
You are a quantitative analyst for a hedge fund which uses predictive models of the economy to try to forecast movements of key economic indicators, and creates trading strategies which attempt to “front run” changes in these indicators. You have been asked by your Director to compare several models of U.S. Retail Sales, evaluate their forecasting ability, and submit a report to be discussed at your next investment committee meeting.

Assignment:
- Within the above context, write a report suitable for sharing with your team on the modelling and predictability of Retail Sales.
- Retail Sales: Non-seasonally adjusted Advance Retail Sales: Retail Trade from FRED from January 1990 to July 2022.

Data

Fetching Data from FRED

# Advance Retail Sales: Retail Trade (RSXFSN)
# Data is not available from 1990, its available from 1992!
# https://fred.stlouisfed.org/series/RSXFSN

retail_sales<-tidyquant::tq_get(c("RSXFSN"), 
                           get="economic.data",
                           from="1992-01-01", 
                           to="2022-07-31") %>% 
  tidyr::fill(price,.direction="down") %>% stats::na.omit()

kable(head(retail_sales)) %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

symbol	date	price
RSXFSN	1992-01-01	130683
RSXFSN	1992-02-01	131244
RSXFSN	1992-03-01	142488
RSXFSN	1992-04-01	147175
RSXFSN	1992-05-01	152420
RSXFSN	1992-06-01	151849

Descriptive Table

A full descriptive table of your time series variable. This includes # of observations, # of missing values, mean, median, standard deviation, skewness, and kurtosis. Discuss any notable statistical features of your variable.

retail_sales %>% 
  dplyr::group_by(symbol) %>% 
  dplyr::summarise(min=min(price), 
                   max=max(price),
                   mean=mean(price),
                   median=median(price), 
                   sd=sd(price), 
                   skew=skewness(price), 
                   kurtosis = kurtosis(price) - 3) %>%
  kable(caption="Advance Retail Sales: Retail Trade (RSXFSN)")  %>% 
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive")) %>% 
  row_spec(1, bold = TRUE, color = "white", background = "#1f77b4")

Advance Retail Sales: Retail Trade (RSXFSN)

symbol	min	max	mean	median	sd	skew	kurtosis
RSXFSN	130683	625145	316065.8	311292	106797.7	0.4703032	-0.2044121

Notable Features

Volatility:

The high standard deviation (106797.7) relative to the mean (316065.8) suggests significant volatility, which could have implications for forecasting and risk assessment.

Potential for Positive Outliers:

The positive skewness (0.4703032) suggests that there could be periods where sales exceed the average significantly, which is important for modeling scenarios that may capture such outliers.

Symmetry in Central Tendency:

The closeness of the mean(316065.8) and median (311292) indicates that while there may be occasional spikes, the general behavior of the series is not heavily affected by extreme values.

Forecasting Implications:

Given the high variability and potential for positive outliers, any forecasting model developed using this data may need to account for periods of rapid growth or spikes in sales.

Time Series & ACF/PACF

Time Series Plot and ACF/PACF Generate a time series plot and an ACF and PACF chart.

- Discuss any interesting features in the time series chart and the ACF and PACF related to potential non-stationarity, AR or MA lag orders, and potential seasonality
- If the data appears non-stationary, take a first difference and re-evaluate the ACF and PACF as above

Time Series Plot and ACF/PACF

# Converting to time series
retail_ts <- ts(retail_sales$price, start = c(1992, 1), frequency = 12)

# Plotting 
p <- autoplot(retail_ts) + ggtitle("U.S. Retail Sales (1992-2022)")
interactive_plot <- ggplotly(p)
interactive_plot

# ACF and PACF plots
acf_plot <- ggAcf(retail_ts) + ggtitle("ACF of Retail Sales") 
interactive_acf <- ggplotly(acf_plot)
interactive_acf

pacf_plot <- ggPacf(retail_ts) + ggtitle("PACF of Retail Sales")
interactive_pacf <- ggplotly(pacf_plot)
interactive_pacf

Time Series Plot and Autocorrelation Analysis Time Series Plot:
Notable Observations:
- Retail sales exhibit a strong upward trend from 1992 to 2022, suggesting potential non-stationarity.
- There are signs of seasonality, particularly in annual cycles, which aligns with economic patterns of consumer behavior.

ACF and PACF Charts (Before Differencing):

plot_grid(acf_plot, pacf_plot, align = "hv", rel_widths = c(1, 1), rel_heights = c(1, 1))

The ACF shows high autocorrelations at low lags that gradually decline, which is a typical sign of non-stationarity. The PACF has significant autocorrelations at the first few lags, suggesting potential AR stationarity.

Differencing to Address Non-Stationarity
ACF and PACF Charts (After First Differencing):

# First difference to remove trend
diff_retail_ts <- diff(retail_ts)
# Plotting differenced time series
autoplot(diff_retail_ts, facets=FALSE, ts.colour="#1f77b4") + 
  ggtitle("Differenced U.S. Retail Sales")

# ACF and PACF for differenced series
diff_retail_ts_acf <- ggAcf(diff_retail_ts) +
  ggtitle("ACF of Differenced Retail Sales")
interactive_diff_acf <- ggplotly(diff_retail_ts_acf)
interactive_diff_acf

diff_retail_ts_pacf <- ggPacf(diff_retail_ts) + ggtitle("PACF of Differenced Retail Sales")
interactive_diff_pacf <- ggplotly(diff_retail_ts_pacf)
interactive_diff_pacf

Given the strong trend and the behavior in the ACF and PACF charts, the data likely exhibits non-stationarity, potentially requiring first-differencing to stabilize the mean.
After taking first differences, the ACF shows significant autocorrelation at seasonal lags (e.g., lag 12), indicating the presence of seasonality. The PACF shows significant partial autocorrelations at lower lags, suggesting AR components.

AR Model & Forecasting

Lag Selection and Unit Root Testing Optimal Lag Order Using BIC: ARIMA Model Selection Using BIC

Auto ARIMA Model Selection:

# Use auto.arima to find the best ARIMA model based on BIC
auto_model <- auto.arima(retail_ts, seasonal = TRUE, ic = "bic")
summary(auto_model)

## Series: retail_ts 
## ARIMA(2,1,0)(0,1,1)[12] 
## 
## Coefficients:
##           ar1      ar2     sma1
##       -0.2897  -0.2934  -0.7146
## s.e.   0.0511   0.0509   0.0365
## 
## sigma^2 = 96894256:  log likelihood = -3760.08
## AIC=7528.16   AICc=7528.27   BIC=7543.63
## 
## Training set error measures:
##                    ME     RMSE    MAE       MPE     MAPE      MASE        ACF1
## Training set 421.7904 9626.525 6265.5 0.0291091 1.894702 0.3587668 0.008958633

Model Specification: ARIMA(2,1,0)(0,1,1)[12]

model_selected <- Arima(retail_ts, order = c(2, 1, 0), seasonal = c(0, 1, 1))
summary(model_selected)

## Series: retail_ts 
## ARIMA(2,1,0)(0,1,1)[12] 
## 
## Coefficients:
##           ar1      ar2     sma1
##       -0.2897  -0.2934  -0.7146
## s.e.   0.0511   0.0509   0.0365
## 
## sigma^2 = 96894256:  log likelihood = -3760.08
## AIC=7528.16   AICc=7528.27   BIC=7543.63
## 
## Training set error measures:
##                    ME     RMSE    MAE       MPE     MAPE      MASE        ACF1
## Training set 421.7904 9626.525 6265.5 0.0291091 1.894702 0.3587668 0.008958633

checkresiduals(model_selected)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,0)(0,1,1)[12]
## Q* = 77.973, df = 21, p-value = 1.757e-08
## 
## Model df: 3.   Total lags used: 24

Overall, this output suggests that the ARIMA(2,1,0)(0,1,1)[12] model provides a good fit for the retail sales data, as evidenced by the low AIC/BIC values and reasonable error measures.
It indicates that the model effectively captures both the trend and seasonality in the data. we can compare this model’s AIC/BIC values with other candidate models to ensure it’s the best choice for forecasting.

Unit Root Test (Augmented Dickey-Fuller):

# Augmented Dickey-Fuller Test
adf_test <- adf.test(diff_retail_ts, alternative = "stationary")
print(adf_test)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff_retail_ts
## Dickey-Fuller = -10.566, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary

Interpretation of the Results Dickey-Fuller Statistic:
The value of the Dickey-Fuller statistic (-10.566) is significantly negative. In general, the more negative this statistic, the stronger the evidence against the null hypothesis of a unit root.
The p-value is 0.01, which is less than the common significance level of 0.05. This indicates that we can reject the null hypothesis.
Since we can reject the null hypothesis, we conclude that the differenced series (diff_retail_ts) is stationary.
This means that the original series (retail_ts) was difference stationary.

Selection & Comparison

Based on your analysis in parts 2 and 3, select 3 potential candidate ARIMA(p,d,q)(P,D,Q) models • Estimate each model and compute information criteria • Compare and contrast the results from each model • Discuss which model appears to be the best fit in-sample

Using the three models proposed in part 4 • Use a training period from January 1990-Dec. 2009 • Generate 1-step ahead recursive forecasts using each model for the out-of-sample period from Jan. 2010 to Jul. 2022 • Plot the forecasts against the actual time series during the out-of-sample period • Use RMSE to compare the forecasting performance of each model • Which model is best?

ARIMA Model Selection and Estimation Candidate ARIMA Models:

# Manually explore alternative models
model_1 <- Arima(retail_ts, order = c(1,1,0), seasonal = c(1,1,0))
model_2 <- Arima(retail_ts, order = c(0,1,1), seasonal = c(0,1,1))
model_3 <- Arima(retail_ts, order = c(1,1,1), seasonal = c(1,1,1))

# Compare Auto Selected model with other models
AIC(model_selected, model_1, model_2, model_3)  

##                df      AIC
## model_selected  4 7528.156
## model_1         3 7619.501
## model_2         3 7540.931
## model_3         5 7534.199

BIC(model_selected, model_1, model_2, model_3)

##                df      BIC
## model_selected  4 7543.633
## model_1         3 7631.109
## model_2         3 7552.539
## model_3         5 7553.546

Forecasting with ARIMA Models

# test periods
train_end <- c(2014, 12)  
test_start <- c(2015, 1) 
test_end <- c(2022, 7)    

train_data <- window(retail_ts, end = train_end)
test_data <- window(retail_ts, start = test_start, end = test_end)

# model on training data
model_selected <- Arima(train_data, order = c(2, 1, 0), seasonal = c(0, 1, 1))

# Generating forecasts for the length of the test data
forecasts <- forecast(model_selected, h = length(test_data))
plot(forecasts)
lines(test_data, col = 'red')  # Add actual test data to the plot

# Calculate accuracy metrics
actuals <- test_data  # Actual values from the test set
predicted <- forecasts$mean  # Predicted values from the model

rmse <- sqrt(mean((actuals - predicted)^2))
mae <- mean(abs(actuals - predicted))
mape <- mean(abs((actuals - predicted) / actuals)) * 100

cat("RMSE:", rmse, "\nMAE:", mae, "\nMAPE:", mape, "%\n")

## RMSE: 35135.85 
## MAE: 24215.66 
## MAPE: 5.009539 %

Training and forecasting ARIMA Models (Manually selected):

# Fit ARIMA models on the training set
train_model_1 <- Arima(train_data, order = c(1,1,0), seasonal = c(1,1,0))
train_model_2 <- Arima(train_data, order = c(0,1,1), seasonal = c(0,1,1))
train_model_3 <- Arima(train_data, order = c(1,1,1), seasonal = c(1,1,1))


# forecasts
fc_1 <- forecast(train_model_1, h = length(test_data))
fc_2 <- forecast(train_model_2, h = length(test_data))
fc_3 <- forecast(train_model_3, h = length(test_data))

plot_fc_1 <- autoplot(fc_1) + ggtitle("Forecast Model 1- ARIMA(1,1,0)(1,1,0)")
plot_fc_2 <- autoplot(fc_2) + ggtitle("Forecast Model 2- ARIMA(0,1,1)(0,1,1)")
plot_fc_3 <- autoplot(fc_3) + ggtitle("Forecast Model 3- ARIMA(1,1,1)(1,1,1)")

# Combination of All plots
plot_grid(plot_fc_1, plot_fc_2, plot_fc_3, ncol = 1, align = "hv", rel_widths = c(1, 1, 1), rel_heights = c(1, 1, 1))

RMSE Comparison

# Calculate RMSE for each model
rmse_1 <- sqrt(mean((test_data - fc_1$mean)^2, na.rm = TRUE))
rmse_2 <- sqrt(mean((test_data - fc_2$mean)^2, na.rm = TRUE))
rmse_3 <- sqrt(mean((test_data - fc_3$mean)^2, na.rm = TRUE))

# Print RMSE for each model
cat("RMSE of ARIMA(1,1,0):", rmse_1, "\n")

## RMSE of ARIMA(1,1,0): 31873.37

cat("RMSE of ARIMA(0,1,1):", rmse_2, "\n")

## RMSE of ARIMA(0,1,1): 35229.9

cat("RMSE of ARIMA(1,1,1):", rmse_3, "\n")

## RMSE of ARIMA(1,1,1): 36625.75

Conclusion and Recommendations
Predictability of U.S. Retail Sales
Based on the analysis conducted on U.S. Retail Sales data from January 1992 to July 2022, it is evident that retail sales can be predicted with a reasonable degree of accuracy using time series modeling techniques. The presence of trends, seasonal patterns, and cyclical fluctuations makes the data amenable to forecasting methodologies, particularly ARIMA models.
Recommended Model
After evaluating various ARIMA models, including ARIMA(2,1,0)(0,1,1)[12], the results indicated that this model provided a suitable balance between fit and complexity, as evidenced by its lower AIC and BIC values compared to other candidate models.

Md Mahmudul Hasan
MQIM 3760573
Faculty of Management
University of New Brunswick
mahmudul.hasan@unb.ca