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.978
Model:                            OLS   Adj. R-squared:                  0.977
Method:                 Least Squares   F-statistic:                     694.9
Date:                Mon, 29 Nov 2021   Prob (F-statistic):           2.65e-38
Time:                        22:05:33   Log-Likelihood:                -6.5079
No. Observations:                  50   AIC:                             21.02
Df Residuals:                      46   BIC:                             28.66
Df Model:                           3
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.0010      0.098     51.063      0.000       4.804       5.198
x1             0.5039      0.015     33.358      0.000       0.473       0.534
x2             0.5111      0.059      8.608      0.000       0.392       0.631
x3            -0.0202      0.001    -15.244      0.000      -0.023      -0.018
==============================================================================
Omnibus:                        6.538   Durbin-Watson:                   1.670
Prob(Omnibus):                  0.038   Jarque-Bera (JB):                2.370
Skew:                           0.100   Prob(JB):                        0.306
Kurtosis:                       1.952   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.49561595  4.98328605  5.43088989  5.81057331  6.10453455  6.30794909
  6.42976223  6.49121945  6.52237583  6.55715807  6.62779046  6.75950048
  6.96637403  7.24904114  7.59457246  7.97860351  8.36933773  8.7327705
  9.03827606  9.26363922  9.39870524  9.44704823  9.42538395  9.36082343
  9.28641873  9.23573257  9.23732407  9.31005655  9.45999839  9.67942595
  9.94809259 10.2365551  10.51101118 10.73885301 10.89402466 10.96130299
 10.93879877 10.8382671  10.68317603 10.50485086 10.33732749 10.21176038
 10.15130402 10.16730864 10.25745545 10.40613761 10.58702436 10.76738711
 10.91347665 10.99606805]

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.98256791 10.83422261 10.57107531 10.23880186  9.89752776  9.60710729
  9.41246918  9.33261641  9.35597366  9.44322136]

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 0x7f83ab8dd9a0>
../../../_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.001029
x1                  0.503858
np.sin(x1)          0.511094
I((x1 - 5) ** 2)   -0.020217
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.982568
1    10.834223
2    10.571075
3    10.238802
4     9.897528
5     9.607107
6     9.412469
7     9.332616
8     9.355974
9     9.443221
dtype: float64