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: 935.7
Date: Fri, 10 Jul 2020 Prob (F-statistic): 3.22e-41
Time: 05:46:36 Log-Likelihood: 1.9020
No. Observations: 50 AIC: 4.196
Df Residuals: 46 BIC: 11.84
Df Model: 3
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 5.0645 0.083 61.184 0.000 4.898 5.231
x1 0.4762 0.013 37.303 0.000 0.451 0.502
x2 0.4994 0.050 9.952 0.000 0.398 0.600
x3 -0.0177 0.001 -15.775 0.000 -0.020 -0.015
==============================================================================
Omnibus: 0.124 Durbin-Watson: 2.203
Prob(Omnibus): 0.940 Jarque-Bera (JB): 0.321
Skew: 0.060 Prob(JB): 0.852
Kurtosis: 2.627 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.6224677 5.0843021 5.50767542 5.86536907 6.13998753 6.32681633
6.43459668 6.48408926 6.50466334 6.52947147 6.59000261 6.71090859
6.9059539 7.17575419 7.50767513 7.87790822 8.25538243 8.6068691
8.90244128 9.12039072 9.25079438 9.29714482 9.27577636 9.21318142
9.14165785 9.09400251 9.09812284 9.17245174 9.32291911 9.54197756
9.80984244 10.0977424 10.37264644 10.60269082 10.76241417 10.83694071
10.82442433 10.73635188 10.59565592 10.43294695 10.2814842 10.17171127
10.12625425 10.15620441 10.25929587 10.42027789 10.61342056 10.80674211
10.96726342 11.06642646]
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)
[11.07678486 10.95641188 10.72489332 10.42686278 10.12107377 9.86601481
9.70558939 9.65836688 9.71303611 9.8311747 ]
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");

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.064507
x1 0.476216
np.sin(x1) 0.499431
I((x1 - 5) ** 2) -0.017682
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 11.076785
1 10.956412
2 10.724893
3 10.426863
4 10.121074
5 9.866015
6 9.705589
7 9.658367
8 9.713036
9 9.831175
dtype: float64