Final-Code

options(stringsAsFactors = FALSE, scipen = 999) set.seed(20260813)

required_packages <- c( “data.table”, “forecast”, “xgboost”, “nnet”, “ggplot2”, “zoo”, “scales”, “patchwork”, “tseries” ) missing_packages <- required_packages[ !vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)] if (length(missing_packages) > 0) { install.packages(missing_packages, repos = “https://cloud.r-project.org”) }

suppressPackageStartupMessages({ library(data.table) library(forecast) library(xgboost) library(nnet) library(ggplot2) library(zoo) library(scales) library(patchwork) library(tseries) })

project_dir <- normalizePath(getwd(), mustWork = TRUE) data_dir <- file.path(project_dir, “data”) output_dir <- file.path(project_dir, “outputs”) figure_dir <- file.path(output_dir, “figures”) dir.create(data_dir, recursive = TRUE, showWarnings = FALSE) dir.create(figure_dir, recursive = TRUE, showWarnings = FALSE)

data_path <- file.path(data_dir, “rossmann.csv”) holdout_days <- 42L validation_days <- 42L

if (!file.exists(data_path)) { zip_path <- file.path(data_dir, “rossmann_data.zip”) data_url <- paste0( “https://raw.githubusercontent.com/raviatkumar/”, “Retail-Sales-Prediction-Capstone-Project/main/”, “DataSet/Rossmann%20Stores%20Data.zip” ) message(“Dataset not found locally; downloading the public mirror…”) download.file(data_url, zip_path, mode = “wb”, quiet = FALSE) unzip(zip_path, exdir = data_dir) extracted <- list.files(data_dir, pattern = “\.csv$”, full.names = TRUE) if (length(extracted) != 1L) stop(“Expected exactly one CSV in the archive.”) file.rename(extracted, data_path) unlink(zip_path) }

raw <- fread(data_path, na.strings = c(““,”NA”), showProgress = FALSE) required_columns <- c( “Store”, “DayOfWeek”, “Date”, “Sales”, “Customers”, “Open”, “Promo”, “StateHoliday”, “SchoolHoliday” ) if (!all(required_columns %in% names(raw))) { stop(“The input file does not contain the expected Rossmann columns.”) }

raw[, Date := as.Date(Date)] raw[, StateHoliday := as.character(StateHoliday)] setorder(raw, Date, Store)

missing_summary <- data.table( Variable = names(raw), Missing = vapply(raw, function(x) sum(is.na(x)), numeric(1)) ) fwrite(missing_summary, file.path(output_dir, “missing_values.csv”))

daily <- raw[, .( Sales = sum(Sales, na.rm = TRUE), OpenStores = sum(Open, na.rm = TRUE), PromoStores = sum(Promo, na.rm = TRUE), SchoolHolidayStores = sum(SchoolHoliday, na.rm = TRUE), StateHolidayStores = sum(StateHoliday != “0”, na.rm = TRUE), StoreRows = .N, DescriptiveCustomers = sum(Customers, na.rm = TRUE) ), by = Date] setorder(daily, Date)

if (anyNA(daily) || anyDuplicated(daily$Date)) { stop("Daily aggregation produced missing or duplicate dates.") } expected_dates <- seq(min(daily$Date), max(daily$Date), by = "day") if (length(daily$Date) != length(expected_dates) || !all(daily$Date == expected_dates)) { stop(“The daily series is not complete; inspect missing dates before modeling.”) }

global_min_date <- min(daily$Date)

add_known_features <- function(dt) { z <- copy(as.data.table(dt)) z[, TimeIndex := as.integer(Date - global_min_date) + 1L] z[, DayOfWeek := as.integer(format(Date, “%u”))] z[, Month := as.integer(format(Date, “%m”))] z[, DayOfMonth := as.integer(format(Date, “%d”))] z[, DayOfYear := as.integer(format(Date, “%j”))] full_span <- as.integer(max(daily$Date) - global_min_date) + 1L z[, Trend := TimeIndex / full_span] z[, OpenShare := OpenStores / pmax(StoreRows, 1)] z[, PromoShare := PromoStores / pmax(StoreRows, 1)] z[, SchoolHolidayShare := SchoolHolidayStores / pmax(StoreRows, 1)] z[, StateHolidayShare := StateHolidayStores / pmax(StoreRows, 1)] z[, LogOpenStores := log1p(OpenStores)]

for (k in 1:3) { z[, paste0(“W_Sin”, k) := sin(2 * pi * k * TimeIndex / 7)] z[, paste0(“W_Cos”, k) := cos(2 * pi * k * TimeIndex / 7)] } for (k in 1:5) { z[, paste0(“A_Sin”, k) := sin(2 * pi * k * TimeIndex / 365.25)] z[, paste0(“A_Cos”, k) := cos(2 * pi * k * TimeIndex / 365.25)] } for (d in 1:7) z[, paste0(“DOW”, d) := as.integer(DayOfWeek == d)] for (m in 1:12) z[, paste0(“Month”, m) := as.integer(Month == m)] z }

daily <- add_known_features(daily)

make_ml_features <- function(dt) { z <- add_known_features(dt) lag_days <- c(1, 7, 14, 28, 365) for (lag_n in lag_days) { z[, paste0(“Lag”, lag_n) := shift(Sales, lag_n)] } previous_sales <- shift(z$Sales, 1) z[, RollMean7 := zoo::rollapplyr(previous_sales, 7, mean, fill = NA)] z[, RollMean14 := zoo::rollapplyr(previous_sales, 14, mean, fill = NA)] z[, RollMean28 := zoo::rollapplyr(previous_sales, 28, mean, fill = NA)] z[, RollSD7 := zoo::rollapplyr(previous_sales, 7, sd, fill = NA)] z[, RollSD28 := zoo::rollapplyr(previous_sales, 28, sd, fill = NA)]

transform_cols <- c( paste0(“Lag”, lag_days), “RollMean7”, “RollMean14”, “RollMean28”, “RollSD7”, “RollSD28” ) for (nm in transform_cols) z[, (nm) := log1p(pmax(get(nm), 0))] z[, LogSales := log1p(Sales)] z }

dynamic_feature_names <- c( “LogOpenStores”, “PromoShare”, “SchoolHolidayShare”, “StateHolidayShare”, “Trend”, paste0(“W_Sin”, 1:3), paste0(“W_Cos”, 1:3), paste0(“A_Sin”, 1:5), paste0(“A_Cos”, 1:5) )

ml_feature_names <- c( “Lag1”, “Lag7”, “Lag14”, “Lag28”, “Lag365”, “RollMean7”, “RollMean14”, “RollMean28”, “RollSD7”, “RollSD28”, “LogOpenStores”, “PromoShare”, “SchoolHolidayShare”, “StateHolidayShare”, “Trend”, paste0(“W_Sin”, 1:3), paste0(“W_Cos”, 1:3), paste0(“A_Sin”, 1:5), paste0(“A_Cos”, 1:5), paste0(“DOW”, 1:7), paste0(“Month”, 1:12) )

max_date <- max(daily$Date) test_start <- max_date - (holdout_days - 1L) validation_start <- test_start - validation_days train_end <- validation_start - 1L

train_daily <- daily[Date <= train_end] validation_daily <- daily[Date >= validation_start & Date < test_start] pretest_daily <- daily[Date < test_start] test_daily <- daily[Date >= test_start]

stopifnot( nrow(validation_daily) == validation_days, nrow(test_daily) == holdout_days )

fit_dynamic <- function(history) { xreg <- as.matrix(history[, ..dynamic_feature_names]) forecast::auto.arima( log1p(history$Sales), xreg = xreg, seasonal = FALSE, stepwise = TRUE, approximation = FALSE, max.p = 5, max.q = 5, max.order = 8, allowdrift = TRUE, allowmean = TRUE ) }

predict_dynamic <- function(model, future_known) { newxreg <- as.matrix(future_known[, ..dynamic_feature_names]) fc <- forecast::forecast(model, xreg = newxreg, h = nrow(future_known), level = 95) data.table( Prediction = pmax(expm1(as.numeric(fc$mean)), 0), Lower95 = pmax(expm1(as.numeric(fc$lower[, 1])), 0), Upper95 = pmax(expm1(as.numeric(fc$upper[, 1])), 0) ) }

complete_ml_rows <- function(history) { f <- make_ml_features(history) f[complete.cases(f[, c(“LogSales”, ml_feature_names), with = FALSE])] }

fit_xgb <- function(history) { f <- complete_ml_rows(history) x_train <- as.matrix(f[, ..ml_feature_names]) storage.mode(x_train) <- “double” dtrain <- xgb.DMatrix( data = x_train, label = f$LogSales ) model <- xgb.train( params = list( objective = “reg:squarederror”, eval_metric = “rmse”, eta = 0.03, max_depth = 4, min_child_weight = 5, subsample = 0.85, colsample_bytree = 0.85, lambda = 1, alpha = 0.05, nthread = 2 ), data = dtrain, nrounds = 650, verbose = 0 ) list(model = model, features = ml_feature_names) }

fit_mlp <- function(history) { f <- complete_ml_rows(history) x <- as.matrix(f[, ..ml_feature_names]) storage.mode(x) <- “double” x_center <- colMeans(x) x_scale <- apply(x, 2, sd) x_scale[!is.finite(x_scale) | x_scale == 0] <- 1 xs <- sweep(sweep(x, 2, x_center, “-”), 2, x_scale, “/”) y <- f$LogSales y_center <- mean(y) y_scale <- sd(y) ys <- (y - y_center) / y_scale set.seed(20260813) model <- nnet::nnet( x = xs, y = ys, size = 12, decay = 0.02, linout = TRUE, maxit = 1500, MaxNWts = 10000, trace = FALSE ) list( model = model, features = ml_feature_names, x_center = x_center, x_scale = x_scale, y_center = y_center, y_scale = y_scale ) }

predict_one_ml <- function(fitted_model, feature_row, model_type) { x <- as.matrix(feature_row[, fitted_model$features, with = FALSE]) storage.mode(x) <- "double" if (model_type == "xgb") { log_prediction <- as.numeric(predict(fitted_model$model, xgb.DMatrix(x))) } else { xs <- sweep(sweep(x, 2, fitted_model$x_center, "-"), 2, fitted_model$x_scale, “/”) scaled_prediction <- as.numeric(predict(fitted_model$model, xs)) log_prediction <- scaled_prediction * fitted_model$y_scale + fitted_model$y_center } pmax(expm1(log_prediction), 0) }

recursive_ml_forecast <- function(fitted_model, history, future_known, model_type) { hist <- copy(history) predictions <- numeric(nrow(future_known)) for (i in seq_len(nrow(future_known))) { next_row <- copy(future_known[i]) next_row[, Sales := NA_real_] next_row[, DescriptiveCustomers := NA_real_] combined <- rbindlist(list(hist, next_row), fill = TRUE, use.names = TRUE) feature_row <- tail(make_ml_features(combined), 1) if (!all(complete.cases(feature_row[, ..ml_feature_names]))) { stop(“Incomplete recursive feature row; check lag history.”) } predictions[i] <- predict_one_ml(fitted_model, feature_row, model_type) next_row[, Sales := predictions[i]] hist <- rbindlist(list(hist, next_row), fill = TRUE, use.names = TRUE) } predictions }

seasonal_naive_forecast <- function(history, horizon) { rep(tail(history$Sales, 7), length.out = horizon) }

metric_table <- function(actual, prediction, model_name, split_name) { data.table( Split = split_name, Model = model_name, RMSE = sqrt(mean((actual - prediction)^2)), MAE = mean(abs(actual - prediction)), MAPE = mean(abs((actual - prediction) / actual)) * 100, R_squared = 1 - sum((actual - prediction)^2) / sum((actual - mean(actual))^2) ) }

message(“Fitting validation models…”) dynamic_validation_model <- fit_dynamic(train_daily) dynamic_validation <- predict_dynamic(dynamic_validation_model, validation_daily)$Prediction

xgb_validation_model <- fit_xgb(train_daily) xgb_validation <- recursive_ml_forecast( xgb_validation_model, train_daily, validation_daily, “xgb” )

mlp_validation_model <- fit_mlp(train_daily) mlp_validation <- recursive_ml_forecast( mlp_validation_model, train_daily, validation_daily, “mlp” )

naive_validation <- seasonal_naive_forecast(train_daily, nrow(validation_daily))

validation_forecasts <- data.table( Date = validation_daily$Date, Actual = validation_daily$Sales, Dynamic = dynamic_validation, XGBoost = xgb_validation, MLP = mlp_validation, SeasonalNaive = naive_validation )

validation_metrics <- rbindlist(list( metric_table(validation_forecasts$Actual, validation_forecasts$Dynamic, “Dynamic regression”, “Validation”), metric_table(validation_forecasts$Actual, validation_forecasts$XGBoost, “XGBoost”, “Validation”), metric_table(validation_forecasts$Actual, validation_forecasts$MLP, “MLP neural network”, “Validation”), metric_table(validation_forecasts$Actual, validation_forecasts$SeasonalNaive, “Seasonal naive”, “Validation”) ))

component_rows <- validation_metrics[Model != “Seasonal naive”] ensemble_weights <- component_rows[, .( Model, Weight = (1 / RMSE) / sum(1 / RMSE) )] weight_lookup <- setNames(ensemble_weights$Weight, ensemble_weights$Model) validation_forecasts[, Ensemble := weight_lookup[[“Dynamic regression”]] * Dynamic + weight_lookup[[“XGBoost”]] * XGBoost + weight_lookup[[“MLP neural network”]] * MLP] validation_metrics <- rbind( validation_metrics, metric_table(validation_forecasts$Actual, validation_forecasts$Ensemble, “Weighted ensemble”, “Validation”) )

message(“Refitting final models and forecasting the six-week holdout…”) dynamic_final_model <- fit_dynamic(pretest_daily) dynamic_test_full <- predict_dynamic(dynamic_final_model, test_daily)

xgb_final_model <- fit_xgb(pretest_daily) xgb_test <- recursive_ml_forecast(xgb_final_model, pretest_daily, test_daily, “xgb”)

mlp_final_model <- fit_mlp(pretest_daily) mlp_test <- recursive_ml_forecast(mlp_final_model, pretest_daily, test_daily, “mlp”)

naive_test <- seasonal_naive_forecast(pretest_daily, nrow(test_daily))

test_forecasts <- data.table( Date = test_daily$Date, Actual = test_daily$Sales, Dynamic = dynamic_test_full$Prediction, XGBoost = xgb_test, MLP = mlp_test, SeasonalNaive = naive_test ) test_forecasts[, Ensemble := weight_lookup[[“Dynamic regression”]] * Dynamic + weight_lookup[[“XGBoost”]] * XGBoost + weight_lookup[[“MLP neural network”]] * MLP]

validation_abs_error <- abs(validation_forecasts$Actual - validation_forecasts$Ensemble) interval_radius <- as.numeric(quantile(validation_abs_error, 0.95, type = 8)) test_forecasts[, :=( EnsembleLower95 = pmax(Ensemble - interval_radius, 0), EnsembleUpper95 = Ensemble + interval_radius )]

test_metrics <- rbindlist(list( metric_table(test_forecasts$Actual, test_forecasts$Dynamic, “Dynamic regression”, “Test”), metric_table(test_forecasts$Actual, test_forecasts$XGBoost, “XGBoost”, “Test”), metric_table(test_forecasts$Actual, test_forecasts$MLP, “MLP neural network”, “Test”), metric_table(test_forecasts$Actual, test_forecasts$Ensemble, “Weighted ensemble”, “Test”), metric_table(test_forecasts$Actual, test_forecasts$SeasonalNaive, “Seasonal naive”, “Test”) )) all_metrics <- rbind(validation_metrics, test_metrics)

dynamic_residuals <- as.numeric(residuals(dynamic_final_model)) dynamic_residuals <- dynamic_residuals[is.finite(dynamic_residuals)] ljung_fitdf <- min(length(coef(dynamic_final_model)), 10) ljung_box <- Box.test( dynamic_residuals, lag = 21, type = “Ljung-Box”, fitdf = ljung_fitdf )

arch_lag <- 7L arch_matrix <- embed(dynamic_residuals^2, arch_lag + 1L) arch_model <- lm(arch_matrix[, 1] ~ arch_matrix[, -1]) arch_lm_stat <- nrow(arch_matrix) * summary(arch_model)$r.squared arch_lm_p <- pchisq(arch_lm_stat, df = arch_lag, lower.tail = FALSE)

garch_fit <- tryCatch( suppressWarnings(tseries::garch(dynamic_residuals, order = c(1, 1), trace = FALSE)), error = function(e) NULL ) garch_status <- if (is.null(garch_fit)) “GARCH fit did not converge” else “GARCH(1,1) fit converged” garch_coefficients <- if (is.null(garch_fit)) { data.table(Parameter = character(), Estimate = numeric()) } else { data.table(Parameter = names(coef(garch_fit)), Estimate = as.numeric(coef(garch_fit))) }

garch_variance_forecast <- data.table() if (!is.null(garch_fit)) { gcoef <- coef(garch_fit) omega <- unname(gcoef[“a0”]) alpha <- unname(gcoef[“a1”]) beta <- unname(gcoef[“b1”]) sigma_history <- as.numeric(predict(garch_fit, newdata = dynamic_residuals)[, 1]) last_sigma2 <- tail(sigma_history[is.finite(sigma_history)], 1)^2 variance_h <- numeric(holdout_days) variance_h[1] <- omega + alpha * tail(dynamic_residuals, 1)^2 + beta * last_sigma2 if (holdout_days > 1) { for (h in 2:holdout_days) { variance_h[h] <- omega + (alpha + beta) * variance_h[h - 1] } } center_log <- log1p(test_forecasts$Dynamic) test_forecasts[, `:=`( DynamicGARCHLower95 = pmax(expm1(center_log - 1.96 * sqrt(variance_h)), 0), DynamicGARCHUpper95 = pmax(expm1(center_log + 1.96 * sqrt(variance_h)), 0) )] garch_variance_forecast <- data.table( Date = test_forecasts$Date, ConditionalVariance = variance_h, ConditionalSD = sqrt(variance_h), Lower95 = test_forecasts$DynamicGARCHLower95, Upper95 = test_forecasts$DynamicGARCHUpper95 ) }

diagnostics <- data.table( Test = c(“Ljung-Box residual autocorrelation”, “ARCH LM volatility clustering”), Statistic = c(as.numeric(ljung_box$statistic), arch_lm_stat), Degrees_of_freedom = c(as.numeric(ljung_box$parameter), arch_lag), P_value = c(ljung_box$p.value, arch_lm_p) )

feature_importance <- as.data.table( xgb.importance(feature_names = ml_feature_names, model = xgb_final_model$model) ) if (nrow(feature_importance) == 0) { feature_importance <- data.table(Feature = character(), Gain = numeric()) }

theme_report <- theme_minimal(base_size = 11, base_family = “DejaVu Sans”) + theme( plot.title = element_text(face = “bold”, color = “#17365D”, size = 14), plot.subtitle = element_text(color = “#555555”, size = 10), panel.grid.minor = element_blank(), legend.position = “bottom”, axis.title = element_text(face = “bold”) )

fig1 <- ggplot(daily, aes(Date, Sales)) + annotate( “rect”, xmin = validation_start, xmax = test_start - 1, ymin = -Inf, ymax = Inf, fill = “#F4B183”, alpha = 0.22 ) + annotate( “rect”, xmin = test_start, xmax = max_date, ymin = -Inf, ymax = Inf, fill = “#5B9BD5”, alpha = 0.18 ) + geom_line(color = “#244A6C”, linewidth = 0.45) + scale_y_continuous(labels = label_number(scale = 1e-6, suffix = “M”)) + labs( title = “Rossmann chain-wide daily sales”, subtitle = “Orange = validation window; blue = final six-week holdout”, x = NULL, y = “Daily sales” ) + theme_report ggsave(file.path(figure_dir, “figure_1_daily_sales.png”), fig1, width = 8.2, height = 4.6, dpi = 300)

profile <- raw[Open == 1, .( MeanSales = mean(Sales), Observations = .N ), by = .(DayOfWeek, Promo)] profile[, Day := factor( DayOfWeek, levels = 1:7, labels = c(“Mon”, “Tue”, “Wed”, “Thu”, “Fri”, “Sat”, “Sun”) )] profile[, Promotion := factor(Promo, levels = c(0, 1), labels = c(“No promotion”, “Promotion”))] fig2 <- ggplot(profile, aes(Day, MeanSales, fill = Promotion)) + geom_col(position = position_dodge(width = 0.78), width = 0.68) + scale_fill_manual(values = c(“#A5A5A5”, “#ED7D31”)) + scale_y_continuous(labels = label_dollar(prefix = ““, accuracy = 1)) + labs( title =”Average open-store sales by weekday and promotion status”, subtitle = “Promotions are associated with higher sales on every trading day”, x = NULL, y = “Mean store sales”, fill = NULL ) + theme_report ggsave(file.path(figure_dir, “figure_2_weekday_promotion.png”), fig2, width = 8.2, height = 4.6, dpi = 300)

forecast_long <- melt( test_forecasts, id.vars = c(“Date”, “Actual”, “EnsembleLower95”, “EnsembleUpper95”), measure.vars = c(“Dynamic”, “XGBoost”, “MLP”, “Ensemble”), variable.name = “Model”, value.name = “Forecast” ) fig3 <- ggplot() + geom_ribbon( data = test_forecasts, aes(Date, ymin = EnsembleLower95, ymax = EnsembleUpper95), fill = “#5B9BD5”, alpha = 0.14 ) + geom_line(data = forecast_long, aes(Date, Forecast, color = Model), linewidth = 0.65) + geom_line(data = test_forecasts, aes(Date, Actual), color = “#111111”, linewidth = 0.9) + scale_color_manual(values = c( Dynamic = “#70AD47”, XGBoost = “#ED7D31”, MLP = “#A64D79”, Ensemble = “#2F5597” )) + scale_y_continuous(labels = label_number(scale = 1e-6, suffix = “M”)) + labs( title = “Six-week holdout: forecasts versus actual sales”, subtitle = “Black = actual; shaded band = validation-calibrated 95% ensemble interval”, x = NULL, y = “Daily sales”, color = NULL ) + theme_report ggsave(file.path(figure_dir, “figure_3_holdout_forecasts.png”), fig3, width = 8.2, height = 4.8, dpi = 300)

metric_plot_data <- melt( test_metrics, id.vars = c(“Split”, “Model”), measure.vars = c(“RMSE”, “MAE”, “MAPE”), variable.name = “Metric”, value.name = “Value” ) fig4 <- ggplot(metric_plot_data, aes(reorder(Model, Value), Value, fill = Model)) + geom_col(width = 0.68, show.legend = FALSE) + coord_flip() + facet_wrap(~ Metric, scales = “free_x”, ncol = 1) + scale_fill_brewer(palette = “Set2”) + scale_y_continuous(labels = function(x) { ifelse(abs(x) >= 100000, paste0(round(x / 1000000, 1), “M”), format(x, trim = TRUE, scientific = FALSE)) }) + labs( title = “Holdout forecast error by model”, subtitle = “Lower values indicate better performance”, x = NULL, y = NULL ) + theme_report + theme(strip.text = element_text(face = “bold”, color = “#17365D”)) ggsave(file.path(figure_dir, “figure_4_model_metrics.png”), fig4, width = 7.4, height = 7.0, dpi = 300)

residual_dt <- data.table( Index = seq_along(dynamic_residuals), Residual = dynamic_residuals ) acf_res <- acf(dynamic_residuals, plot = FALSE, lag.max = 35) acf_sq <- acf(dynamic_residuals^2, plot = FALSE, lag.max = 35) acf_dt <- rbindlist(list( data.table(Lag = as.numeric(acf_res$lag)[-1], ACF = as.numeric(acf_res$acf)[-1], Series = “Residuals”), data.table(Lag = as.numeric(acf_sq$lag)[-1], ACF = as.numeric(acf_sq$acf)[-1], Series = “Squared residuals”) )) p_resid <- ggplot(residual_dt, aes(Index, Residual)) + geom_hline(yintercept = 0, color = “#777777”, linewidth = 0.3) + geom_line(color = “#2F5597”, linewidth = 0.4) + labs(title = “Dynamic-model residuals”, x = “Training day”, y = “Log-scale residual”) + theme_report p_acf <- ggplot(acf_dt, aes(Lag, ACF, color = Series)) + geom_hline(yintercept = 0, color = “#777777”, linewidth = 0.3) + geom_segment(aes(xend = Lag, y = 0, yend = ACF), linewidth = 0.55) + scale_color_manual(values = c(“Residuals” = “#70AD47”, “Squared residuals” = “#C00000”)) + labs(title = “Residual dependence and volatility”, x = “Lag (days)”, y = “Autocorrelation”, color = NULL) + theme_report fig5 <- p_resid / p_acf + plot_annotation( title = “Residual diagnostics for the final dynamic regression”, theme = theme(plot.title = element_text(face = “bold”, color = “#17365D”, size = 14)) ) ggsave(file.path(figure_dir, “figure_5_residual_diagnostics.png”), fig5, width = 8.2, height = 7.1, dpi = 300)

top_importance <- head(feature_importance, 15) if (nrow(top_importance) > 0) { top_importance[, Feature := factor(Feature, levels = rev(Feature))] fig6 <- ggplot(top_importance, aes(Feature, Gain)) + geom_col(fill = “#4472C4”, width = 0.68) + coord_flip() + scale_y_continuous(labels = percent_format(accuracy = 1)) + labs( title = “XGBoost feature importance”, x = NULL, y = “Gain” ) + theme_report + theme(plot.title = element_text(face = “bold”, color = “#17365D”, size = 12)) ggsave(file.path(figure_dir, “figure_6_feature_importance.png”), fig6, width = 7.6, height = 5.4, dpi = 300) fig5_complete <- p_resid / (p_acf | fig6) + plot_annotation( title = “Residual diagnostics and XGBoost feature importance”, theme = theme(plot.title = element_text(face = “bold”, color = “#17365D”, size = 14)) ) ggsave( file.path(figure_dir, “figure_5_diagnostics_importance.png”), fig5_complete, width = 8.2, height = 8.2, dpi = 300 ) }

fwrite( daily[, setdiff(names(daily), “DescriptiveCustomers”), with = FALSE], file.path(output_dir, “daily_processed.csv”) ) fwrite(validation_forecasts, file.path(output_dir, “validation_forecasts.csv”)) fwrite(test_forecasts, file.path(output_dir, “test_forecasts.csv”)) fwrite(all_metrics, file.path(output_dir, “model_metrics.csv”)) fwrite(ensemble_weights, file.path(output_dir, “ensemble_weights.csv”)) fwrite(feature_importance, file.path(output_dir, “xgboost_feature_importance.csv”)) fwrite(diagnostics, file.path(output_dir, “residual_diagnostics.csv”)) fwrite(garch_coefficients, file.path(output_dir, “garch_coefficients.csv”)) fwrite(garch_variance_forecast, file.path(output_dir, “garch_variance_forecast.csv”))

data_summary <- data.table( Item = c( “Raw observations”, “Stores”, “Daily observations”, “First date”, “Last date”, “Training end”, “Validation start”, “Test start”, “Missing values”, “Customers used as predictor” ), Value = c( format(nrow(raw), big.mark = “,”), format(uniqueN(raw$Store), big.mark = ","), format(nrow(daily), big.mark = ","), as.character(min(daily$Date)), as.character(max(daily$Date)), as.character(train_end), as.character(validation_start), as.character(test_start), format(sum(missing_summary$Missing), big.mark = “,”), “No” ) ) fwrite(data_summary, file.path(output_dir, “data_summary.csv”))

best_test_model <- test_metrics[which.min(RMSE)] total_missing <- sum(missing_summary$Missing) open_sales_mean <- raw[Open == 1, mean(Sales)] promo_uplift <- raw[Open == 1, mean(Sales), by = Promo][order(Promo)] promo_uplift_pct <- if (nrow(promo_uplift) == 2) { (promo_uplift$V1[2] / promo_uplift$V1[1] - 1) * 100 } else NA_real_ best_vs_naive_rmse_pct <- ( 1 - best_test_model$RMSE / test_metrics[Model == “Seasonal naive”, RMSE] ) * 100 summary_lines <- c( “ROSSMANN SALES FORECASTING - EMPIRICAL SUMMARY”, paste(“Generated:”, format(Sys.time(), “%Y-%m-%d %H:%M:%S”)), paste(“Data:”, min(daily$Date), "through", max(daily$Date)), paste(“Raw rows:”, format(nrow(raw), big.mark = “,”)), paste(“Stores:”, uniqueN(raw$Store)), paste("Daily observations:", nrow(daily)), paste("Missing values:", total_missing), paste("Mean sales per open store-day:", round(open_sales_mean, 2)), paste("Descriptive promotion uplift among open stores (%):", round(promo_uplift_pct, 2)), paste("Training ends:", train_end), paste("Validation:", validation_start, "through", test_start - 1L), paste("Final holdout:", test_start, "through", max_date), "", "Ensemble weights:", paste(ensemble_weights$Model, sprintf(“%.4f”, ensemble_weights$Weight), sep = ": "), "", "Final holdout metrics:", apply(test_metrics, 1, function(x) paste(names(test_metrics), x, sep = "=", collapse = "; ")), "", paste("Best model by holdout RMSE:", best_test_model$Model), paste(“Best-model RMSE improvement over seasonal naive (%):”, round(best_vs_naive_rmse_pct, 2)), paste(“Ljung-Box p-value:”, signif(ljung_box$p.value, 5)), paste(“ARCH LM p-value:”, signif(arch_lm_p, 5)), paste(“GARCH status:”, garch_status) ) writeLines(summary_lines, file.path(output_dir, “analysis_summary.txt”))

saveRDS( list( dynamic = dynamic_final_model, xgboost = xgb_final_model, mlp = mlp_final_model, ensemble_weights = ensemble_weights, split_dates = list( train_end = train_end, validation_start = validation_start, test_start = test_start, test_end = max_date ) ), file.path(output_dir, “fitted_models.rds”) )

message(“Analysis complete. Outputs saved to:”, output_dir) print(test_metrics)