Session 21: Advanced Python for Insurance Analytics

1. Time Series in Insurance

Time series analysis is one of the most important quantitative tools in insurance analytics. Claims volumes and amounts vary systematically over time — driven by seasonality (monsoon claims spike, holiday travel accidents, winter health claims), trends (growing portfolio, medical inflation), and external shocks (regulatory changes, COVID-19). Understanding and forecasting these patterns is essential for reserving, pricing, operational planning, and capital management.

1.1 Why Time Series Matters in Insurance

1.2 Loading and Preparing Time Series Data

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller, acf, pacf
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from sklearn.linear_model import LinearRegression, Ridge, Lasso
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
import warnings
warnings.filterwarnings('ignore')

# Load the claims dataset
df = pd.read_csv('data/insurance_cleaned.csv')
df['claim_date'] = pd.to_datetime(df['claim_date'])

# Aggregate claims by month
monthly_claims = df.set_index('claim_date').resample('M').agg(
    claim_count=('claim_id', 'count'),
    total_claim_amount=('claim_amount', 'sum'),
    avg_claim_amount=('claim_amount', 'mean'),
    premium_amount=('premium', 'sum')
).reset_index()

monthly_claims = monthly_claims.sort_values('claim_date')
monthly_claims = monthly_claims.set_index('claim_date')

print(f"Monthly claims data: {len(monthly_claims):,} months")
print(f"Date range: {monthly_claims.index.min().date()} to {monthly_claims.index.max().date()}")
print(f"\nSummary statistics:")
print(monthly_claims[['claim_count', 'total_claim_amount']].describe())

# Plot the series
fig, axes = plt.subplots(2, 1, figsize=(14, 8))

axes[0].plot(monthly_claims.index, monthly_claims['claim_count'],
             color='#6c5ce7', linewidth=1.5, marker='o', markersize=3)
axes[0].set_title('Monthly Insurance Claims Count — Time Series', fontweight='bold')
axes[0].set_ylabel('Number of Claims')
axes[0].spines['top'].set_visible(False)
axes[0].spines['right'].set_visible(False)

axes[1].plot(monthly_claims.index, monthly_claims['total_claim_amount'] / 1e7,
             color='#e17055', linewidth=1.5, marker='s', markersize=3)
axes[1].set_title('Monthly Total Claim Amount (₹ Crores)', fontweight='bold')
axes[1].set_ylabel('Total Claim Amount (₹ Cr)')
axes[1].set_xlabel('Date')
axes[1].spines['top'].set_visible(False)
axes[1].spines['right'].set_visible(False)

plt.tight_layout()
plt.show()

print("Observe: Is there a visible trend? (upward/downward/stable)")
print("          Is there seasonality? (regular annual pattern)")
print("          Are there outliers? (one-off spikes or dips)")

2. Time Series Decomposition

Decomposition separates a time series into its constituent components — trend (long-term direction), seasonality (regular periodic patterns), and residual (random noise/irregular). Decomposition is the first step in understanding what drives the series, which components are predictable, and where the forecasting challenge lies.

2.1 Additive vs. Multiplicative Decomposition

# Decompose the claims count time series
# We use the additive model: Y(t) = Trend(t) + Seasonal(t) + Residual(t)
# For multiplicative: Y(t) = Trend(t) × Seasonal(t) × Residual(t)
# Multiplicative is appropriate when the seasonal amplitude grows with the trend level

from statsmodels.tsa.seasonal import seasonal_decompose

# Ensure we have at least 2 full periods of data (24 months for 12-month seasonality)
result_add = seasonal_decompose(monthly_claims['claim_count'], model='additive', period=12)
result_mul = seasonal_decompose(monthly_claims['claim_count'], model='multiplicative', period=12)

fig, axes = plt.subplots(4, 2, figsize=(16, 10), sharex=True)

# Original
axes[0, 0].plot(result_add.observed, color='#6c5ce7', linewidth=1.5)
axes[0, 0].set_ylabel('Original')
axes[0, 0].set_title('Additive Decomposition')
axes[0, 0].spines['top'].set_visible(False)
axes[0, 0].spines['right'].set_visible(False)

axes[0, 1].plot(result_mul.observed, color='#6c5ce7', linewidth=1.5)
axes[0, 1].set_ylabel('Original')
axes[0, 1].set_title('Multiplicative Decomposition')
axes[0, 1].spines['top'].set_visible(False)
axes[0, 1].spines['right'].set_visible(False)

# Trend
axes[1, 0].plot(result_add.trend, color='#e17055', linewidth=2)
axes[1, 0].set_ylabel('Trend')
axes[1, 0].spines['top'].set_visible(False)
axes[1, 0].spines['right'].set_visible(False)

axes[1, 1].plot(result_mul.trend, color='#e17055', linewidth=2)
axes[1, 1].set_ylabel('Trend')
axes[1, 1].spines['top'].set_visible(False)
axes[1, 1].spines['right'].set_visible(False)

# Seasonal
axes[2, 0].plot(result_add.seasonal, color='#00d2d3', linewidth=1.5)
axes[2, 0].set_ylabel('Seasonal')
axes[2, 0].spines['top'].set_visible(False)
axes[2, 0].spines['right'].set_visible(False)

axes[2, 1].plot(result_mul.seasonal, color='#00d2d3', linewidth=1.5)
axes[2, 1].set_ylabel('Seasonal')
axes[2, 1].spines['top'].set_visible(False)
axes[2, 1].spines['right'].set_visible(False)

# Residual
axes[3, 0].plot(result_add.resid, color='#74b9ff', linewidth=1, marker='o', markersize=2)
axes[3, 0].axhline(y=0, color='gray', linestyle='--', linewidth=0.5)
axes[3, 0].set_ylabel('Residual')
axes[3, 0].set_xlabel('Date')
axes[3, 0].spines['top'].set_visible(False)
axes[3, 0].spines['right'].set_visible(False)

axes[3, 1].plot(result_mul.resid, color='#74b9ff', linewidth=1, marker='o', markersize=2)
axes[3, 1].axhline(y=1, color='gray', linestyle='--', linewidth=0.5)
axes[3, 1].set_ylabel('Residual')
axes[3, 1].set_xlabel('Date')
axes[3, 1].spines['top'].set_visible(False)
axes[3, 1].spines['right'].set_visible(False)

plt.tight_layout()
plt.show()

# Check which model has more stable residuals (lower variance, fewer extreme values)
add_resid_var = result_add.resid.var()
mul_resid_var = result_mul.resid.var()

print(f"Residual variance — Additive: {add_resid_var:.2f}")
print(f"Residual variance — Multiplicative: {mul_resid_var:.2f}")
better_model = 'Additive' if add_resid_var < mul_resid_var else 'Multiplicative'
print(f"✓ The {better_model} model has lower residual variance — better fit.")
print(f"  (Lower variance in residuals means more of the pattern is captured by trend + seasonality)")

# Extract seasonal pattern
if result_add.seasonal is not None:
    seasonal_pattern = result_add.seasonal.groupby(result_add.seasonal.index.month).mean()
    print(f"\nAverage seasonal impact on claims by month:")
    for month in range(1, 13):
        impact = seasonal_pattern[month]
        print(f"  Month {month:>2d}: {'+' if impact > 0 else ''}{impact:.0f} claims (vs. trend)")
    peak_month = seasonal_pattern.idxmax()
    trough_month = seasonal_pattern.idxmin()
    print(f"\nPeak month:   Month {peak_month} ({seasonal_pattern[peak_month]:.0f} claims above trend)")
    print(f"Trough month: Month {trough_month} ({seasonal_pattern[trough_month]:.0f} claims below trend)")

3. ARIMA Modelling

ARIMA (AutoRegressive Integrated Moving Average) is the most widely used statistical method for time series forecasting. It models the time series as a function of: its own past values (AR — auto-regressive component), past forecast errors (MA — moving average component), and differencing to make the series stationary (I — integrated component). An ARIMA model is specified by three parameters: p (AR order), d (differencing order), and q (MA order).

3.1 Testing for Stationarity

ARIMA requires the time series to be stationary — meaning the statistical properties (mean, variance, autocorrelation) do not change over time. Most insurance time series are not stationary — they have trends (growing portfolio, inflation) and seasonality. Differencing removes these non-stationary components.

# Augmented Dickey-Fuller test for stationarity
def test_stationarity(series, name):
    result = adfuller(series.dropna())
    print(f"ADF Test — {name}")
    print(f"  ADF Statistic: {result[0]:.4f}")
    print(f"  p-value:       {result[1]:.4f}")
    print(f"  Critical values:")
    for key, val in result[4].items():
        print(f"    {key}: {val:.4f}")
    if result[1] < 0.05:
        print(f"  ✓ p < 0.05 — Series IS stationary. No differencing needed (d=0).")
        return 0
    else:
        print(f"  ✗ p ≥ 0.05 — Series is NOT stationary. Differencing needed (d ≥ 1).")
        return 1

# Test the original series
ts = monthly_claims['claim_count'].dropna()
d_order = test_stationarity(ts, 'Original Claims Count')

# If non-stationary, difference once and test again
if d_order > 0:
    ts_diff1 = ts.diff().dropna()
    d_order = test_stationarity(ts_diff1, 'Claims Count — First Difference')
    if d_order > 0:
        ts_diff2 = ts_diff1.diff().dropna()
        d_order = test_stationarity(ts_diff2, 'Claims Count — Second Difference')
        final_d = 2
    else:
        final_d = 1
else:
    final_d = 0

print(f"\nRecommended differencing order (d) for ARIMA: {final_d}")

3.2 Selecting p and q — ACF and PACF Plots

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

# ACF and PACF of original series
plot_acf(ts.dropna(), lags=24, ax=axes[0, 0])
axes[0, 0].set_title('ACF — Original Series')
plot_pacf(ts.dropna(), lags=24, ax=axes[0, 1], method='ywm')
axes[0, 1].set_title('PACF — Original Series')

# ACF and PACF of differenced series
differenced = ts.diff().dropna()
plot_acf(differenced.dropna(), lags=24, ax=axes[1, 0])
axes[1, 0].set_title(f'ACF — Differenced (d={final_d}) Series')
plot_pacf(differenced.dropna(), lags=24, ax=axes[1, 1], method='ywm')
axes[1, 1].set_title(f'PACF — Differenced (d={final_d}) Series')

plt.tight_layout()
plt.show()

print("How to read ACF/PACF plots for ARIMA parameter selection:")
print("  p (AR order): Look at the PACF — the number of significant lags before")
print("                they cut off (become insignificant) is the suggested p.")
print("  q (MA order): Look at the ACF — the number of significant lags before")
print("                they cut off is the suggested q.")
print("  d (differencing): The number of differencing steps needed to achieve stationarity.")
print()
print("Suggested starting parameters from visual inspection:")
print("  Based on the plots above, try ARIMA(p=2, d=1, q=2) as a starting point.")

3.3 Fitting the ARIMA Model

# Split into training and test sets (last 12 months for testing)
train = ts.iloc[:-12]
test = ts.iloc[-12:]

print(f"Training data: {train.index[0].date()} to {train.index[-1].date()} ({len(train)} months)")
print(f"Test data:     {test.index[0].date()} to {test.index[-1].date()} ({len(test)} months)")

# Fit ARIMA model
# Order = (p, d, q) — starting with p=2, d=1, q=2 based on ACF/PACF analysis
arima_order = (2, 1, 2)
model = ARIMA(train, order=arima_order)
fitted_model = model.fit()

print(f"\nARIMA{arima_order} Model Results:")
print(f"Log Likelihood: {fitted_model.llf:.2f}")
print(f"AIC: {fitted_model.aic:.2f}")
print(f"BIC: {fitted_model.bic:.2f}")
print(f"\nCoefficients:")
for param, value, pval in zip(fitted_model.params.index, fitted_model.params.values, fitted_model.pvalues):
    print(f"  {param:15s}: {value:>8.4f} (p={pval:.4f})")

# Forecast
forecast_result = fitted_model.get_forecast(steps=12)
forecast_mean = forecast_result.predicted_mean
forecast_ci = forecast_result.conf_int()

# Assess accuracy
if len(test) == 12:
    mae = mean_absolute_error(test, forecast_mean)
    rmse = np.sqrt(mean_squared_error(test, forecast_mean))
    mape = np.mean(np.abs((test.values - forecast_mean.values) / test.values)) * 100
    print(f"\nForecast Accuracy (next 12 months):")
    print(f"  MAE:  {mae:.0f} claims")
    print(f"  RMSE: {rmse:.0f} claims")
    print(f"  MAPE: {mape:.1f}%")
    print(f"  Mean monthly claims: {test.mean():.0f}")
    print(f"  Error as % of mean:  {mae/test.mean()*100:.1f}%")

3.4 Visualising the Forecast

# Plot historical data and forecast
fig, ax = plt.subplots(figsize=(14, 6))

# Historical
ax.plot(train.index[-36:], train[-36:], color='#6c5ce7', linewidth=2, label='Training Data (last 36 months)')
ax.plot(test.index, test.values, color='#00b894', linewidth=2, marker='o',
        markersize=6, label='Actual (test period)')

# Forecast
ax.plot(test.index, forecast_mean.values, color='#e17055', linewidth=2, linestyle='--',
        marker='s', markersize=6, label='ARIMA Forecast')
ax.fill_between(test.index,
                forecast_ci.iloc[:, 0].values,
                forecast_ci.iloc[:, 1].values,
                alpha=0.2, color='#e17055', label='95% Confidence Interval')

ax.set_xlabel('Date')
ax.set_ylabel('Monthly Claims Count')
ax.set_title(f'ARIMA{arima_order} — Claims Forecast (12 Months)', fontweight='bold')
ax.legend(loc='upper left')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.tight_layout()
plt.show()

print("Key questions after reviewing the forecast:")
print("1. Does the forecast capture the seasonal pattern? (Should show peaks and troughs)")
print("2. Are the confidence intervals widening? (Wider = more uncertainty further out)")
print("3. Is the forecast systematically over/under predicting? (Bias check)")
print("4. Does the actual test data fall within the confidence intervals?")

4. Model Evaluation — Beyond Basic Metrics

Forecast accuracy metrics (MAE, RMSE, MAPE) tell part of the story — but in insurance, we care about more than how close the forecast is to the actual. We also care about whether the forecast is unbiased, whether it is useful compared to a simple benchmark, and whether the directional accuracy is good.

4.1 Comparing Against a Naive Benchmark

# Naive forecast: assume next 12 months = same as last 12 months (seasonal naive)
naive_forecast = train[-12:].values  # Use last 12 months of training data as forecast

# This is actually a "seasonal naive" model — it uses the most recent complete seasonal cycle
# as the forecast for the next cycle. It's surprisingly hard to beat for strongly seasonal series.

# Or simpler: last value repeated (flat forecast)
naive_flat = np.full(12, train[-1])

# Evaluate naive models
mae_seasonal_naive = mean_absolute_error(test, naive_forecast)
mae_arima = mean_absolute_error(test, forecast_mean)
rmse_seasonal_naive = np.sqrt(mean_squared_error(test, naive_forecast))
rmse_arima = np.sqrt(mean_squared_error(test, forecast_mean))

print("=" * 60)
print("MODEL COMPARISON — ARIMA vs. NAIVE BENCHMARKS")
print("=" * 60)
print(f"{'Model':30s} {'MAE':>12s} {'RMSE':>12s}")
print("-" * 54)
print(f"{'Seasonal Naive (last 12 months)':30s} {mae_seasonal_naive:>8.0f}     {rmse_seasonal_naive:>8.0f}")
print(f"{'ARIMA' + str(arima_order)}{'':30s} {mae_arima:>8.0f}     {rmse_arima:>8.0f}")

if mae_arima < mae_seasonal_naive:
    improvement = (mae_seasonal_naive - mae_arima) / mae_seasonal_naive * 100
    print(f"\n✓ ARIMA outperforms seasonal naive by {improvement:.1f}%")
else:
    print(f"\n✗ Seasonal naive outperforms ARIMA — the series may be too noisy for ARIMA.")
    print(f"  Consider SARIMA (Seasonal ARIMA) which explicitly models seasonality.")

4.2 Directional Accuracy and Bias

# Directional accuracy: does the model correctly predict whether claims will go up or down?
actual_changes = np.sign(np.diff(test.values))
forecast_changes = np.sign(np.diff(forecast_mean.values))
directional_accuracy = (actual_changes == forecast_changes[:len(actual_changes)]).mean()

print(f"\nDirectional Accuracy: {directional_accuracy*100:.0f}%")
print(f"(The model correctly predicts whether claims will increase or decrease")
print(f" month-over-month {directional_accuracy*100:.0f}% of the time)")

# Bias: is the model systematically over or under forecasting?
bias = (forecast_mean.values - test.values).mean()
bias_pct = bias / test.mean() * 100
print(f"\nForecast Bias: {bias:+.0f} claims per month ({bias_pct:+.1f}% of mean)")

if abs(bias_pct) < 5:
    print(f"✓ Bias is within acceptable range (< 5%) — model is not systematically over/under forecasting.")
elif bias_pct > 0:
    print(f"△ Model systematically OVER-forecasts by {bias_pct:.1f}% — may lead to over-reserving.")
else:
    print(f"△ Model systematically UNDER-forecasts by {bias_pct:.1f}% — may lead to under-reserving.")

# Track where the forecast is within the confidence interval
in_ci = ((test.values >= forecast_ci.iloc[:, 0].values) &
         (test.values <= forecast_ci.iloc[:, 1].values)).mean()
print(f"\nActual values within 95% CI: {in_ci*100:.0f}%")
print(f"(If the model is well-calibrated, ~95% of actual values should be within the CI.)")

5. Premium Prediction with Multiple Regression

Beyond time series, insurers need to predict the relationship between policyholder characteristics and premium amounts — for pricing optimisation, underwriting, and customer segmentation. Multiple linear regression with regularisation (Ridge and Lasso) is the standard approach.

5.1 Feature Preparation

# Build a dataset for premium prediction — one row per policy
policy_df = df.groupby('policy_id').agg(
    premium=('premium', 'first'),
    sum_assured=('sum_assured', 'first'),
    age=('age', 'first'),
    income=('income', 'first'),
    credit_score=('credit_score', 'first'),
    policy_type=('policy_type', 'first'),
    channel=('channel', 'first'),
    tenure_months=('tenure_months', 'first'),
    claim_count=('claim_id', 'count'),
    total_claim_amount=('claim_amount', 'sum')
).reset_index()

# Create derived features
policy_df['claim_flag'] = (policy_df['claim_count'] > 0).astype(int)
policy_df['age_squared'] = policy_df['age'] ** 2
policy_df['risk_income_ratio'] = policy_df['sum_assured'] / (policy_df['income'] + 1)
policy_df['claim_severity'] = policy_df['total_claim_amount'] / (policy_df['claim_count'] + 1)

# Drop missing
model_data = policy_df.dropna(subset=['premium']).copy()

print(f"Premium prediction dataset: {len(model_data):,} policies, {model_data.shape[1]} features")
print(f"Premium range: ₹{model_data['premium'].min():,.0f} — ₹{model_data['premium'].max():,.0f}")
print(f"Average premium: ₹{model_data['premium'].mean():,.0f}")

5.2 OLS, Ridge, and Lasso

from sklearn.preprocessing import OneHotEncoder, StandardScaler
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.model_selection import cross_val_score

# Define features
num_features = ['age', 'income', 'credit_score', 'sum_assured', 'tenure_months',
                'age_squared', 'risk_income_ratio', 'claim_count']
cat_features = ['policy_type', 'channel']

# Limit features to available ones
num_features = [c for c in num_features if c in model_data.columns]
cat_features = [c for c in cat_features if c in model_data.columns]

X = model_data[num_features + cat_features]
y = model_data['premium']

# Preprocessor
preprocessor = ColumnTransformer(
    transformers=[
        ('num', StandardScaler(), num_features),
        ('cat', OneHotEncoder(drop='first', handle_unknown='ignore', sparse_output=False), cat_features)
    ])

# Split
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

# Transform
X_train_t = preprocessor.fit_transform(X_train)
X_test_t = preprocessor.transform(X_test)

# OLS (Linear Regression)
ols = LinearRegression()
ols.fit(X_train_t, y_train)
y_pred_ols = ols.predict(X_test_t)

# Ridge
ridge = Ridge(alpha=10)
ridge.fit(X_train_t, y_train)
y_pred_ridge = ridge.predict(X_test_t)

# Lasso
lasso = Lasso(alpha=1)
lasso.fit(X_train_t, y_train)
y_pred_lasso = lasso.predict(X_test_t)

# Compare
for name, pred in [('OLS', y_pred_ols), ('Ridge', y_pred_ridge), ('Lasso', y_pred_lasso)]:
    mae = mean_absolute_error(y_test, pred)
    rmse = np.sqrt(mean_squared_error(y_test, pred))
    r2 = r2_score(y_test, pred)
    non_zero = "all" if name != 'Lasso' else f"{np.sum(lasso.coef_ != 0)}"
    print(f"{'─' * 50}")
    print(f"{name} Regression:")
    print(f"  MAE:  ₹{mae:,.0f} ({mae/y_test.mean()*100:.1f}% of mean premium)")
    print(f"  RMSE: ₹{rmse:,.0f}")
    print(f"  R²:   {r2:.3f}")
    print(f"  Features used: {non_zero}")

5.3 Cross-Validation for Regularisation Strength

from sklearn.linear_model import RidgeCV, LassoCV

# Ridge with cross-validated alpha
ridge_cv = RidgeCV(alist=[0.1, 1, 10, 50, 100], cv=5)
ridge_cv.fit(X_train_t, y_train)
print(f"Best Ridge alpha (CV): {ridge_cv.alpha_}")
print(f"Ridge CV R²: {r2_score(y_test, ridge_cv.predict(X_test_t)):.3f}")

# Lasso with cross-validated alist
lasso_cv = LassoCV(alphas=[0.1, 0.5, 1, 5, 10], cv=5, max_iter=10000)
lasso_cv.fit(X_train_t, y_train)
print(f"Best Lasso alpha (CV): {lasso_cv.alpha_}")
print(f"Lasso CV R²: {r2_score(y_test, lasso_cv.predict(X_test_t)):.3f}")
print(f"Features retained by Lasso: {np.sum(lasso_cv.coef_ != 0)} of {X_train_t.shape[1]}")

# Which features survived Lasso regularisation?
feature_names = num_features + list(
    preprocessor.named_transformers_['cat']
    .get_feature_names_out(cat_features)
) if len(cat_features) > 0 else num_features

retained_features = [fn for fn, coef in zip(feature_names, lasso_cv.coef_)
                     if abs(coef) > 0]
print(f"\nFeatures retained by Lasso ({len(retained_features)} of {len(feature_names)}):")
for fn in retained_features:
    print(f"  ✓ {fn}")

6. Integrated Claims Forecasting Pipeline

The most powerful approach combines time series and regression techniques into an integrated forecasting pipeline. The time series model captures the temporal dynamics (trend, seasonality), while the regression model captures the relationship between claims and external drivers (policy growth, economic indicators, weather). The combined forecast is more accurate and more interpretable than either model alone.

6.1 Building the Combined Pipeline

# Create additional features for the time series model
ts_data = monthly_claims[['claim_count']].copy()
ts_data['month'] = ts_data.index.month
ts_data['quarter'] = ts_data.index.quarter
ts_data['year'] = ts_data.index.year
ts_data['t'] = np.arange(len(ts_data))

# Add lag features
ts_data['lag_1'] = ts_data['claim_count'].shift(1)
ts_data['lag_12'] = ts_data['claim_count'].shift(12)
ts_data['rolling_3m_mean'] = ts_data['claim_count'].rolling(3).mean()
ts_data['rolling_12m_mean'] = ts_data['claim_count'].rolling(12).mean()

# Drop NaN from lags
ts_data = ts_data.dropna()

# Define features for the regression-on-time-series model
X_ts = ts_data[['month', 'quarter', 't', 'lag_1', 'lag_12', 'rolling_3m_mean']]
y_ts = ts_data['claim_count']

# Re-split by time
split_idx = len(X_ts) - 12
X_train_ts = X_ts.iloc[:split_idx]
X_test_ts = X_ts.iloc[split_idx:]
y_train_ts = y_ts.iloc[:split_idx]
y_test_ts = y_ts.iloc[split_idx:]

# Train Ridge on the time series features
ts_reg = Ridge(alpha=1)
ts_reg.fit(X_train_ts, y_train_ts)
y_pred_ts = ts_reg.predict(X_test_ts)

# Compare with ARIMA
mae_ts_reg = mean_absolute_error(y_test_ts, y_pred_ts)
print("=" * 55)
print("INTEGRATED FORECASTING PIPELINE — COMPARISON")
print("=" * 55)
print(f"{'Method':35s} {'MAE':>10s} {'RMSE':>10s}")
print("-" * 55)
print(f"{'ARIMA' + str(arima_order):35s} {mae_arima:>8.0f}   {rmse_arima:>8.0f}")
print(f"{'Seasonal Naive':35s} {mae_seasonal_naive:>8.0f}   {rmse_seasonal_naive:>8.0f}")
print(f"{'Ridge on Lag Features':35s} {mae_ts_reg:>8.0f}   {np.sqrt(mean_squared_error(y_test_ts, y_pred_ts)):>8.0f}")

# Final recommendation
best_model = min([
    ('ARIMA', mae_arima),
    ('Seasonal Naive', mae_seasonal_naive),
    ('Ridge on Lags', mae_ts_reg)
], key=lambda x: x[1])

print(f"\nBest model: {best_model[0]} (MAE = {best_model[1]:.0f})")