Fppy section 7.7 Adding Spline model forecast visualization for Boston marathon winning times

[presmerats]

Hi, I've been working on a version of the Spline model forecast for the Boston Marathon winning times in Python trying to replicate the R version of it. It uses the libraries presented in "An introduction to statistical learning” by James et al. from its python version, available on the ISLP package.

The code can be added and run just after the one in section 7.7, as it uses different elements defined there (boston_men, futr_df). Here it is:

import numpy as np, pandas as pd
from matplotlib.pyplot import subplots
import statsmodels.api as sm
from ISLP import load_data
from ISLP.models import (
    summarize ,
    poly ,
    ModelSpec as MS)
from pygam import (s as s_gam ,
    l as l_gam,
    f as f_gam,
    LinearGAM ,
    LogisticGAM)
from ISLP.transforms import (
    BSpline ,
    NaturalSpline)
from ISLP.models import bs, ns
from ISLP.pygam import (
    approx_lam ,
    degrees_of_freedom ,
    plot as plot_gam ,
    anova as anova_gam)

X_train = boston_men[["trend"]]
y_train = boston_men["y"].tolist()
X_futr =  futr_df[["trend"]]


y = boston_men["y"]
X_futr =  futr_df[["trend"]]

ns_trend = MS([bs('trend', df=3, degree=1)]).fit(boston_men)
X = ns_trend.transform(boston_men)
Xnew = ns_trend.transform(X_futr)

M = sm.OLS(y, X).fit()
#summarize(M)

# in sample
insample_preds =  M.get_prediction(X)
insample_preds_mean = insample_preds.predicted_mean

# forecast
preds = M.get_prediction(Xnew)
preds_mean = preds.predicted_mean
bands = preds.conf_int(alpha=0.05)

# add the spline  in sample predictions to the df to plot 
boston_men2 = boston_men.copy()
boston_men2['linear'] =  mf_linear.forecast_fitted_values()["LinearRegression"]
boston_men2['exponential'] = [np.exp(v) for v in mf_log.forecast_fitted_values()["LinearRegression"]]
boston_men2['spline'] = insample_preds_mean

# add the forecast part 
spline_fcast = futr_df.iloc[0:4,].copy() 
spline_fcast['spline'] = preds_mean[0:4,]
spline_fcast['spline-lo-95'] = bands[0:4,0]
spline_fcast['spline-hi-95'] = bands[0:4,1]
spline_fcast = spline_fcast[['unique_id','ds','spline','spline-lo-95','spline-hi-95']]
#spline_fcast

boston_forecast = forecasts_linear.merge(forecasts_log).merge(spline_fcast)

fig, ax = plt.subplots()

ax.set_prop_cycle(
    color=["black", "darkorange", "blue", "teal", "darkorange", "blue",   "teal"]
)
ax.plot(boston_men["ds"], boston_men["y"],label="y")

ax.plot(boston_men["ds"], boston_men2["linear"])
ax.plot(boston_men["ds"], boston_men2["exponential"])
ax.plot(boston_men["ds"], boston_men2["spline"])

ax.plot(boston_forecast["ds"], boston_forecast["linear"], label="linear")
ax.fill_between(
    boston_forecast["ds"],
    boston_forecast["linear-lo-95"],
    boston_forecast["linear-hi-95"], facecolor='darkorange', alpha=.1, label="linear_95")


ax.plot(boston_forecast["ds"], boston_forecast["exponential"], label="exponential")
ax.fill_between(
    boston_forecast["ds"],
    boston_forecast["exponential-lo-95"],
    boston_forecast["exponential-hi-95"], 
    facecolor='blue', alpha=.1, label="exponential_95")

ax.plot(boston_forecast["ds"], boston_forecast["spline"], label="spline")
ax.fill_between(
    boston_forecast["ds"],
    boston_forecast["spline-lo-95"],
    boston_forecast["spline-hi-95"], 
    facecolor='teal', alpha=.1, label="spline_level_95")

ax.set_title("Boston marathon winning times")
ax.set_xlabel("Year")
ax.set_ylabel("Minutes")

ax.legend(bbox_to_anchor=(1.04, 1), loc="upper left")

plt.show()

[robjhyndman]

We've now added piecewise linear trends to this example in https://otexts.com/fpppy/nbs/07-regression.html#example-boston-marathon-winning-times