pacman::p_load(pacman, tidyverse, gridExtra, zoo, forecast)

# Function to read, process, and predict future sales for a given data file
process_and_predict <- function(file_path, currency_col) {
  # Read data, specifying column types
  data <- read_csv(file_path, col_types = cols(
    Year = col_character(),
    !!currency_col := col_double()
  ))
  
  # Process data: separate year and quarter, create date, filter columns
  data <- data %>%
    separate(Year, into = c("Year", "Quarter"), sep = " ") %>%
    mutate(
      Quarter = recode(Quarter, Q1 = "01", Q2 = "04", Q3 = "07", Q4 = "10"),
      Date = as.Date(paste(Year, Quarter, "01", sep = "-"))
    ) %>%
    select(Date, all_of(currency_col))
  
  # Create time series object
  ts_data <- ts(data[[currency_col]], start = c(2013, 1), frequency = 4)
  
  # Custom Plot Function
  customize_plot <- function(forecast_obj, title, actual_data = data) {
    forecast_years <- as.numeric(time(forecast_obj$mean))
    forecast_dates <- as.Date(zoo::as.yearqtr(forecast_years), frac = 1)
    
    plot_data <- data.frame(
      Date = forecast_dates,
      Forecast = as.vector(forecast_obj$mean),
      Upper = as.vector(forecast_obj$upper[, 2]),
      Lower = as.vector(forecast_obj$lower[, 2])
    )
    
    ggplot() +
      geom_line(data = actual_data, aes(x = Date, y = .data[[currency_col]]), color = "blue", linetype = "solid") +
      geom_line(data = plot_data, aes(x = Date, y = Forecast), color = "red", linetype = "dashed") +
      geom_ribbon(data = plot_data, aes(x = Date, ymin = Lower, ymax = Upper), fill = "grey", alpha = 0.2) +
      labs(title = title, y = currency_col, x = "Year") +
      theme_minimal() +
      scale_x_date(date_breaks = "1 year", date_labels = "%Y")
  }
  
  # Seasonal Decomposition of Time Series (STL)
  stl_model <- stl(ts_data, s.window = "periodic")
  stl_forecast <- forecast(stl_model, h = 40)
  stl_plot <- customize_plot(stl_forecast, paste(str_extract(file_path, "^[^\\.]+"), "STL Forecast"))
  
  # Exponential Smoothing with Seasonality (ETS)
  ets_model <- ets(ts_data)
  ets_forecast <- forecast(ets_model, h = 40)
  ets_plot <- customize_plot(ets_forecast, paste(str_extract(file_path, "^[^\\.]+"), "ETS Forecast"))
  
  # Autoregressive Integrated Moving Average (ARIMA) with Seasonality (SARIMA)
  sarima_model <- auto.arima(ts_data, seasonal = TRUE)
  sarima_forecast <- forecast(sarima_model, h = 40)
  sarima_plot <- customize_plot(sarima_forecast, paste(str_extract(file_path, "^[^\\.]+"), "SARIMA Forecast"))
  
  # Combine all plots for one company
  combined_plot <- grid.arrange(stl_plot, ets_plot, sarima_plot, ncol = 1)
  
  return(combined_plot)
}
# Generate plots for each company
amazon_plot  <- process_and_predict("Amazon.csv",  "Billion USD")

alibaba_plot <- process_and_predict("Alibaba.csv", "Billion CYD")

ebay_plot    <- process_and_predict("eBay.csv",    "Billion USD")