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.982
Model: OLS Adj. R-squared: 0.981
Method: Least Squares F-statistic: 851.2
Date: Sun, 06 Feb 2022 Prob (F-statistic): 2.73e-40
Time: 19:08:08 Log-Likelihood: -1.7968
No. Observations: 50 AIC: 11.59
Df Residuals: 46 BIC: 19.24
Df Model: 3
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 4.9213 0.089 55.214 0.000 4.742 5.101
x1 0.5156 0.014 37.510 0.000 0.488 0.543
x2 0.4931 0.054 9.125 0.000 0.384 0.602
x3 -0.0213 0.001 -17.687 0.000 -0.024 -0.019
==============================================================================
Omnibus: 0.315 Durbin-Watson: 2.170
Prob(Omnibus): 0.854 Jarque-Bera (JB): 0.055
Skew: 0.078 Prob(JB): 0.973
Kurtosis: 3.050 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.38760694 4.87735555 5.3278358 5.71217509 6.01319901 6.22625304
6.35996726 6.43483846 6.47986257 6.52777072 6.6096516 6.74984368
6.96193648 7.24653781 7.59117396 7.97233907 8.35935706 8.71942166
9.02298637 9.24861884 9.3865222 9.4401446 9.42561274 9.36908226
9.30244045 9.2580673 9.26351547 9.33698358 9.48432612 9.69809166
9.95874703 10.23788663 10.50289932 10.7223264 10.87103015 10.93432389
10.91038499 10.81055427 10.65747277 10.48136189 10.31505811 10.188618
10.12438008 10.1332946 10.21312448 10.34881287 10.51495691 10.67998137
10.81132635 10.88079683]
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.85442785 10.69704421 10.42798273 10.09130966 9.74503165 9.44689374
9.24024135 9.14340779 9.14522567 9.20776103]
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 0x7f1d6b24eee0>
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.921282
x1 0.515622
np.sin(x1) 0.493083
I((x1 - 5) ** 2) -0.021347
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.854428
1 10.697044
2 10.427983
3 10.091310
4 9.745032
5 9.446894
6 9.240241
7 9.143408
8 9.145226
9 9.207761
dtype: float64
