Time Series Analysis in R

library(forecast)
library(ggplot2)
library(dplyr)

# Create ts object
# ts(data, start, frequency)
# frequency: 12=monthly, 4=quarterly, 52=weekly, 365=daily

monthly_sales <- ts(
  data      = c(120, 135, 128, 142, 156, 168, 145, 162, 175, 188, 195, 220,
                 130, 148, 142, 158, 172, 184, 160, 178, 192, 205, 212, 245),
  start     = c(2023, 1),   # Year 2023, January
  frequency = 12            # Monthly
)

# Properties
start(monthly_sales)    # 2023 1
end(monthly_sales)      # 2024 12
frequency(monthly_sales) # 12
length(monthly_sales)   # 24

# Plot
autoplot(monthly_sales) +
  labs(title="Monthly Sales 2023-2024", x="Time", y="Units Sold") +
  theme_minimal()

# Seasonal plots
ggseasonplot(monthly_sales) +
  labs(title="Seasonal Pattern by Month")

ggsubseriesplot(monthly_sales) +
  labs(title="Monthly Subseries Plot")

# Decompose = Trend + Seasonal + Residual

# Classical decomposition
decomp <- decompose(monthly_sales, type="additive")  # or "multiplicative"
plot(decomp)
# Components:
decomp$trend     # Smooth trend
decomp$seasonal  # Seasonal pattern
decomp$random    # Residuals

# STL decomposition (better — handles changing seasonality)
stl_decomp <- stl(monthly_sales, s.window="periodic")
autoplot(stl_decomp) + theme_minimal()

# Moving averages (manual smoothing)
ma_3  <- ma(monthly_sales, order=3)   # 3-month moving average
ma_12 <- ma(monthly_sales, order=12)  # 12-month (removes seasonality)

autoplot(monthly_sales, series="Original", alpha=0.5) +
  autolayer(ma_3,  series="3-Month MA") +
  autolayer(ma_12, series="12-Month MA") +
  labs(title="Sales with Moving Averages") +
  theme_minimal()

# AUTO-ARIMA — automatically selects best model
model <- auto.arima(monthly_sales, seasonal=TRUE)
summary(model)

# Forecast 12 months ahead
forecast_result <- forecast(model, h=12)
print(forecast_result)

# Plot forecast with confidence intervals
autoplot(forecast_result) +
  labs(title="12-Month Sales Forecast",
       subtitle="80% and 95% prediction intervals",
       x="Time", y="Units Sold") +
  theme_minimal()

# Exponential Smoothing (ETS — handles trend + seasonality)
ets_model    <- ets(monthly_sales)
ets_forecast <- forecast(ets_model, h=12)
autoplot(ets_forecast) + theme_minimal()

# Model accuracy
accuracy(model)
# Shows: ME, RMSE, MAE, MAPE, MASE
cat(sprintf("MAPE: %.2f%% (Mean Absolute Percentage Error)\n",
             accuracy(model)[5]))

# Forecast specific values
cat("\n12-Month Forecast:\n")
for (i in 1:12) {
  cat(sprintf("  Month %d: %.0f (80%% CI: %.0f - %.0f)\n",
               i, forecast_result$mean[i],
               forecast_result$lower[i,1], forecast_result$upper[i,1]))
}