Consider the following dataset (from FRED) on Annual Estimates of the Resident Population of California (units are thousands of persons) from 1900 to 2024.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as smcapop = pd.read_csv('CAPOP_29Sept2026.csv')
print(capop.head(10))
print(capop.tail(10))
tme = np.arange(1900, 2026)
plt.figure(figsize=(6, 4))
plt.plot(tme, capop['CAPOP'], label='California Population')
plt.xlabel('Year')
plt.ylabel('Population (in thousands)')
plt.title('California Annual Resident Population (in thousands)')
plt.legend()
plt.show() observation_date CAPOP
0 1900-01-01 1490.0
1 1901-01-01 1550.0
2 1902-01-01 1623.0
3 1903-01-01 1702.0
4 1904-01-01 1792.0
5 1905-01-01 1893.0
6 1906-01-01 1976.0
7 1907-01-01 2054.0
8 1908-01-01 2161.0
9 1909-01-01 2282.0
observation_date CAPOP
116 2016-01-01 39149.186
117 2017-01-01 39337.785
118 2018-01-01 39437.463
119 2019-01-01 39437.610
120 2020-01-01 39527.808
121 2021-01-01 39152.927
122 2022-01-01 39125.347
123 2023-01-01 39181.667
124 2024-01-01 39364.774
125 2025-01-01 39355.309
We work with the logarithms of the population data as this will lead to models with better interpretability.
y = np.log(capop['CAPOP'])
n = len(y)
plt.plot(tme, y)
plt.xlabel('Year')
plt.ylabel('Log(Population in thousands)')
plt.title('Logarithm of California Annual Resident Population (in thousands)')
plt.show()
The first model we fit is linear regression:
with being i.i.d . here can be interpreted in terms of the overall growth rate.
x = np.arange(1, n+1)
X = np.column_stack((np.ones(n), x))
linreg = sm.OLS(y, X).fit()
print(linreg.summary())
plt.plot(tme, y, label='Log(Population in thousands)')
plt.plot(tme, linreg.fittedvalues, color='red', label='Fitted least squares line')
plt.xlabel('Year')
plt.ylabel('Log(Population in thousands)')
plt.title('Log Data with fitted least squares line')
plt.show() OLS Regression Results
==============================================================================
Dep. Variable: CAPOP R-squared: 0.948
Model: OLS Adj. R-squared: 0.948
Method: Least Squares F-statistic: 2278.
Date: Tue, 29 Sep 2026 Prob (F-statistic): 1.16e-81
Time: 19:00:09 Log-Likelihood: 8.2264
No. Observations: 126 AIC: -12.45
Df Residuals: 124 BIC: -6.780
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 7.7387 0.041 188.952 0.000 7.658 7.820
x1 0.0267 0.001 47.726 0.000 0.026 0.028
==============================================================================
Omnibus: 11.681 Durbin-Watson: 0.006
Prob(Omnibus): 0.003 Jarque-Bera (JB): 8.717
Skew: -0.529 Prob(JB): 0.0128
Kurtosis: 2.265 Cond. No. 147.
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
The fitted slope coefficient here is 0.0267. The interpretation is that the population increases by 2.67% each year.
Broken-stick Regression or Change of Slope Model¶
The simple linear regression model obviously does not provide a good fit to the data. For an improved model, let us consider:
This model uses for the slope before , and for the slope after . We estimated by the following:
For each fixed value of , calculate RSS
Use the value of with the smallest RSS as the estimate
def rss(c):
x = np.arange(1, n+1)
x_c = np.maximum(0, x - c)
X = np.column_stack((np.ones(n), x, x_c))
md = sm.OLS(y, X).fit()
rss = np.sum(md.resid ** 2)
return rssWe compute for each in a grid of values of as follows.
num_c_vals = 1000 #this is the number of different values of c we will try
allcvals = np.linspace(1, n, num_c_vals)
rssvals = np.array([rss(c) for c in allcvals])
plt.plot(allcvals, rssvals)
plt.show()
The estimate is obtained by minimizing as follows.
c_hat = allcvals[np.argmin(rssvals)]
print(c_hat)
print(c_hat - 1 + tme[0]) #this is the estimated year when the slope changes
rss_smallest = np.min(rssvals)
print(rss_smallest)
66.56556556556558
1965.5655655655655
0.3677408605674205
The fitted values will now look much better than before.
c = c_hat
x = np.arange(1, n+1)
x_c = np.maximum(0, x - c)
X = np.column_stack((np.ones(n), x, x_c))
md = sm.OLS(y, X).fit()
plt.plot(tme, y)
plt.plot(tme, md.fittedvalues, color = 'red')
plt.axvline(c_hat - 1 + tme[0], color='green', linestyle='--')
plt.show()
Point estimates of the other parameters are obtained as follows.
#Estimates of other parameters:
c = c_hat
x = np.arange(1, n+1)
x_c = np.maximum(0, x - c)
X = np.column_stack((np.ones(n), x, x_c))
md = sm.OLS(y, X).fit()
print(md.params) #this gives estimates of beta_0, beta_1, beta_2 const 7.369872
x1 0.037974
x2 -0.024299
dtype: float64
The estimate of is 0.038 and the estimate of is -0.024. This means that the growth rate before 1965 was 3.8% while the growth rate after 1965 is %.
Models with more changes of Slope¶
We can seek to obtain further improved models by using more points of change of slope. The following model has two points where the slope changes:
The unknown parameters now are as well as (also ).
For parameter estimation, we can define the RSS function:
and then we can do some sort of grid optimization for and . The code for this is given below.
def rss(c): #c is now a vector of change of slope points
n = len(y)
x = np.arange(1, n+1)
X = np.column_stack((np.ones(n), x))
if np.isscalar(c):
c = [c]
for j in range(len(c)):
x_c = np.maximum(0, x - c[j])
X = np.column_stack((X, x_c))
md = sm.OLS(y, X).fit()
rss = np.sum(md.resid ** 2)
return rssnum_c_vals = 200 #this is the number of different values of c we will try
c1_gr = np.linspace(1, n, num_c_vals)
c2_gr = np.linspace(1, n, num_c_vals)
C1, C2 = np.meshgrid(c1_gr, c2_gr)
g = pd.DataFrame({'c1': C1.ravel(), 'c2': C2.ravel()})
g['rss'] = g.apply(lambda row: rss([row['c1'], row['c2']]), axis=1)
min_row = g.loc[g['rss'].idxmin()]
print(min_row)
c_opt = min_row[:-1]
rss_opt = min_row.iloc[-1]
print(c_opt - 1 + tme[0]) #these are the estimated years when the slope changes
print(c_hat) #this was the best c in a single change of slope model
print(rss_smallest) #this was the smallest RSS in a single change of slope model
print(rss_opt)c1 102.130653
c2 62.557789
rss 0.199761
Name: 19761, dtype: float64
c1 2001.130653
c2 1961.557789
Name: 19761, dtype: float64
66.56556556556558
0.3677408605674205
0.19976107958305042
c = np.array(c_opt)
x = np.arange(1, n+1)
x_c1 = np.maximum(0, x - c[0])
x_c2 = np.maximum(0, x - c[1])
X = np.column_stack((np.ones(n), x, x_c1, x_c2))
md_2 = sm.OLS(y, X).fit()
print(md_2.summary())
#plt.plot(tme, y)
#plt.plot(tme, md.fittedvalues, color = 'orange', label='1 change of slope')
plt.plot(tme, md_2.fittedvalues, color = 'red')
plt.axvline(c[0] - 1 + tme[0], color='green', linestyle='--')
plt.axvline(c[1] - 1 + tme[0], color='green', linestyle='--')
plt.show() OLS Regression Results
==============================================================================
Dep. Variable: CAPOP R-squared: 0.998
Model: OLS Adj. R-squared: 0.998
Method: Least Squares F-statistic: 2.549e+04
Date: Tue, 29 Sep 2026 Prob (F-statistic): 1.92e-170
Time: 19:03:40 Log-Likelihood: 227.37
No. Observations: 126 AIC: -446.7
Df Residuals: 122 BIC: -435.4
Df Model: 3
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 7.3659 0.010 749.118 0.000 7.346 7.385
x1 0.0382 0.000 160.182 0.000 0.038 0.039
x2 -0.0126 0.001 -11.621 0.000 -0.015 -0.010
x3 -0.0200 0.001 -37.521 0.000 -0.021 -0.019
==============================================================================
Omnibus: 14.458 Durbin-Watson: 0.084
Prob(Omnibus): 0.001 Jarque-Bera (JB): 21.254
Skew: 0.585 Prob(JB): 2.43e-05
Kurtosis: 4.638 Cond. No. 210.
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
The interpretation of this model is that: before 1961, the growth rate is 3.82%. Between 1961 and 2000, the growth fell to 3.82 - 1.26 = 2.56%. After 2000, the growth rate further dropped to 2.56 - 2 = 0.56 %.
There is a different way of fitting this model that can be computationally more efficient. This involves jointly minimizing least squares over all the parameters. Specifically, we aim to optimize:
over all the parameters simultaneously. This optimization is nonlinear and nontrivial but we can use black box optimization libraries. For example, we can use the pytorch library to solve this optimization.
First import the torch libraries.
import torch
import torch.nn as nn
import torch.optim as optimThe following is the code for representing the broken stick regression (or change of slope) model in PyTorch. PyTorch organizes models as classes.
class BrokenStickRegression(nn.Module):
def __init__(self, knots_init, beta_init):
super().__init__()
self.num_knots = len(knots_init)
self.beta = nn.Parameter(torch.tensor(beta_init, dtype=torch.float32))
self.knots = nn.Parameter(torch.tensor(knots_init, dtype=torch.float32))
def forward(self, x):
knots_sorted,_ = torch.sort(self.knots)
out = self.beta[0] + self.beta[1] * x
for j in range(self.num_knots):
out = out + self.beta[j+2] * torch.relu(x - knots_sorted[j])
return outThe class BrokenStickRegression inherits from nn.Module, which is the base class for all models in PyTorch (linear regression, neural networks, etc.). By subclassing nn.Module, we get access to a lot of PyTorch functionality automatically: for example, it keeps track of parameters, makes them trainable, and allows us to use optimizers like gradient descent.
The term knots refers to .
In the constructor (__init__), the code takes two inputs: a list of initial knot positions (knots_init) and initial coefficients (beta_init). These are wrapped in nn.Parameter, which tells PyTorch: “this is a parameter we want to optimize.” So, during training, PyTorch will update both the regression coefficients and the knot locations.
The forward method defines the model equation (how does the model convert input to output ). When you pass in an input x, it:
Sorts the knot locations (so they are in order).
Starts with the base linear trend,
For each knot, adds a term
The effect is exactly the broken-stick regression function: linear to start, with slope changes at each knot. The nice thing about using nn.Module is that all of this can be trained with PyTorch’s optimization tools.
The raw data here is given by with covariate given by time . We do not apply the model to this raw data but we scale the raw data first. The algorithm for parameter fitting will work better with scaled data as opposed to raw data.
The code below first standardizes the raw input arrays y_raw and x_raw by subtracting their means and dividing by their standard deviations, resulting in y_scaled and x_scaled, respectively. This scaling step ensures that both inputs have mean 0 and standard deviation 1, which helps neural network training converge faster and more reliably. After scaling, the arrays are converted into PyTorch tensors (y_torch and x_torch) with dtype=torch.float32, and an additional singleton dimension is added using unsqueeze(1) to ensure each data point is treated as a one-dimensional feature vector. Conversion to PyTorch tensors is necessary because PyTorch models and optimization routines (like gradient computation and parameter updates) operate on tensors, not on NumPy arrays.
y_raw = y
n = len(y_raw)
x_raw = np.arange(1, n+1)
y_scaled = (y_raw - np.mean(y_raw)) / np.std(y_raw)
x_scaled = (x_raw - np.mean(x_raw)) / np.std(x_raw)
y_torch = torch.tensor(y_scaled, dtype=torch.float32).unsqueeze(1)
x_torch = torch.tensor(x_scaled, dtype=torch.float32).unsqueeze(1)For using PyTorch, the first step is to obtain suitable initial estimates for the parameters. These are deduced as follows. We take to be quantiles of the scaled covariate at equal levels. Then we fit the model with fixed at these initial values, and obtain the initial values of the coefficients.
k = 2 #this is the number of knots (changes of slope)
quantile_levels = np.linspace(1/(k+1), k/(k+1), k)
knots_init = np.quantile(x_scaled, quantile_levels)
n = len(y_scaled)
X = np.column_stack((np.ones(n), x_scaled))
for j in range(k):
x_c = np.maximum(0, x_scaled - knots_init[j])
X = np.column_stack((X, x_c))
md_init = sm.OLS(y_scaled, X).fit()
beta_init = md_init.params.values
print(knots_init)
print(knots_init * np.std(x_raw) + np.mean(x_raw)) #these are the initial knots in original scale
print(beta_init)[-0.57278616 0.57278616]
[42.66666667 84.33333333]
[ 0.48135629 1.48815313 -0.45562361 -0.68045908]
The next code line md_nn = BrokenStickRegression(knots_init=knots_init, beta_init=beta_init) creates an instance of the custom neural network model BrokenStickRegression. It initializes the learnable parameters of the model: the knot locations are set to the values provided in knots_init, and the coefficients are set to beta_init. This prepares the model for training by defining its initial piecewise linear structure.
md_nn = BrokenStickRegression(knots_init, beta_init)
#This code creates an instance of the BrokenStickRegression class with initial knots and beta values.The next block of code sets up and runs the training loop for the model BrokenStickRegression. The Adam optimizer is initialized with the model’s parameters and a learning rate of 0.01, and the loss function is set to mean squared error (MSELoss). For 10,000 epochs, the code repeatedly performs one training step: it clears previous gradients with optimizer.zero_grad(), computes predictions y_pred by passing x_torch through the model, evaluates the loss between predictions and true values, backpropagates the loss with loss.backward(), and updates the model parameters using optimizer.step(). Every 100 epochs, the current epoch and loss value are printed to monitor training progress. Running the code multiple times may be necessary to ensure good convergence, especially because this is a non-convex optimization problem.
optimizer = optim.Adam(md_nn.parameters(), lr=0.01)
loss_fn = nn.MSELoss()
for epoch in range(10000):
optimizer.zero_grad()
y_pred = md_nn(x_torch)
loss = loss_fn(y_pred, y_torch)
loss.backward()
optimizer.step()
if epoch % 500 == 0:
print(f'Epoch {epoch}, Loss: {loss.item()}')Epoch 0, Loss: 0.0015925031621009111
Epoch 500, Loss: 0.0015920239966362715
Epoch 1000, Loss: 0.001592024345882237
Epoch 1500, Loss: 0.0015920238802209496
Epoch 2000, Loss: 0.0015925345942378044
Epoch 2500, Loss: 0.0015925955958664417
Epoch 3000, Loss: 0.0015921215526759624
Epoch 3500, Loss: 0.0015965248458087444
Epoch 4000, Loss: 0.0015920241130515933
Epoch 4500, Loss: 0.001592024927958846
Epoch 5000, Loss: 0.0015920244622975588
Epoch 5500, Loss: 0.0015920245787128806
Epoch 6000, Loss: 0.0015920244622975588
Epoch 6500, Loss: 0.0015920235309749842
Epoch 7000, Loss: 0.0015920302830636501
Epoch 7500, Loss: 0.0015933793038129807
Epoch 8000, Loss: 0.0016028147656470537
Epoch 8500, Loss: 0.0015920246951282024
Epoch 9000, Loss: 0.0015920241130515933
Epoch 9500, Loss: 0.0015920239966362715
The next code prints out the current loss value, the estimated model coefficients (beta), and the estimated knot locations after training. The detach().numpy() calls are used to move each tensor from the computation graph to a regular NumPy array so they can be easily printed or further processed without tracking gradients. Finally, it prints rss_opt, which represents the smallest residual sum of squares (RSS) achieved with two knots using our previous grid-based search method. This provides a measure of the quality of the PyTorch parameter estimates.
print(loss.detach().numpy())
print(md_nn.beta.detach().numpy())
print(md_nn.knots.detach().numpy())
print(rss_opt)0.0015920238
[ 0.35588792 1.3921012 -0.72420883 -0.45621735]
[-0.03157811 1.0487218 ]
0.19976107958305042
rss_opt seems to be on a different scale to the loss. This is because of the change of scale when we went from y_raw and x_raw to y_scaled and x_scaled, as well as because the loss here is Mean-Squared Error (not just Squared Error). We can bring the loss to the same scale as RSS as follows.
#converting back to original scale
print(loss.detach().numpy() * (np.std(y_raw)**2) * n ) #this rescales the loss to the rss scale
print(md_nn.beta.detach().numpy())
print(md_nn.knots.detach().numpy() * np.std(x_raw) + np.mean(x_raw)) #these are the knots in original scale
print(rss_opt)0.1996371551927383
[ 0.35588792 1.3921012 -0.72420883 -0.45621735]
[ 62.35144338 101.64402653]
0.19976107958305042
The loss achieved by PyTorch seems slightly better than the RSS we got from our grid search method. This suggests that PyTorch is working well here. Next let us look at the estimates of obtained from PyTorch (converted to the original scale). We compare them with the estimates from our previous grid-search method.
#drop intercept
print(md_nn.beta.detach().numpy()[1:] * (np.std(y_raw)/np.std(x_raw)))
print(md_2.params.values[1:])[ 0.03818257 -0.01986361 -0.01251314]
[ 0.03815355 -0.01255336 -0.019984 ]
This is essentially the same model that we found before (in the grid search method, and seem to be swapped though).
This PyTorch method scales easily to higher values of the number of knots i.e., to fit models of the form:
On the other hand, the grid-search method will take very long to run. Below we fit the model with .
k = 6 #this is the number of knots (changes of slope)
quantile_levels = np.linspace(1/(k+1), k/(k+1), k)
knots_init = np.quantile(x_scaled, quantile_levels)
n = len(y_scaled)
X = np.column_stack((np.ones(n), x_scaled))
for j in range(k):
x_c = np.maximum(0, x_scaled - knots_init[j])
X = np.column_stack((X, x_c))
md_init = sm.OLS(y_scaled, X).fit()
beta_init = md_init.params.values
print(knots_init)
print(knots_init * np.std(x_raw) + np.mean(x_raw)) #these are the initial knots in original scale
print(beta_init)[-1.22739892 -0.73643935 -0.24547978 0.24547978 0.73643935 1.22739892]
[ 18.85714286 36.71428571 54.57142857 72.42857143 90.28571429
108.14285714]
[ 0.86633443 1.74925369 -0.50458691 0.16019777 -0.39630201 -0.34868515
-0.14694489 -0.36487881]
md_nn_k = BrokenStickRegression(knots_init, beta_init)
#This code creates an instance of the BrokenStickRegression class with initial knots and beta values.optimizer = optim.Adam(md_nn_k.parameters(), lr=0.001)
loss_fn = nn.MSELoss()
for epoch in range(10000):
optimizer.zero_grad()
y_pred = md_nn_k(x_torch)
loss = loss_fn(y_pred, y_torch)
loss.backward()
optimizer.step()
if epoch % 500 == 0:
print(f'Epoch {epoch}, Loss: {loss.item()}')Epoch 0, Loss: 0.0011327910469844937
Epoch 500, Loss: 0.0004229986807331443
Epoch 1000, Loss: 0.00035374544677324593
Epoch 1500, Loss: 0.0003272918111179024
Epoch 2000, Loss: 0.00031581372604705393
Epoch 2500, Loss: 0.00031160947401076555
Epoch 3000, Loss: 0.0003096174623351544
Epoch 3500, Loss: 0.00030928454361855984
Epoch 4000, Loss: 0.00030810970929451287
Epoch 4500, Loss: 0.0003078346198890358
Epoch 5000, Loss: 0.0003078820009250194
Epoch 5500, Loss: 0.0003075877029914409
Epoch 6000, Loss: 0.00030752806924283504
Epoch 6500, Loss: 0.00030754992621950805
Epoch 7000, Loss: 0.0003074596752412617
Epoch 7500, Loss: 0.0003074312990065664
Epoch 8000, Loss: 0.0003074094420298934
Epoch 8500, Loss: 0.000307392911054194
Epoch 9000, Loss: 0.00030737632187083364
Epoch 9500, Loss: 0.0003073618281632662
The estimates of the knots are given below.
c_est = md_nn_k.knots.detach().numpy() * np.std(x_raw) + np.mean(x_raw)
print(c_est)[ 29.95563226 38.27802863 59.08876163 63.65471383 94.58710829
113.89943503]
The rss value and the fitted values are computed as below.
c = c_est
x = np.arange(1, n+1)
X = np.column_stack((np.ones(n), x))
if np.isscalar(c):
c = [c]
for j in range(len(c)):
x_c = np.maximum(0, x - c[j])
X = np.column_stack((X, x_c))
md_k = sm.OLS(y, X).fit()
rss_k = np.sum(md_k.resid ** 2)
print(rss_k)0.038539956283362334
plt.plot(tme, y, color= 'black')
#plt.plot(tme, md_2.fittedvalues, color = 'blue', label='2 changes of slope')
plt.plot(tme, md_k.fittedvalues, color = 'red', label = 'k changes of slope')
#plt.axvline(c[0] - 1 + tme[0], color='green', linestyle='--')
#plt.axvline(c[1] - 1 + tme[0], color='green', linestyle='--')
plt.show() 
The fitted values from the model seem to be quite close to the fitted values for the model. Go back and run the method for . This should lead to overfitting and the fitted values will look a lot like the data.
Importance of Initialization¶
This method based on PyTorch requires careful initialization. This can be seen from the following simple simulated example.
rng = np.random.default_rng(1)
n = 400
x = np.arange(1, n + 1)
y_true = 0.6 * (
np.maximum(0, x - 40)
- 2 * np.maximum(0, x - 50)
+ np.maximum(0, x - 60)
)
sig_noise = 0.4
y = y_true + sig_noise * rng.standard_normal(n)
n = len(y)
tme = np.arange(1, n + 1)
plt.plot(tme, y)
plt.xlabel('Time')
plt.ylabel('Data')
plt.title('A simulated dataset')
plt.show()
y_raw = y
n = len(y_raw)
x_raw = np.arange(1, n + 1)
y_scaled = (y_raw - np.mean(y_raw)) / np.std(y_raw)
x_scaled = (x_raw - np.mean(x_raw)) / np.std(x_raw)
y_torch = torch.tensor(y_scaled, dtype=torch.float32).unsqueeze(1)
x_torch = torch.tensor(x_scaled, dtype=torch.float32).unsqueeze(1)k = 3
quantile_levels = np.linspace(1 / (k + 1), k / (k + 1), k)
knots_init = np.quantile(x_scaled, quantile_levels)
#Better initialization:
c_init = np.array([38., 55., 62.])
knots_init = (c_init - np.mean(x_raw)) / np.std(x_raw)
n = len(y_scaled)
X = np.column_stack((np.ones(n), x_scaled))
for j in range(k):
x_c = np.maximum(0, x_scaled - knots_init[j])
X = np.column_stack((X, x_c))
md_init = sm.OLS(y_scaled, X).fit()
beta_init = np.asarray(md_init.params)
print(knots_init)
print(knots_init * np.std(x_raw) + np.mean(x_raw))
print(beta_init)
[-1.40729568 -1.2600709 -1.19944893]
[38. 55. 62.]
[ 1.67032957 1.10106531 36.13272113 -133.93530714 96.74733646]
md_nn_k = BrokenStickRegression(knots_init, beta_init)
optimizer = optim.Adam(md_nn_k.parameters(), lr=0.001)
loss_fn = nn.MSELoss()
for epoch in range(100000):
optimizer.zero_grad()
y_pred = md_nn_k(x_torch)
loss = loss_fn(y_pred, y_torch)
loss.backward()
optimizer.step()
if epoch % 500 == 0:
print(f'Epoch {epoch}, Loss: {loss.item()}')
Epoch 0, Loss: 0.3337884843349457
Epoch 500, Loss: 0.2502412497997284
Epoch 1000, Loss: 0.2500342130661011
Epoch 1500, Loss: 0.25066688656806946
Epoch 2000, Loss: 0.2502841055393219
Epoch 2500, Loss: 0.24942947924137115
Epoch 3000, Loss: 0.24923084676265717
Epoch 3500, Loss: 0.24903494119644165
Epoch 4000, Loss: 0.24883756041526794
Epoch 4500, Loss: 0.24873390793800354
Epoch 5000, Loss: 0.24844516813755035
Epoch 5500, Loss: 0.24862927198410034
Epoch 6000, Loss: 0.2480580359697342
Epoch 6500, Loss: 0.2478664070367813
Epoch 7000, Loss: 0.247695654630661
Epoch 7500, Loss: 0.24748843908309937
Epoch 8000, Loss: 0.2473001331090927
Epoch 8500, Loss: 0.24711519479751587
Epoch 9000, Loss: 0.24692924320697784
Epoch 9500, Loss: 0.2467452883720398
Epoch 10000, Loss: 0.24655698239803314
Epoch 10500, Loss: 0.24696698784828186
Epoch 11000, Loss: 0.2461887001991272
Epoch 11500, Loss: 0.2460019290447235
Epoch 12000, Loss: 0.24582156538963318
Epoch 12500, Loss: 0.24564018845558167
Epoch 13000, Loss: 0.24545890092849731
Epoch 13500, Loss: 0.2452784925699234
Epoch 14000, Loss: 0.24509967863559723
Epoch 14500, Loss: 0.24492177367210388
Epoch 15000, Loss: 0.2447425127029419
Epoch 15500, Loss: 0.24457122385501862
Epoch 16000, Loss: 0.24438896775245667
Epoch 16500, Loss: 0.2442915290594101
Epoch 17000, Loss: 0.24403856694698334
Epoch 17500, Loss: 0.24386468529701233
Epoch 18000, Loss: 0.2436923235654831
Epoch 18500, Loss: 0.24351993203163147
Epoch 19000, Loss: 0.24335099756717682
Epoch 19500, Loss: 0.2431982010602951
Epoch 20000, Loss: 0.24304541945457458
Epoch 20500, Loss: 0.24292001128196716
Epoch 21000, Loss: 0.24268215894699097
Epoch 21500, Loss: 0.2425023317337036
Epoch 22000, Loss: 0.2423669993877411
Epoch 22500, Loss: 0.24217264354228973
Epoch 23000, Loss: 0.2420063465833664
Epoch 23500, Loss: 0.24184300005435944
Epoch 24000, Loss: 0.24170631170272827
Epoch 24500, Loss: 0.24181854724884033
Epoch 25000, Loss: 0.2413560152053833
Epoch 25500, Loss: 0.24119283258914948
Epoch 26000, Loss: 0.24103349447250366
Epoch 26500, Loss: 0.24087320268154144
Epoch 27000, Loss: 0.24071462452411652
Epoch 27500, Loss: 0.240559920668602
Epoch 28000, Loss: 0.24042747914791107
Epoch 28500, Loss: 0.24024488031864166
Epoch 29000, Loss: 0.24018265306949615
Epoch 29500, Loss: 0.23998668789863586
Epoch 30000, Loss: 0.23980534076690674
Epoch 30500, Loss: 0.23967331647872925
Epoch 31000, Loss: 0.2394835501909256
Epoch 31500, Loss: 0.23933224380016327
Epoch 32000, Loss: 0.23929759860038757
Epoch 32500, Loss: 0.23903590440750122
Epoch 33000, Loss: 0.23888865113258362
Epoch 33500, Loss: 0.23879866302013397
Epoch 34000, Loss: 0.2385980188846588
Epoch 34500, Loss: 0.23845332860946655
Epoch 35000, Loss: 0.23831051588058472
Epoch 35500, Loss: 0.23816941678524017
Epoch 36000, Loss: 0.23826450109481812
Epoch 36500, Loss: 0.23788872361183167
Epoch 37000, Loss: 0.23775054514408112
Epoch 37500, Loss: 0.23760537803173065
Epoch 38000, Loss: 0.2374647706747055
Epoch 38500, Loss: 0.23732320964336395
Epoch 39000, Loss: 0.23718249797821045
Epoch 39500, Loss: 0.23704391717910767
Epoch 40000, Loss: 0.23690450191497803
Epoch 40500, Loss: 0.23676681518554688
Epoch 41000, Loss: 0.2366296350955963
Epoch 41500, Loss: 0.23649120330810547
Epoch 42000, Loss: 0.23635481297969818
Epoch 42500, Loss: 0.23639366030693054
Epoch 43000, Loss: 0.23608079552650452
Epoch 43500, Loss: 0.23594680428504944
Epoch 44000, Loss: 0.2360156625509262
Epoch 44500, Loss: 0.23567607998847961
Epoch 45000, Loss: 0.2355482429265976
Epoch 45500, Loss: 0.2354096919298172
Epoch 46000, Loss: 0.23538152873516083
Epoch 46500, Loss: 0.23514504730701447
Epoch 47000, Loss: 0.23501411080360413
Epoch 47500, Loss: 0.2348833978176117
Epoch 48000, Loss: 0.23475231230258942
Epoch 48500, Loss: 0.23462317883968353
Epoch 49000, Loss: 0.23449382185935974
Epoch 49500, Loss: 0.23437954485416412
Epoch 50000, Loss: 0.23423728346824646
Epoch 50500, Loss: 0.23410965502262115
Epoch 51000, Loss: 0.23398250341415405
Epoch 51500, Loss: 0.23385700583457947
Epoch 52000, Loss: 0.2338639497756958
Epoch 52500, Loss: 0.2336178570985794
Epoch 53000, Loss: 0.23352116346359253
Epoch 53500, Loss: 0.23336033523082733
Epoch 54000, Loss: 0.23323959112167358
Epoch 54500, Loss: 0.23316413164138794
Epoch 55000, Loss: 0.233122780919075
Epoch 55500, Loss: 0.23290196061134338
Epoch 56000, Loss: 0.23281234502792358
Epoch 56500, Loss: 0.2327202409505844
Epoch 57000, Loss: 0.2325470894575119
Epoch 57500, Loss: 0.23245932161808014
Epoch 58000, Loss: 0.23233044147491455
Epoch 58500, Loss: 0.23218590021133423
Epoch 59000, Loss: 0.2320631742477417
Epoch 59500, Loss: 0.2319961041212082
Epoch 60000, Loss: 0.23184096813201904
Epoch 60500, Loss: 0.2317439615726471
Epoch 61000, Loss: 0.23161433637142181
Epoch 61500, Loss: 0.23151201009750366
Epoch 62000, Loss: 0.2314615249633789
Epoch 62500, Loss: 0.23129260540008545
Epoch 63000, Loss: 0.2311883121728897
Epoch 63500, Loss: 0.23107726871967316
Epoch 64000, Loss: 0.23098143935203552
Epoch 64500, Loss: 0.23090913891792297
Epoch 65000, Loss: 0.2309906929731369
Epoch 65500, Loss: 0.23065637052059174
Epoch 66000, Loss: 0.23057611286640167
Epoch 66500, Loss: 0.23045001924037933
Epoch 67000, Loss: 0.23073549568653107
Epoch 67500, Loss: 0.23025599122047424
Epoch 68000, Loss: 0.23014375567436218
Epoch 68500, Loss: 0.23004059493541718
Epoch 69000, Loss: 0.22994939982891083
Epoch 69500, Loss: 0.2300117462873459
Epoch 70000, Loss: 0.2301529049873352
Epoch 70500, Loss: 0.22963611781597137
Epoch 71000, Loss: 0.22953739762306213
Epoch 71500, Loss: 0.2294856309890747
Epoch 72000, Loss: 0.22933930158615112
Epoch 72500, Loss: 0.2292517125606537
Epoch 73000, Loss: 0.22930678725242615
Epoch 73500, Loss: 0.22905199229717255
Epoch 74000, Loss: 0.22894662618637085
Epoch 74500, Loss: 0.22884993255138397
Epoch 75000, Loss: 0.22875282168388367
Epoch 75500, Loss: 0.22866447269916534
Epoch 76000, Loss: 0.228560671210289
Epoch 76500, Loss: 0.22846391797065735
Epoch 77000, Loss: 0.22836750745773315
Epoch 77500, Loss: 0.2282736748456955
Epoch 78000, Loss: 0.2281828671693802
Epoch 78500, Loss: 0.2281033843755722
Epoch 79000, Loss: 0.2279931604862213
Epoch 79500, Loss: 0.22789736092090607
Epoch 80000, Loss: 0.22780433297157288
Epoch 80500, Loss: 0.22771085798740387
Epoch 81000, Loss: 0.22761978209018707
Epoch 81500, Loss: 0.22752809524536133
Epoch 82000, Loss: 0.22722014784812927
Epoch 82500, Loss: 0.22711309790611267
Epoch 83000, Loss: 0.226999431848526
Epoch 83500, Loss: 0.22693751752376556
Epoch 84000, Loss: 0.2267604023218155
Epoch 84500, Loss: 0.22683244943618774
Epoch 85000, Loss: 0.22653193771839142
Epoch 85500, Loss: 0.2264166623353958
Epoch 86000, Loss: 0.22630451619625092
Epoch 86500, Loss: 0.22619356215000153
Epoch 87000, Loss: 0.22609089314937592
Epoch 87500, Loss: 0.226101353764534
Epoch 88000, Loss: 0.22589576244354248
Epoch 88500, Loss: 0.2257642298936844
Epoch 89000, Loss: 0.2256300002336502
Epoch 89500, Loss: 0.225519597530365
Epoch 90000, Loss: 0.22541120648384094
Epoch 90500, Loss: 0.2253202497959137
Epoch 91000, Loss: 0.22567005455493927
Epoch 91500, Loss: 0.22508090734481812
Epoch 92000, Loss: 0.22496801614761353
Epoch 92500, Loss: 0.2248597890138626
Epoch 93000, Loss: 0.22475849092006683
Epoch 93500, Loss: 0.2247171401977539
Epoch 94000, Loss: 0.22458332777023315
Epoch 94500, Loss: 0.22442626953125
Epoch 95000, Loss: 0.2243170142173767
Epoch 95500, Loss: 0.22421088814735413
Epoch 96000, Loss: 0.22410598397254944
Epoch 96500, Loss: 0.2239968329668045
Epoch 97000, Loss: 0.2238941639661789
Epoch 97500, Loss: 0.22378449141979218
Epoch 98000, Loss: 0.22367841005325317
Epoch 98500, Loss: 0.22357258200645447
Epoch 99000, Loss: 0.223481222987175
Epoch 99500, Loss: 0.22349247336387634
c_est = md_nn_k.knots.detach().numpy() * np.std(x_raw) + np.mean(x_raw)
print(c_est)[37.2653834 51.78968288 59.48162598]
c = c_est
x = np.arange(1, n + 1)
X = np.column_stack((np.ones(n), x))
if np.isscalar(c):
c = [c]
for j in range(len(c)):
x_c = np.maximum(0, x - c[j])
X = np.column_stack((X, x_c))
md_k = sm.OLS(y, X).fit()
rss_k = np.sum(md_k.resid ** 2)
print(rss_k)
plt.plot(x, y, color='lightgray', label='Data')
plt.plot(x, y_true, color='blue', label='True function')
plt.plot(x, md_k.fittedvalues, color='red', label='Fitted function')
plt.legend()
plt.show()60.99424395355864
High-dimensional Linear Regression and Regularization¶
We now consider the model:
We can write this as
where is constructed as below.
y = np.log(capop['CAPOP'])
n = len(y)
tme = np.arange(1, n+1)
plt.plot(tme, y)
plt.xlabel('Year')
plt.ylabel('Log(Population in thousands)')
plt.title('Logarithm of California Annual Resident Population (in thousands)')
plt.show()
x = np.arange(1, n+1)
Xfull = np.column_stack([np.ones(n), x-1])
for i in range(n-2):
c = i+2
xc = ((x > c).astype(float))*(x-c)
Xfull = np.column_stack([Xfull, xc])
print(Xfull)[[ 1. 0. -0. ... -0. -0. -0.]
[ 1. 1. 0. ... -0. -0. -0.]
[ 1. 2. 1. ... -0. -0. -0.]
...
[ 1. 123. 122. ... 1. 0. -0.]
[ 1. 124. 123. ... 2. 1. 0.]
[ 1. 125. 124. ... 3. 2. 1.]]
If we fit this model without any regularization, we get the following.
mdfull = sm.OLS(y, Xfull).fit()
print(mdfull.summary()) OLS Regression Results
==============================================================================
Dep. Variable: CAPOP R-squared: 1.000
Model: OLS Adj. R-squared: nan
Method: Least Squares F-statistic: nan
Date: Tue, 29 Sep 2026 Prob (F-statistic): nan
Time: 19:48:05 Log-Likelihood: 3207.8
No. Observations: 126 AIC: -6164.
Df Residuals: 0 BIC: -5806.
Df Model: 125
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 7.3065 inf 0 nan nan nan
x1 0.0395 inf 0 nan nan nan
x2 0.0065 inf 0 nan nan nan
x3 0.0015 inf 0 nan nan nan
x4 0.0040 inf 0 nan nan nan
x5 0.0033 inf 0 nan nan nan
x6 -0.0119 inf -0 nan nan nan
x7 -0.0042 inf -0 nan nan nan
x8 0.0121 inf 0 nan nan nan
x9 0.0037 inf 0 nan nan nan
x10 -0.0016 inf -0 nan nan nan
x11 -0.0011 inf -0 nan nan nan
x12 -0.0003 inf -0 nan nan nan
x13 0.0007 inf 0 nan nan nan
x14 -0.0094 inf -0 nan nan nan
x15 -0.0179 inf -0 nan nan nan
x16 -0.0042 inf -0 nan nan nan
x17 0.0113 inf 0 nan nan nan
x18 -0.0038 inf -0 nan nan nan
x19 -0.0050 inf -0 nan nan nan
x20 0.0391 inf 0 nan nan nan
x21 0.0032 inf 0 nan nan nan
x22 -0.0153 inf -0 nan nan nan
x23 0.0172 inf 0 nan nan nan
x24 -0.0060 inf -0 nan nan nan
x25 -0.0208 inf -0 nan nan nan
x26 0.0004 inf 0 nan nan nan
x27 0.0021 inf 0 nan nan nan
x28 -0.0057 inf -0 nan nan nan
x29 -0.0032 inf -0 nan nan nan
x30 -0.0024 inf -0 nan nan nan
x31 -0.0124 inf -0 nan nan nan
x32 -0.0076 inf -0 nan nan nan
x33 -0.0003 inf -0 nan nan nan
x34 0.0045 inf 0 nan nan nan
x35 0.0027 inf 0 nan nan nan
x36 0.0077 inf 0 nan nan nan
x37 0.0025 inf 0 nan nan nan
x38 -0.0096 inf -0 nan nan nan
x39 -0.0002 inf -0 nan nan nan
x40 0.0048 inf 0 nan nan nan
x41 0.0164 inf 0 nan nan nan
x42 0.0261 inf 0 nan nan nan
x43 0.0285 inf 0 nan nan nan
x44 -0.0447 inf -0 nan nan nan
x45 -0.0067 inf -0 nan nan nan
x46 -0.0209 inf -0 nan nan nan
x47 0.0054 inf 0 nan nan nan
x48 -0.0048 inf -0 nan nan nan
x49 0.0034 inf 0 nan nan nan
x50 0.0056 inf 0 nan nan nan
x51 0.0095 inf 0 nan nan nan
x52 0.0021 inf 0 nan nan nan
x53 0.0076 inf 0 nan nan nan
x54 -0.0120 inf -0 nan nan nan
x55 -0.0097 inf -0 nan nan nan
x56 0.0133 inf 0 nan nan nan
x57 -0.0038 inf -0 nan nan nan
x58 0.0029 inf 0 nan nan nan
x59 -0.0036 inf -0 nan nan nan
x60 -0.0130 inf -0 nan nan nan
x61 0.0130 inf 0 nan nan nan
x62 -0.0045 inf -0 nan nan nan
x63 5.425e-05 inf 0 nan nan nan
x64 -0.0073 inf -0 nan nan nan
x65 -0.0033 inf -0 nan nan nan
x66 -0.0090 inf -0 nan nan nan
x67 0.0021 inf 0 nan nan nan
x68 -0.0054 inf -0 nan nan nan
x69 0.0049 inf 0 nan nan nan
x70 -0.0031 inf -0 nan nan nan
x71 0.0055 inf 0 nan nan nan
x72 -0.0069 inf -0 nan nan nan
x73 0.0020 inf 0 nan nan nan
x74 0.0008 inf 0 nan nan nan
x75 0.0025 inf 0 nan nan nan
x76 0.0013 inf 0 nan nan nan
x77 0.0005 inf 0 nan nan nan
x78 0.0026 inf 0 nan nan nan
x79 -0.0031 inf -0 nan nan nan
x80 0.0049 inf 0 nan nan nan
x81 -0.0029 inf -0 nan nan nan
x82 0.0016 inf 0 nan nan nan
x83 -0.0002 inf -0 nan nan nan
x84 -0.0026 inf -0 nan nan nan
x85 0.0039 inf 0 nan nan nan
x86 0.0019 inf 0 nan nan nan
x87 -9.855e-05 inf -0 nan nan nan
x88 -0.0002 inf -0 nan nan nan
x89 0.0017 inf 0 nan nan nan
x90 -0.0014 inf -0 nan nan nan
x91 -0.0094 inf -0 nan nan nan
x92 -0.0003 inf -0 nan nan nan
x93 -0.0063 inf -0 nan nan nan
x94 -0.0033 inf -0 nan nan nan
x95 0.0002 inf 0 nan nan nan
x96 0.0035 inf 0 nan nan nan
x97 0.0046 inf 0 nan nan nan
x98 0.0007 inf 0 nan nan nan
x99 -0.0003 inf -0 nan nan nan
x100 0.0111 inf 0 nan nan nan
x101 -0.0108 inf -0 nan nan nan
x102 -0.0030 inf -0 nan nan nan
x103 -0.0004 inf -0 nan nan nan
x104 -0.0018 inf -0 nan nan nan
x105 -0.0020 inf -0 nan nan nan
x106 -0.0017 inf -0 nan nan nan
x107 0.0010 inf 0 nan nan nan
x108 0.0034 inf 0 nan nan nan
x109 -1.6e-05 inf -0 nan nan nan
x110 -5.495e-05 inf -0 nan nan nan
x111 -0.0012 inf -0 nan nan nan
x112 -0.0003 inf -0 nan nan nan
x113 -4.045e-05 inf -0 nan nan nan
x114 0.0005 inf 0 nan nan nan
x115 -0.0005 inf -0 nan nan nan
x116 -0.0019 inf -0 nan nan nan
x117 -0.0015 inf -0 nan nan nan
x118 -0.0023 inf -0 nan nan nan
x119 -0.0025 inf -0 nan nan nan
x120 0.0023 inf 0 nan nan nan
x121 -0.0118 inf -0 nan nan nan
x122 0.0088 inf 0 nan nan nan
x123 0.0021 inf 0 nan nan nan
x124 0.0032 inf 0 nan nan nan
x125 -0.0049 inf -0 nan nan nan
==============================================================================
Omnibus: 3889.456 Durbin-Watson: 0.002
Prob(Omnibus): 0.000 Jarque-Bera (JB): 13.618
Skew: -0.149 Prob(JB): 0.00110
Kurtosis: 1.417 Cond. No. 1.81e+04
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.81e+04. This might indicate that there are
strong multicollinearity or other numerical problems.
/Users/aditya/mambaforge/envs/stat153fall2026/lib/python3.14/site-packages/statsmodels/regression/linear_model.py:1795: RuntimeWarning: divide by zero encountered in divide
return 1 - (np.divide(self.nobs - self.k_constant, self.df_resid)
/Users/aditya/mambaforge/envs/stat153fall2026/lib/python3.14/site-packages/statsmodels/regression/linear_model.py:1795: RuntimeWarning: invalid value encountered in scalar multiply
return 1 - (np.divide(self.nobs - self.k_constant, self.df_resid)
/Users/aditya/mambaforge/envs/stat153fall2026/lib/python3.14/site-packages/statsmodels/regression/linear_model.py:1717: RuntimeWarning: divide by zero encountered in scalar divide
return np.dot(wresid, wresid) / self.df_resid
In this unregularized estimation, the estimate of equals , the estimate of is , and the estimate of is for . Let us check this.
print(y[0])
print(y[1] - y[0])
print(y[2] - y[1] - y[1] + y[0])
print(mdfull.params)7.306531398939505
0.03947881097378758
0.006542546627510859
const 7.306531
x1 0.039479
x2 0.006543
x3 0.001506
x4 0.004001
...
x121 -0.011814
x122 0.008825
x123 0.002143
x124 0.003224
x125 -0.004903
Length: 126, dtype: float64
We shall study regularization in the next lecture.