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.985
Model:                            OLS   Adj. R-squared:                  0.984
Method:                 Least Squares   F-statistic:                     1000.
Date:                Fri, 04 Aug 2023   Prob (F-statistic):           7.14e-42
Time:                        17:42:09   Log-Likelihood:                 3.7584
No. Observations:                  50   AIC:                            0.4831
Df Residuals:                      46   BIC:                             8.131
Df Model:                           3
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.1011      0.080     63.956      0.000       4.941       5.262
x1             0.4930      0.012     40.077      0.000       0.468       0.518
x2             0.4562      0.048      9.434      0.000       0.359       0.554
x3            -0.0198      0.001    -18.364      0.000      -0.022      -0.018
==============================================================================
Omnibus:                        2.512   Durbin-Watson:                   2.663
Prob(Omnibus):                  0.285   Jarque-Bera (JB):                2.282
Skew:                          -0.513   Prob(JB):                        0.320
Kurtosis:                       2.789   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.60522507  5.06516918  5.48875385  5.85111604  6.13636566  6.34019625
  6.47059253  6.54651852  6.59480183  6.64572595  6.72805472  6.86430659
  7.06705483  7.33686181  7.66218661  8.02128113  8.38576329  8.72528008
  9.01249449  9.22757691  9.36146306  9.41734326  9.41013824  9.36404769
  9.30857439  9.2736769   9.28484752  9.35892392  9.50132283  9.70515001
  9.95233302 10.21659042 10.46774975 10.67670483 10.82019772 10.88464005
 10.86834545 10.78180669 10.6459717  10.48880191 10.34067826 10.2294097
 10.17566423 10.18957326 10.26906701 10.40021456 10.5595127  10.71774741
 10.84479364 10.9145648 ]

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.89686239 10.75790128 10.51557229 10.21064636  9.89679238  9.62743718
  9.44268479  9.35949757  9.3675432   9.43172446]

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 0x7f31e9500850>
../../../_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.101060
x1                  0.492976
np.sin(x1)          0.456210
I((x1 - 5) ** 2)   -0.019833
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.896862
1    10.757901
2    10.515572
3    10.210646
4     9.896792
5     9.627437
6     9.442685
7     9.359498
8     9.367543
9     9.431724
dtype: float64