Prediction (out of sample)

[1]:
%matplotlib inline
[2]:
import numpy as np
import statsmodels.api as sm

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.984
Model:                            OLS   Adj. R-squared:                  0.983
Method:                 Least Squares   F-statistic:                     950.4
Date:                Mon, 25 Apr 2022   Prob (F-statistic):           2.26e-41
Time:                        14:03:04   Log-Likelihood:                 3.0632
No. Observations:                  50   AIC:                             1.874
Df Residuals:                      46   BIC:                             9.522
Df Model:                           3
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.1691      0.081     63.914      0.000       5.006       5.332
x1             0.4801      0.012     38.493      0.000       0.455       0.505
x2             0.4910      0.049     10.013      0.000       0.392       0.590
x3            -0.0187      0.001    -17.115      0.000      -0.021      -0.017
==============================================================================
Omnibus:                        0.338   Durbin-Watson:                   1.657
Prob(Omnibus):                  0.844   Jarque-Bera (JB):                0.468
Skew:                          -0.171   Prob(JB):                        0.791
Kurtosis:                       2.672   Cond. No.                         221.
==============================================================================

Warnings:
[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.70052172  5.16475425  5.59072261  5.9516683   6.23048981  6.42255237
  6.53644942  6.59259076  6.61984923  6.65081704  6.7164509   6.84098599
  7.03795412  7.30796055  7.6385845   8.00641989  8.38092085  8.72942037
  9.02249734  9.23881042  9.36860429  9.4153125   9.39499327  9.33369123
  9.26315831  9.21563705  9.21856319  9.29005806  9.43595133  9.64882317
  9.90922336 10.18886683 10.4552809  10.6771404  10.82941429 10.89747826
 10.87951739 10.78682437 10.64194414 10.47496983 10.3185985  10.20275918
 10.14969555 10.17031138 10.2623791  10.41090575 10.59059631 10.77000957
 10.91672385 11.0026638 ]

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.99808449 10.86430345 10.62057544 10.31077966  9.99267653  9.72376598
  9.54720944  9.48126241  9.5148049   9.61006397]

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.169090
x1                  0.480122
np.sin(x1)          0.490990
I((x1 - 5) ** 2)   -0.018743
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.998084
1    10.864303
2    10.620575
3    10.310780
4     9.992677
5     9.723766
6     9.547209
7     9.481262
8     9.514805
9     9.610064
dtype: float64