library(fpp3)
library(fredr)
library(dplyr)
library(tseries)
library(fabletools)
library(ggplot2)
library(readr)
library(tidyverse)
library(tsibble)
library(feasts)
library(scales)
library(patchwork)
library(ggtime)
library(tseries)
fredr_has_key()[1] TRUE#Retail Sales: Retail Trade (MRTSSM44000USN)
Retail_Sales <- fredr(series_id = "MRTSSM44000USN") |>
transmute(month = yearmonth(date),value) |>
as_tsibble(index = month)
autoplot(Retail_Sales, value) +
labs(title = "Retail Sales: Retail Trade Millions of DOllars, U.S.",x = "Month",y = "Millions of Dollars") #Not Stationary, Upward Trend With Extreme Seasonality
n <- nrow(Retail_Sales)
train <- Retail_Sales |> slice(1:(n-12))
test <- Retail_Sales |> slice((n-11):n)ets_fit <- train |> model(ETS(value))
arima_fit <- train |> model(ARIMA(value))
ets_fc <- forecast(ets_fit, h = 12)
arima_fc <- forecast(arima_fit, h = 12)
#ETS Plot
autoplot(Retail_Sales |> filter_index("2022 Jan" ~ .), value) +
autolayer(test, value, series = "Actual", color = "black", size = 1) +
autolayer(ets_fc, .mean, series = "ETS Forecast", color = "blue", size = 1) +
autolayer(ets_fc, level = 95, alpha = 0.25, fill = "blue") +
labs(title = "Zoomed ETS Forecast with 95% Prediction Interval",x = "Month",y = "Millions of Dollars",color = "Series")
#ARIMA Plot
autoplot(Retail_Sales |> filter_index("2022 Jan" ~ .), value) +
autolayer(test, value, series = "Actual", color = "black", size = 1) +
autolayer(arima_fc, .mean, series = "ARIMA Forecast", color = "red", size = 1) +
autolayer(arima_fc, level = 95, alpha = 0.25, fill = "red") +
labs(title = "Zoomed ARIMA Forecast with 95% Prediction Interval",x = "Month",y = "Millions of Dollars",color = "Series")
accuracy(ets_fc, test)# A tibble: 1 × 10
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ETS(value) Test -982. 8717. 7593. -0.193 1.21 NaN NaN -0.0350accuracy(arima_fc, test)# A tibble: 1 × 10
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ARIMA(value) Test -5727. 10476. 9203. -0.975 1.49 NaN NaN -0.461#Add Personal Disposable Income Variable - FRED Series DSPIC96
income <- fredr(series_id = "DSPIC96") |>
transmute(month = yearmonth(date), income = value) |>
as_tsibble(index = month)
#Combine With Retail Sales
data <- Retail_Sales |>
inner_join(income, by = "month")
#Test/Train Split
n <- nrow(data)
train <- data |> slice(1:(n-12))
test <- data |> slice((n-11):n)
#ARIMAX
arimax_fit <- train |> model(ARIMA(value ~ income))
#Forecasts
arima_fc <- forecast(arima_fit, h = 12)
arimax_fc <- forecast(arimax_fit, new_data = test)
#Accuracy Merics
accuracy(arima_fc, test)# A tibble: 1 × 10
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ARIMA(value) Test -5727. 10476. 9203. -0.975 1.49 NaN NaN -0.461accuracy(arimax_fc, test)# A tibble: 1 × 10
.model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ARIMA(value ~ income) Test -726. 9323. 7841. -0.183 1.25 NaN NaN -0.303arimax_fc <- forecast(arimax_fit, new_data = test)
#ARIMA vs Actual Plot
autoplot(data |> filter_index("2022 Jan" ~ .), value) +
autolayer(test, value, series = "Actual", color = "black", size = 1) +
autolayer(arima_fc, .mean, series = "ARIMA", color = "red", size = 1) +
labs(title = "Zoomed Forecast: ARIMA vs Actual",x = "Month",y = "Retail Sales",color = "Series")
#ARIMAX vs. Actual Plot
autoplot(data |> filter_index("2022 Jan" ~ .), value) +
autolayer(test, value, series = "Actual", color = "black", size = 1) +
autolayer(arimax_fc, .mean, series = "ARIMAX", color = "blue", size = 1) +
labs(title = "Zoomed Forecast: ARIMAX vs Actual",x = "Month",y = "Retail Sales",color = "Series")