Prediction (out of sample)

%matplotlib inline

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

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

olsmod = sm.OLS(y, X)
olsres = olsmod.fit()
print(olsres.summary())

                            OLS Regression Results
==============================================================================
Dep. Variable:                      y   R-squared:                       0.986
Model:                            OLS   Adj. R-squared:                  0.985
Method:                 Least Squares   F-statistic:                     1056.
Date:                Tue, 08 Feb 2022   Prob (F-statistic):           2.09e-42
Time:                        18:18:22   Log-Likelihood:                 3.5963
No. Observations:                  50   AIC:                            0.8074
Df Residuals:                      46   BIC:                             8.455
Df Model:                           3
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          4.9445      0.080     61.793      0.000       4.783       5.106
x1             0.5051      0.012     40.926      0.000       0.480       0.530
x2             0.5786      0.049     11.927      0.000       0.481       0.676
x3            -0.0201      0.001    -18.548      0.000      -0.022      -0.018
==============================================================================
Omnibus:                        5.168   Durbin-Watson:                   2.286
Prob(Omnibus):                  0.075   Jarque-Bera (JB):                4.064
Skew:                          -0.646   Prob(JB):                        0.131
Kurtosis:                       3.530   Cond. No.                         221.
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

In-sample prediction

ypred = olsres.predict(X)
print(ypred)

[ 4.44209751  4.95658678  5.4266468   5.82074414  6.1187256   6.31512937
  6.42008236  6.45763636  6.46181645  6.47103071  6.52175996  6.64256397
  6.84938926  7.14294909  7.50860625  7.91877796  8.33746763  8.72617897
  9.05024085  9.28450366  9.41747131  9.45319003  9.41058356  9.32034405
  9.2198894   9.14721566  9.13465432  9.20356039  9.36080379  9.59764068
  9.89115002 10.20799947 10.50992173 10.76000182 10.92874206 10.99890851
 10.96836237 10.85041146 10.67162376 10.46746233 10.27645888 10.1338831
 10.06594792 10.08550286 10.18992285 10.3615395  10.57054403 10.77988471
 10.9513544  11.05186784]

Create a new sample of explanatory variables Xnew, predict and plot

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)

[11.04644133 10.88975255 10.60449217 10.24236938  9.87145168  9.55949953
  9.35737629  9.28659508  9.33405158  9.45523248]

Plot comparison

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")

<matplotlib.legend.Legend at 0x7fd1f43624c0>

Predicting with Formulas

Using formulas can make both estimation and prediction a lot easier

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

res.params

Intercept           4.944539
x1                  0.505053
np.sin(x1)          0.578604
I((x1 - 5) ** 2)   -0.020098
dtype: float64

Now we only have to pass the single variable and we get the transformed right-hand side variables automatically

res.predict(exog=dict(x1=x1n))

0    11.046441
1    10.889753
2    10.604492
3    10.242369
4     9.871452
5     9.559500
6     9.357376
7     9.286595
8     9.334052
9     9.455232
dtype: float64