required_packages <- c( “fpp3”, “readr”, “scales” )

missing_packages <- required_packages[ !required_packages %in% rownames(installed.packages())]

if (length(missing_packages) > 0) { install.packages(missing_packages) }

library(fpp3) library(readr) library(scales)

fred_url <- paste0( “https://fred.stlouisfed.org/graph/”, “fredgraph.csv?id=RSXFSN” )

retail_sales_raw <- read_csv( fred_url, na = “.”, show_col_types = FALSE )

names(retail_sales_raw)[1:2] <- c( “Date”, “Sales” )

retail_sales <- retail_sales_raw %>% transmute( Month = yearmonth(as.Date(Date)), Sales = as.numeric(Sales) ) %>% filter( !is.na(Month), !is.na(Sales) ) %>% arrange(Month) %>% as_tsibble( index = Month )

print(retail_sales)

train <- retail_sales %>% filter( Month < yearmonth(“2023 Jan”) )

test <- retail_sales %>% filter( Month >= yearmonth(“2023 Jan”), Month <= yearmonth(“2025 Dec”) )

if (nrow(train) == 0) { stop(“The training dataset is empty.”) }

if (nrow(test) == 0) { stop( paste( “The test dataset does not contain observations”, “from January 2023 through December 2025.” ) ) }

cat( “period:”, as.character(min(train\(Month)), "through", as.character(max(train\)Month)), “” )

cat( “Test period:”, as.character(min(test\(Month)), "through", as.character(max(test\)Month)), “” )

cat( “Training observations:”, nrow(train), “” )

cat( “Test observations:”, nrow(test), “” )

models <- train %>% model( ETS = ETS(Sales), NAIVE = NAIVE(Sales), SNAIVE = SNAIVE( Sales ~ lag(“year”) ) )

print(models)

base_forecasts <- models %>% forecast( new_data = test ) %>% as_tibble() %>% transmute( Month, .model, Forecast = .mean )

ensemble_forecast <- base_forecasts %>% group_by(Month) %>% summarise( Forecast = mean( Forecast, na.rm = TRUE ), .groups = “drop” ) %>% mutate( .model = “Ensemble” ) %>% select( Month, .model, Forecast )

all_forecasts <- bind_rows( base_forecasts, ensemble_forecast ) %>% arrange( .model, Month )

print(head(all_forecasts))

evaluation_data <- all_forecasts %>% left_join( test %>% as_tibble() %>% select( Month, Actual = Sales ), by = “Month” ) %>% mutate( Error = Actual - Forecast, Absolute_Error = abs(Error), Squared_Error = Error^2, Absolute_Percentage_Error = if_else( Actual == 0, NA_real_, abs(Error / Actual) * 100 ) )

print(head(evaluation_data))

accuracy_results_raw <- evaluation_data %>% group_by(.model) %>% summarise( RMSE = sqrt( mean( Squared_Error, na.rm = TRUE ) ), MAE = mean( Absolute_Error, na.rm = TRUE ), MAPE = mean( Absolute_Percentage_Error, na.rm = TRUE ), ME = mean( Error, na.rm = TRUE ), R_squared = { total_variation <- sum( ( Actual - mean( Actual, na.rm = TRUE ) )^2, na.rm = TRUE )

  if (total_variation == 0) {
    NA_real_
  } else {
    1 -
      sum(
        Squared_Error,
        na.rm = TRUE
      ) /
      total_variation
  }
},
.groups = "drop"

) %>% arrange(RMSE)

accuracy_results <- accuracy_results_raw %>% mutate( across( c( RMSE, MAE, MAPE, ME, R_squared ), function(x) round(x, 3) ) )

print(accuracy_results)

write_csv( accuracy_results, “RSXFSN_forecast_accuracy_results.csv” )

best_model <- accuracy_results %>% slice_min( order_by = RMSE, n = 1, with_ties = FALSE )

cat( “model based on lowest RMSE:” )

print(best_model)

figure1 <- ggplot() + geom_line( data = test, aes( x = Month, y = Sales ), linewidth = 1.2 ) + geom_line( data = all_forecasts, aes( x = Month, y = Forecast, linetype = .model ), linewidth = 0.8 ) + labs( title = “Actual and Forecast U.S. Retail Trade Sales”, subtitle = “RSXFSN test period: January 2023–December 2025”, x = “Month”, y = “Retail sales, millions of dollars”, linetype = “Model” ) + scale_y_continuous( labels = scales::label_comma() ) + theme_minimal()

print(figure1)

figure2 <- evaluation_data %>% ggplot( aes( x = Month ) ) + geom_line( aes( y = Actual ), linewidth = 1.1 ) + geom_line( aes( y = Forecast ), linewidth = 0.8, linetype = “dashed” ) + facet_wrap( ~ .model, ncol = 2 ) + labs( title = “Actual and Forecast Retail Trade Sales by Model”, subtitle = “Solid line represents actual sales; dashed line represents forecasts”, x = “Month”, y = “Retail sales, millions of dollars” ) + scale_y_continuous( labels = scales::label_comma() ) + theme_minimal()

print(figure2)

accuracy_long <- accuracy_results_raw %>% select( .model, RMSE, MAE, MAPE ) %>% pivot_longer( cols = c( RMSE, MAE, MAPE ), names_to = “Metric”, values_to = “Value” )

figure3 <- accuracy_long %>% ggplot( aes( x = .model, y = Value ) ) + geom_col() + facet_wrap( ~ Metric, scales = “free_y” ) + labs( title = “Forecast Accuracy Comparison”, subtitle = “Lower values indicate better performance”, x = “Forecasting model”, y = “Metric value” ) + theme_minimal()

print(figure3)

figure4 <- evaluation_data %>% ggplot( aes( x = Month, y = Error, linetype = .model ) ) + geom_hline( yintercept = 0, linetype = “dashed” ) + geom_line( linewidth = 0.8 ) + labs( title = “Forecast Errors During the Test Period”, subtitle = “Error equals actual sales minus forecast sales”, x = “Month”, y = “Forecast error, millions of dollars”, linetype = “Model” ) + scale_y_continuous( labels = scales::label_comma() ) + theme_minimal()

print(figure4)

ets_residual_plot <- models %>% select(ETS) %>% gg_tsresiduals()

print(ets_residual_plot)

ets_residuals <- models %>% select(ETS) %>% augment()

ljung_box_results <- ets_residuals %>% features( .innov, ljung_box, lag = 24, dof = 0 )

print(ljung_box_results)

ggsave( filename = “RSXFSN_figure1_actual_vs_forecasts.png”, plot = figure1, width = 11, height = 7, units = “in”, dpi = 300 )

ggsave( filename = “RSXFSN_figure2_forecasts_by_model.png”, plot = figure2, width = 11, height = 8, units = “in”, dpi = 300 )

ggsave( filename = “RSXFSN_figure3_accuracy_comparison.png”, plot = figure3, width = 10, height = 7, units = “in”, dpi = 300 )

ggsave( filename = “RSXFSN_figure4_forecast_errors.png”, plot = figure4, width = 11, height = 7, units = “in”, dpi = 300 )

ggsave( filename = “RSXFSN_figure5_ets_residuals.png”, plot = ets_residual_plot, width = 10, height = 8, units = “in”, dpi = 300 )

write_csv( retail_sales %>% as_tibble(), “RSXFSN_retail_sales_data.csv” )

write_csv( evaluation_data, “RSXFSN_forecast_evaluation_data.csv” )

write_csv( all_forecasts, “RSXFSN_all_model_forecasts.csv” )

saveRDS( train, “RSXFSN_training_data.rds” )

saveRDS( test, “RSXFSN_test_data.rds” )

saveRDS( models, “RSXFSN_forecast_models.rds” )

saveRDS( all_forecasts, “RSXFSN_all_forecasts.rds” )

cat( “analysis completed successfully.” )

cat( “The data, fitted models, forecasts, accuracy results, and five figures were saved to the working directory.” )