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.983
Model: OLS Adj. R-squared: 0.982
Method: Least Squares F-statistic: 893.9
Date: Wed, 26 Nov 2025 Prob (F-statistic): 9.06e-41
Time: 17:00:24 Log-Likelihood: 0.081726
No. Observations: 50 AIC: 7.837
Df Residuals: 46 BIC: 15.48
Df Model: 3
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 5.0030 0.086 58.280 0.000 4.830 5.176
x1 0.5025 0.013 37.953 0.000 0.476 0.529
x2 0.5095 0.052 9.790 0.000 0.405 0.614
x3 -0.0203 0.001 -17.454 0.000 -0.023 -0.018
==============================================================================
Omnibus: 2.872 Durbin-Watson: 2.017
Prob(Omnibus): 0.238 Jarque-Bera (JB): 2.615
Skew: -0.478 Prob(JB): 0.271
Kurtosis: 2.415 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.49578477 4.98255505 5.42933672 5.80836054 6.10187901 6.30508231
6.42688847 6.48847822 6.51981501 6.55472213 6.62532562 6.75677606
6.96311641 7.24497479 7.58946125 7.9722856 8.36174819 8.72394807
9.02835285 9.25281504 9.38721086 9.43510354 9.41315787 9.34840223
9.27378794 9.2227754 9.22383666 9.29577756 9.44464796 9.66274777
9.92989207 10.21672742 10.48955453 10.71586501 10.86968235 10.93582982
10.91242396 10.811184 10.65550648 10.4766212 10.30846031 10.18208329
10.12057405 10.13524821 10.22379407 10.37065234 10.54957254 10.72792586
10.87206595 10.952857 ]
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.93744036 10.78752286 10.52308658 10.18966811 9.84720966 9.55538307
9.35898009 9.27694578 9.29774011 9.38216336]
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 0x7fbdcb337a60>
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.002995
x1 0.502479
np.sin(x1) 0.509536
I((x1 - 5) ** 2) -0.020288
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.937440
1 10.787523
2 10.523087
3 10.189668
4 9.847210
5 9.555383
6 9.358980
7 9.276946
8 9.297740
9 9.382163
dtype: float64