model.py

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score, mean_squared_error

# expects df with columns: date, prime
df = pd.read_csv('FRB_H15.csv')
df = df.copy()
df['date'] = pd.to_datetime(df['date'])
df['year'] = df['date'].dt.year

# yearly average prime rate
yearly = (
    df.groupby('year', as_index=False)['prime']
      .mean()
      .rename(columns={'prime': 'avg_prime'})
      .sort_values('year')
      .reset_index(drop=True)
)

# training window (best bang-for-buck for regime changes)
window_years = 30
last_year = int(yearly['year'].max())
train_start_year = max(int(yearly['year'].min()), last_year - window_years + 1)

train = yearly[yearly['year'] >= train_start_year].copy()

# model year-over-year change (delta), not level
train['delta'] = train['avg_prime'].diff()
train = train.dropna().reset_index(drop=True)

# features
t0 = int(train['year'].min())
train['t'] = train['year'] - t0

X = train[['t']]
y = train['delta']

degree = 2
delta_model = Pipeline([
    ('poly', PolynomialFeatures(degree=degree, include_bias=False)),
    ('lin', LinearRegression())
])
delta_model.fit(X, y)

# quick sanity metrics (in-sample on training deltas)
y_hat = delta_model.predict(X)
rmse = mean_squared_error(y, y_hat) ** 0.5
r2 = r2_score(y, y_hat)
print(f"Training window: {train_start_year}–{last_year} ({window_years} yrs)")
print(f"Delta model: polynomial degree={degree} | R^2={r2:.3f} | RMSE={rmse:.3f}")

# forecast years (include 2026–2036)
future_years = pd.DataFrame({'year': list(range(2026, 2036))})
future_years['t'] = future_years['year'] - t0

# predict deltas and accumulate from last observed level
future_years['delta_pred'] = delta_model.predict(future_years[['t']])

last_level = float(yearly.loc[yearly['year'] == last_year, 'avg_prime'].iloc[0])
future_years['avg_prime_pred'] = last_level + future_years['delta_pred'].cumsum()

# optional: keep predictions within a sane range
future_years['avg_prime_pred'] = future_years['avg_prime_pred'].clip(0, 25)

# plot: historical yearly averages + dotted forecast after 2025
plt.figure(figsize=(10, 6))
plt.plot(yearly['year'], yearly['avg_prime'], linewidth=2, label='Historical (yearly avg)')
plt.plot(future_years['year'], future_years['avg_prime_pred'], linestyle='--', linewidth=2, label='Forecast (2026–2036)')
plt.axvline(x=2025, linestyle=':', alpha=0.6)
plt.title('Average Prime Rate by Year with 10-Year Forecast (delta + polynomial)')
plt.xlabel('Year')
plt.ylabel('Average Prime Rate')
plt.grid(True)
plt.legend()
plt.show()

print("\nForecast (yearly avg prime):")
print(future_years[['year', 'avg_prime_pred']].to_string(index=False))