Prediction (out of sample)#
[1]:
%matplotlib inline
[2]:
import matplotlib.pyplot as plt
import numpy as np
import statsmodels.api as sm
plt.rc("figure", figsize=(16, 8))
plt.rc("font", size=14)
Artificial data#
[3]:
nsample = 50
sig = 0.25
x1 = np.linspace(0, 20, nsample)
X = np.column_stack((x1, np.sin(x1), (x1 - 5) ** 2))
X = sm.add_constant(X)
beta = [5.0, 0.5, 0.5, -0.02]
y_true = np.dot(X, beta)
y = y_true + sig * np.random.normal(size=nsample)
Estimation#
[4]:
olsmod = sm.OLS(y, X)
olsres = olsmod.fit()
print(olsres.summary())
OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.981
Model: OLS Adj. R-squared: 0.980
Method: Least Squares F-statistic: 803.2
Date: Wed, 28 Oct 2026 Prob (F-statistic): 1.01e-39
Time: 15:13:35 Log-Likelihood: -3.4845
No. Observations: 50 AIC: 14.97
Df Residuals: 46 BIC: 22.62
Df Model: 3
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 4.9256 0.092 53.428 0.000 4.740 5.111
x1 0.5006 0.014 35.206 0.000 0.472 0.529
x2 0.4256 0.056 7.614 0.000 0.313 0.538
x3 -0.0193 0.001 -15.488 0.000 -0.022 -0.017
==============================================================================
Omnibus: 0.585 Durbin-Watson: 2.082
Prob(Omnibus): 0.746 Jarque-Bera (JB): 0.716
Skew: -0.195 Prob(JB): 0.699
Kurtosis: 2.563 Cond. No. 221.
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
In-sample prediction#
[5]:
ypred = olsres.predict(X)
print(ypred)
[ 4.44219718 4.89112022 5.30584899 5.66319111 5.94832418 6.15723111
6.2973601 6.38640092 6.44937853 6.51454154 6.60872109 6.75292249
6.95887419 7.22710078 7.5468368 7.89779549 8.25350186 8.58564243
8.86871705 9.08422855 9.22372171 9.29017248 9.297499 9.26827491
9.23002058 9.21068156 9.23403722 9.31579386 9.46100422 9.66323734
9.90563509 10.16368201 10.40923303 10.61513769 10.75970087 10.83024723
10.82520373 10.75435807 10.63725064 10.49996403 10.37083777 10.27581222
10.23416664 10.25535168 10.3374366 10.46742659 10.62339769 10.77809886
10.90342928 10.9750552 ]
Create a new sample of explanatory variables Xnew, predict and plot#
[6]:
x1n = np.linspace(20.5, 25, 10)
Xnew = np.column_stack((x1n, np.sin(x1n), (x1n - 5) ** 2))
Xnew = sm.add_constant(Xnew)
ynewpred = olsres.predict(Xnew) # predict out of sample
print(ynewpred)
[10.9661786 10.843774 10.6245301 10.34647828 10.05968123 9.81397584
9.64677141 9.57389052 9.58569514 9.64944634]
Plot comparison#
[7]:
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.plot(x1, y, "o", label="Data")
ax.plot(x1, y_true, "b-", label="True")
ax.plot(np.hstack((x1, x1n)), np.hstack((ypred, ynewpred)), "r", label="OLS prediction")
ax.legend(loc="best")
[7]:
<matplotlib.legend.Legend at 0x7f4b51aa6ba0>
Predicting with Formulas#
Using formulas can make both estimation and prediction a lot easier
[8]:
from statsmodels.formula.api import ols
data = {"x1": x1, "y": y}
res = ols("y ~ x1 + np.sin(x1) + I((x1-5)**2)", data=data).fit()
We use the I to indicate use of the Identity transform. Ie., we do not want any expansion magic from using **2
[9]:
res.params
[9]:
Intercept 4.925570
x1 0.500567
np.sin(x1) 0.425556
I((x1 - 5) ** 2) -0.019335
dtype: float64
Now we only have to pass the single variable and we get the transformed right-hand side variables automatically
[10]:
res.predict(exog=dict(x1=x1n))
[10]:
0 10.966179
1 10.843774
2 10.624530
3 10.346478
4 10.059681
5 9.813976
6 9.646771
7 9.573891
8 9.585695
9 9.649446
dtype: float64