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.985
Model:                            OLS   Adj. R-squared:                  0.984
Method:                 Least Squares   F-statistic:                     1000.
Date:                Sat, 19 Aug 2023   Prob (F-statistic):           7.08e-42
Time:                        20:59:16   Log-Likelihood:                 2.5007
No. Observations:                  50   AIC:                             2.999
Df Residuals:                      46   BIC:                             10.65
Df Model:                           3
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          4.9979      0.082     61.106      0.000       4.833       5.163
x1             0.5028      0.013     39.857      0.000       0.477       0.528
x2             0.6145      0.050     12.393      0.000       0.515       0.714
x3            -0.0201      0.001    -18.106      0.000      -0.022      -0.018
==============================================================================
Omnibus:                        2.166   Durbin-Watson:                   2.203
Prob(Omnibus):                  0.339   Jarque-Bera (JB):                1.651
Skew:                          -0.445   Prob(JB):                        0.438
Kurtosis:                       3.043   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.49655741  5.02419218  5.50507012  5.90570016  6.20467798  6.39620268
  6.49102986  6.5147045   6.50336396  6.49780078  6.53676052  6.65057602
  6.85618375  7.15434124  7.52950262  7.95237272  8.38472022  8.78565894
  9.11836533  9.35612846  9.48673854  9.514493    9.45949039  9.35432827
  9.23874749  9.15310293  9.13173328  9.19731927  9.35715735  9.60196082
  9.9073857  10.2380303  10.55325167 10.81384309 10.98847563 11.05884542
 11.02268091 10.89411609 10.70136828 10.48210216 10.27724154 10.12424568
 10.05095485 10.07101629 10.18164194 10.36406657 10.58663081 10.80998272
 10.9935432  11.10217301]

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)
[11.09929507 10.93641054 10.63761877 10.25783911  9.86936477  9.54416288
  9.33625447  9.2684882   9.32694623  9.46435179]

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 0x7f69e0fb1650>
../../../_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.997879
x1                  0.502758
np.sin(x1)          0.614525
I((x1 - 5) ** 2)   -0.020053
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    11.099295
1    10.936411
2    10.637619
3    10.257839
4     9.869365
5     9.544163
6     9.336254
7     9.268488
8     9.326946
9     9.464352
dtype: float64