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.983
Model:                            OLS   Adj. R-squared:                  0.982
Method:                 Least Squares   F-statistic:                     902.9
Date:                Fri, 17 Dec 2021   Prob (F-statistic):           7.22e-41
Time:                        22:32:59   Log-Likelihood:                 2.4641
No. Observations:                  50   AIC:                             3.072
Df Residuals:                      46   BIC:                             10.72
Df Model:                           3
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.2118      0.082     63.674      0.000       5.047       5.377
x1             0.4768      0.013     37.770      0.000       0.451       0.502
x2             0.4920      0.050      9.916      0.000       0.392       0.592
x3            -0.0189      0.001    -17.027      0.000      -0.021      -0.017
==============================================================================
Omnibus:                        0.349   Durbin-Watson:                   2.108
Prob(Omnibus):                  0.840   Jarque-Bera (JB):                0.522
Skew:                          -0.074   Prob(JB):                        0.770
Kurtosis:                       2.522   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.7399829   5.20377305  5.6291874   5.98940991  6.26730232  6.45821989
  6.57077456  6.62541997  6.65109102  6.68044998  6.74452009  6.86758861
  7.06321645  7.33201015  7.66152225  8.02829636  8.40172095  8.74905881
  9.0408256   9.2556341   9.38370803  9.42848831  9.40606768  9.34254674
  9.26974569  9.21997634  9.2207334   9.29017716  9.43414981  9.64521559
  9.90388245 10.18180456 10.44643941 10.66639452 10.81658507 10.88235528
 10.86188635 10.76649525 10.61877536 10.44888449 10.28958999 10.1708852
 10.11506161 10.13304644 10.22260732 10.36871906 10.54603229 10.72303821
 10.86724547 10.95051836]

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.94239571 10.80474951 10.5568759  10.24274839  9.92025154  9.64700887
  9.46627458  9.3963434   9.42607143  9.51760456]

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 0x7f9be25b23d0>
../../../_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.211765
x1                  0.476779
np.sin(x1)          0.492045
I((x1 - 5) ** 2)   -0.018871
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.942396
1    10.804750
2    10.556876
3    10.242748
4     9.920252
5     9.647009
6     9.466275
7     9.396343
8     9.426071
9     9.517605
dtype: float64