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.986
Model: OLS Adj. R-squared: 0.985
Method: Least Squares F-statistic: 1056.
Date: Tue, 08 Feb 2022 Prob (F-statistic): 2.09e-42
Time: 18:18:22 Log-Likelihood: 3.5963
No. Observations: 50 AIC: 0.8074
Df Residuals: 46 BIC: 8.455
Df Model: 3
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 4.9445 0.080 61.793 0.000 4.783 5.106
x1 0.5051 0.012 40.926 0.000 0.480 0.530
x2 0.5786 0.049 11.927 0.000 0.481 0.676
x3 -0.0201 0.001 -18.548 0.000 -0.022 -0.018
==============================================================================
Omnibus: 5.168 Durbin-Watson: 2.286
Prob(Omnibus): 0.075 Jarque-Bera (JB): 4.064
Skew: -0.646 Prob(JB): 0.131
Kurtosis: 3.530 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.44209751 4.95658678 5.4266468 5.82074414 6.1187256 6.31512937
6.42008236 6.45763636 6.46181645 6.47103071 6.52175996 6.64256397
6.84938926 7.14294909 7.50860625 7.91877796 8.33746763 8.72617897
9.05024085 9.28450366 9.41747131 9.45319003 9.41058356 9.32034405
9.2198894 9.14721566 9.13465432 9.20356039 9.36080379 9.59764068
9.89115002 10.20799947 10.50992173 10.76000182 10.92874206 10.99890851
10.96836237 10.85041146 10.67162376 10.46746233 10.27645888 10.1338831
10.06594792 10.08550286 10.18992285 10.3615395 10.57054403 10.77988471
10.9513544 11.05186784]
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.04644133 10.88975255 10.60449217 10.24236938 9.87145168 9.55949953
9.35737629 9.28659508 9.33405158 9.45523248]
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 0x7fd1f43624c0>

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.944539
x1 0.505053
np.sin(x1) 0.578604
I((x1 - 5) ** 2) -0.020098
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.046441
1 10.889753
2 10.604492
3 10.242369
4 9.871452
5 9.559500
6 9.357376
7 9.286595
8 9.334052
9 9.455232
dtype: float64