Prediction (out of sample)

[1]:
%matplotlib inline
[2]:
import numpy as np
import statsmodels.api as sm

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.987
Model:                            OLS   Adj. R-squared:                  0.986
Method:                 Least Squares   F-statistic:                     1159.
Date:                Wed, 22 Jan 2020   Prob (F-statistic):           2.52e-43
Time:                        18:52:17   Log-Likelihood:                 5.3942
No. Observations:                  50   AIC:                            -2.788
Df Residuals:                      46   BIC:                             4.860
Df Model:                           3
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          4.9063      0.077     63.561      0.000       4.751       5.062
x1             0.5143      0.012     43.200      0.000       0.490       0.538
x2             0.5408      0.047     11.556      0.000       0.447       0.635
x3            -0.0208      0.001    -19.855      0.000      -0.023      -0.019
==============================================================================
Omnibus:                       14.714   Durbin-Watson:                   2.165
Prob(Omnibus):                  0.001   Jarque-Bera (JB):               21.111
Skew:                          -0.945   Prob(JB):                     2.61e-05
Kurtosis:                       5.562   Cond. No.                         221.
==============================================================================

Warnings:
[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.3875047   4.89333499  5.35698133  5.74896927  6.05046156  6.25635301
  6.37610929  6.4322117   6.45646362  6.48476527  6.55121538  6.68250892
  6.89355116  7.18500877  7.54320031  7.94234404  8.34879379  8.72656692
  9.04325617  9.27535412  9.41211531  9.45732177  9.4286616   9.35482287
  9.27078027  9.21204884  9.20884898  9.28114128  9.43534711  9.66329356
  9.94355623 10.24497924 10.53179399 10.76949591 10.93051339 10.99873777
 10.97216993 10.86324875 10.6968074  10.50599338 10.32682243 10.19226127
 10.12681124 10.14248284 10.23682239 10.3933152  10.58409886 10.77454098
 10.92892917 11.01633851]

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.00235995 10.84601188 10.56850337 10.21816722  9.85862638  9.55321685
  9.34948084  9.26752672  9.29510584  9.39061143]

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           4.906340
x1                  0.514289
np.sin(x1)          0.540824
I((x1 - 5) ** 2)   -0.020753
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.002360
1    10.846012
2    10.568503
3    10.218167
4     9.858626
5     9.553217
6     9.349481
7     9.267527
8     9.295106
9     9.390611
dtype: float64