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.980
Model: OLS Adj. R-squared: 0.979
Method: Least Squares F-statistic: 764.5
Date: Thu, 12 Mar 2026 Prob (F-statistic): 3.09e-39
Time: 16:11:17 Log-Likelihood: -4.0099
No. Observations: 50 AIC: 16.02
Df Residuals: 46 BIC: 23.67
Df Model: 3
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 4.9305 0.093 52.923 0.000 4.743 5.118
x1 0.5143 0.014 35.793 0.000 0.485 0.543
x2 0.5335 0.056 9.446 0.000 0.420 0.647
x3 -0.0216 0.001 -17.112 0.000 -0.024 -0.019
==============================================================================
Omnibus: 0.523 Durbin-Watson: 1.922
Prob(Omnibus): 0.770 Jarque-Bera (JB): 0.582
Skew: 0.225 Prob(JB): 0.747
Kurtosis: 2.722 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.39083819 4.8970421 5.36125956 5.7544136 6.05792096 6.26674532
6.39022468 6.45053721 6.47905741 6.51120137 6.58060803 6.71361236
6.92491864 7.21518455 7.57091315 7.96667035 8.3692637 8.7431958
9.05649639 9.28597503 9.42103119 9.46539589 9.43651869 9.36270084
9.27844536 9.21878839 9.2135428 9.28240003 9.43169483 9.65336453
9.92627362 10.21968636 10.4983165 10.72812474 10.88191128 10.94378469
10.91177281 10.79814676 10.62740487 10.43224763 10.2482051 10.10779946
10.03520165 10.0422602 10.12655453 10.27179263 10.45048784 10.62847477
10.7705226 10.84612345]
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.81880584 10.65048132 10.36207295 10.00126177 9.63081278 9.3132079
9.09534822 8.99707096 9.00629266 9.08196758]
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 0x7fe88f75e020>
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.930540
x1 0.514291
np.sin(x1) 0.533531
I((x1 - 5) ** 2) -0.021588
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.818806
1 10.650481
2 10.362073
3 10.001262
4 9.630813
5 9.313208
6 9.095348
7 8.997071
8 9.006293
9 9.081968
dtype: float64