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.984
Model:                            OLS   Adj. R-squared:                  0.983
Method:                 Least Squares   F-statistic:                     932.3
Date:                Sun, 28 Nov 2021   Prob (F-statistic):           3.50e-41
Time:                        16:30:27   Log-Likelihood:                 1.7799
No. Observations:                  50   AIC:                             4.440
Df Residuals:                      46   BIC:                             12.09
Df Model:                           3
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.0294      0.083     60.611      0.000       4.862       5.196
x1             0.5006      0.013     39.119      0.000       0.475       0.526
x2             0.4773      0.050      9.488      0.000       0.376       0.579
x3            -0.0206      0.001    -18.319      0.000      -0.023      -0.018
==============================================================================
Omnibus:                        3.288   Durbin-Watson:                   2.127
Prob(Omnibus):                  0.193   Jarque-Bera (JB):                1.758
Skew:                           0.146   Prob(JB):                        0.415
Kurtosis:                       2.129   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.51480037  4.98917759  5.42556862  5.79795994  6.08972616  6.29636154
  6.42622029  6.49914396  6.54320157  6.59007791  6.66986786  6.80613186
  7.01202498  7.28813562  7.62238883  7.99203015  8.36736415  8.71663328
  9.01123555  9.23042378  9.36471424  9.41744488  9.40422686  9.35037969
  9.28677132  9.24474663  9.25097751  9.32308079  9.46672384  9.6746934
  9.92808075 10.19938832 10.45704748 10.67060518 10.81572727 10.87819658
 10.85624873 10.76086201 10.61395369 10.44477908 10.28512495 10.16408705
 10.10328978 10.11333336 10.19205223 10.32487085 10.48719842 10.64846909
 10.77716383 10.84598891]

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.82242939 10.67208374 10.41367066 10.0898477   9.75676721  9.47032825
  9.27249062  9.18100157  9.1850507   9.24791667]

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 0x7faf013c9c40>
../../../_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.029398
x1                  0.500610
np.sin(x1)          0.477320
I((x1 - 5) ** 2)   -0.020584
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.822429
1    10.672084
2    10.413671
3    10.089848
4     9.756767
5     9.470328
6     9.272491
7     9.181002
8     9.185051
9     9.247917
dtype: float64