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.974
Model:                            OLS   Adj. R-squared:                  0.972
Method:                 Least Squares   F-statistic:                     578.4
Date:                Thu, 03 Aug 2023   Prob (F-statistic):           1.63e-36
Time:                        21:28:34   Log-Likelihood:                -9.4278
No. Observations:                  50   AIC:                             26.86
Df Residuals:                      46   BIC:                             34.50
Df Model:                           3
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.1255      0.104     49.366      0.000       4.917       5.335
x1             0.4860      0.016     30.350      0.000       0.454       0.518
x2             0.5127      0.063      8.144      0.000       0.386       0.639
x3            -0.0194      0.001    -13.796      0.000      -0.022      -0.017
==============================================================================
Omnibus:                        0.357   Durbin-Watson:                   2.065
Prob(Omnibus):                  0.836   Jarque-Bera (JB):                0.115
Skew:                          -0.117   Prob(JB):                        0.944
Kurtosis:                       3.027   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.64061756  5.11840817  5.55630284  5.92636128  6.21072672  6.40455968
  6.51683315  6.56885828  6.59078308  6.61663916  6.6787504   6.80242231
  7.00178433  7.27746844  7.61650531  7.99445514  8.37942304  8.73729929
  9.03736342  9.25733156  9.38701757  9.43000662  9.40306615  9.333391
  9.25413543  9.19896593  9.1965299   9.26574901  9.41271037  9.62966648
  9.8963081  10.18310099 10.45613833 10.68271149 10.83668399 10.90278561
 10.8791213  10.77748248 10.62140976 10.44232531 10.27437024 10.14879536
 10.08882671 10.10584963 10.19753818 10.34823723 10.53153464 10.71460062
 10.86358162 10.94916236]

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.93937053 10.7947581  10.53543021 10.20720392  9.87039062  9.58502965
  9.39618865  9.32292935  9.35364065  9.4488814 ]

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 0x7f650562a510>
../../../_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.125519
x1                  0.485986
np.sin(x1)          0.512674
I((x1 - 5) ** 2)   -0.019396
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.939371
1    10.794758
2    10.535430
3    10.207204
4     9.870391
5     9.585030
6     9.396189
7     9.322929
8     9.353641
9     9.448881
dtype: float64