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:                     1018.
Date:                Thu, 14 Dec 2023   Prob (F-statistic):           4.80e-42
Time:                        14:39:40   Log-Likelihood:                 3.0746
No. Observations:                  50   AIC:                             1.851
Df Residuals:                      46   BIC:                             9.499
Df Model:                           3
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          4.9317      0.081     60.993      0.000       4.769       5.094
x1             0.5066      0.012     40.625      0.000       0.482       0.532
x2             0.5504      0.049     11.227      0.000       0.452       0.649
x3            -0.0206      0.001    -18.810      0.000      -0.023      -0.018
==============================================================================
Omnibus:                        2.839   Durbin-Watson:                   2.081
Prob(Omnibus):                  0.242   Jarque-Bera (JB):                2.094
Skew:                          -0.492   Prob(JB):                        0.351
Kurtosis:                       3.190   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.41682151  4.92268515  5.38579485  5.77615583  6.07459828  6.27592688
  6.38977441  6.43901906  6.45602567  6.47732824  6.53762749  6.66408968
  6.8718832   7.1616865   7.51957661  7.91931663  8.32666638  8.705008
  9.02136211  9.25180625  9.3854053   9.42600804  9.39161484  9.31142024
  9.22101636  9.15654525  9.14876081  9.21797603  9.37072562  9.59869241
  9.88007377 10.18316377 10.47156237 10.71015582 10.87088578 10.9373592
 10.90754167 10.79409166 10.62228076 10.42584165 10.2414259  10.10258223
 10.03424454 10.04863535 10.14325767 10.30130525 10.49442381 10.68736953
 10.84379957 10.93224304]

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.91763232 10.75840089 10.47613227 10.12001251  9.75478779  9.44491225
  9.23876741  9.15681664  9.18659505  9.28576134]

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 0x7f2614e341c0>
../../../_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.931706
x1                  0.506602
np.sin(x1)          0.550372
I((x1 - 5) ** 2)   -0.020595
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.917632
1    10.758401
2    10.476132
3    10.120013
4     9.754788
5     9.444912
6     9.238767
7     9.156817
8     9.186595
9     9.285761
dtype: float64

Last update: Dec 14, 2023