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.977
Model:                            OLS   Adj. R-squared:                  0.976
Method:                 Least Squares   F-statistic:                     658.4
Date:                Sun, 06 Feb 2022   Prob (F-statistic):           8.91e-38
Time:                        19:08:08   Log-Likelihood:                -7.1077
No. Observations:                  50   AIC:                             22.22
Df Residuals:                      46   BIC:                             29.86
Df Model:                           3
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.0450      0.099     50.898      0.000       4.845       5.245
x1             0.5011      0.015     32.779      0.000       0.470       0.532
x2             0.4107      0.060      6.834      0.000       0.290       0.532
x3            -0.0205      0.001    -15.270      0.000      -0.023      -0.018
==============================================================================
Omnibus:                        0.434   Durbin-Watson:                   2.493
Prob(Omnibus):                  0.805   Jarque-Bera (JB):                0.584
Skew:                          -0.080   Prob(JB):                        0.747
Kurtosis:                       2.495   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.53262493  4.9804055   5.394574    5.75274767  6.04062159  6.25431889
  6.40102776  6.49782041  6.56884817  6.64137349  6.74129083  6.8888722
  7.09543661  7.36149031  7.67664369  8.02131821  8.36996312  8.69525336
  8.97257902  9.18408865  9.32162221  9.38805175  9.39680955  9.36968118
  9.33322627  9.31441478  9.33619594  9.41372774  9.55188653  9.74446584
  9.97519609 10.22041763 10.45296797 10.64664446 10.78050938 10.84232986
 10.8305878  10.75472921 10.6336123  10.49240913 10.35846988 10.25682933
 10.20609369 10.21538368 10.28283599 10.39590944 10.53344543 10.66914429
 10.77588619 10.83018688]

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.80262641 10.66458354 10.43216434 10.14207251  9.84262305  9.58191303
  9.39604577  9.30029149  9.28534863  9.31962117]

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 0x7f2ff3eb48e0>
../../../_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.045008
x1                  0.501084
np.sin(x1)          0.410699
I((x1 - 5) ** 2)   -0.020495
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.802626
1    10.664584
2    10.432164
3    10.142073
4     9.842623
5     9.581913
6     9.396046
7     9.300291
8     9.285349
9     9.319621
dtype: float64