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.984
Model: OLS Adj. R-squared: 0.983
Method: Least Squares F-statistic: 923.2
Date: Wed, 09 Sep 2026 Prob (F-statistic): 4.37e-41
Time: 05:30:58 Log-Likelihood: 1.0739
No. Observations: 50 AIC: 5.852
Df Residuals: 46 BIC: 13.50
Df Model: 3
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 5.0381 0.084 59.865 0.000 4.869 5.207
x1 0.4870 0.013 37.518 0.000 0.461 0.513
x2 0.4937 0.051 9.675 0.000 0.391 0.596
x3 -0.0186 0.001 -16.295 0.000 -0.021 -0.016
==============================================================================
Omnibus: 0.272 Durbin-Watson: 1.752
Prob(Omnibus): 0.873 Jarque-Bera (JB): 0.105
Skew: 0.112 Prob(JB): 0.949
Kurtosis: 2.975 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.57385062 5.04125498 5.47027837 5.83401646 6.11527457 6.30939267
6.42501103 6.48265065 6.51134172 6.54385402 6.61131298 6.7380858
6.93777805 7.21099828 7.54525803 7.91702367 8.2955829 8.6480907
8.94496565 9.1647501 9.2976356 9.34707438 9.32921219 9.27023555
9.20206944 9.15713207 9.16300867 9.23791934 9.38772549 9.60496689
9.8700873 10.15464755 10.425998 10.65264278 10.80941449 10.88160909
10.8674018 10.77814683 10.63651198 10.47275434 10.31974916 10.20758854
10.15863736 10.18385882 10.28101321 10.43502587 10.62046394 10.80571489
10.95818038 11.04963207]
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)
[11.05133099 10.92326175 10.68478404 10.38001614 10.06703322 9.80364853
9.63325869 9.57421867 9.61534789 9.71866757]
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 0x7f4994979010>
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 5.038093
x1 0.486951
np.sin(x1) 0.493665
I((x1 - 5) ** 2) -0.018570
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 11.051331
1 10.923262
2 10.684784
3 10.380016
4 10.067033
5 9.803649
6 9.633259
7 9.574219
8 9.615348
9 9.718668
dtype: float64