Module 9: Time Series Lab - Build Your Forecasting Engine

⚡ URGENT MISSION BRIEF

Situation: MegaRetail's inventory system is hemorrhaging $75M annually. The board has given you 90 minutes to build and deploy a forecasting engine that will transform the company's supply chain.

Your Objective: Implement a complete time series forecasting system with ARIMA modeling, dynamic safety stock calculation, and multi-SKU optimization.

Success Criteria: Achieve <10% MAPE, >95% service level, and demonstrate $50M+ in annual savings.

Stakes: Success = Chief Innovation Officer promotion. Failure = Company faces hostile takeover.

📊 Task 1: Load and Prepare Time Series Data

Estimated: 15 min

First, we need to load three years of sales history and prepare it for analysis. This data contains hidden patterns worth millions.

data_preparation.py

import pandas as pd
import numpy as np
from datetime import datetime, timedelta

# Generate realistic sales data with patterns
def generate_sales_data():
    """
    Create 3 years of daily sales data with:
    - Linear growth trend (12% annually)
    - Monthly seasonality
    - Holiday spikes
    - Random noise
    """
    dates = pd.date_range(start='2022-01-01', end='2024-12-31', freq='D')
    n_days = len(dates)
    
    # TODO: Create trend component (12% annual growth)
    trend = # Your code here
    
    # TODO: Add monthly seasonality (peaks in November/December)
    seasonal = # Your code here
    
    # TODO: Add holiday spikes (Black Friday, Christmas)
    holidays = # Your code here
    
    # Combine components
    sales = trend + seasonal + holidays + np.random.normal(0, 50, n_days)
    
    # Create DataFrame
    df = pd.DataFrame({
        'date': dates,
        'sales': sales.clip(min=0),  # No negative sales
        'units': (sales / 50).clip(min=0).astype(int)
    })
    
    return df

# Load and inspect data
sales_df = generate_sales_data()
print(f"Data shape: {sales_df.shape}")
print(f"Date range: {sales_df['date'].min()} to {sales_df['date'].max()}")
print(f"Average daily sales: ${sales_df['sales'].mean():,.2f}")
print(f"Peak sales day: ${sales_df['sales'].max():,.2f}")
print(f"Total revenue: ${sales_df['sales'].sum():,.2f}")

Hint: For the trend, use: 100000 + np.arange(n_days) * 100. For seasonality, use sine waves with period 365.25. For holidays, create spikes on specific dates using conditionals.

🔍 Task 2: Decompose Time Series Components

Estimated: 10 min

Identify the hidden patterns in your sales data using STL decomposition. Each component reveals a different business insight.

decomposition.py

from statsmodels.tsa.seasonal import STL
import matplotlib.pyplot as plt

def decompose_time_series(df):
    """
    Perform STL decomposition to extract:
    - Trend: Long-term direction
    - Seasonal: Repeating patterns
    - Residual: Random variations
    """
    # Set date as index
    df_indexed = df.set_index('date')
    
    # TODO: Perform STL decomposition
    # Use period=365 for yearly seasonality
    stl = # Your code here
    result = # Fit the model
    
    # Extract components
    trend = result.trend
    seasonal = result.seasonal
    residual = result.resid
    
    # Calculate component strengths
    trend_strength = 1 - np.var(residual) / np.var(trend + residual)
    seasonal_strength = 1 - np.var(residual) / np.var(seasonal + residual)
    
    # Business insights
    insights = {
        'trend_direction': 'increasing' if trend.iloc[-1] > trend.iloc[0] else 'decreasing',
        'trend_strength': trend_strength,
        'seasonal_strength': seasonal_strength,
        'peak_season': df_indexed.groupby(df_indexed.index.month)['sales'].mean().idxmax(),
        'volatility': residual.std()
    }
    
    return result, insights

# Run decomposition
decomposition, insights = decompose_time_series(sales_df)

print("📈 Decomposition Insights:")
print(f"Trend: {insights['trend_direction']} (strength: {insights['trend_strength']:.2%})")
print(f"Seasonality strength: {insights['seasonal_strength']:.2%}")
print(f"Peak sales month: {insights['peak_season']}")
print(f"Daily volatility: ${insights['volatility']:,.2f}")

# Calculate financial impact
seasonal_impact = insights['seasonal_strength'] * sales_df['sales'].mean() * 365
print(f"\n💰 Financial Impact:")
print(f"Seasonal variation worth: ${seasonal_impact:,.0f} annually")

🤖 Task 3: Build ARIMA Forecasting Model

Estimated: 20 min

Implement an auto-ARIMA model that will predict future demand with precision. This is where the $50M savings begin.

arima_model.py

from statsmodels.tsa.arima.model import ARIMA
from sklearn.metrics import mean_absolute_percentage_error
import warnings
warnings.filterwarnings('ignore')

def build_arima_model(df, test_days=30):
    """
    Build and optimize ARIMA model
    Target: MAPE < 10% for board approval
    """
    # Split data
    train_size = len(df) - test_days
    train = df.iloc[:train_size]
    test = df.iloc[train_size:]
    
    # TODO: Grid search for best ARIMA parameters
    best_aic = np.inf
    best_params = None
    best_model = None
    
    for p in range(0, 3):
        for d in range(0, 2):
            for q in range(0, 3):
                try:
                    # TODO: Fit ARIMA model
                    model = # Your code here
                    fitted = # Fit the model
                    
                    if fitted.aic < best_aic:
                        best_aic = fitted.aic
                        best_params = (p, d, q)
                        best_model = fitted
                except:
                    continue
    
    # Generate forecast
    forecast = best_model.forecast(steps=test_days)
    
    # Calculate metrics
    mape = mean_absolute_percentage_error(test['sales'], forecast)
    rmse = np.sqrt(np.mean((test['sales'] - forecast) ** 2))
    
    # Calculate confidence intervals
    forecast_df = best_model.get_forecast(steps=test_days)
    confidence = forecast_df.conf_int()
    
    return {
        'model': best_model,
        'params': best_params,
        'mape': mape,
        'rmse': rmse,
        'forecast': forecast,
        'confidence': confidence,
        'test_actual': test['sales'].values
    }

# Build and evaluate model
results = build_arima_model(sales_df)

print(f"🎯 Model Performance:")
print(f"Best ARIMA parameters: {results['params']}")
print(f"MAPE: {results['mape']:.2%}")
print(f"RMSE: ${results['rmse']:,.2f}")

# Financial impact
avg_daily_sales = sales_df['sales'].mean()
accuracy_improvement = 0.15 - results['mape']  # vs 15% baseline
daily_savings = avg_daily_sales * 0.25 * accuracy_improvement  # 25% holding cost
annual_savings = daily_savings * 365

print(f"\n💰 Business Impact:")
print(f"Accuracy improvement: {accuracy_improvement:.2%}")
print(f"Daily inventory savings: ${daily_savings:,.2f}")
print(f"Annual savings: ${annual_savings:,.2f}")

# Check if we meet board requirements
if results['mape'] < 0.10:
    print("\n✅ SUCCESS: Model meets board requirements!")
else:
    print(f"\n⚠️ WARNING: MAPE {results['mape']:.2%} > 10% target")

Model Validation Tests:

MAPE < 10% PENDING

Residuals White Noise PENDING

No Autocorrelation PENDING

Forecast Within CI PENDING

⚡ URGENT CHALLENGE: Black Friday in 30 Days!

The CEO just called - Black Friday is approaching and last year we lost $10M in sales due to stockouts.

Modify your forecast to account for the 500% demand spike!

🛡️ Task 4: Calculate Dynamic Safety Stock

Estimated: 15 min

Replace static safety stock with dynamic calculations based on forecast uncertainty. This is where we save $35M annually.

safety_stock.py

🔬 Time Series Lab: Operation Forecast

⚡ URGENT MISSION BRIEF

Model Validation Tests:

⚡ URGENT CHALLENGE: Black Friday in 30 Days!