Time Series Forecasting EDA

import pandas as pd
import numpy as np
from scipy import stats

np.random.seed(42)   # reproducible results

# ── 1. Build the date index ──────────────────────────────────────────────────
weeks = pd.date_range(start="2022-01-03", periods=104, freq="W-MON")

# ── 2. Decompose into signal components ─────────────────────────────────────
t         = np.arange(104)
trend     = 140_000 + t * 600                      # upward trend: +$600/week
amplitude = 22_000                                 # seasonal swing
seasonality = amplitude * np.sin(2 * np.pi * t / 52)  # 52-week annual cycle
noise     = np.random.normal(0, 5_000, 104)        # week-to-week random noise
promo     = np.random.randint(7_000, 15_000, 104)  # promo spend

sales    = trend + seasonality + noise
units    = (sales / 100 + np.random.normal(0, 20, 104)).astype(int)

df = pd.DataFrame({
    "sales_revenue": sales.round(0),
    "units_sold":    units,
    "promo_spend":   promo.astype(float)
}, index=weeks)
df.index.name = "week"

# ── 3. Extract time features from DatetimeIndex ──────────────────────────────
df["month"]   = df.index.month
df["quarter"] = df.index.quarter.map({1:"Q1",2:"Q2",3:"Q3",4:"Q4"})
df["year"]    = df.index.year
df["week_num"]= df.index.isocalendar().week.astype(int)

print(f"Series: {len(df)} weeks | {df.index.min().date()} to {df.index.max().date()}")
print(f"\nSales revenue summary:")
print(df["sales_revenue"].describe().round(0).to_string())
print(f"\nFirst 3 rows:")
print(df[["sales_revenue","units_sold","promo_spend","month","quarter"]].head(3).to_string())

# Rolling means — different window lengths reveal different structure
df["roll4_mean"]  = df["sales_revenue"].rolling(window=4,  min_periods=1).mean().round(0)
df["roll13_mean"] = df["sales_revenue"].rolling(window=13, min_periods=1).mean().round(0)
df["roll52_mean"] = df["sales_revenue"].rolling(window=52, min_periods=1).mean().round(0)

# Rolling std — volatility over time
df["roll13_std"]  = df["sales_revenue"].rolling(window=13, min_periods=1).std().round(0)

# Trend strength: ratio of rolling-mean variance to raw variance
trend_var  = df["roll13_mean"].var()
total_var  = df["sales_revenue"].var()
trend_strength = trend_var / total_var

print("Rolling statistics (first 5 rows shown, full 104-row series computed):")
print(df[["sales_revenue","roll4_mean","roll13_mean","roll13_std"]].head(5).round(0).to_string())

print(f"\nTrend strength (roll13 variance / total variance): {trend_strength:.3f}")
print(f"Interpretation: {trend_strength*100:.1f}% of total variance is explained by the 13-week trend")

# Year-over-year growth
yr1_mean = df[df["year"]==2022]["sales_revenue"].mean()
yr2_mean = df[df["year"]==2023]["sales_revenue"].mean()
yoy_growth = (yr2_mean - yr1_mean) / yr1_mean * 100
print(f"\nYear-over-year comparison:")
print(f"  2022 mean weekly sales: ${yr1_mean:,.0f}")
print(f"  2023 mean weekly sales: ${yr2_mean:,.0f}")
print(f"  YoY growth:             {yoy_growth:.1f}%")

overall_mean = df["sales_revenue"].mean()

# Monthly seasonal profile
monthly = df.groupby("month")["sales_revenue"].agg(
    weeks         = "count",
    mean_sales    = "mean",
    std_sales     = "std",
    min_sales     = "min",
    max_sales     = "max"
).round(0).reset_index()
monthly["seasonal_index"] = (monthly["mean_sales"] / overall_mean).round(3)
monthly["vs_avg_pct"]     = ((monthly["seasonal_index"] - 1) * 100).round(1)
monthly["month_name"]     = pd.to_datetime(monthly["month"], format="%m").dt.strftime("%b")

print("Monthly seasonal profile:")
print(monthly[["month_name","weeks","mean_sales","seasonal_index","vs_avg_pct"]].to_string(index=False))

# Quarterly seasonal index
quarterly = df.groupby("quarter")["sales_revenue"].agg(
    weeks      = "count",
    mean_sales = "mean"
).round(0).reset_index()
quarterly["seasonal_index"] = (quarterly["mean_sales"] / overall_mean).round(3)
quarterly["vs_avg_pct"]     = ((quarterly["seasonal_index"] - 1) * 100).round(1)

# Reorder Q1→Q4
qorder = ["Q1","Q2","Q3","Q4"]
quarterly["quarter"] = pd.Categorical(quarterly["quarter"], categories=qorder, ordered=True)
quarterly = quarterly.sort_values("quarter")
print("\nQuarterly seasonal index:")
print(quarterly[["quarter","weeks","mean_sales","seasonal_index","vs_avg_pct"]].to_string(index=False))

# Peak vs trough ratio
peak_month   = monthly.loc[monthly["seasonal_index"].idxmax(), "month_name"]
trough_month = monthly.loc[monthly["seasonal_index"].idxmin(), "month_name"]
peak_idx     = monthly["seasonal_index"].max()
trough_idx   = monthly["seasonal_index"].min()
print(f"\nPeak month: {peak_month} (index {peak_idx:.3f})")
print(f"Trough month: {trough_month} (index {trough_idx:.3f})")
print(f"Peak-to-trough ratio: {peak_idx/trough_idx:.2f}x")

# Create lag features using .shift()
df["lag1"]  = df["sales_revenue"].shift(1)    # 1 week ago
df["lag4"]  = df["sales_revenue"].shift(4)    # 4 weeks ago (~1 month)
df["lag13"] = df["sales_revenue"].shift(13)   # 13 weeks ago (~1 quarter)
df["lag52"] = df["sales_revenue"].shift(52)   # 52 weeks ago (same week last year)

# Pearson correlations with current sales (dropna to align lengths)
lags = {"lag1": 1, "lag4": 4, "lag13": 13, "lag52": 52}
print("Autocorrelation — sales_revenue vs lag features:")
for col, lag in lags.items():
    valid = df[["sales_revenue", col]].dropna()
    r, p  = stats.pearsonr(valid["sales_revenue"], valid[col])
    sig   = "***" if p < 0.001 else "**" if p < 0.01 else "*" if p < 0.05 else "ns"
    print(f"  lag-{lag:>2}:  r = {r:+.3f}  p = {p:.4f}  {sig}")

# First-order differencing: removes trend, tests if differenced series is stationary
df["diff1"] = df["sales_revenue"].diff(1)    # week-to-week change
df["diff52"]= df["sales_revenue"].diff(52)   # year-over-year change

print(f"\nFirst-difference (week-over-week) stats:")
print(df["diff1"].dropna().describe().round(1).to_string())

print(f"\nSeasonal-difference (52-week) stats:")
print(df["diff52"].dropna().describe().round(1).to_string())

# Promotional lift: correlation of promo spend with sales
promo_r, promo_p = stats.pearsonr(df["sales_revenue"], df["promo_spend"])
print(f"\nPromo spend vs sales: r = {promo_r:+.3f}  p = {promo_p:.4f}")

from scipy.stats import shapiro

def adf_test_manual(series, label):
    """
    Manual ADF approximation using OLS regression on lagged differences.
    scipy-only implementation — no statsmodels required.
    Returns ADF-like statistic and stationarity conclusion.
    """
    s     = series.dropna().values
    dy    = np.diff(s)            # first differences
    y_lag = s[:-1]                # lagged level (the 'unit root' term)

    # OLS: regress dy on y_lag and a constant
    X      = np.column_stack([y_lag, np.ones(len(y_lag))])
    beta, res, _, _ = np.linalg.lstsq(X, dy, rcond=None)
    rho    = beta[0]              # coefficient on lagged level
    ss_res = np.sum((dy - X @ beta) ** 2)
    se_rho = np.sqrt(ss_res / (len(dy) - 2) / np.sum((y_lag - y_lag.mean())**2))
    t_stat = rho / se_rho         # ADF t-statistic

    # Critical value approximation (n=100, 5% level ≈ -2.89)
    stationary = t_stat < -2.89
    print(f"  {label:<30} t-stat: {t_stat:7.3f}  {'STATIONARY ✓' if stationary else 'NON-STATIONARY ✗'}")
    return t_stat, stationary

print("ADF Stationarity Tests (5% critical value ≈ -2.89):")
adf_test_manual(df["sales_revenue"],          "Raw sales")
adf_test_manual(np.log(df["sales_revenue"]), "Log(sales)")
adf_test_manual(df["diff1"].dropna(),        "First-differenced sales")
adf_test_manual(df["diff52"].dropna(),       "52-week-differenced sales")

# Normality of raw vs differenced series
_, p_raw  = shapiro(df["sales_revenue"].dropna())
_, p_diff = shapiro(df["diff1"].dropna())
print(f"\nShapiro-Wilk normality test:")
print(f"  Raw sales:          p = {p_raw:.4f}  {'Normal' if p_raw > 0.05 else 'Non-normal'}")
print(f"  First-differenced:  p = {p_diff:.4f}  {'Normal' if p_diff > 0.05 else 'Non-normal'}")

# Baseline 1: lag-52 (same week last year)
baseline_df = df[df["lag52"].notna()].copy()   # only 2023 rows have lag52
baseline_df["forecast_lag52"]   = baseline_df["lag52"]
baseline_df["resid_lag52"]      = baseline_df["sales_revenue"] - baseline_df["forecast_lag52"]
baseline_df["pct_error_lag52"]  = (baseline_df["resid_lag52"] / baseline_df["sales_revenue"] * 100)

# Baseline 2: seasonal mean — average sales for that week_num across all years
week_means = df.groupby("week_num")["sales_revenue"].mean()
baseline_df["forecast_seasonal"] = baseline_df["week_num"].map(week_means)
baseline_df["resid_seasonal"]    = baseline_df["sales_revenue"] - baseline_df["forecast_seasonal"]
baseline_df["pct_error_seasonal"]= (baseline_df["resid_seasonal"] / baseline_df["sales_revenue"] * 100)

# Accuracy metrics
def mae(resid):  return resid.abs().mean()
def rmse(resid): return np.sqrt((resid**2).mean())
def mape(pct):   return pct.abs().mean()

print("Forecast baseline accuracy (2023 hold-out — 52 weeks):")
print(f"  {'Metric':<8}  {'Lag-52':>12}  {'Seasonal Mean':>14}")
print(f"  {'MAE':<8}  ${mae(baseline_df['resid_lag52']):>10,.0f}  ${mae(baseline_df['resid_seasonal']):>12,.0f}")
print(f"  {'RMSE':<8}  ${rmse(baseline_df['resid_lag52']):>10,.0f}  ${rmse(baseline_df['resid_seasonal']):>12,.0f}")
print(f"  {'MAPE':<8}  {mape(baseline_df['pct_error_lag52']):>11.1f}%  {mape(baseline_df['pct_error_seasonal']):>13.1f}%")

# Where does lag-52 fail most? — worst 5 weeks by absolute error
worst5 = baseline_df.nlargest(5, "resid_lag52")[
    ["sales_revenue","forecast_lag52","resid_lag52","pct_error_lag52","month"]
].round(0)
print("\nWorst 5 weeks for lag-52 baseline (largest positive residuals):")
print(worst5.to_string())