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.983
Method:                 Least Squares   F-statistic:                     970.0
Date:                Thu, 05 Nov 2020   Prob (F-statistic):           1.43e-41
Time:                        07:28:38   Log-Likelihood:                 1.5386
No. Observations:                  50   AIC:                             4.923
Df Residuals:                      46   BIC:                             12.57
Df Model:                           3
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          4.9649      0.083     59.546      0.000       4.797       5.133
x1             0.4954      0.013     38.528      0.000       0.470       0.521
x2             0.4615      0.051      9.129      0.000       0.360       0.563
x3            -0.0190      0.001    -16.796      0.000      -0.021      -0.017
==============================================================================
Omnibus:                        4.447   Durbin-Watson:                   2.049
Prob(Omnibus):                  0.108   Jarque-Bera (JB):                3.703
Skew:                           0.662   Prob(JB):                        0.157
Kurtosis:                       3.163   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.49083841  4.95047481  5.3736972   5.73535453  6.01937263  6.22139516
  6.34949933  6.42286885  6.46864209  6.51745333  6.59839963  6.73426028
  6.93775419  7.20945026  7.53767386  7.900425    8.26899278  8.61267261
  8.9038108   9.12234804  9.25911506  9.31733912  9.31211379  9.26791911
  9.21459959  9.18246084  9.19729045  9.27612116  9.42443242  9.63525015
  9.89029254 10.16297372 10.42277173 10.64024318 10.79186077 10.86387868
 10.85459102 10.77461212 10.64513255 10.49443765 10.35326029 10.24973155
 10.20475886 10.22859099 10.31913416 10.46229611 10.63430134 10.80559744
 10.94571039 11.0282514 ]

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.02549655 10.90062849 10.67174529 10.38009018 10.07995374  9.8253817
  9.65694264  9.59179541  9.61948802  9.70451667]

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.964910
x1                  0.495433
np.sin(x1)          0.461495
I((x1 - 5) ** 2)   -0.018963
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.025497
1    10.900628
2    10.671745
3    10.380090
4    10.079954
5     9.825382
6     9.656943
7     9.591795
8     9.619488
9     9.704517
dtype: float64