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.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.984
Model:                            OLS   Adj. R-squared:                  0.982
Method:                 Least Squares   F-statistic:                     915.8
Date:                Sun, 29 Aug 2021   Prob (F-statistic):           5.24e-41
Time:                        19:09:28   Log-Likelihood:                 1.5726
No. Observations:                  50   AIC:                             4.855
Df Residuals:                      46   BIC:                             12.50
Df Model:                           3
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.1183      0.083     61.428      0.000       4.951       5.286
x1             0.4856      0.013     37.792      0.000       0.460       0.512
x2             0.4874      0.051      9.649      0.000       0.386       0.589
x3            -0.0190      0.001    -16.808      0.000      -0.021      -0.017
==============================================================================
Omnibus:                        1.610   Durbin-Watson:                   2.437
Prob(Omnibus):                  0.447   Jarque-Bera (JB):                1.236
Skew:                          -0.385   Prob(JB):                        0.539
Kurtosis:                       2.975   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.64425117  5.11019528  5.53803293  5.90119891  6.18271521  6.37798056
  6.49552632  6.55561486  6.58691048  6.62176998  6.69092646  6.81844004
  7.01774481  7.28944164  7.62119947  7.98978128  8.36486185  8.7140101
  9.00801747  9.22569673  9.35736305  9.40642516  9.3888253   9.33041997
  9.26273177  9.21777034  9.22277325  9.29573094  9.44243097  9.65550711
  9.91564962 10.19477775 10.46065336 10.68217753 10.83450006 10.90310223
 10.88618226 10.79495125 10.65179132 10.48657821 10.33177284 10.21708802
 10.16460679 10.18515415 10.27651863 10.42381556 10.60193262 10.77965572
 10.92479735 11.00948495]

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.00400327 10.87007501 10.62681585 10.31778804 10.0003348   9.73154074
  9.55425552  9.4866028   9.51754288  9.60957567]

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");
../../../_images/examples_notebooks_generated_predict_12_0.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.118340
x1                  0.485647
np.sin(x1)          0.487443
I((x1 - 5) ** 2)   -0.018964
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.004003
1    10.870075
2    10.626816
3    10.317788
4    10.000335
5     9.731541
6     9.554256
7     9.486603
8     9.517543
9     9.609576
dtype: float64