Linear Changepoint Model with LASSO

I was experimenting with time series data specifically using fourier basis functions and fitting with a LASSO. I then decided to try just connecting the first point with another point and that point to the final point and see what that gave me. So the X matrix is built like this (definitely not efficient):

import quandl
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

data = quandl.get("BITSTAMP/USD")

y = data['Low']
y = y[-730:]

n_changepoints = 25
array_splits = np.array_split(np.array(y),n_changepoints)
y = np.array(y)
initial_point = y[0]
final_point = y[-1]
changepoints = np.zeros(shape=(len(y),n_changepoints))
for i in range(n_changepoints):
moving_point = array_splits[i][-1]
if i == 0:
len_splits = len(array_splits[i])
else:
len_splits = len(np.concatenate(array_splits[:i+1]).ravel())
slope = (moving_point - initial_point)/(len_splits)
if i != n_changepoints - 1:
reverse_slope = (final_point - moving_point)/(len(y) - len_splits)
reverse_slope = reverse_slope#*.5
changepoints[0:len_splits, i] = slope * (1+np.array(list(range(len_splits))))
changepoints[len_splits:, i] = changepoints[len_splits-1, i] + reverse_slope * (1+np.array(list(range((len(y) - len_splits)))))

plt.plot(y - y[0])
plt.plot(changepoints)
plt.show()


Then if I fit it with the lasso and some random alpha I get something like this:

from sklearn.linear_model import Lasso
lasso_model = Lasso(alpha = 10000).fit(changepoints, y)
fitted = lasso_model.predict(changepoints)
plt.plot(y)
plt.plot(fitted)
plt.show()


It now looks like a trend changepoint model reminiscent of the output from fbprophet.

And when I then combine the fourier terms with these and fit the lasso and do some tuning I end up with some fairly reasonable results that are in line with something like fbprophet. So my main question is, what is this idea of connecting the points (i.e. using the derivative) related to? From some research it looks like it is just some spline stuff but I am obviously new to this area and would love some direction!

It can be found here along with a quick example. Really going for something that at its core is feature engineering that is passed off to a simple fit for a linear regression, here it is a GLM with lasso regularization with the goal of resembling something from fbProphet. So is this approach of the 'weighted' matrix a consistent one or is it definitely going to be worse than using another set of linear piecewise functions?

Example output with those fourier basis functions:

Nice question.

I have more to say tomorrow, but for now the main point here is that your predictor matrix, changepoints is just a weighted matrix of piecewise linear basis functions, where the weights are set to fit the 25 equally-spaced points in your data.

It is common enough to use lasso to penalise regression splines, although other methods are more popular. If you haven't already read it, chapter 5 of Elements of Statistical Learning covers this.

Let $$X$$ be standard matrix of basis functions (middle panel), $$w$$ be the weights applied to each basis (right panel), and $$X_w$$ be the weighted basis matrix (left panel). If you were to use OLS instead of lasso, it wouldn't matter if you used the $$X$$ or $$X_w$$, since the regression weights could capture the difference: $$\hat y = X_w \times \beta = X \times \beta_w$$, where $$\beta_w = w \times \beta$$.

With lasso, however, you end up with a slightly different solution, since penalising $$\beta$$ isn't the same as penalising $$\beta_w$$. I need to do a little more poking around to be sure, but using $$X_w$$, which is already tailored to the data, should just lead to a tighter fit to the data for a given L1 penalty term (alpha), and as a result there exists some penalty, $$\alpha > 10,000$$, that gives exactly the same solution when using $$X$$.

Note: I've skipped over some details on centring the basis functions and subtracting the linear trend.

Code:

def plot_weights(weights, colour=True):
n = weights.shape[0]
if colour:
for i in range(n):
plt.plot([i, i], [0, weights[i]])
else:
plt.vlines(x=np.arange(n), ymin=0, ymax=weights)
plt.hlines(0, 0, n, linestyle='dashed')
plt.ylabel('Weight')

slope = changepoints[:, 0]
weighted_bases = changepoints[:, 1:] - slope.reshape(-1, 1)

ix_argmax = np.abs(weighted_bases).argmax(0)
basis_weights = np.array([weighted_bases[ix, i] for i, ix in enumerate(ix_argmax)])
raw_bases = weighted_bases / basis_weights

fig, axes = plt.subplots(1, 4, figsize=(14, 2))
plt.sca(axes[0])
plt.plot(changepoints)
plt.title('Changepoint Matrix =')

plt.sca(axes[1])
plt.plot(slope)
plt.title('Slope +')

plt.sca(axes[2])
plt.plot(bases)
plt.title('(Basis functions')

plt.sca(axes[3])
plot_weights(basis_weights)
plt.title('× Weights)')
plt.tight_layout()

• hey awesome answer. Definitely connects the dots for me getting from the piecewise linear to our setup. Like you were saying with OLS I get the same results whether I use the weighted matrix or the standard (except for some endpoint problems with the standard basis functions). And as you alluded to, LASSO gives me pretty uninspiring results with the standard basis functions so I guess a good follow up question would be: Is this weighted matrix an acceptable way to be fitting for this trend changepoint? I edited my question with a more complete implementation for what I'm going for! Commented Sep 22, 2020 at 0:01
• I haven't worked this out in full, but it's possible (but maybe not in sklearn) to apply a different penalty to each term in lasso. Rescaling the predictors and applying the same penalty to each is a different way of doing this. I suspect this means your approach works by adaptively rescaling (and thus changing the penalty applied to) each of the basis functions, which seems cool.
– Eoin
Commented Sep 22, 2020 at 9:08
• Have you tried using lasso with the standard basis functions, but decreasing the penalty term (increasing alpha)? Two other minor points: a) this kind of piecewise linear basis function isn't commonly used because it's difficult to fit linear models with such strongly correlated predictors, and b) if you want to know whether you can use your model for prediction, the only real answer is that you can if it performs well in cross-validation, and you can't otherwise!
– Eoin
Commented Sep 22, 2020 at 9:12
• yeah I moved away from sklearn and to statmodels which takes an array for the regularization. In terms of the standard basis functions it seems to always 'look' like a worse fit no matter the alpha. Another concern for the standard basis functions would be how I would use them to forecast forward since if I continue their trajectory they would all go negative, with the wonky basis stuff I just continue their trajectory then average them and pass that along. Looks like it does better overall and beats prophet on average for my 4k time series in a rolling CV. Interesting stuff! Commented Sep 22, 2020 at 13:59