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.985
Model:                            OLS   Adj. R-squared:                  0.984
Method:                 Least Squares   F-statistic:                     1034.
Date:                Sun, 24 Jul 2022   Prob (F-statistic):           3.37e-42
Time:                        16:51:07   Log-Likelihood:                 5.1690
No. Observations:                  50   AIC:                            -2.338
Df Residuals:                      46   BIC:                             5.310
Df Model:                           3
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.1114      0.078     65.919      0.000       4.955       5.267
x1             0.4797      0.012     40.113      0.000       0.456       0.504
x2             0.5009      0.047     10.654      0.000       0.406       0.595
x3            -0.0187      0.001    -17.806      0.000      -0.021      -0.017
==============================================================================
Omnibus:                        2.264   Durbin-Watson:                   1.954
Prob(Omnibus):                  0.322   Jarque-Bera (JB):                1.407
Skew:                           0.365   Prob(JB):                        0.495
Kurtosis:                       3.377   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.64398191  5.11177769  5.54068099  5.90339513  6.18247467  6.37319163
  6.48431228  6.53665687  6.55967896  6.5866262   6.64907792  6.7717567
  6.9684665   7.23982473  7.57316084  7.94459799  8.32297601  8.67497087
  8.97056979  9.18800219  9.31731658  9.36201554  9.33848027  9.27327918
  9.1988028   9.1479421   9.14868461  9.21951611  9.36638351  9.58171777
  9.84567746 10.12940869 10.39978569 10.62485326 10.77907699 10.84753842
 10.82838615 10.73313997 10.58479824 10.4140593  10.25427784 10.13598476
 10.08187113 10.10306    10.197279   10.34923364 10.53312018 10.71686507
 10.86739458 10.95606872]

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.95282189 10.81798679 10.57120546 10.25723955  9.93501115  9.66317661
  9.4857654   9.42140017  9.45873728  9.55924412]

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 0x7fcfa1d68d30>
../../../_images/examples_notebooks_generated_predict_12_1.png

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.111378
x1                  0.479700
np.sin(x1)          0.500864
I((x1 - 5) ** 2)   -0.018696
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.952822
1    10.817987
2    10.571205
3    10.257240
4     9.935011
5     9.663177
6     9.485765
7     9.421400
8     9.458737
9     9.559244
dtype: float64