---
title: "A Comparative Analysis of U.S. Retail Sales Forecasting Models: Implications for Investment Decisions"
output:
  html_document: 
    code_download: true
    highlight: espresso
date: "`r Sys.Date()`"
editor_options: 
  markdown: 
    wrap: 72
---

#### **Md Mahmudul Hasan**
<style type="text/css"> body, td {font-size: 14px;} code.r{font-size: 12px;} pre {font-size: 14px} </style>
<script src="https://code.jquery.com/jquery-3.7.1.slim.min.js"></script>
<script type="text/javascript">
  $(document).ready(function() {
  $('a').attr('target', '_blank');
  });
</script>
<style>
.boxed {
  background: #232023;
  color: white;
  border: 0px solid #535353;
  margin: 0px auto;
  width: auto;
  padding: 10px;
  border-radius: 0px;
}
</style>

#  {.tabset}

```{r setup, include=FALSE}
knitr::opts_chunk$set(
	fig.align = "center",
	message = FALSE,
	warning = FALSE
)
library(fpp2)
library(tidyquant)
library(tidyverse)
library(quantmod)
library(dplyr)
library(ggthemes)
library(curl)
library(readxl)
library(cowplot)
library(ggpubr)
library(plotly)
library(ggplot2)
library(reshape2)
library(gt)
library(gtExtras)
library(knitr)
library(ggfortify)
library(moments)
library(skimr)
library(psych)
library(kableExtra)
library(forecast)
library(tseries)
library(tibble)
```

## Assignment: 
<br>
*Context:* 
<br>
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.<br>
<br>
*Assignment:* 
<br>
- Within the above context, write a report suitable for sharing with your team on the modelling and predictability of Retail Sales.
<br>
- Retail Sales: Non-seasonally adjusted Advance Retail Sales: Retail Trade from FRED from January 1990 to July 2022.

## Data
Fetching Data from FRED

```{r message=FALSE, warning=FALSE}
# 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"))

```
## Descriptive Table
<div class="boxed">
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.
</div>

```{r message=FALSE, warning=FALSE}
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")
```

**Notable Features** 
<br>
<br>
*Volatility:* 
<br>
<br>
The high standard deviation (106797.7) relative to the mean (316065.8) suggests significant volatility, which could have implications for forecasting and risk assessment.
<br>
<br>
*Potential for Positive Outliers:*
<br>
<br>
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.
<br>
<br>
*Symmetry in Central Tendency:* 
<br>
<br>
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.
<br>
<br>
*Forecasting Implications:* 
<br>
<br>
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.
<br>
<br>

## Time Series & ACF/PACF

<div class="boxed">
Time Series Plot and ACF/PACF
Generate a time series plot and an ACF and PACF chart.\
<br>
- 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\
</div>
<br>
**Time Series Plot and ACF/PACF**
<br>
```{r message=FALSE, warning=FALSE}
# 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

```
<br>
Time Series Plot and Autocorrelation Analysis Time Series Plot:
<br>
**Notable Observations:** 
<br>
- 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. \
<br>
**ACF and PACF Charts (Before Differencing):**
<br>
```{r message=FALSE, warning=FALSE}
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.*

<br>
<br>

**Differencing to Address Non-Stationarity**
<br>
ACF and PACF Charts (After First Differencing):
<br>

```{r message=FALSE, warning=FALSE}
# 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

```
<br>
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. 
<br>
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:**
```{r message=FALSE, warning=FALSE}
# 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)

```

<!-- Model Specification: ARIMA(2,1,0)(0,1,1)[12] -->
<!-- ARIMA(2,1,0): -->

<!-- AR (AutoRegressive) order = 2: This indicates that the model uses the past 2 values of the time series to predict future values.\ -->
<!-- I (Integrated) order = 1: This suggests that the data was differenced once to achieve stationarity.\ -->
<!-- MA (Moving Average) order = 0: This means that there is no moving average component in this part of the model.\ -->

<!-- <br> -->
<!-- <br> -->
<!-- Seasonal part (0,1,1):\ -->
<!-- P (Seasonal AR) order = 0: No seasonal autoregressive terms.\ -->
<!-- D (Seasonal differencing) order = 1: Seasonal differencing was applied once.\ -->
<!-- Q (Seasonal MA) order = 1: A seasonal moving average component was included.\ -->
<!-- <br> -->
<!-- 12: Indicates the seasonality period (12 months for monthly data).\ -->


<!-- Coefficients\ -->
<!-- ar1: -0.2897 (with a standard error of 0.0511)\ -->
<!-- ar2: -0.2934 (with a standard error of 0.0509)\ -->
<!-- sma1: -0.7146 (with a standard error of 0.0365)\ -->
<!-- These coefficients represent the weights applied to the previous values and the seasonal moving average:\ -->

<!-- The ar1 and ar2 coefficients are negative, indicating that the current value is inversely related to its two previous values.\ -->
<!-- The sma1 coefficient being negative suggests that the current value is negatively influenced by the previous seasonal error term.\ -->
<!-- 3. Model Fit Statistics\ -->
<!-- sigma^2 = 96894256: This is the estimated variance of the residuals (the errors from the model). A lower value indicates a better fit, though it should be compared relative to other models.\ -->
<!-- log likelihood = -3760.08: This statistic helps evaluate the model fit; higher values (less negative) indicate a better fit.\ -->
<!-- 4. Information Criteria\ -->
<!-- AIC = 7528.16: Akaike Information Criterion. Lower values indicate a better fit, taking into account the number of parameters used in the model. It’s useful for model comparison.\ -->
<!-- AICc = 7528.27: AIC corrected for small sample sizes. This is slightly higher than AIC but has a similar interpretation.\ -->
<!-- BIC = 7543.63: Bayesian Information Criterion. Like AIC, but with a greater penalty for the number of parameters. Lower values are better, and it can be used for model selection.\ -->
<!-- 5. Training Set Error Measures\ -->
<!-- ME (Mean Error) = 421.7904: The average of the residuals. Ideally, it should be close to zero, indicating no systematic bias in predictions.\ -->
<!-- RMSE (Root Mean Square Error) = 9626.525: This measures the standard deviation of the residuals. A lower value indicates a better fit.\ -->
<!-- MAE (Mean Absolute Error) = 6265.5: The average absolute error. It provides a straightforward measure of prediction accuracy.\ -->
<!-- MPE (Mean Percentage Error) = 0.0291091: The average of the percentage errors. A value near zero indicates good accuracy without systematic bias.\ -->
<!-- MAPE (Mean Absolute Percentage Error) = 1.894702: Average percentage error; values below 5% are typically considered excellent.\ -->
<!-- MASE (Mean Absolute Scaled Error) = 0.3587668: A measure that compares the forecast error to a naive forecast. Values less than 1 indicate better performance than a naive model.\ -->
<!-- ACF1 = 0.008958633: The autocorrelation of the residuals at lag 1. A value near zero suggests no significant autocorrelation, indicating the model has captured the relationship in the data well.\ -->
**Model Specification: ARIMA(2,1,0)(0,1,1)[12]**
<br>
```{r message=FALSE, warning=FALSE}
model_selected <- Arima(retail_ts, order = c(2, 1, 0), seasonal = c(0, 1, 1))
summary(model_selected)

checkresiduals(model_selected)
```
<br>
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.
<br>
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.
<br>
<br>
**Unit Root Test (Augmented Dickey-Fuller):**
<br>
```{r message=FALSE, warning=FALSE}
# Augmented Dickey-Fuller Test
adf_test <- adf.test(diff_retail_ts, alternative = "stationary")
print(adf_test)

```
<br>
**Interpretation of the Results Dickey-Fuller Statistic:**
<br>
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.
<br>
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.
<br>
**Since we can reject the null hypothesis, we conclude that the differenced series (diff_retail_ts) is stationary.** 
<br>
**This means that the original series (retail_ts) was difference stationary.**


## Selection & Comparison
<div class="boxed">
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 
</div>
<br>
<div class="boxed">
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?
</div>
<br>
**ARIMA Model Selection and Estimation Candidate ARIMA Models:**
<br>
```{r message=FALSE, warning=FALSE}
# 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)  
BIC(model_selected, model_1, model_2, model_3)
```

<br>
Forecasting with ARIMA Models
<br>
```{r message=FALSE, warning=FALSE}
# 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")

```
<br>
**Training and forecasting ARIMA Models (Manually selected):**
<br>
```{r message=FALSE, warning=FALSE}

# 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))


```
<br>
**RMSE Comparison**
<br>
```{r message=FALSE, warning=FALSE}
# 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")
cat("RMSE of ARIMA(0,1,1):", rmse_2, "\n")
cat("RMSE of ARIMA(1,1,1):", rmse_3, "\n")

```
<br>
**Conclusion and Recommendations**
<br>
*Predictability of U.S. Retail Sales*
<br>
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.
<br>
**Recommended Model**
<br>
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.


<br>

<span>  <span style="color: Steelblue;">**Md Mahmudul Hasan**</span> \
MQIM 3760573\
Faculty of Management \
University of New Brunswick \
mahmudul.hasan\@unb.ca \

<br>