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:                     959.2
Date:                Sat, 05 Feb 2022   Prob (F-statistic):           1.84e-41
Time:                        17:46:33   Log-Likelihood:                 1.5235
No. Observations:                  50   AIC:                             4.953
Df Residuals:                      46   BIC:                             12.60
Df Model:                           3
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          4.9810      0.083     59.721      0.000       4.813       5.149
x1             0.5254      0.013     40.843      0.000       0.499       0.551
x2             0.4484      0.051      8.868      0.000       0.347       0.550
x3            -0.0228      0.001    -20.226      0.000      -0.025      -0.021
==============================================================================
Omnibus:                        0.849   Durbin-Watson:                   2.239
Prob(Omnibus):                  0.654   Jarque-Bera (JB):                0.886
Skew:                          -0.285   Prob(JB):                        0.642
Kurtosis:                       2.684   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.40990193  4.89176252  5.33676768  5.72047784  6.02727357  6.25292176
  6.40527112  6.50296268  6.57236709  6.64325202  6.74389141  6.89642014
  7.11319733  7.39477578  7.72981112  8.09692572  8.46822088  8.8138606
  9.10697349  9.3280677   9.46823333  9.53060637  9.5298534   9.48976193
  9.43933222  9.40801257  9.42086106  9.49442858  9.6340396   9.83291731
 10.07329708 10.32934512 10.57140289 10.77085996 10.90485455 10.96002971
 10.93472782 10.83926292 10.69422614 10.52710265 10.36775593 10.24352125
 10.1747146  10.17129481 10.23122735 10.3408186  10.47696567 10.61095211
 10.71316601 10.75796509]

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.70990923 10.54098639 10.26878266  9.93337454  9.5875168   9.28372626
  9.06142388  8.93728309  8.90114757  8.9185179 ]

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 0x7f585f3ac970>
../../../_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           4.980974
x1                  0.525362
np.sin(x1)          0.448439
I((x1 - 5) ** 2)   -0.022843
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.709909
1    10.540986
2    10.268783
3     9.933375
4     9.587517
5     9.283726
6     9.061424
7     8.937283
8     8.901148
9     8.918518
dtype: float64