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.980
Method:                 Least Squares   F-statistic:                     783.3
Date:                Sun, 17 Jan 2021   Prob (F-statistic):           1.78e-39
Time:                        14:02:15   Log-Likelihood:                -2.4471
No. Observations:                  50   AIC:                             12.89
Df Residuals:                      46   BIC:                             20.54
Df Model:                           3
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          4.9938      0.090     55.304      0.000       4.812       5.176
x1             0.5017      0.014     36.028      0.000       0.474       0.530
x2             0.4618      0.055      8.435      0.000       0.352       0.572
x3            -0.0208      0.001    -17.030      0.000      -0.023      -0.018
==============================================================================
Omnibus:                       11.511   Durbin-Watson:                   2.007
Prob(Omnibus):                  0.003   Jarque-Bera (JB):               11.423
Skew:                          -1.054   Prob(JB):                      0.00331
Kurtosis:                       4.020   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.47322233  4.94282036  5.37536716  5.74569718  6.03772697  6.24709769
  6.38189119  6.46130226  6.51248492  6.56609111  6.65123469  6.79070806
  6.9972375   7.27139221  7.60149085  7.96552071  8.33475437  8.6784695
  8.96899648  9.18626441  9.32109874  9.37672859  9.36825625  9.32017586
  9.26234916  9.22509922  9.23422852  9.30677962  9.44823513  9.65161695
  9.89863272 10.16268123 10.41322292 10.62079744 10.76186389 10.82266853
 10.80150451 10.70899223 10.56633454 10.40183306 10.2462384  10.12769785
 10.06713091 10.07479226 10.14858691 10.27441447 10.42848574 10.58123125
 10.70215941 10.76486621]

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.73691476 10.58585495 10.32979531 10.01000284  9.68079944  9.39626201
  9.19698254  9.10012976  9.09524566  9.14680596]

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.993803
x1                  0.501736
np.sin(x1)          0.461761
I((x1 - 5) ** 2)   -0.020823
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.736915
1    10.585855
2    10.329795
3    10.010003
4     9.680799
5     9.396262
6     9.196983
7     9.100130
8     9.095246
9     9.146806
dtype: float64