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.981
Model:                            OLS   Adj. R-squared:                  0.979
Method:                 Least Squares   F-statistic:                     773.1
Date:                Mon, 26 Oct 2020   Prob (F-statistic):           2.40e-39
Time:                        17:34:14   Log-Likelihood:                -2.7725
No. Observations:                  50   AIC:                             13.54
Df Residuals:                      46   BIC:                             21.19
Df Model:                           3
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.0829      0.091     55.925      0.000       4.900       5.266
x1             0.4980      0.014     35.528      0.000       0.470       0.526
x2             0.5271      0.055      9.566      0.000       0.416       0.638
x3            -0.0204      0.001    -16.559      0.000      -0.023      -0.018
==============================================================================
Omnibus:                        6.936   Durbin-Watson:                   1.540
Prob(Omnibus):                  0.031   Jarque-Bera (JB):                2.372
Skew:                           0.019   Prob(JB):                        0.305
Kurtosis:                       1.934   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.57338364  5.06564798  5.51674852  5.89795929  6.19092135  6.39065916
  6.50639802  6.56004837  6.58260599  6.60905947  6.6726417   6.79936975
  7.00377042  7.28649344  7.63420478  8.02177739  8.41641947  8.7830621
  9.09012096  9.31468552  9.4462832   9.4886      9.45887486  9.38506736
  9.30126409  9.24207823  9.2369628   9.30537153  9.45356272  9.67357101
  9.94451603 10.23603303 10.51326172 10.74257367 10.89709742 10.96113337
 10.93273315 10.82401978 10.65919584 10.47056711 10.29323478 10.15932842
 10.09272724 10.10513673 10.19416552 10.34371829 10.52664013 10.70917809
 10.85652662 10.93854608]

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.92096936 10.76456178 10.48999384 10.14437102  9.78970068  9.48771054
  9.28473553  9.20037344  9.22268663  9.31112486]

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.082872
x1                  0.497993
np.sin(x1)          0.527091
I((x1 - 5) ** 2)   -0.020380
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.920969
1    10.764562
2    10.489994
3    10.144371
4     9.789701
5     9.487711
6     9.284736
7     9.200373
8     9.222687
9     9.311125
dtype: float64