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.984
Model:                            OLS   Adj. R-squared:                  0.983
Method:                 Least Squares   F-statistic:                     923.4
Date:                Wed, 30 Nov 2022   Prob (F-statistic):           4.34e-41
Time:                        21:29:37   Log-Likelihood:                 1.1397
No. Observations:                  50   AIC:                             5.721
Df Residuals:                      46   BIC:                             13.37
Df Model:                           3
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.0182      0.084     59.707      0.000       4.849       5.187
x1             0.4945      0.013     38.153      0.000       0.468       0.521
x2             0.4081      0.051      8.009      0.000       0.306       0.511
x3            -0.0195      0.001    -17.163      0.000      -0.022      -0.017
==============================================================================
Omnibus:                        0.829   Durbin-Watson:                   1.771
Prob(Omnibus):                  0.661   Jarque-Bera (JB):                0.865
Skew:                          -0.282   Prob(JB):                        0.649
Kurtosis:                       2.689   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.52991477  4.9702252   5.37741426  5.72924135  6.01149236  6.22031506
  6.36285198  6.45606693  6.52395788  6.59361416  6.6907658   6.83555609
  7.03923216  7.30229705  7.61442711  7.95616814  8.30213162  8.62516602
  8.90081769  9.11134854  9.24865036  9.31557694  9.32547527  9.2999926
  9.26551988  9.24885584  9.27280396  9.35242598  9.4925673   9.68706086
  9.91974037 10.16709645 10.40213905 10.59883195 10.73637034 10.80259907
 10.79600966 10.72598813 10.61127295 10.47687634 10.34997481 10.25544416
 10.21177247 10.22802252 10.30234284 10.42227188 10.56678542 10.70975077
 10.82422032 10.88685894]

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.87048771 10.74475174 10.52565482 10.24966759  9.96479824  9.71883843
  9.54766222  9.46644288  9.46593789  9.51475189]

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 0x7f3b50e96c20>
../../../_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.018240
x1                  0.494549
np.sin(x1)          0.408091
I((x1 - 5) ** 2)   -0.019533
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.870488
1    10.744752
2    10.525655
3    10.249668
4     9.964798
5     9.718838
6     9.547662
7     9.466443
8     9.465938
9     9.514752
dtype: float64