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.985
Model:                            OLS   Adj. R-squared:                  0.984
Method:                 Least Squares   F-statistic:                     1004.
Date:                Fri, 10 Jul 2020   Prob (F-statistic):           6.49e-42
Time:                        05:46:36   Log-Likelihood:                 4.1649
No. Observations:                  50   AIC:                           -0.3297
Df Residuals:                      46   BIC:                             7.318
Df Model:                           3
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          5.0480      0.079     63.808      0.000       4.889       5.207
x1             0.4969      0.012     40.728      0.000       0.472       0.521
x2             0.4630      0.048      9.653      0.000       0.366       0.560
x3            -0.0206      0.001    -19.190      0.000      -0.023      -0.018
==============================================================================
Omnibus:                        2.802   Durbin-Watson:                   2.088
Prob(Omnibus):                  0.246   Jarque-Bera (JB):                2.620
Skew:                           0.543   Prob(JB):                        0.270
Kurtosis:                       2.721   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.53407269  5.0011611   5.43120627  5.79897547  6.08834231  6.29493629
  6.42686088  6.50336208  6.55166641  6.60250768  6.68507745  6.82222899
  7.02672255  7.29912885  7.62773537  7.99047062  8.35853048  8.7011106
  8.99046747  9.20647663  9.33993901  9.39409238  9.38407931  9.33445937
  9.27517401  9.23662728  9.24469053  9.31645182  9.45740842  9.66056358
  9.90757611 10.17177377 10.42253548 10.63032221 10.77153022 10.83236906
 10.81112729 10.7184537  10.57560773 10.41096657 10.25536293 10.13701925
 10.07691089 10.08531989 10.16014569 10.28725024 10.44278082 10.59708902
 10.71960184 10.78384478]

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.75764109 10.60816081 10.35356079 10.03521823  9.70760006  9.42492756
  9.22790109  9.13373523  9.13194404  9.1869084 ]

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.048021
x1                  0.496934
np.sin(x1)          0.462994
I((x1 - 5) ** 2)   -0.020558
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.757641
1    10.608161
2    10.353561
3    10.035218
4     9.707600
5     9.424928
6     9.227901
7     9.133735
8     9.131944
9     9.186908
dtype: float64