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.983
Model: OLS Adj. R-squared: 0.982
Method: Least Squares F-statistic: 891.0
Date: Tue, 27 Oct 2026 Prob (F-statistic): 9.74e-41
Time: 21:59:23 Log-Likelihood: 1.2159
No. Observations: 50 AIC: 5.568
Df Residuals: 46 BIC: 13.22
Df Model: 3
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 5.0610 0.084 60.308 0.000 4.892 5.230
x1 0.4928 0.013 38.074 0.000 0.467 0.519
x2 0.4890 0.051 9.610 0.000 0.387 0.591
x3 -0.0201 0.001 -17.676 0.000 -0.022 -0.018
==============================================================================
Omnibus: 4.387 Durbin-Watson: 2.148
Prob(Omnibus): 0.112 Jarque-Bera (JB): 3.357
Skew: -0.599 Prob(JB): 0.187
Kurtosis: 3.417 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.55885307 5.03270294 5.46797333 5.83801624 6.12580079 6.32671136
6.4493059 6.5139099 6.54927698 6.5878647 6.66050195 6.79132379
6.99380616 7.2685518 7.60319125 7.97441509 8.35180361 8.70282472
8.99817883 9.21661273 9.34841141 9.39699421 9.37835305 9.31842498
9.2488308 9.20167989 9.20429464 9.27472158 9.41876628 9.62903957
9.88617154 10.16199396 10.42416826 10.64149864 10.78905746 10.85228079
10.82936122 10.73154479 10.58128334 10.40854569 10.24589388 10.12313332
10.0624158 10.07460009 10.1574678 10.2960879 10.46526983 10.63370226
10.7690974 10.84349585]
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.82458682 10.67629935 10.41780868 10.09281276 9.75883347 9.47313326
9.27869529 9.19369966 9.20707231 9.2811963 ]
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 0x720c44086a50>
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.061003
x1 0.492772
np.sin(x1) 0.488962
I((x1 - 5) ** 2) -0.020086
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.824587
1 10.676299
2 10.417809
3 10.092813
4 9.758833
5 9.473133
6 9.278695
7 9.193700
8 9.207072
9 9.281196
dtype: float64