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.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.988
Model:                            OLS   Adj. R-squared:                  0.987
Method:                 Least Squares   F-statistic:                     1248.
Date:                Mon, 07 Dec 2020   Prob (F-statistic):           4.71e-44
Time:                        17:22:22   Log-Likelihood:                 8.8151
No. Observations:                  50   AIC:                            -9.630
Df Residuals:                      46   BIC:                            -1.982
Df Model:                           3
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.0615      0.072     70.213      0.000       4.916       5.207
x1             0.4840      0.011     43.534      0.000       0.462       0.506
x2             0.4200      0.044      9.610      0.000       0.332       0.508
x3            -0.0184      0.001    -18.830      0.000      -0.020      -0.016
==============================================================================
Omnibus:                       13.108   Durbin-Watson:                   1.819
Prob(Omnibus):                  0.001   Jarque-Bera (JB):               17.300
Skew:                           0.880   Prob(JB):                     0.000175
Kurtosis:                       5.282   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.60194513  5.03816284  5.44086647  5.7871662   6.06243306  6.26270236
  6.39532512  6.4777603   6.53470635  6.59404332  6.68225225  6.82006441
  7.0190552   7.27974236  7.59150091  7.93430879  8.28203642  8.60673974
  8.88325133  9.09331534  9.22858689  9.29200308  9.29730047  9.26675817
  9.22753747  9.20721931  9.22927278  9.30919928  9.45198562  9.65128472
  9.89045845 10.14531117 10.38806506 10.5919239  10.73547564 10.80621031
 10.80257539 10.73423089 10.62046224 10.48701172 10.36184911 10.2705763
 10.23222101 10.25611057 10.34033946 10.4720824  10.62970147 10.78630122
 10.91414745 10.98922392]

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.98615378 10.87138808 10.66139777 10.39371801 10.11775831  9.8827053
  9.72548019  9.66169913  9.68184982  9.7536205 ]

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");
../../../_images/examples_notebooks_generated_predict_12_0.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.061452
x1                  0.483995
np.sin(x1)          0.420003
I((x1 - 5) ** 2)   -0.018380
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.986154
1    10.871388
2    10.661398
3    10.393718
4    10.117758
5     9.882705
6     9.725480
7     9.661699
8     9.681850
9     9.753621
dtype: float64