Prediction (out of sample)¶
[1]:
%matplotlib inline
[2]:
import numpy as np
import matplotlib.pyplot as plt
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.979
Method: Least Squares F-statistic: 776.1
Date: Tue, 02 Sep 2025 Prob (F-statistic): 2.20e-39
Time: 08:44:57 Log-Likelihood: -1.8200
No. Observations: 50 AIC: 11.64
Df Residuals: 46 BIC: 19.29
Df Model: 3
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 5.1041 0.089 57.239 0.000 4.925 5.284
x1 0.4896 0.014 35.602 0.000 0.462 0.517
x2 0.3586 0.054 6.632 0.000 0.250 0.467
x3 -0.0201 0.001 -16.614 0.000 -0.022 -0.018
==============================================================================
Omnibus: 1.777 Durbin-Watson: 2.201
Prob(Omnibus): 0.411 Jarque-Bera (JB): 0.983
Skew: -0.285 Prob(JB): 0.612
Kurtosis: 3.384 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.60259219 5.02330623 5.4139524 5.75498903 6.03392695 6.24738141
6.40162818 6.51157245 6.59829995 6.68561267 6.79611818 6.94751534
7.1496864 7.40307336 7.69860538 8.01918907 8.34251691 8.64473227
8.90434905 9.10578172 9.24190602 9.31522943 9.33747925 9.32767575
9.30900728 9.30502044 9.33575144 9.41443424 9.54532623 9.72300874
9.93327731 10.15547549 10.36588858 10.54163991 10.66444942 10.72363701
10.71787726 10.65541719 10.55272127 10.43176641 10.31643118 10.22857262
10.18443522 10.19198197 10.24958607 10.34629824 10.46364558 10.57866672
10.6676842 10.71019469]
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.67914857 10.55055999 10.33849066 10.07498537 9.80222637 9.56220565
9.38644386 9.28827292 9.25957195 9.27275548]
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 0x70a5b77e1610>
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.104130
x1 0.489628
np.sin(x1) 0.358568
I((x1 - 5) ** 2) -0.020062
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.679149
1 10.550560
2 10.338491
3 10.074985
4 9.802226
5 9.562206
6 9.386444
7 9.288273
8 9.259572
9 9.272755
dtype: float64
