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.987
Model:                            OLS   Adj. R-squared:                  0.986
Method:                 Least Squares   F-statistic:                     1170.
Date:                Sat, 09 Jul 2022   Prob (F-statistic):           2.03e-43
Time:                        10:15:30   Log-Likelihood:                 7.7837
No. Observations:                  50   AIC:                            -7.567
Df Residuals:                      46   BIC:                           0.08071
Df Model:                           3
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.1052      0.074     69.374      0.000       4.957       5.253
x1             0.4843      0.011     42.670      0.000       0.461       0.507
x2             0.4330      0.045      9.705      0.000       0.343       0.523
x3            -0.0189      0.001    -18.929      0.000      -0.021      -0.017
==============================================================================
Omnibus:                        1.123   Durbin-Watson:                   2.568
Prob(Omnibus):                  0.570   Jarque-Bera (JB):                0.429
Skew:                           0.094   Prob(JB):                        0.807
Kurtosis:                       3.413   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.63363988  5.07701837  5.4858746   5.83661057  6.11414468  6.31438958
  6.44492373  6.52374625  6.57631976  6.63138689  6.71624801  6.85227578
  7.05140376  7.3141657   7.62960788  7.97708878  8.32967037  8.65854404
  8.93776368  9.14850863  9.2821758   9.34179319  9.34152234  9.30433168
  9.25822285  9.23163012  9.24874873  9.32555965  9.467204    9.66713833
  9.90820971 10.16547393 10.41029371 10.61504379 10.75764963 10.82521421
 10.81613699 10.74037684 10.61781587 10.47499283 10.34074284 10.24145993
 10.19676051 10.21626051 10.29799543 10.4287432  10.58619643 10.74262753
 10.86944424 10.94188724]

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.93268982 10.80838416 10.58595081 10.30408625 10.01372868  9.76558657
  9.59772343  9.52623856  9.54132537  9.60967239]

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 0x7ff3f04c7e20>
../../../_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.105208
x1                  0.484274
np.sin(x1)          0.432998
I((x1 - 5) ** 2)   -0.018863
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.932690
1    10.808384
2    10.585951
3    10.304086
4    10.013729
5     9.765587
6     9.597723
7     9.526239
8     9.541325
9     9.609672
dtype: float64