---
title: Supply Chain Data Analytics
subtitle: Analyzing and Forcasting Supermarket Sales
authors:
  - name: Stan Brouwer (2671939)
    orchid: 0009-0006-3447-0096
    affiliation: Vrije Universiteit
    corresponding: true
  - name: Liz Chen (2840693)
    affiliation: Master TSCM
  - name: Maaike Lamberts (2854979)
    affiliation: Supply Chain Data analysis
  - name: Niek Schroor (2786837)
    affiliation: Group 10

date: last-modified
bibliography: references.bib
number-sections: false

output:
  html: default
  pdf: default
prefer-html: true
---

::: column-margin

please visit 
https://sjbrou.github.io/Supply_Chain_Data_Analysis/
for an interactive version with better visualizations!

:::

## Data selection

We analyze, forecast and interpret the [Superstore sales](https://public.tableau.com/app/sample-data/sample_-_superstore.xls) provided by [Tableau](https://public.tableau.com/app/learn/sample-data) using different statistical and machine learning methods.

The dataset provided contains information about products, sales and profits of a fictitious US company. The dataset contains about 10,000 rows with 1,850 unique product names and 17 product subcategories, covering four consecutive years on a daily basis.

We describe our work in the PDF version. However, we would like to recommend reading our quarto manuscript *here* as it contains the **relevant** R code in the Article Notebook.

## Data Pre-processing

The superstore data set we selected is of high quality: At first glance (which needs to be verified during the visualization), the data appears to have been recorded regularly and without interruptions. There is no sign of a sudden structural change. Since the data are consumer products, it should contain both trends and seasonality.
Nevertheless, we have included hypothetical steps to demonstrate our understanding of the data preprocessing procedure. In detail, we did:

```{r message=FALSE, warning=FALSE, error=FALSE}
#| label: setup

# Clear workspace
rm(list = ls())
# Function to load (and install if necessary) dependencies
install_and_load <- function(packages) {
  install.packages(setdiff(packages, rownames(installed.packages())), dependencies = TRUE)
  invisible(lapply(packages, require, character.only = TRUE))
}
install_and_load(c("tidyverse", "readxl", "ggplot2", "lubridate", "stats", "Amelia","forecast", "tseries", "plotly", "stringr", "knitr", "kableExtra"))
```

-   Remove whitespaces from column names
-   Remove the <code>Row_ID</code> column as it can be inferred by it's index
-   Remove all columns with a single unique value, as storing these would be [redundant](https://few.vu.nl/~molenaar/courses/StatR/chapters/B-06-raw_data.html)
-   Ensure machine-readable date formats in <code>yyyy-mm-dd</code> as these usually differ per locale.
-   Ensure proper decimal separators
-   Calculate the number of missing values (both NA and empty string "") per column.

```{r}
#| label: data_preprocessing

# Load the data
suppressWarnings({data <- read_excel("data/sample_-_superstore.xls")}) # The Postal code column is stored as 'text' but coerced to numeric, causing warnings which we suppress

# Improve column names
colnames(data) <- str_replace_all(colnames(data), " ", "_")
colnames(data) <- str_replace_all(colnames(data), "-", "_")

# Remove the 'Row_ID' column as it can be inferred by it's index
data <- subset(data, select = -`Row_ID`)

# Remove all columns that have only one unique value, as storing these would be redundant
data <- data[, sapply(data, function(col) length(unique(col)) > 1)]

# Ensure a machine-readable date format as these are usually horrible in excel files
data$Order_Date <- as.Date(data$Order_Date, format = "%Y-%m-%d")
data$Ship_Date <- as.Date(data$Ship_Date, format = "%Y-%m-%d")

# The readxl package by default uses the correct decimal separator (as opposed to base R)

# Calculate the number of missing values per column.
# Origional dates and R date objects are in unix time, which return NA when compared to text (empty string). These dates are stored as 'double' datatype, Thus we check character columns for empty strings, and all columns for NA values. 
missing_values <- sapply(data, function(col) {
  if (inherits(col, "Date")) {
    sum(is.na(col))
  } else if (is.character(col)) {
    sum(is.na(col) | col == "")
  } else {
    sum(is.na(col))
  }
})

# sum(missing_values) returns 0!

# Optionally, print the missing values as a nice table
missing_values_table <- data.frame(
  Column = names(missing_values),
  Missing_or_Empty = missing_values
)
# Note that there are no missing values, thus we do not print them
# kable(missing_values_table, caption = "Missing or Empty Values in Columns", format = "pipe")


rm(missing_values, missing_values_table)
```

After these steps (and transposing the table for better document formatting), the data looks as follows:
```{r}
#| label: data_table1
kable(t(head(data, 3)), caption = "First 3 Rows of the Data (Transposed)", format = "markdown")
```

We did not find any missing values, confirming the quality of the data set. There is some more processing to do, for instance the removal of outliers. However, by doing so we impose our own assumptions on the data. Let's start by evaluating the descriptive statistics of our data and check if further processing is required.

```{r}
#| label: descriptive statistics

descriptive_statistics <- function(column) {
  if (is.numeric(column)) {
    stats <- list(
      Min = min(column, na.rm = TRUE), # Note that handling NA values increases robustness (and I copied the funciton from some of my earlier work)
      Max = max(column, na.rm = TRUE),
      Mean = mean(column, na.rm = TRUE),
      Median = median(column, na.rm = TRUE),
      StdDev = sd(column, na.rm = TRUE)
    )
  } else if (inherits(column, "Date")) {
    stats <- list(
      Earliest = format(min(column, na.rm = TRUE), "%Y-%m-%d"),
      Latest = format(max(column, na.rm = TRUE), "%Y-%m-%d")
    )
  } else if (is.character(column)) {
    stats <- list(
      Unique = length(unique(column)),
      Mode = names(sort(table(column), decreasing = TRUE)[1])
    )
  } else {
    stats <- NULL
  }
  return(stats)
}

# Call function on dataframe
descriptive_stats <- lapply(data, descriptive_statistics)

# Separate to tables dependent on data type
numeric_stats <- as.data.frame(do.call(rbind, lapply(names(data), function(col_name) {
  if (is.numeric(data[[col_name]])) {
    c(Column = col_name, descriptive_stats[[col_name]])
  }
})), stringsAsFactors = FALSE)
date_stats <- as.data.frame(do.call(rbind, lapply(names(data), function(col_name) {
  if (inherits(data[[col_name]], "Date")) {
    c(Column = col_name, descriptive_stats[[col_name]])
  }
})), stringsAsFactors = FALSE)
character_stats <- as.data.frame(do.call(rbind, lapply(names(data), function(col_name) {
  if (is.character(data[[col_name]])) {
    c(Column = col_name, descriptive_stats[[col_name]])
  }
})), stringsAsFactors = FALSE)
```

```{r}
#| Label: Descriptive_Output_tables
kable(
  numeric_stats,
  caption = "Descriptive Statistics for Numeric Columns",
  format = "pipe")

kable(
  date_stats,
  caption = "Descriptive Statistics for Date Columns",
  format = "pipe")
```
We inspect the orders with the lowest and highest Sales amount (in USD). The most expensive orders were professional printers, cameras and teleconferencing units with high unit prices. The orders with the lowest sales amount were often binders and had a high Discount rate.

Interestingly there are orders with a negative profit. They typically have high Discount rates and often concern the same item, such as the “Cubify CubeX 3D Printer Triple Head Print”. The orders with a negative Profit were often part of a larger order (for instance CA-2016-108196), and placed by customers with multiple orders. We suspect these negative Profit's to be caused by items of lower quality that receive discounts, general discount codes, or volume discounts. However, due to the high discounts especially on orders with negative profit, we assume these to be valid orders. 

** Some negative profit products **

In figure x we plotted the quantities of the most sold products. Unfortunately, the sold quantities of individual products were too low to determine any meaningful trends.

```{r message=FALSE, warnng=FALSE, error=FALSE}
#| label: Quantity_top_products
# Optionally: print top 10 sale quantity barplot
# # Sum of Quantity for top products
# top_products <- data %>%
#   group_by(Product_Name) %>%
#   summarize(total_quantity = sum(Quantity, na.rm = TRUE)) %>%
#   arrange(desc(total_quantity)) %>%
#   slice_head(n = 10) %>% 
#   mutate(ProdName8 = substr(Product_Name, 1, 8)) # Truncate product names to the first 8 characters. Long names mess up formatting
# 
# # Plot
# ggplot(top_products, aes(x = reorder(ProdName8, -total_quantity), y = total_quantity)) +
#   geom_bar(stat = "identity", fill = "steelblue") +
#   labs(title = "Top 20 Most Sold Products",
#        x = "Product ID",
#        y = "Total Quantity") +
#   theme_minimal() +
#   coord_flip()

# Aggregate quantity by Product Name and Order Date
time_series_data <- data %>%
  group_by(Product_Name, Order_Date) %>%
  summarize(total_quantity = sum(Quantity, na.rm = TRUE)) %>%
  ungroup()
# Filter by total quantity sold 
top_products <- time_series_data %>%
  group_by(Product_Name) %>%
  summarize(total_quantity = sum(total_quantity)) %>%
  arrange(desc(total_quantity)) %>%
  slice_head(n = 10)

# Filter the time-series data
filtered_time_series_data <- time_series_data %>%
  filter(Product_Name %in% top_products$Product_Name) %>%
  mutate(ProdName10 = substr(Product_Name, 1, 10)) # Product names can be quite long and mess up layouts. Lets only plot the first 10 chars.

# Here we do some special plotting. We want to show the plot with only one selected line by default, but make sure that the other 9 top sold products can be selected. We first create the ggplotly object, and than modify the visibility of the traces
```

```{r}
#| label: Top_Products_Quantity
#| fig-cap: Figure X Sale quantity of the most popular products

p_ly <- ggplotly(ggplot(filtered_time_series_data, aes(x = Order_Date, y = total_quantity, color = ProdName10)) +
  geom_line(size = 1) +
  labs(title = "Quantity Sold Over Time per Product",
       x = "Order Date",
       y = "Quantity Sold") +
  theme_minimal() +
  theme(legend.position = "bottom") +
  scale_color_discrete(name = "Product Name"))

for (i in seq_along(p_ly$x$data)) {
  if (i == 1) {
    p_ly$x$data[[i]]$visible <- TRUE
  } else {
    p_ly$x$data[[i]]$visible <- "legendonly"
  }
}

p_ly
```

Our proposed workaround is to aggregate <code>Product_Name</code> by <code>Sub_Category</code>, and treat it as a single product for the rest of the assignment, which we plotted in figure X.

```{r}
#| label: Aggregated_Sub_Category_sales
#| fig-caption: Aggregated Sub_Category sales (toggle )

top_categories <- data %>%
  count(Sub_Category, sort = TRUE)

top_10_categories <- top_categories$Sub_Category[0:10]

# Filter the data for  top 10 products
top_10_data <- data %>% filter(Sub_Category %in% top_10_categories)

# calculate sales per month
top_10_data <- top_10_data %>%
  mutate(Month = floor_date(Order_Date, unit = "month"))

# Aggregate data by month for each sub-category
top_10_data_aggregated <- top_10_data %>%
  group_by(Month, Sub_Category) %>%
  summarise(Sales_Count = n(), .groups = 'drop')

# Some special interactive plot formatting (see previous plot)
p_ly <- ggplotly(ggplot(top_10_data_aggregated, aes(x = Month, y = Sales_Count, color = Sub_Category, group = Sub_Category)) +
    geom_line(size = 1) +
    geom_point(size = 2) +
    labs(title = "Monthly Sales for the Top 3 Most Sold Sub Categories",
         x = "Month",
         y = "Sales Count",
         color = "Sub Category") +
    theme_minimal())

for (i in seq_along(p_ly$x$data)) {
  if (i == 1) {
    p_ly$x$data[[i]]$visible <- TRUE
  } else {
    p_ly$x$data[[i]]$visible <- "legendonly"
  }
}

p_ly
```

This aggregated Quantity starts to show trends and seasonality, and is much more useful to base predictions on! We will use these aggregated sub-categories for the rest of the assignment.

To properly finish our data preprocessing we ran some statistics on Quantity aggregated by Sub_Category. Table x contains some descriptive statistics.

```{r}
#| label: sub_category_descriptive_statistics_table
#| 
library(dplyr)
library(kableExtra)

# Summarize the data
outlier_summary <- data %>%
  group_by(Sub_Category) %>%
  summarize(
    Min = round(min(Quantity), 2),
    Mean = round(mean(Quantity), 2),
    Max = round(max(Quantity), 2),
    Sd = round(sd(Quantity), 2),
    CI_lower = round(Mean - 1.96 * (Sd / sqrt(n())), 2),
    CI_upper = round(Mean + 1.96 * (Sd / sqrt(n())), 2),
    .groups = "drop"
  )

# Output tables
kable(
  outlier_summary,
  caption = "Statistics for Sub_Category quantity",
  format = "pipe")
```

The statistics for <code>Quantity</code> aggregated by <code>Sub_Category</code> looks valid. We can visualize it as histogram and check for anomalies. Figure y contains histograms of <code>Quantity</code> per <code>Sub_Category</code>.

```{r}
#| label: sub_category_histograms

sub_categories <- unique(data$Sub_Category)

p <- plot_ly()
for (i in seq_along(sub_categories)) {
  sub <- sub_categories[i]
  subset_data <- data %>% filter(Sub_Category == sub)
  p <- add_trace(
    p,
    x = subset_data$Quantity,
    type = "histogram",
    name = sub,
    visible = ifelse(i == 1, TRUE, FALSE)
  )
}

# We add a drop down menu for Sub_Category as toggling visibility in default ggplot2 adds the histograms up. Instead we want to be able to show each histogram seperately. 
dropdown_buttons <- lapply(seq_along(sub_categories), function(i) {
  list(
    method = "update",
    args = list(
      list(visible = lapply(seq_along(sub_categories), function(j) j == i)),
      list(xaxis = list(title = "Quantity", autorange = TRUE), 
           yaxis = list(title = "Frequency", autorange = TRUE))
    ),
    label = sub_categories[i]
  )
})

# Style drop down layout
p <- p %>%
  layout(
    title = "Distribution of Quantity Sold per Order by Sub Category",
    xaxis = list(title = "Quantity"),
    yaxis = list(title = "Frequency"),
    showlegend = FALSE,  # Drop down instead of legend
    updatemenus = list(
      list(
        type = "dropdown",
        buttons = dropdown_buttons,
        direction = "down",
        x = 0.99,
        y = 0.99,
        showactive = TRUE,
        xanchor = "left",
        yanchor = "top"
      )
    )
  )
p
```

The histograms show that the quantities a right-skewed distributed. This is to be expected since most orders contain only a small number of items. We will not remove the outliers with large quantities since they appear valid..


## Forecasting Method Evaluation

### Forecasting top 3 product categories (4a)

Let's forecast sold quantities for the three most sold sub-categories:

The steps taken for data preparation were:

-   Identifying Top Subcategories: The top three subcategories are selected from our dataset based on their sold quantities. The top three were: Binders, furnishing and paper.
-   The sold quantities are aggregated monthly to create a time series object which we can use in the forecasting.
-   A KPSS showed that the data is non stationary. First-order differencing is applied to transform the data from non-stationary to stationary. The KPSS results in a p-value >0.05 showing the stationarity. 

```{r message=FALSE, error=FALSE, warning=FALSE}
#| label: plots1

top_categories <- data %>%
  group_by(Sub_Category) %>%
  summarise(Total_Quantity = sum(Quantity)) %>%
  arrange(desc(Total_Quantity))
top_3_subcategories <- top_categories$Sub_Category[0:3]

top_3_data <- data %>% filter(Sub_Category %in% top_3_subcategories)

# calculate sales per month
top_3_data <- top_3_data %>%
  mutate(Month = floor_date(Order_Date, unit = "month"))

# Aggregate data by month for each product
top_3_data_aggregated <- top_3_data %>%
  group_by(Month, Sub_Category) %>%
  summarise(Sales_Count = n(), .groups = 'drop')

# Create a time series object for each product
ts_data <- top_3_data_aggregated %>%
  pivot_wider(names_from = Sub_Category, values_from = Sales_Count, values_fill = 0) %>%
  select(-Month) %>%
  as.matrix()

ts_data <- ts(ts_data, start = c(2014, 1), end = c(2017, 12), frequency = 12) # @Stan This appears to be duplicate code?
ts_list <- list()
for (subcategory in top_3_subcategories) {
  # Filter data for the subcategory
  subcategory_data <- top_3_data_aggregated %>% filter(Sub_Category == subcategory)
  ts_list[[subcategory]] <- ts(subcategory_data$Sales_Count,
                                start = c(2014, 1),
                                end = c(2017, 12),
                                frequency = 12)
}


# forecasting methods on top 3 sub-categories
forecast_results <- list()

for (subcategory in names(ts_list)) {
  ts_current <- ts_list[[subcategory]]
  train_size <- floor(0.7 * length(ts_current))
  train_ts <- window(ts_current, end = c(2014 + (train_size - 1) %/% 12, (train_size - 1) %% 12 + 1))
  test_ts <- window(ts_current, start = c(2014 + train_size %/% 12, train_size %% 12 + 1))

  arima_model <- auto.arima(train_ts)
  arima_forecast <- forecast(arima_model, h = length(test_ts))
  arima_accuracy <- accuracy(arima_forecast, test_ts)

  hw_model <- HoltWinters(train_ts)
  hw_forecast <- forecast(hw_model, h = length(test_ts))
  hw_accuracy <- accuracy(hw_forecast, test_ts)

  ets_model <- ets(train_ts)
  ets_forecast <- forecast(ets_model, h = length(test_ts))
  ets_accuracy <- accuracy(ets_forecast, test_ts)

  forecast_results[[subcategory]] <- list(
    ARIMA = list(Model = arima_model, Forecast = arima_forecast, Accuracy = arima_accuracy),
    HoltWinters = list(Model = hw_model, Forecast = hw_forecast, Accuracy = hw_accuracy),
    ETS = list(Model = ets_model, Forecast = ets_forecast, Accuracy = ets_accuracy)
  )
}


# For formatting, we ommitted almost all output. You can uncomment the code and check your output if you like.

# # Step 5: Print results
#  for (subcategory in names(forecast_results)) {
#   cat("\n\nResults for Sub_Category:", subcategory, "\n")
# 
#   cat("\nARIMA Accuracy:\n")
#   print(forecast_results[[subcategory]]$ARIMA$Accuracy)
# 
#   cat("\nHolt-Winters Accuracy:\n")
#   print(forecast_results[[subcategory]]$HoltWinters$Accuracy)
# 
#   cat("\nETS Accuracy:\n")
#   print(forecast_results[[subcategory]]$ETS$Accuracy)
# }
```

Three models are applied to each subcategory to forecast it. The models we use are: ARIMA, Holt-Winters and ETS. We have chosen these models because of their level of suitability for discrete time series data with all different levels of trend and seasonality. To evaluate the methods and its effectiveness , the data is split into a training set (70%) and testing set (30%). 

```{r}
#| label: Results1

# Lets plot the results
for (subcategory in names(forecast_results)) {
  ts_current <- ts_list[[subcategory]]
  train_size <- floor(0.7 * length(ts_current))
  test_ts <- window(ts_current, start = c(2014 + train_size %/% 12, train_size %% 12 + 1))
  arima_forecast <- forecast_results[[subcategory]]$ARIMA$Forecast
  hw_forecast <- forecast_results[[subcategory]]$HoltWinters$Forecast
  ets_forecast <- forecast_results[[subcategory]]$ETS$Forecast
  
  # Combined plot
  plot(arima_forecast$mean, col = "blue", lwd = 2, 
       ylim = range(c(arima_forecast$mean, hw_forecast$mean, ets_forecast$mean, test_ts)),
       main = paste("Combined Forecasts for", subcategory, "before differencing"),
       xlab = "Time", ylab = "Forecasted Values")
  lines(test_ts, col = "red", lty = 2, lwd = 2)
  lines(hw_forecast$mean, col = "green", lwd = 2)
  lines(ets_forecast$mean, col = "purple", lwd = 2)
  legend("topleft", legend = c("ARIMA", "Holt-Winters", "ETS", "Test Data"),
         col = c("blue", "green", "purple", "red"), lty = c(1, 1, 1, 2), lwd = 2)
}

```
To assess the results, we use the following performance metrics: ME, RMSE, MAE and MAPE. They are calculated for the training and testing phases of the forecast. 


```{r message=FALSE, warning=FALSE, error=FALSE}

# Lets plot it as an table
# As the relevant data is stored in lists nested in lists nested in lists nested in lists, this was some work
results_table <- data.frame(
  Sub_Category = character(),
  Method = character(),
  Dataset = character(),
  ME = numeric(),
  RMSE = numeric(),
  MAE = numeric(),
  MPE = numeric(),
  MAPE = numeric(),
  MASE = numeric(),
  ACF1 = numeric(),
  Theil_U = numeric(),
  stringsAsFactors = FALSE
)

for (subcategory in names(forecast_results)) {
  subcategory_results <- forecast_results[[subcategory]]
  for (method in c("ARIMA", "HoltWinters", "ETS")) {
    if (!is.null(subcategory_results[[method]])) {
      accuracy_metrics <- subcategory_results[[method]]$Accuracy
      for (dataset in rownames(accuracy_metrics)) {
        row <- data.frame(
          Sub_Category = subcategory,
          Method = method,
          Dataset = dataset,
          ME = round(accuracy_metrics[dataset, "ME"], 3),
          RMSE = round(accuracy_metrics[dataset, "RMSE"], 3),
          MAE = round(accuracy_metrics[dataset, "MAE"], 3),
          MPE = round(accuracy_metrics[dataset, "MPE"], 3),
          MAPE = round(accuracy_metrics[dataset, "MAPE"], 3),
          MASE = round(accuracy_metrics[dataset, "MASE"], 3),
          ACF1 = round(accuracy_metrics[dataset, "ACF1"], 3),
          Theil_U = ifelse("Theil's U" %in% colnames(accuracy_metrics), round(accuracy_metrics[dataset, "Theil's U"], 3), NA)
        )
        results_table <- rbind(results_table, row)
      }
    }
  }
}
filtered_results <- results_table[results_table$Dataset == "Training set", ]
knitr::kable(filtered_results, format = "html", caption = "Training Set Forecast Accuracy")

```



As we can see on the forecasting results ARIMA performed well for binders. We can state this because of the lowest RMSE. 

For the subcategory furnishings we can see that the ETS forecasting method is the most stable across the training and testing phase. 

For the last subcategory and product paper the ETS model is again the most consistent, comparing the statistics for training and test set. The high variability in the test data leads to larger forecasting errors in all the 3 models.

Concerning the residual diagnostics, the checks show no real autocorrelation for ARIMA models. Which indicates a good fitting forecast. 

### Conclusion (4a)

The most effective model is not the same in all the subcategories. Each model was validated based on its ability to capture seasonality and trend. ARIMA performed better for Binders, while ETS performed better for Furnishings and Paper. 

### Clustering Subcategories and forecasting (4b)

Steps Taken:

- Trend strength, random variation and seasonal strength were calculated for the subcategory using the time series we have.


Clustering:

The clustering method used is the hierarchical clustering method. To group subcategories into three different clusters based on features which are normalized. The hierarchical clustering gave the following results:

- Cluster 1: Stronger seasonality

- Cluster 2: moderate trend and seasonality

- Cluster 3: Lower trend and seasonal strength


For each different cluster we also used the models introduced in 4A (ARIMA, Holt-Winters, ETS). These results where all aggregated at the level of each cluster so we can assess mean RMSE and MAPE for each different model.

```{r message=FALSE, warning=FALSE, error=FALSE}
filtered_results <- results_table[results_table$Dataset == "Test set", ]
knitr::kable(filtered_results, format = "html", caption = "Forecast Accuracy Results")

```



```{r message=FALSE, warning=FALSE, error=FALSE}
top_3_subcategories <- top_categories$Sub_Category[1:3]
kpss_results <- data.frame(
  Sub_Category = character(),
  KPSS_Statistic = numeric(),
  P_Value = numeric(),
  Null_Hypothesis = character(),
  stringsAsFactors = FALSE
)

# Perform KPSS Test for each subcategory
for (subcategory in top_3_subcategories) {
  if (subcategory %in% names(ts_list)) {
    ts_current <- ts_list[[subcategory]]
    test_result <- kpss.test(ts_current)
    
    kpss_results <- kpss_results %>%
      add_row(
        Sub_Category = subcategory,
        KPSS_Statistic = test_result$statistic,
        P_Value = test_result$p.value,
        Null_Hypothesis = ifelse(test_result$p.value < 0.05, 
                                 "Rejected (Non-Stationary)", 
                                 "Not Rejected (Stationary)")
      )
  } else {
    cat("\nSub-Category not found in ts_list:", subcategory, "\n")
  }
}

# Print nicely formattted table
kable(kpss_results, caption = "KPSS Test Results for Top 3 Subcategories", digits = 3)
```

```{r message=FALSE, warning=FALSE, error=FALSE}
#| label: Results2
options(warn = -1)
# # Step 6: Visualization of Forecasts
# for (subcategory in names(forecast_results)) {
#   plot(forecast_results[[subcategory]]$ARIMA$Forecast, main = paste("ARIMA Forecast for", subcategory))
#   lines(test_ts, col = "red", lty = 2)
# 
#   plot(forecast_results[[subcategory]]$HoltWinters$Forecast, main = paste("Holt-Winters Forecast for", subcategory))
#   lines(test_ts, col = "red", lty = 2)
# 
#   plot(forecast_results[[subcategory]]$ETS$Forecast, main = paste("ETS Forecast for", subcategory))
#   lines(test_ts, col = "red", lty = 2)
# }

#more stationary tests
# Perform KPSS Test for the top 3 subcategories
# top_3_subcategories <- top_categories$Sub_Category[0:3]
# 
# for (subcategory in top_3_subcategories) {
#   if (subcategory %in% names(ts_list)) {
#     ts_current <- ts_list[[subcategory]]
#     cat("\nKPSS Test for Sub-Category:", subcategory, "\n")
#     print(kpss.test(ts_current))
#   } else {
#     cat("\nSub-Category not found in ts_list:", subcategory, "\n")
#   }
# }

#Because the all the 3 subcategory are non stationary because of a P value which is <=0.05 we need to use differencing
# Apply differencing to each of the top 3 subcategories
differenced_series <- list()

for (subcategory in top_3_subcategories) {
  if (subcategory %in% names(ts_list)) {
    ts_current <- ts_list[[subcategory]]  # Get the time series for the subcategory
    ts_diff <- diff(ts_current, differences = 1)  # Apply first-order differencing
    differenced_series[[subcategory]] <- ts_diff  # Store the differenced series

    # Recheck stationarity with KPSS test
    #cat("\nKPSS Test for Differenced Sub-Category:", subcategory, "\n")
    print(kpss.test(ts_diff))
  } else {
    #cat("\nSub-Category not found in ts_list:", subcategory, "\n")
  }
}

```

```{r}
#| label: results3
# Time differenced plots
# # Now P value is larger then 0.05 so we have stationary data
# # Plot the differenced series for each subcategory
# for (subcategory in top_3_subcategories) {
#   if (subcategory %in% names(differenced_series)) {
#     ts_diff <- differenced_series[[subcategory]]
#     cat("\nPlotting Differenced Series for Sub-Category:", subcategory, "\n")
#     plot(ts_diff, main = paste("Differenced Series for Sub-Category:", subcategory),
#          ylab = "Differenced Values", xlab = "Time")
#   }
# }

# Combine differenced series plots for all top subcategories
combined_diff_plot <- function(differenced_series, top_3_subcategories) {
  plot(NULL, xlim = range(time(differenced_series[[top_3_subcategories[1]]])), 
       ylim = range(sapply(differenced_series[top_3_subcategories], range, na.rm = TRUE)),
       xlab = "Time", ylab = "Differenced Values",
       main = "Differenced Series for Top Subcategories")
  colors <- c("blue", "green", "purple")
  for (i in seq_along(top_3_subcategories)) {
    subcategory <- top_3_subcategories[i]
    if (subcategory %in% names(differenced_series)) {
      ts_diff <- differenced_series[[subcategory]]
      lines(ts_diff, col = colors[i], lwd = 2, lty = i)
    }
  }
  legend("topright", legend = top_3_subcategories, 
         col = colors, lty = 1:length(top_3_subcategories), lwd = 2)
}
combined_diff_plot(differenced_series, top_3_subcategories)
```
 
 
```{r}
#| label: Results4

#NEW FORECASTING FOR 4A with stationary data
# Step 4: Apply forecasting methods to the differenced top 3 sub-categories
forecast_results <- list()  # Store results

for (subcategory in names(differenced_series)) {
  ts_current <- differenced_series[[subcategory]]

  train_size <- floor(0.7 * length(ts_current))
  train_ts <- window(ts_current, end = c(2014 + (train_size - 1) %/% 12, (train_size - 1) %% 12 + 1))
  test_ts <- window(ts_current, start = c(2014 + train_size %/% 12, train_size %% 12 + 1))

  arima_model <- auto.arima(train_ts)
  arima_forecast <- forecast(arima_model, h = length(test_ts))
  arima_accuracy <- accuracy(arima_forecast, test_ts)

  hw_model <- HoltWinters(train_ts)
  hw_forecast <- forecast(hw_model, h = length(test_ts))
  hw_accuracy <- accuracy(hw_forecast, test_ts)

  ets_model <- ets(train_ts)
  ets_forecast <- forecast(ets_model, h = length(test_ts))
  ets_accuracy <- accuracy(ets_forecast, test_ts)

  forecast_results[[subcategory]] <- list(
    ARIMA = list(Model = arima_model, Forecast = arima_forecast, Accuracy = arima_accuracy),
    HoltWinters = list(Model = hw_model, Forecast = hw_forecast, Accuracy = hw_accuracy),
    ETS = list(Model = ets_model, Forecast = ets_forecast, Accuracy = ets_accuracy)
  )
}

# # Step 5: Print results
# for (subcategory in names(forecast_results)) {
#   cat("\n\nResults for Sub_Category:", subcategory, "\n")
# 
#   cat("\nARIMA Accuracy:\n")
#   print(forecast_results[[subcategory]]$ARIMA$Accuracy)
# 
#   cat("\nHolt-Winters Accuracy:\n")
#   print(forecast_results[[subcategory]]$HoltWinters$Accuracy)
# 
#   cat("\nETS Accuracy:\n")
#   print(forecast_results[[subcategory]]$ETS$Accuracy)
# }

# # Step 6: Visualization of Forecasts
# for (subcategory in names(forecast_results)) {
#   plot(forecast_results[[subcategory]]$ARIMA$Forecast,
#        main = paste("ARIMA Forecast for", subcategory),
#        ylab = "Differenced Values", xlab = "Time")
#   lines(test_ts, col = "red", lty = 2)
# 
#   plot(forecast_results[[subcategory]]$HoltWinters$Forecast,
#        main = paste("Holt-Winters Forecast for", subcategory),
#        ylab = "Differenced Values", xlab = "Time")
#   lines(test_ts, col = "red", lty = 2)
# 
#   plot(forecast_results[[subcategory]]$ETS$Forecast,
#        main = paste("ETS Forecast for", subcategory),
#        ylab = "Differenced Values", xlab = "Time")
#   lines(test_ts, col = "red", lty = 2)
# }

# Lets plot the results
for (subcategory in names(forecast_results)) {
  ts_current <- ts_list[[subcategory]]
  train_size <- floor(0.7 * length(ts_current))
  test_ts <- window(ts_current, start = c(2014 + train_size %/% 12, train_size %% 12 + 1))
  arima_forecast <- forecast_results[[subcategory]]$ARIMA$Forecast
  hw_forecast <- forecast_results[[subcategory]]$HoltWinters$Forecast
  ets_forecast <- forecast_results[[subcategory]]$ETS$Forecast
  
  # Combined plot
  plot(arima_forecast$mean, col = "blue", lwd = 2, 
       ylim = range(c(arima_forecast$mean, hw_forecast$mean, ets_forecast$mean, test_ts)),
       main = paste("Combined Forecasts for", subcategory, "before differencing"),
       xlab = "Time", ylab = "Forecasted Values")
  lines(test_ts, col = "red", lty = 2, lwd = 2)
  lines(hw_forecast$mean, col = "green", lwd = 2)
  lines(ets_forecast$mean, col = "purple", lwd = 2)
  legend("topleft", legend = c("ARIMA", "Holt-Winters", "ETS", "Test Data"),
         col = c("blue", "green", "purple", "red"), lty = c(1, 1, 1, 2), lwd = 2)
}

```



```{r}
#| label: KPSStest
# final KPSS test

# Perform KPSS Test for the differenced series in the top 3 subcategories
for (subcategory in top_3_subcategories) {
  if (subcategory %in% names(differenced_series)) {
    ts_current <- differenced_series[[subcategory]]  # Get the differenced series
    cat("\nKPSS Test for Differenced Sub-Category:", subcategory, "\n")
    print(kpss.test(ts_current))
  } else {
    cat("\nSub-Category not found in differenced_series:", subcategory, "\n")
  }
}
# now they are all 0.1
```



### Clustering (4b)

```{r}
#| label: somelable
# 4B
# 4B: Group Products into Clusters and Apply Forecasting Techniques
time_series_features <- data.frame()

for (subcategory in names(ts_list)) {
  ts_current <- ts_list[[subcategory]]

  # Decompose the time series to extract features
  decomposition <- decompose(ts_current)
  trend_strength <- var(decomposition$trend, na.rm = TRUE) / var(ts_current, na.rm = TRUE)
  seasonal_strength <- var(decomposition$seasonal, na.rm = TRUE) / var(ts_current, na.rm = TRUE)
  random_strength <- var(decomposition$random, na.rm = TRUE) / var(ts_current, na.rm = TRUE)

  # Store extracted features
  time_series_features <- rbind(time_series_features,
                                data.frame(SubCategory = subcategory,
                                           TrendStrength = trend_strength,
                                           SeasonalStrength = seasonal_strength,
                                           RandomStrength = random_strength))
}

# Step 2: Normalize the Features for Clustering
time_series_features_scaled <- time_series_features %>%
  select(-SubCategory) %>%
  scale()

# verify rows
# nrow(time_series_features_scaled)

# Step 3: Perform K-Means Clustering
# Determine the number of clusters dynamically
k <- min(3, nrow(time_series_features_scaled))  # Set k to the smaller of 3 or the number of rows
# Hierarchical Clustering
distance_matrix <- dist(time_series_features_scaled)  # Calculate distance matrix
hc <- hclust(distance_matrix)  # Perform hierarchical clustering
time_series_features$Cluster <- cutree(hc, k = k)  # Cut tree into 'k' clusters
# Add cluster information to the original data
time_series_features$Cluster <- cutree(hc, k = k)

# Step 4: Apply Forecasting Techniques to Each Cluster
forecast_results_by_cluster <- list()

for (cluster_id in unique(time_series_features$Cluster)) {
  # Get subcategories in the current cluster
  subcategories_in_cluster <- time_series_features$SubCategory[time_series_features$Cluster == cluster_id]

  # Initialize storage for cluster results
  cluster_forecast_results <- list()

  for (subcategory in subcategories_in_cluster) {
    if (subcategory %in% names(ts_list)) {
      ts_current <- ts_list[[subcategory]]  # Access the time series

      # Split the data into training and validation sets (70% training, 30% testing)
      train_size <- floor(0.7 * length(ts_current))
      train_ts <- window(ts_current, end = c(2014 + (train_size - 1) %/% 12, (train_size - 1) %% 12 + 1))
      test_ts <- window(ts_current, start = c(2014 + train_size %/% 12, train_size %% 12 + 1))

      # 1. ARIMA
      arima_model <- auto.arima(train_ts)
      arima_forecast <- forecast(arima_model, h = length(test_ts))
      arima_accuracy <- accuracy(arima_forecast, test_ts)

      # 2. Holt-Winters
      hw_model <- HoltWinters(train_ts)
      hw_forecast <- forecast(hw_model, h = length(test_ts))
      hw_accuracy <- accuracy(hw_forecast, test_ts)

      # 3. ETS
      ets_model <- ets(train_ts)
      ets_forecast <- forecast(ets_model, h = length(test_ts))
      ets_accuracy <- accuracy(ets_forecast, test_ts)

      # Store results for the subcategory
      cluster_forecast_results[[subcategory]] <- list(
        ARIMA = list(Model = arima_model, Forecast = arima_forecast, Accuracy = arima_accuracy),
        HoltWinters = list(Model = hw_model, Forecast = hw_forecast, Accuracy = hw_accuracy),
        ETS = list(Model = ets_model, Forecast = ets_forecast, Accuracy = ets_accuracy)
      )
    } else {
      cat("\nSub-Category not found in ts_list:", subcategory, "\n")
    }
  }

  # Store results for the cluster
  forecast_results_by_cluster[[paste0("Cluster_", cluster_id)]] <- cluster_forecast_results
}

# Step 5: Print Forecasting Accuracy for Each Cluster
for (cluster_id in names(forecast_results_by_cluster)) {
  cat("\n\nResults for", cluster_id, "\n")
  cluster_results <- forecast_results_by_cluster[[cluster_id]]

  for (subcategory in names(cluster_results)) {
    cat("\nSub-Category:", subcategory, "\n")

    cat("\nARIMA Accuracy:\n")
    print(cluster_results[[subcategory]]$ARIMA$Accuracy)

    cat("\nHolt-Winters Accuracy:\n")
    print(cluster_results[[subcategory]]$HoltWinters$Accuracy)

    cat("\nETS Accuracy:\n")
    print(cluster_results[[subcategory]]$ETS$Accuracy)
  }
}

## Step 6: Visualize the Clusters
library(ggplot2)

ggplot(time_series_features, aes(x = TrendStrength, y = SeasonalStrength, color = as.factor(Cluster))) +
  geom_point(size = 3) +
  labs(title = "Clusters of Subcategories Based on Time-Series Features",
       x = "Trend Strength", y = "Seasonal Strength", color = "Cluster") +
  theme_minimal()
#
#check
#residual diagnostic
for (cluster_id in names(forecast_results_by_cluster)) {
  cluster_results <- forecast_results_by_cluster[[cluster_id]]
  for (subcategory in names(cluster_results)) {
    cat("\nResidual Diagnostics for Sub-Category:", subcategory, "\n")
    checkresiduals(cluster_results[[subcategory]]$ARIMA$Model)
  }
}
# P-value is higher then 0.1 so we have stationary data, this is good
#cluster level metrics
cluster_metrics <- data.frame()
for (cluster_id in names(forecast_results_by_cluster)) {
  cluster_results <- forecast_results_by_cluster[[cluster_id]]
  cluster_rmse <- sapply(cluster_results, function(x) x$ARIMA$Accuracy["Test set", "RMSE"])
  cluster_mape <- sapply(cluster_results, function(x) x$ARIMA$Accuracy["Test set", "MAPE"])
  cluster_metrics <- rbind(cluster_metrics, data.frame(Cluster = cluster_id, MeanRMSE = mean(cluster_rmse), MeanMAPE = mean(cluster_mape)))
}
print(cluster_metrics)

```
Cluster 1 (e.g., Binders): ARIMA outperformed other methods due to significant autocorrelation and trend components.

Cluster 2 (e.g., Furnishings): ETS was the most accurate method, effectively balancing trend and seasonality.

Cluster 3 (e.g., Paper): ETS also performed best, with ARIMA showing higher error rates due to variability in random components.

Residual diagnostics were performed for all ARIMA models, confirming no significant autocorrelation (p > 0.05).

Cluster-Level Metrics based on mean RMSE and MAPE show:
-   Cluster 1 had the lowest RMSE using ARIMA.
-   Cluster 2 and 3 were better modeled with ETS

### Conclusion (4b)

Clustering allows for tailored forecasting strategies. We conclude that for the given data set ARIMA is more effective for clusters with strong trends, while ETS is preferable for clusters with mixed seasonal and trend characteristics. The approach aligns with lecture notes, emphasizing the importance of adapting models based on time series characteristics.

## Forecasting future values

### Forecasting 3 products (5a)

In this session, we focused on evaluating different forecasting models (ARIMA, Holt-Winters, and ETS) for multiple sub-categories by analyzing their accuracy metrics, such as RMSE, MAPE, and residual diagnostics. Based on the evaluation results, we selected the best-performing model for each sub-category. We then used these models to forecast the future outcomes for each sub-category, projecting the data for the next year. 

```{r}
#Binders->choose ARIMA
binders_ts <- ts_list[["Binders"]]
arima_model <- auto.arima(binders_ts)
summary(arima_model)
arima_forecast <- forecast(arima_model, h = 12)
print(arima_forecast)
plot(arima_forecast, main = "ARIMA Forecast for Binders (Next 12 Months)")

#Paper->choose ETS
paper_ts <- ts_list[["Paper"]]
ets_model <- ets(paper_ts)
summary(ets_model)
ets_forecast <- forecast(ets_model, h = 12)
print(ets_forecast)
plot(ets_forecast, main = "ETS Forecast for Paper (Next 12 Months)")

#Furnishings->choose ETS
furnishings_ts <- ts_list[["Furnishings"]]
ets_model <- ets(furnishings_ts)
summary(ets_model)
ets_forecast <- forecast(ets_model, h = 12)
print(ets_forecast)
plot(ets_forecast, main = "ETS Forecast for Furnishings (Next 12 Months)")

```

### Applying to all data (5b)

In this session, we first grouped the sub-categories into clusters based on key time-series features, including trend strength, seasonal strength, and random strength, using hierarchical clustering. Once the clusters were formed, we applied and evaluated multiple forecasting models—ARIMA, Holt-Winters, and ETS—on each sub-category within its respective cluster, comparing their accuracy metrics such as RMSE and MAPE. Based on the evaluation results, we selected the best-performing model for each sub-category and used it to forecast the future outcomes within a year, leveraging the clustering to enhance the accuracy and relevance of our predictions.

```{r}
#Cluster_Binders->Holt-Winters
cluster_id <- 1
subcategory <- "Binders"
hw_model <- forecast_results_by_cluster[[paste0("Cluster_", cluster_id)]][[subcategory]]$HoltWinters$Model
hw_forecast <- forecast(hw_model, h = 12)
print(hw_forecast)
plot(hw_forecast, main = "Holt-Winters Forecast for Binders (Next 12 Months)", xlab = "Time", ylab = "Forecasted Values")

#Cluster_paper->ETS
cluster_id <- 2
subcategory <- "Paper"
ets_model <- forecast_results_by_cluster[[paste0("Cluster_", cluster_id)]][[subcategory]]$ETS$Model
ets_forecast <- forecast(ets_model, h = 12)
print(ets_forecast)
plot(ets_forecast, main = "ETS Forecast for Paper (Next 12 Months)", xlab = "Time", ylab = "Forecasted Values")

#Cluster_Furnishings->Holt-Winters
cluster_id <- 3
subcategory <- "Furnishings"
hw_model <- forecast_results_by_cluster[[paste0("Cluster_", cluster_id)]][[subcategory]]$HoltWinters$Model
hw_forecast <- forecast(hw_model, h = 12)
print(hw_forecast)
plot(hw_forecast, main = "Holt-Winters Forecast for Furnishings (Next 12 Months)", xlab = "Time", ylab = "Forecasted Values")

```
## Forecast interpretation





```{r eval=FALSE}

# There was some missing pieces of code that I did not want to remove, but added no value either, 


# Check for missing values
missing_values <- colSums(is.na(data))
print(missing_values)  # Print missing values for reference
# heat map
library(Amelia)
missmap(data, main = "Missing Data Pattern")
#distribution of key variables
#plot Quantity
ggplot(data, aes(x = Quantity)) +
  geom_histogram(binwidth = 1, fill = "steelblue") +
  labs(title = "Distribution of Quantity", x = "Quantity", y = "Frequency") +
  theme_minimal()
#plot sales
ggplot(data, aes(x = Sales)) +
  geom_histogram(binwidth = 50, fill = "tomato") +
  labs(title = "Distribution of Sales", x = "Sales", y = "Frequency") +
  theme_minimal()
# plot profit
ggplot(data, aes(x = Profit)) +
  geom_histogram(binwidth = 10, fill = "purple") +
  labs(title = "Distribution of Profit", x = "Profit", y = "Frequency") +
  theme_minimal()
# time based trends
data$Order_Date <- as.Date(data$Order_Date, format = "%Y-%m-%d")  # Ensure date format
time_series <- data %>%
  group_by(Order_Date) %>%
  summarize(total_sales = sum(Sales), total_profit = sum(Profit), total_quantity = sum(Quantity))

ggplot(time_series, aes(x = Order_Date)) +
  geom_line(aes(y = total_sales, color = "Sales")) +
  geom_line(aes(y = total_profit, color = "Profit")) +
  geom_line(aes(y = total_quantity, color = "Quantity")) +
  labs(title = "Sales, Profit, and Quantity Over Time", x = "Date", y = "Value") +
  theme_minimal() +
  scale_color_manual(name = "Metrics", values = c("Sales" = "blue", "Profit" = "green", "Quantity" = "red"))

#sales by category and sub category
category_sales <- data %>%
  group_by(Category, Sub_Category) %>%
  summarize(total_sales = sum(Sales))

ggplot(category_sales, aes(x = reorder(Sub_Category, -total_sales), y = total_sales, fill = Category)) +
  geom_bar(stat = "identity") +
  labs(title = "Sales by Category and Sub-Category", x = "Sub-Category", y = "Total Sales") +
  theme_minimal() +
  coord_flip()

#Outliers detection
#Quantity
ggplot(data, aes(x = Category, y = Quantity)) +
  geom_boxplot() +
  labs(title = "Boxplot of Quantity by Category", x = "Category", y = "Quantity")
#sales
ggplot(data, aes(x = Category, y = Sales)) +
  geom_boxplot() +
  labs(title = "Boxplot of Sales by Category", x = "Category", y = "Sales")

#profit
ggplot(data, aes(x = Category, y = Profit)) +
  geom_boxplot() +
  labs(title = "Boxplot of Profit by Category", x = "Category", y = "Profit")
#Geo visualization

us_map <- map_data("state")
if("State" %in% colnames(data)) {
  state_sales <- data %>%
    group_by(State) %>%
    summarize(total_sales = sum(Sales))

  # Convert state names to lowercase to match map data
  state_sales$State <- tolower(state_sales$State)

  # Merge state sales data with map data
  state_sales_map <- merge(us_map, state_sales, by.x = "region", by.y = "State", all.x = TRUE)

  # Plot sales by state
  ggplot(state_sales_map, aes(long, lat, group = group, fill = total_sales)) +
    geom_polygon(color = "white") +
    scale_fill_continuous(low = "lightblue", high = "darkblue", na.value = "gray90") +
    labs(title = "Sales by State", fill = "Total Sales") +
    theme_void() +
    coord_fixed(1.3)
}

# correlation matrix
numerical_data <- data %>% select(where(is.numeric))

cor_matrix <- cor(numerical_data, use = "complete.obs")

# Convert the correlation matrix to a long format
cor_data <- as.data.frame(as.table(cor_matrix))

# Plot the correlation matrix using ggplot2
ggplot(cor_data, aes(Var1, Var2, fill = Freq)) +
  geom_tile(color = "white") +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white",
                       midpoint = 0, limit = c(-1, 1), space = "Lab",
                       name="Correlation") +
  geom_text(aes(label = round(Freq, 2)), color = "black", size = 4) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1,
                                   size = 12, hjust = 1)) +
  coord_fixed() +
  labs(title = "Correlation Matrix of Key Variables", x = "", y = "")
```


```{r eval=FALSE}

#Aggregate sales per month
monthly_sales <- data %>%
  mutate(Month = floor_date(Order_Date, "month")) %>%
  group_by(Month) %>%
  summarize(total_sales = sum(Sales))
#Convert to time series
sales_ts <- ts(monthly_sales$total_sales, frequency = 12, start = c(year(min(monthly_sales$Month)), month(min(monthly_sales$Month))))
#Arima model
arima_model <- auto.arima(sales_ts)
arima_forecast <- forecast(arima_model, h = 12)
autoplot(arima_forecast) + labs(title = "ARIMA Forecast for Monthly Sales")
#Holts winter model
hw_model <- HoltWinters(sales_ts)
hw_forecast <- forecast(hw_model, h = 12)
autoplot(hw_forecast) + labs(title = "Holt-Winters Forecast for Monthly Sales")
# clustering for segmentation
library(cluster)
#data clustering
clustering_data <- data %>%
  select(Sales, Quantity, Discount, Profit) %>%
  na.omit()
set.seed(123)
kmeans_model <- kmeans(clustering_data, centers = 3)
data$Cluster <- as.factor(kmeans_model$cluster)
# visualize clustering result
ggplot(data, aes(x = Sales, y = Profit, color = Cluster)) +
  geom_point(alpha = 0.6) +
  labs(title = "K-Means Clustering of Sales and Profit", x = "Sales", y = "Profit") +
  theme_minimal()

```