Week 11: Time Series Analysis and Forecasting

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from datetime import datetime, timedelta

print("="*50)
print("TIME SERIES DATA FUNDAMENTALS")
print("="*50)

print("""
WHAT IS TIME SERIES DATA?
- Measurements taken at regular intervals over time
- Data points are DEPENDENT (not independent)
- Order matters: can't shuffle randomly
- Sequential patterns: patterns repeat, trends emerge

EXAMPLES:
Daily:
- Stock closing prices
- Temperature readings
- Website page views

Weekly:
- Sales figures
- Customer counts

Monthly:
- Unemployment rate
- Retail sales

Hourly:
- Server load
- Heart rate monitoring
- Power consumption
""")

# Create sample time series
print("\nEXAMPLE: SALES DATA OVER TIME")
print("-"*50)

# Create date range
dates = pd.date_range(start='2023-01-01', periods=365, freq='D')

# Create synthetic sales data with trend and seasonality
np.random.seed(42)
trend = np.linspace(100, 150, 365) # Increasing trend
seasonality = 20 * np.sin(np.arange(365) * 2 * np.pi / 365) # Yearly cycle
noise = np.random.normal(0, 5, 365) # Random noise
sales = trend + seasonality + noise

# Create DataFrame
ts_data = pd.DataFrame({
'Date': dates,
'Sales': sales
})

print(f"\nDataset: {len(ts_data)} daily sales observations")
print(f"Date range: {ts_data['Date'].min()} to {ts_data['Date'].max()}")
print(f"Sales range: ${ts_data['Sales'].min():.2f} to ${ts_data['Sales'].max():.2f}")
print(f"\nFirst few rows:")
print(ts_data.head(10))

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Plot 1: Full time series
axes[0, 0].plot(ts_data['Date'], ts_data['Sales'], linewidth=1.5)
axes[0, 0].set_xlabel('Date')
axes[0, 0].set_ylabel('Sales ($)')
axes[0, 0].set_title('Sales Over Time (Full Year)')
axes[0, 0].grid(True, alpha=0.3)

# Plot 2: Components (trend, seasonality, noise)
axes[0, 1].plot(dates, trend, label='Trend', linewidth=2)
axes[0, 1].plot(dates, seasonality, label='Seasonality', linewidth=1.5, alpha=0.7)
axes[0, 1].plot(dates, noise, label='Noise', linewidth=1, alpha=0.5)
axes[0, 1].set_xlabel('Date')
axes[0, 1].set_ylabel('Component Value')
axes[0, 1].set_title('Time Series Components')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# Plot 3: First 90 days (zoom in)
axes[1, 0].plot(ts_data['Date'][:90], ts_data['Sales'][:90], marker='o', linewidth=2, markersize=4)
axes[1, 0].set_xlabel('Date')
axes[1, 0].set_ylabel('Sales ($)')
axes[1, 0].set_title('First 90 Days (Zoomed In)')
axes[1, 0].grid(True, alpha=0.3)

# Plot 4: Month comparison
ts_data['Month'] = ts_data['Date'].dt.month
ts_data.boxplot(column='Sales', by='Month', ax=axes[1, 1])
axes[1, 1].set_xlabel('Month')
axes[1, 1].set_ylabel('Sales ($)')
axes[1, 1].set_title('Sales Distribution by Month')
plt.sca(axes[1, 1])
plt.xticks(range(1, 13), ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun',
'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'])

plt.tight_layout()
plt.show()

# Time series components
print("\n" + "="*50)
print("TIME SERIES COMPONENTS")
print("="*50)

print("""
EVERY TIME SERIES CAN BE DECOMPOSED INTO:

1. TREND (Long-term direction)
- Is data going up, down, or staying flat?
- Often due to fundamental changes
- Example: Sales growing due to business expansion
- Can be linear or non-linear
- Duration: Months to years

2. SEASONALITY (Regular repeating patterns)
- Patterns that repeat at fixed intervals
- Can be: Daily, Weekly, Monthly, Yearly
- Due to systematic factors
- Examples:
* Retail: Higher sales before holidays
* Temperature: Hot in summer, cold in winter
* Traffic: More cars during rush hours
- Seasonal period: Time for pattern to repeat

3. CYCLICAL (Long-term oscillations)
- Rises and falls that are NOT fixed intervals
- Longer than seasonality
- Often tied to economic cycles
- Difficult to forecast
- Example: Business cycle, stock market cycles

4. NOISE/IRREGULAR (Random variation)
- Unexplained random fluctuations
- Due to: Measurement error, unexpected events
- No predictable pattern
- Cannot be forecast

TIME SERIES DECOMPOSITION FORMULA:
Y_t = Trend_t + Seasonal_t + Cyclic_t + Noise_t
or (multiplicative):
Y_t = Trend_t × Seasonal_t × Noise_t

ADDITIVE: Use when seasonal variation constant
MULTIPLICATIVE: Use when seasonal variation increases with level

IMPORTANCE FOR FORECASTING:
- Understand each component
- Forecast each separately
- Combine forecasts
- Better results than treating as one series
""")

# Autocorrelation
print("\n" + "="*50)
print("AUTOCORRELATION - KEY CONCEPT FOR TIME SERIES")
print("="*50)

print("""
AUTOCORRELATION (ACF):
- Correlation of time series with itself at different lags
- Lag: How many time steps back
- ACF(1): Correlation with 1 step ago
- ACF(2): Correlation with 2 steps ago
- Range: -1 to 1 (like correlation)

WHAT IT TELLS US:
- How dependent today is on yesterday
- Whether data has memory
- Helps identify patterns

PATTERNS:
- High ACF at many lags: Strong trend (data not stationary)
- ACF decays quickly: Stationary data (good for forecasting)
- Spikes at specific lags: Seasonal pattern
- Random low values: White noise (can't predict)

PARTIAL AUTOCORRELATION (PACF):
- Direct correlation at each lag
- Removes intermediate lag effects
- Used for ARIMA model identification

IMPORTANCE:
- First check when analyzing time series
- Tells if data is stationary (needed for ARIMA)
- Helps choose forecasting method
""")

# Calculate ACF manually
from pandas.plotting import autocorrelation_plot

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# ACF plot
from statsmodels.graphics.tsaplots import plot_acf
plot_acf(ts_data['Sales'], lags=40, ax=axes[0])
axes[0].set_title('Autocorrelation Function (ACF)')
axes[0].set_xlabel('Lag (days)')

# Manual lagged plot
for lag in [1, 7, 30]:
lagged_sales = ts_data['Sales'].shift(lag)
correlation = ts_data['Sales'].corr(lagged_sales)
print(f"Correlation at lag {lag}: {correlation:.3f}")

# Seasonal subseries plot
ts_data['DayOfWeek'] = ts_data['Date'].dt.dayofweek
ts_data['Week'] = ts_data['Date'].dt.isocalendar().week

axes[1].scatter(ts_data['DayOfWeek'], ts_data['Sales'], alpha=0.5)
axes[1].set_xlabel('Day of Week (0=Monday)')
axes[1].set_ylabel('Sales ($)')
axes[1].set_title('Seasonal Pattern by Day of Week')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

import numpy as np
import pandas as pd
from statsmodels.tsa.stattools import adfuller, kpss
import matplotlib.pyplot as plt

print("="*50)
print("STATIONARITY - CRITICAL FOR FORECASTING")
print("="*50)

print("""
WHAT IS STATIONARY DATA?
Stationary time series has:
✓ Constant mean over time
✓ Constant variance over time
✓ No trend (flat over long term)
✓ No systematic patterns

Non-stationary data has:
✗ Changing mean (upward/downward trend)
✗ Changing variance (becoming more/less volatile)
✗ Clear trend or drift
✗ Seasonal patterns

WHY STATIONARITY MATTERS:
- Most forecasting methods (ARIMA, etc.) REQUIRE stationary data
- Non-stationary → unreliable forecasts
- Violations: Forecasts will be poor, confidence intervals wrong
- Solution: Transform data to make stationary

EXAMPLES:
Stationary: Random walk around fixed value
Non-stationary: Stock price (has trend)

TESTING FOR STATIONARITY:
1. Visual inspection: Plot and look for trend/changes
2. Statistical tests:
- Augmented Dickey-Fuller (ADF) test
- KPSS test
- Phillips-Perron test
""")

# Create stationary vs non-stationary data
np.random.seed(42)
n = 200

# Non-stationary (random walk with drift)
non_stationary = np.cumsum(np.random.normal(0.5, 1, n))

# Stationary (mean-reverting)
stationary = np.random.normal(0, 1, n)

# Seasonal non-stationary
t = np.arange(n)
seasonal_trend = 0.1 * t + 10 * np.sin(2 * np.pi * t / 30) + np.random.normal(0, 1, n)

# Visualization
fig, axes = plt.subplots(3, 2, figsize=(14, 12))

# Row 1: Non-stationary (random walk)
axes[0, 0].plot(non_stationary, linewidth=1.5)
axes[0, 0].set_title('Non-Stationary: Random Walk with Drift')
axes[0, 0].set_ylabel('Value')
axes[0, 0].grid(True, alpha=0.3)

# Differences
diff_non_stat = np.diff(non_stationary)
axes[0, 1].plot(diff_non_stat, linewidth=1.5)
axes[0, 1].set_title('Differenced (Random Walk) → Stationary')
axes[0, 1].set_ylabel('Change')
axes[0, 1].grid(True, alpha=0.3)

# Row 2: Stationary
axes[1, 0].plot(stationary, linewidth=1.5)
axes[1, 0].set_title('Stationary: White Noise')
axes[1, 0].set_ylabel('Value')
axes[1, 0].grid(True, alpha=0.3)

# Differences (shouldn't need)
diff_stat = np.diff(stationary)
axes[1, 1].plot(diff_stat, linewidth=1.5)
axes[1, 1].set_title('Differenced (Already Stationary)')
axes[1, 1].set_ylabel('Change')
axes[1, 1].grid(True, alpha=0.3)

# Row 3: Seasonal + trend
axes[2, 0].plot(seasonal_trend, linewidth=1.5)
axes[2, 0].set_title('Non-Stationary: Trend + Seasonality')
axes[2, 0].set_ylabel('Value')
axes[2, 0].grid(True, alpha=0.3)

diff_seasonal = np.diff(seasonal_trend)
axes[2, 1].plot(diff_seasonal, linewidth=1.5)
axes[2, 1].set_title('Differenced: Removed Trend')
axes[2, 1].set_ylabel('Change')
axes[2, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# ADF test
print("\n" + "="*50)
print("AUGMENTED DICKEY-FULLER (ADF) TEST")
print("="*50)

print("""
HYPOTHESIS TEST FOR STATIONARITY:
- H₀: Data is non-stationary (has unit root)
- H₁: Data is stationary
- If p < 0.05: Reject H₀ → Likely stationary
- If p ≥ 0.05: Fail to reject → Likely non-stationary

INTERPRETATION:
- p < 0.05: Stationary (good for ARIMA)
- p ≥ 0.05: Non-stationary (need to difference)

WHEN TO USE:
- First test on raw time series
- If non-stationary: Difference once
- Test again until stationary
- Record number of differences needed (this is 'd' in ARIMA)
""")

print("\nTesting example series:")

# Test non-stationary
result_ns = adfuller(non_stationary)
print(f"\nNon-stationary series (random walk):")
print(f" ADF Statistic: {result_ns[0]:.4f}")
print(f" P-value: {result_ns[1]:.4f}")
print(f" Verdict: {'Stationary' if result_ns[1] < 0.05 else 'NON-STATIONARY'}")

# Test differences
result_diff_ns = adfuller(diff_non_stat)
print(f"\nAfter differencing:")
print(f" ADF Statistic: {result_diff_ns[0]:.4f}")
print(f" P-value: {result_diff_ns[1]:.4f}")
print(f" Verdict: {'Stationary' if result_diff_ns[1] < 0.05 else 'NON-STATIONARY'}")

# Test stationary
result_s = adfuller(stationary)
print(f"\nStationary series (white noise):")
print(f" ADF Statistic: {result_s[0]:.4f}")
print(f" P-value: {result_s[1]:.4f}")
print(f" Verdict: {'Stationary' if result_s[1] < 0.05 else 'NON-STATIONARY'}")

# Differencing guide
print("\n" + "="*50)
print("DIFFERENCING STRATEGY")
print("="*50)

print("""
DIFFERENCING:
- Take difference between consecutive values
- Y_diff = Y_t - Y_(t-1)
- Removes trend from data
- Usually makes non-stationary data stationary

HOW TO KNOW HOW MANY DIFFERENCES NEEDED:
1. Plot data: Is there trend?
- If yes → Need differencing
2. Apply 1st difference, test with ADF
- If stationary (p < 0.05) → Done
- If not → Apply 2nd difference
3. Usually need 0-2 differences
- 0: Already stationary
- 1: Has trend
- 2: Has trend in the trend (rare)

WARNING: Over-differencing!
- Too many differences → Creates false patterns
- Lose information
- Rule: Use least differences needed to achieve stationarity

SEASONAL DIFFERENCING:
- For seasonal data: Difference at seasonal lag
- Example: Monthly data with yearly seasonality
* Regular difference: Y_t - Y_(t-1)
* Seasonal difference: Y_t - Y_(t-12)
- Can combine both types
- Usually need both for seasonal ARIMA
""")

print("\nDifferencing example:")
print(f"Original mean: {non_stationary.mean():.2f}")
print(f"Original variance: {non_stationary.var():.2f}")
print(f"\nAfter 1st difference:")
print(f"Differenced mean: {diff_non_stat.mean():.2f}")
print(f"Differenced variance: {diff_non_stat.var():.2f}")
print("→ Mean now around 0 (drift removed)")
print("→ Likely stationary")

import numpy as np
import pandas as pd
from statsmodels.tsa.holtwinters import ExponentialSmoothing, SimpleExpSmoothing
from sklearn.metrics import mean_squared_error, mean_absolute_error
import matplotlib.pyplot as plt

print("="*50)
print("EXPONENTIAL SMOOTHING METHODS")
print("="*50)

print("""
EXPONENTIAL SMOOTHING:
- Family of forecasting methods
- Give more weight to recent observations
- Recent data trusted more than old data
- Intuitive and easy to implement
- Good for data without strong trend/seasonality

TYPES:

1. SIMPLE EXPONENTIAL SMOOTHING (SES)
- For stationary data (no trend, no seasonality)
- Works: Weighted average of past observations
- Weights decay exponentially into past
- Parameter α (0-1): smoothing parameter
* α=0: Ignore new data
* α=0.5: Equal weight to recent and past
* α=1: Only use latest value (no smoothing)
- Good when: Noise is main issue
- Example: Stationary sales with noise

2. HOLT'S LINEAR TREND METHOD
- For data with trend (but no seasonality)
- Tracks trend separately from level
- Parameters: α (level), β (trend)
- Forecasts into future: Includes trend
- Good when: Clear upward or downward trend
- Example: Growing sales

3. HOLT-WINTERS SEASONAL METHOD
- For data with both trend AND seasonality
- Tracks: Level, Trend, Seasonal components
- Two variants:
* Additive: For constant seasonal variation
* Multiplicative: For increasing seasonal variation
- Parameters: α (level), β (trend), γ (seasonal)
- Best overall for real-world data
- Example: Retail sales (trend + holiday seasonality)
""")

# Create synthetic data with different characteristics
np.random.seed(42)
n = 200

# 1. Stationary data
stationary_data = 50 + np.random.normal(0, 3, n)

# 2. Data with trend
trend_data = 50 + 0.3 * np.arange(n) + np.random.normal(0, 3, n)

# 3. Data with trend and seasonality
t = np.arange(n)
seasonal_data = (50 + 0.3 * t +
10 * np.sin(2 * np.pi * t / 30) +
np.random.normal(0, 2, n))

# Fit models
print("\nFITTING MODELS TO DIFFERENT DATA TYPES:")
print("-"*50)

# Split into train/test
train_size = int(0.8 * n)
train_stat, test_stat = stationary_data[:train_size], stationary_data[train_size:]
train_trend, test_trend = trend_data[:train_size], trend_data[train_size:]
train_seas, test_seas = seasonal_data[:train_size], seasonal_data[train_size:]

# 1. Simple Exponential Smoothing
print("\n1. Simple Exponential Smoothing (Stationary data):")
model_ses = SimpleExpSmoothing(train_stat).fit(smoothing_level=0.2)
pred_ses = model_ses.forecast(steps=len(test_stat))
rmse_ses = np.sqrt(mean_squared_error(test_stat, pred_ses))
print(f" RMSE: {rmse_ses:.2f}")
print(f" Smoothing parameter (α): {model_ses.params['smoothing_level']:.3f}")

# 2. Holt's linear trend
print("\n2. Holt's Linear Trend (Trend data):")
from statsmodels.tsa.holtwinters import Holt
model_holt = Holt(train_trend).fit()
pred_holt = model_holt.forecast(steps=len(test_trend))
rmse_holt = np.sqrt(mean_squared_error(test_trend, pred_holt))
print(f" RMSE: {rmse_holt:.2f}")

# 3. Holt-Winters
print("\n3. Holt-Winters (Trend + Seasonality):")
model_hw = ExponentialSmoothing(train_seas, seasonal_periods=30,
trend='add', seasonal='add').fit()
pred_hw = model_hw.forecast(steps=len(test_seas))
rmse_hw = np.sqrt(mean_squared_error(test_seas, pred_hw))
print(f" RMSE: {rmse_hw:.2f}")

# Visualization
fig, axes = plt.subplots(3, 1, figsize=(14, 12))

# Plot 1: Stationary
axes[0].plot(range(n), stationary_data, label='Actual', linewidth=1.5, alpha=0.7)
axes[0].plot(range(train_size, n), pred_ses, label='SES Forecast',
linewidth=2, color='red')
axes[0].axvline(x=train_size, color='gray', linestyle='--', alpha=0.5)
axes[0].set_title(f'Simple Exponential Smoothing (RMSE: {rmse_ses:.2f})')
axes[0].set_ylabel('Value')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Plot 2: Trend
axes[1].plot(range(n), trend_data, label='Actual', linewidth=1.5, alpha=0.7)
axes[1].plot(range(train_size, n), pred_holt, label='Holt Forecast',
linewidth=2, color='red')
axes[1].axvline(x=train_size, color='gray', linestyle='--', alpha=0.5)
axes[1].set_title(f"Holt's Linear Trend (RMSE: {rmse_holt:.2f})")
axes[1].set_ylabel('Value')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

# Plot 3: Seasonal
axes[2].plot(range(n), seasonal_data, label='Actual', linewidth=1.5, alpha=0.7)
axes[2].plot(range(train_size, n), pred_hw, label='Holt-Winters Forecast',
linewidth=2, color='red')
axes[2].axvline(x=train_size, color='gray', linestyle='--', alpha=0.5)
axes[2].set_title(f'Holt-Winters (RMSE: {rmse_hw:.2f})')
axes[2].set_ylabel('Value')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print("\n" + "="*50)
print("CHOOSING EXPONENTIAL SMOOTHING METHOD")
print("="*50)

print("""
DECISION GUIDE:

1. Plot your data
2. Check for trend:
- No trend + no seasonality → Simple Exponential Smoothing
- Trend present → Holt's Linear Trend
- Trend + seasonality → Holt-Winters

3. If using Holt-Winters:
- Additive: Seasonal variation constant
- Multiplicative: Seasonal variation increases with level

ADVANTAGES:
✓ Fast to implement
✓ Good for short-term forecasts
✓ Intuitive interpretation
✓ Requires little data
✓ Works well in practice

DISADVANTAGES:
✗ Limited to short-term (trend line extrapolated)
✗ No uncertainty intervals (in simple form)
✗ Can overfit with wrong parameters
✗ Doesn't handle complex patterns

WHEN TO USE:
- Need quick forecast
- Short time horizon (days/weeks)
- Limited data
- Trend and seasonality are stable
""")

import numpy as np
import pandas as pd
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt

print("="*50)
print("ARIMA MODELS")
print("="*50)

print("""
ARIMA = AutoRegressive Integrated Moving Average

COMPONENTS:

1. AR (AutoRegressive): p
- Use past values to predict future
- Y_t = c + φ₁*Y_(t-1) + φ₂*Y_(t-2) + ... + ε
- "Regress on itself"
- p = how many past values to use

2. I (Integrated): d
- Number of differences to make stationary
- d=0: Already stationary
- d=1: One difference needed
- d=2: Two differences needed (rare)

3. MA (Moving Average): q
- Use past forecast errors
- Y_t = μ + ε_t + θ₁*ε_(t-1) + θ₂*ε_(t-2) + ...
- Smooth out noise
- q = how many past errors to use

ARIMA(p,d,q) NOTATION:
- (0,0,0): White noise (can't predict)
- (1,0,0): AR(1) - depends on previous value
- (0,1,0): Random walk - difference once to get constant
- (0,0,1): MA(1) - depends on previous error
- (1,1,1): Typical time series model

HOW TO CHOOSE p, d, q:

d (differencing):
1. Plot data
2. If trend: d=1
3. If still non-stationary: d=2
4. ADF test: Until p < 0.05

p and q (from ACF/PACF plots):
- ACF plot: Look for MA pattern (q)
* Cuts off after lag q → MA(q)
* Tails off → AR or ARMA

- PACF plot: Look for AR pattern (p)
* Cuts off after lag p → AR(p)
* Tails off → MA or ARMA

PRACTICAL APPROACH:
1. Make data stationary (find d)
2. Plot ACF/PACF on differenced data
3. Try several (p,q) combinations
4. Use AIC/BIC to choose best
5. Check residuals (should be white noise)
""")

# Create example data
np.random.seed(42)
n = 200
t = np.arange(n)

# Non-stationary with trend
data = 50 + 0.2*t + np.random.normal(0, 2, n)

# Split
train_size = int(0.8 * n)
train, test = data[:train_size], data[train_size:]

print("\nFITTING ARIMA MODELS:")
print("-"*50)

# Try different ARIMA orders
models_to_try = [(1,0,0), (1,1,0), (0,1,0), (1,1,1), (2,1,1)]

results = []
for order in models_to_try:
try:
model = ARIMA(train, order=order)
fitted_model = model.fit()
pred = fitted_model.forecast(steps=len(test))
rmse = np.sqrt(mean_squared_error(test, pred))
aic = fitted_model.aic

results.append({
'order': order,
'rmse': rmse,
'aic': aic
})
print(f"ARIMA{order}: RMSE={rmse:.2f}, AIC={aic:.1f}")
except:
print(f"ARIMA{order}: Failed to fit")

# Find best model
best_order = min(results, key=lambda x: x['aic'])['order']
print(f"\nBest model: ARIMA{best_order}")

# Fit best model
best_model = ARIMA(train, order=best_order).fit()
best_pred = best_model.forecast(steps=len(test))
best_rmse = np.sqrt(mean_squared_error(test, best_pred))

print(f"Best model RMSE: {best_rmse:.2f}")
print(f"\nModel Summary:")
print(best_model.summary())

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Plot 1: Original data and forecast
axes[0, 0].plot(range(n), data, label='Actual', linewidth=1.5, alpha=0.7)
axes[0, 0].plot(range(train_size, n), best_pred, label=f'ARIMA{best_order} Forecast',
linewidth=2, color='red')
axes[0, 0].axvline(x=train_size, color='gray', linestyle='--', alpha=0.5)
axes[0, 0].set_title(f'ARIMA Model Forecast (RMSE: {best_rmse:.2f})')
axes[0, 0].set_ylabel('Value')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# Plot 2: ACF of original data
plot_acf(train, lags=30, ax=axes[0, 1])
axes[0, 1].set_title('ACF of Original Data')

# Plot 3: Residuals
residuals = best_model.resid
axes[1, 0].plot(residuals, linewidth=1)
axes[1, 0].axhline(y=0, color='red', linestyle='--')
axes[1, 0].set_title('Residuals (Should be White Noise)')
axes[1, 0].set_ylabel('Residual')
axes[1, 0].grid(True, alpha=0.3)

# Plot 4: Residual distribution
axes[1, 1].hist(residuals, bins=20, edgecolor='black')
axes[1, 1].set_title('Residual Distribution')
axes[1, 1].set_xlabel('Residual')
axes[1, 1].grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.show()

print("\n" + "="*50)
print("CHECKING MODEL ADEQUACY")
print("="*50)

print("""
RESIDUAL DIAGNOSTICS:
After fitting ARIMA, check residuals:

1. SHOULD LOOK LIKE WHITE NOISE:
- Mean = 0
- Constant variance
- No patterns
- Normally distributed

2. ACF OF RESIDUALS:
- Should not have significant spikes
- If spikes → Model missing patterns
- Solution: Try different (p,d,q)

3. LJUNG-BOX TEST:
- Tests if residuals are independent
- H₀: Residuals are white noise
- p > 0.05: Good (residuals are white noise)
- p < 0.05: Bad (model missing something)

4. FORECAST ERRORS:
- Directional accuracy: % of correct up/down
- Magnitude errors: RMSE, MAE
- Forecast bias: Mean of errors

IF DIAGNOSTICS BAD:
→ Try different (p,d,q)
→ Check for outliers
→ Consider other methods
→ May have structural breaks

import numpy as np
import pandas as pd
from statsmodels.tsa.statespace.sarimax import SARIMAX
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt

print("="*50)
print("SARIMA - SEASONAL ARIMA")
print("="*50)

print("""
SARIMA(p,d,q)(P,D,Q,s):
- ARIMA + SEASONAL components
- For data with clear seasonal patterns

PARAMETERS:
(p,d,q): Non-seasonal ARIMA (same as ARIMA)
(P,D,Q,s): Seasonal ARIMA
- P: Seasonal AR order
- D: Seasonal differencing
- Q: Seasonal MA order
- s: Seasonal period (frequency of pattern)

Examples:
- Monthly data, yearly seasonality: s=12
- Daily data, weekly seasonality: s=7
- Quarterly data, yearly seasonality: s=4

WHEN TO USE:
- Clear seasonal patterns visible
- ACF/PACF show spikes at seasonal lags
- Holt-Winters worked but want ARIMA
- Need probabilistic forecasts

CHOOSING PARAMETERS:
1. Determine seasonal period (s)
- Look at data
- Auto-detect with autocorrelation

2. Regular differencing (d)
- Until trend removed

3. Seasonal differencing (D)
- Usually 0 or 1

4. Seasonal AR/MA (P, Q)
- From ACF/PACF at seasonal lags

5. Non-seasonal AR/MA (p, q)
- From ACF/PACF of differenced data

EXAMPLE: Electricity usage
- Daily data (365 days/year)
- Strong weekly seasonality (7 days)
- Gradual growth trend
→ SARIMA(1,1,1)(1,0,1,7)
→ Try auto_arima()
""")

# Create seasonal data
np.random.seed(42)
n = 365
t = np.arange(n)

# Seasonal pattern (weekly)
seasonal = 5 * np.sin(2 * np.pi * t / 7)

# Trend
trend = 0.01 * t

# Noise
noise = np.random.normal(0, 1, n)

# Combined
data_seasonal = 50 + trend + seasonal + noise

# Split
train_size = int(0.8 * n)
train_s, test_s = data_seasonal[:train_size], data_seasonal[train_size:]

print("\nFITTING SARIMA MODELS:")
print("-"*50)

# Fit SARIMA
model_sarima = SARIMAX(train_s, order=(1,1,1), seasonal_order=(1,0,1,7))
fitted_sarima = model_sarima.fit(disp=False)

pred_sarima = fitted_sarima.forecast(steps=len(test_s))
rmse_sarima = np.sqrt(mean_squared_error(test_s, pred_sarima))

print(f"SARIMA(1,1,1)(1,0,1,7):")
print(f" RMSE: {rmse_sarima:.2f}")

# Get confidence intervals
pred_sarima_ci = fitted_sarima.get_forecast(steps=len(test_s))
pred_ci = pred_sarima_ci.conf_int()

# Visualization
fig, axes = plt.subplots(2, 1, figsize=(14, 10))

# Plot 1: Forecast with confidence intervals
axes[0].plot(range(n), data_seasonal, label='Actual', linewidth=1.5, alpha=0.7)
axes[0].plot(range(train_size, n), pred_sarima, label='SARIMA Forecast',
linewidth=2, color='red')
axes[0].fill_between(range(train_size, n),
pred_ci.iloc[:, 0], pred_ci.iloc[:, 1],
alpha=0.2, color='red', label='95% CI')
axes[0].axvline(x=train_size, color='gray', linestyle='--', alpha=0.5)
axes[0].set_title(f'SARIMA Forecast with Confidence Intervals (RMSE: {rmse_sarima:.2f})')
axes[0].set_ylabel('Value')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Plot 2: Residuals
residuals_s = fitted_sarima.resid
axes[1].plot(residuals_s, linewidth=1)
axes[1].axhline(y=0, color='red', linestyle='--')
axes[1].set_title('SARIMA Residuals')
axes[1].set_xlabel('Time')
axes[1].set_ylabel('Residual')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print("\nAUTO-ARIMA:")
print("-"*50)

print("""
AUTO_ARIMA: Automated parameter selection

Advantages:
✓ No need to manually choose p,d,q,P,D,Q
✓ Tests stationarity
✓ Searches parameter space
✓ Returns best model by AIC/BIC

Disadvantages:
✗ Black box (you don't understand choice)
✗ May be slow with large data
✗ Not always best (automatic choices)

Usage:
from statsmodels.tsa.arima.auto_arima import auto_arima
auto_model = auto_arima(data, seasonal=True, m=7)

Common approach:
1. Try auto_arima for quick baseline
2. Manually tune based on ACF/PACF
3. Compare models
4. Choose best by validation
""")

import numpy as np
import pandas as pd
from sklearn.metrics import mean_squared_error, mean_absolute_error, mean_absolute_percentage_error
import matplotlib.pyplot as plt

print("="*50)
print("EVALUATING TIME SERIES FORECASTS")
print("="*50)

print("""
EVALUATION METRICS:

1. RMSE (Root Mean Squared Error)
- √(Σ(actual - predicted)² / n)
- Units: Same as original data
- Penalizes large errors more
- Range: 0 to ∞ (lower better)
- Use: When all errors equally important

2. MAE (Mean Absolute Error)
- Average absolute error
- More robust to outliers
- Easy to interpret (average error)
- Use: When few large errors shouldn't dominate

3. MAPE (Mean Absolute Percentage Error)
- Average percentage error
- |actual - predicted| / |actual|
- Scale-independent
- Good for comparing across different scales
- Use: For relative accuracy

4. THEIL'S U STATISTIC
- Compares naive forecast to model
- U < 1: Model better than naive
- U = 1: Model = naive
- U > 1: Naive better

5. DIRECTIONAL ACCURACY
- % of time model predicts correct direction
- Did price go up/down correctly?
- Range: 0-100%
- Business metric (catches trends)

CHOOSING METRIC:
- RMSE: Default, penalizes large errors
- MAE: Robust, easy to explain
- MAPE: Compare across different products
- Directional: For trading/investment
""")

# Create comparison data
np.random.seed(42)
n = 100
actual = 50 + 0.1*np.arange(n) + 10*np.sin(np.arange(n)/10) + np.random.normal(0,2,n)

# Three forecasts (naive, simple model, good model)
naive_forecast = np.concatenate([[actual[0]], actual[:-1]]) # Previous value

simple_forecast = actual.mean() + 10*np.sin(np.arange(n)/10) # Only seasonal

good_forecast = actual + np.random.normal(0,0.5,n) # Near perfect

# Calculate metrics
def calc_metrics(actual, pred):
rmse = np.sqrt(mean_squared_error(actual, pred))
mae = mean_absolute_error(actual, pred)
mape = mean_absolute_percentage_error(actual, pred)

# Directional accuracy
actual_direction = np.sign(np.diff(actual))
pred_direction = np.sign(np.diff(pred))
dir_acc = np.mean(actual_direction == pred_direction) * 100

return {'RMSE': rmse, 'MAE': mae, 'MAPE': mape, 'Dir. Acc': dir_acc}

print("\nMETRIC COMPARISON:")
print("-"*50)

metrics_naive = calc_metrics(actual, naive_forecast)
metrics_simple = calc_metrics(actual, simple_forecast)
metrics_good = calc_metrics(actual, good_forecast)

print("\nNaive (previous value):")
for k, v in metrics_naive.items():
print(f" {k}: {v:.2f}")

print("\nSimple (seasonal only):")
for k, v in metrics_simple.items():
print(f" {k}: {v:.2f}")

print("\nGood (near actual):")
for k, v in metrics_good.items():
print(f" {k}: {v:.2f}")

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Plot 1: Forecasts comparison
axes[0, 0].plot(actual, label='Actual', linewidth=2)
axes[0, 0].plot(naive_forecast, label='Naive', alpha=0.7, linewidth=1.5)
axes[0, 0].plot(simple_forecast, label='Simple', alpha=0.7, linewidth=1.5)
axes[0, 0].plot(good_forecast, label='Good', alpha=0.7, linewidth=1.5)
axes[0, 0].set_title('Forecast Comparison')
axes[0, 0].set_ylabel('Value')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# Plot 2: Errors by method
errors = {
'Naive': np.abs(actual - naive_forecast),
'Simple': np.abs(actual - simple_forecast),
'Good': np.abs(actual - good_forecast)
}

axes[0, 1].boxplot([errors['Naive'], errors['Simple'], errors['Good']],
labels=['Naive', 'Simple', 'Good'])
axes[0, 1].set_title('Forecast Errors Distribution')
axes[0, 1].set_ylabel('Absolute Error')
axes[0, 1].grid(True, alpha=0.3, axis='y')

# Plot 3: Error over time (Naive)
axes[1, 0].plot(np.abs(actual - naive_forecast), label='Naive', linewidth=1.5)
axes[1, 0].plot(np.abs(actual - good_forecast), label='Good', linewidth=1.5)
axes[1, 0].set_title('Forecast Error Over Time')
axes[1, 0].set_ylabel('Absolute Error')
axes[1, 0].set_xlabel('Time')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

# Plot 4: Scatter of actual vs predicted
axes[1, 1].scatter(actual, good_forecast, alpha=0.6)
axes[1, 1].plot([actual.min(), actual.max()],
[actual.min(), actual.max()], 'r--', label='Perfect')
axes[1, 1].set_xlabel('Actual')
axes[1, 1].set_ylabel('Predicted')
axes[1, 1].set_title('Actual vs Predicted (Good Model)')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print("\n" + "="*50)
print("TIME SERIES FORECASTING BEST PRACTICES")
print("="*50)

print("""
1. PROPER DATA SPLITTING
✓ Use time-based split (not random)
✓ Train set: Earlier data
✓ Test set: Later data
✓ Validate on future data
✗ Don't: Random shuffle, look ahead

2. HANDLE MISSING DATA
✓ Identify missing values
✓ Interpolate if small gaps
✓ Forward/backward fill if needed
✓ Document approach
✗ Don't: Ignore missing values

3. OUTLIERS
✓ Identify extreme values
✓ Investigate: Errors or real events?
✓ If errors: Correct or remove
✓ If real: Keep (important for model)
✗ Don't: Blindly remove outliers

4. EVALUATE ON MULTIPLE METRICS
✓ Use RMSE, MAE, MAPE, directional
✓ Plot actual vs predicted
✓ Check residuals (white noise)
✓ Validate on unseen test data
✗ Don't: Only look at one metric

5. FORECAST UNCERTAINTY
✓ Report confidence intervals
✓ Further future = wider intervals
✓ Consider prediction intervals (not just CI)
✗ Don't: Report single point forecast

6. MODEL COMPARISON
✓ Try multiple methods
✓ Simple baseline (naive, mean)
✓ Linear models (exponential smoothing)
✓ ARIMA for complex patterns
✓ Machine learning if good predictors

7. AVOID COMMON PITFALLS
✗ DATA LEAKAGE: Using future info
✗ OVERFITTING: Too complex for amount of data
✗ NON-STATIONARITY: Not differencing
✗ SEASONAL IGNORANCE: Ignoring seasonality
✗ HORIZON: Forecasting too far into future

8. MONITORING IN PRODUCTION
✓ Track forecast accuracy over time
✓ Detect performance degradation
✓ Retrain with new data regularly
✓ Adapt to structural changes
✓ Update when accuracy drops

WORKFLOW:
1. Explore data → Visualize trends, seasonality
2. Make stationary → Differencing, transformations
3. Identify ARIMA parameters → ACF/PACF, ADF test
4. Fit model → Train on historical data
5. Validate → Evaluate on test data
6. Check residuals → Diagnostics
7. Forecast future → With uncertainty intervals
8. Monitor → Track performance, retrain as needed
""")