Time Series Forecasting

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime, timedelta
import warnings
warnings.filterwarnings('ignore')
 
# Set seed for reproducible results
np.random.seed(42)
 
# Generate 2 years of daily sales data
date_range = pd.date_range(start='2022-01-01', end='2023-12-31', freq='D')
n_days = len(date_range)
 
# Base sales (trend + noise)
base_sales = 10000 + np.linspace(0, 2000, n_days)  # Upward trend
 
# Weekly seasonality (weekend sales increase)
weekly_pattern = np.array([1.0, 0.95, 0.9, 0.95, 1.1, 1.3, 1.2])
weekly_seasonality = np.tile(weekly_pattern, n_days // 7 + 1)[:n_days]
 
# Monthly seasonality (end of month sales increase)
monthly_seasonality = 1 + 0.1 * np.sin(2 * np.pi * np.arange(n_days) / 30)
 
# Yearly seasonality (Nov-Dec holiday season)
yearly_seasonality = np.ones(n_days)
for i, date in enumerate(date_range):
    if date.month in [11, 12]:
        yearly_seasonality[i] = 1.3
    elif date.month in [1, 2]:
        yearly_seasonality[i] = 0.8
 
# Final sales
daily_sales = pd.DataFrame({
    'ds': date_range,
    'y': (base_sales * weekly_seasonality * monthly_seasonality *
          yearly_seasonality * (1 + np.random.normal(0, 0.05, n_days)))
})
 
print(f"Data period: {daily_sales['ds'].min().date()} ~ {daily_sales['ds'].max().date()}")
print(f"Total days: {len(daily_sales)}")
print(f"Average daily sales: ${daily_sales['y'].mean():,.0f}")
print(daily_sales.head())

# Full time series visualization
fig, axes = plt.subplots(3, 1, figsize=(14, 10))
 
# Overall sales trend
axes[0].plot(daily_sales['ds'], daily_sales['y'], linewidth=0.8)
axes[0].set_title('Daily Sales Trend', fontsize=14, fontweight='bold')
axes[0].set_xlabel('Date')
axes[0].set_ylabel('Sales ($)')
axes[0].grid(True, alpha=0.3)
 
# Monthly aggregation
monthly_sales = daily_sales.set_index('ds').resample('M').sum()
axes[1].bar(monthly_sales.index, monthly_sales['y'], width=20, color='steelblue')
axes[1].set_title('Monthly Sales', fontsize=14, fontweight='bold')
axes[1].set_xlabel('Month')
axes[1].set_ylabel('Sales ($)')
axes[1].grid(True, alpha=0.3)
 
# Day of week average
daily_sales['dayofweek'] = daily_sales['ds'].dt.day_name()
dow_order = ['Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday', 'Saturday', 'Sunday']
dow_sales = daily_sales.groupby('dayofweek')['y'].mean().reindex(dow_order)
axes[2].bar(dow_sales.index, dow_sales.values, color='coral')
axes[2].set_title('Average Sales by Day of Week', fontsize=14, fontweight='bold')
axes[2].set_xlabel('Day of Week')
axes[2].set_ylabel('Average Sales ($)')
axes[2].tick_params(axis='x', rotation=45)
 
plt.tight_layout()
plt.show()

from statsmodels.tsa.seasonal import seasonal_decompose
 
# Time series decomposition (multiplicative model)
ts_data = daily_sales.set_index('ds')['y']
decomposition = seasonal_decompose(ts_data, model='multiplicative', period=7)
 
fig, axes = plt.subplots(4, 1, figsize=(14, 12))
 
# Original data
axes[0].plot(decomposition.observed)
axes[0].set_title('Observed', fontsize=12)
axes[0].set_ylabel('Sales')
 
# Trend
axes[1].plot(decomposition.trend, color='orange')
axes[1].set_title('Trend', fontsize=12)
axes[1].set_ylabel('Sales')
 
# Seasonality
axes[2].plot(decomposition.seasonal, color='green')
axes[2].set_title('Seasonal', fontsize=12)
axes[2].set_ylabel('Seasonal Index')
 
# Residual
axes[3].plot(decomposition.resid, color='red')
axes[3].set_title('Residual', fontsize=12)
axes[3].set_ylabel('Residual')
 
plt.tight_layout()
plt.show()
 
print("=== Time Series Decomposition Summary ===")
print(f"Trend range: ${decomposition.trend.min():,.0f} ~ ${decomposition.trend.max():,.0f}")
print(f"Seasonality range: {decomposition.seasonal.min():.2f} ~ {decomposition.seasonal.max():.2f}")

# Calculate moving averages
daily_sales['ma_7'] = daily_sales['y'].rolling(window=7).mean()
daily_sales['ma_30'] = daily_sales['y'].rolling(window=30).mean()
 
# Visualization
plt.figure(figsize=(14, 6))
plt.plot(daily_sales['ds'][-180:], daily_sales['y'][-180:],
         alpha=0.5, label='Actual Sales')
plt.plot(daily_sales['ds'][-180:], daily_sales['ma_7'][-180:],
         linewidth=2, label='7-day Moving Average')
plt.plot(daily_sales['ds'][-180:], daily_sales['ma_30'][-180:],
         linewidth=2, label='30-day Moving Average')
plt.title('Moving Averages (Last 180 days)', fontsize=14, fontweight='bold')
plt.xlabel('Date')
plt.ylabel('Sales ($)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

from statsmodels.tsa.holtwinters import SimpleExpSmoothing, ExponentialSmoothing
 
# Train/test split
train_size = int(len(daily_sales) * 0.8)
train = daily_sales['y'][:train_size]
test = daily_sales['y'][train_size:]
 
print(f"Training data: {train_size} days")
print(f"Test data: {len(test)} days")
 
# Simple exponential smoothing
ses_model = SimpleExpSmoothing(train).fit(smoothing_level=0.2, optimized=False)
ses_forecast = ses_model.forecast(len(test))
 
# Evaluation
from sklearn.metrics import mean_absolute_error, mean_squared_error
 
mae_ses = mean_absolute_error(test, ses_forecast)
rmse_ses = np.sqrt(mean_squared_error(test, ses_forecast))
print(f"\n=== Simple Exponential Smoothing Results ===")
print(f"MAE: ${mae_ses:,.2f}")
print(f"RMSE: ${rmse_ses:,.2f}")

# Holt-Winters model (trend + seasonality)
hw_model = ExponentialSmoothing(
    train,
    trend='add',           # Additive trend
    seasonal='mul',        # Multiplicative seasonality
    seasonal_periods=7     # Weekly seasonality
).fit()
 
hw_forecast = hw_model.forecast(len(test))
 
# Evaluation
mae_hw = mean_absolute_error(test, hw_forecast)
rmse_hw = np.sqrt(mean_squared_error(test, hw_forecast))
print("=== Holt-Winters Results ===")
print(f"MAE: ${mae_hw:,.2f}")
print(f"RMSE: ${rmse_hw:,.2f}")
 
# Forecast visualization
plt.figure(figsize=(14, 6))
plt.plot(daily_sales['ds'][:train_size], train, label='Training Data', alpha=0.7)
plt.plot(daily_sales['ds'][train_size:], test, label='Actual Values', alpha=0.7)
plt.plot(daily_sales['ds'][train_size:], hw_forecast,
         label='Holt-Winters Forecast', linewidth=2, color='red')
plt.axvline(x=daily_sales['ds'].iloc[train_size], color='gray',
            linestyle='--', label='Forecast Start')
plt.title('Holt-Winters Forecast Results', fontsize=14, fontweight='bold')
plt.xlabel('Date')
plt.ylabel('Sales ($)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.stattools import adfuller
 
# Stationarity test (ADF test)
adf_result = adfuller(daily_sales['y'])
print("=== ADF Test (Stationarity Test) ===")
print(f"ADF Statistic: {adf_result[0]:.4f}")
print(f"p-value: {adf_result[1]:.4f}")
print(f"Stationarity: {'Stationary' if adf_result[1] < 0.05 else 'Non-stationary (differencing needed)'}")
 
# ACF, PACF visualization
fig, axes = plt.subplots(1, 2, figsize=(14, 4))
plot_acf(daily_sales['y'], lags=30, ax=axes[0])
axes[0].set_title('Autocorrelation Function (ACF)')
plot_pacf(daily_sales['y'], lags=30, ax=axes[1])
axes[1].set_title('Partial Autocorrelation Function (PACF)')
plt.tight_layout()
plt.show()

from statsmodels.tsa.arima.model import ARIMA
 
# ARIMA model (p=1, d=1, q=1)
arima_model = ARIMA(train, order=(1, 1, 1)).fit()
 
print("=== ARIMA Model Summary ===")
print(f"AIC: {arima_model.aic:.2f}")
print(f"BIC: {arima_model.bic:.2f}")
 
# Forecast
arima_forecast = arima_model.forecast(len(test))
 
# Evaluation
mae_arima = mean_absolute_error(test, arima_forecast)
rmse_arima = np.sqrt(mean_squared_error(test, arima_forecast))
print(f"\nMAE: ${mae_arima:,.2f}")
print(f"RMSE: ${rmse_arima:,.2f}")

from statsmodels.tsa.statespace.sarimax import SARIMAX
 
# SARIMA model
sarima_model = SARIMAX(
    train,
    order=(1, 1, 1),
    seasonal_order=(1, 1, 1, 7)  # Weekly seasonality
).fit(disp=False)
 
sarima_forecast = sarima_model.forecast(len(test))
 
# Evaluation
mae_sarima = mean_absolute_error(test, sarima_forecast)
rmse_sarima = np.sqrt(mean_squared_error(test, sarima_forecast))
print("=== SARIMA Results ===")
print(f"MAE: ${mae_sarima:,.2f}")
print(f"RMSE: ${rmse_sarima:,.2f}")

# Model performance comparison
models_comparison = pd.DataFrame({
    'Model': ['Simple Exp. Smoothing', 'Holt-Winters', 'ARIMA', 'SARIMA'],
    'MAE': [mae_ses, mae_hw, mae_arima, mae_sarima],
    'RMSE': [rmse_ses, rmse_hw, rmse_arima, rmse_sarima]
}).round(2)
 
models_comparison['MAPE'] = (models_comparison['MAE'] / test.mean() * 100).round(2)
 
print("=== Time Series Model Performance Comparison ===")
print(models_comparison.to_string(index=False))
 
# Visualization
fig, ax = plt.subplots(figsize=(10, 6))
x = np.arange(len(models_comparison))
width = 0.35
 
bars1 = ax.bar(x - width/2, models_comparison['MAE'], width, label='MAE', color='steelblue')
bars2 = ax.bar(x + width/2, models_comparison['RMSE'], width, label='RMSE', color='coral')
 
ax.set_xlabel('Model')
ax.set_ylabel('Error ($)')
ax.set_title('Time Series Model Performance Comparison', fontsize=14, fontweight='bold')
ax.set_xticks(x)
ax.set_xticklabels(models_comparison['Model'])
ax.legend()
ax.grid(True, alpha=0.3)
 
plt.tight_layout()
plt.show()

# Train final model on full data
final_model = ExponentialSmoothing(
    daily_sales['y'],
    trend='add',
    seasonal='mul',
    seasonal_periods=7
).fit()
 
# 90-day forecast
forecast_days = 90
future_forecast = final_model.forecast(forecast_days)
future_dates = pd.date_range(
    start=daily_sales['ds'].max() + timedelta(days=1),
    periods=forecast_days
)
 
# Confidence interval (simple estimation)
forecast_std = daily_sales['y'].std() * 0.1
upper_bound = future_forecast + 1.96 * forecast_std
lower_bound = future_forecast - 1.96 * forecast_std
 
# Visualization
plt.figure(figsize=(14, 6))
plt.plot(daily_sales['ds'][-90:], daily_sales['y'][-90:],
         label='Actual Data', alpha=0.7)
plt.plot(future_dates, future_forecast,
         label='Forecast', color='red', linewidth=2)
plt.fill_between(future_dates, lower_bound, upper_bound,
                 alpha=0.3, color='red', label='95% Confidence Interval')
plt.axvline(x=daily_sales['ds'].max(), color='gray',
            linestyle='--', label='Forecast Start')
plt.title('90-Day Sales Forecast', fontsize=14, fontweight='bold')
plt.xlabel('Date')
plt.ylabel('Sales ($)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Monthly forecast totals
forecast_df = pd.DataFrame({
    'ds': future_dates,
    'yhat': future_forecast,
    'yhat_lower': lower_bound,
    'yhat_upper': upper_bound
})
 
forecast_df['month'] = forecast_df['ds'].dt.to_period('M')
monthly_forecast = forecast_df.groupby('month').agg(
    Forecast_Sales=('yhat', 'sum'),
    Lower=('yhat_lower', 'sum'),
    Upper=('yhat_upper', 'sum')
).round(0)
 
print("=== Monthly Sales Forecast ===")
print(monthly_forecast)
 
# Day of week average forecast
forecast_df['dayofweek'] = forecast_df['ds'].dt.day_name()
daily_pattern = forecast_df.groupby('dayofweek')['yhat'].mean().reindex(dow_order)
print("\n=== Average Forecast Sales by Day of Week ===")
for day, value in daily_pattern.items():
    print(f"{day}: ${value:,.0f}")

# Define major holidays
holidays = pd.DataFrame({
    'holiday': ['New Year', 'Christmas', 'Black Friday', 'Thanksgiving'],
    'ds': pd.to_datetime(['2023-01-01', '2023-12-25', '2023-11-24', '2023-11-23'])
})
 
# Analyze sales around holidays
def analyze_holiday_effect(data, holiday_date, window=7):
    mask = (data['ds'] >= holiday_date - timedelta(days=window)) & \
           (data['ds'] <= holiday_date + timedelta(days=window))
    holiday_sales = data.loc[mask, 'y'].mean()
    normal_sales = data['y'].mean()
    effect = (holiday_sales - normal_sales) / normal_sales * 100
    return effect
 
print("=== Holiday Effect Analysis ===")
for _, row in holidays.iterrows():
    effect = analyze_holiday_effect(daily_sales, row['ds'])
    print(f"{row['holiday']}: {effect:+.1f}%")