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.983
Model:                            OLS   Adj. R-squared:                  0.982
Method:                 Least Squares   F-statistic:                     892.3
Date:                Sun, 29 Aug 2021   Prob (F-statistic):           9.44e-41
Time:                        19:09:28   Log-Likelihood:                0.13260
No. Observations:                  50   AIC:                             7.735
Df Residuals:                      46   BIC:                             15.38
Df Model:                           3
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.0656      0.086     59.069      0.000       4.893       5.238
x1             0.5009      0.013     37.870      0.000       0.474       0.527
x2             0.3695      0.052      7.107      0.000       0.265       0.474
x3            -0.0202      0.001    -17.391      0.000      -0.023      -0.018
==============================================================================
Omnibus:                        0.882   Durbin-Watson:                   2.305
Prob(Omnibus):                  0.644   Jarque-Bera (JB):                0.353
Skew:                          -0.175   Prob(JB):                        0.838
Kurtosis:                       3.217   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.5607078   4.99086943  5.39020562  5.73857891  6.02311935  6.240339
  6.39670495  6.50757685  6.5946833   6.68255198  6.79447974  6.94870506
  7.15541146  7.41505445  7.71828662  8.04749342  8.3796872   8.69028399
  8.95714258  9.16420228  9.3041216   9.37948449  9.40237579  9.39239567
  9.37343947  9.36977169  9.40203943  9.48388012  9.61968097  9.80385807
 10.02177393 10.25214228 10.4705254  10.65334914 10.78177609 10.84480055
 10.84105667 10.77904291 10.67572575 10.55375212 10.43772862 10.35017869
 10.30784209 10.3189246  10.38174981 10.48503461 10.6097429  10.73221322
 10.82804624 10.8761135 ]

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.84955707 10.7227231  10.510102   10.24471558  9.97003218  9.72932406
  9.55507289  9.46101707  9.43878811  9.45995953]

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.065595
x1                  0.500858
np.sin(x1)          0.369500
I((x1 - 5) ** 2)   -0.020195
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.849557
1    10.722723
2    10.510102
3    10.244716
4     9.970032
5     9.729324
6     9.555073
7     9.461017
8     9.438788
9     9.459960
dtype: float64