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.987
Model: OLS Adj. R-squared: 0.986
Method: Least Squares F-statistic: 1130.
Date: Sat, 12 Sep 2026 Prob (F-statistic): 4.52e-43
Time: 23:49:14 Log-Likelihood: 4.2737
No. Observations: 50 AIC: -0.5473
Df Residuals: 46 BIC: 7.101
Df Model: 3
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 4.9113 0.079 62.215 0.000 4.752 5.070
x1 0.5198 0.012 42.695 0.000 0.495 0.544
x2 0.5309 0.048 11.093 0.000 0.435 0.627
x3 -0.0210 0.001 -19.663 0.000 -0.023 -0.019
==============================================================================
Omnibus: 1.579 Durbin-Watson: 1.702
Prob(Omnibus): 0.454 Jarque-Bera (JB): 1.416
Skew: 0.274 Prob(JB): 0.493
Kurtosis: 2.383 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.38583627 4.89102108 5.35457924 5.74757578 6.05151819 6.26139473
6.38649781 6.4488973 6.47981449 6.51449252 6.58640606 6.72176149
6.93519128 7.22734989 7.5848063 7.98225063 8.38665249 8.76268784
9.07854254 9.31113939 9.44992963 9.49862605 9.47459302 9.40599355
9.32716221 9.27296404 9.27306632 9.34706434 9.50126191 9.72763572
10.0051535 10.30322958 10.58675002 10.82184174 10.98143778 11.0497243
11.02473877 10.9186924 10.75596387 10.569094 10.39343934 10.26136325
10.19691865 10.21189621 10.30388716 10.45667906 10.64291946 10.82860982
10.97869137 11.06280517]
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.04671252 10.89053379 10.61508987 10.26782897 9.91120958 9.60740851
9.40309783 9.31801788 9.3401442 9.42863151]
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 0x7fe08c6dcec0>
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.911289
x1 0.519794
np.sin(x1) 0.530926
I((x1 - 5) ** 2) -0.021018
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.046713
1 10.890534
2 10.615090
3 10.267829
4 9.911210
5 9.607409
6 9.403098
7 9.318018
8 9.340144
9 9.428632
dtype: float64