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.987
Model: OLS Adj. R-squared: 0.987
Method: Least Squares F-statistic: 1208.
Date: Mon, 24 Feb 2020 Prob (F-statistic): 9.91e-44
Time: 22:49:06 Log-Likelihood: 7.1613
No. Observations: 50 AIC: -6.323
Df Residuals: 46 BIC: 1.325
Df Model: 3
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 4.9615 0.075 66.588 0.000 4.812 5.112
x1 0.4948 0.011 43.061 0.000 0.472 0.518
x2 0.5677 0.045 12.567 0.000 0.477 0.659
x3 -0.0190 0.001 -18.864 0.000 -0.021 -0.017
==============================================================================
Omnibus: 0.123 Durbin-Watson: 1.994
Prob(Omnibus): 0.940 Jarque-Bera (JB): 0.021
Skew: -0.041 Prob(JB): 0.990
Kurtosis: 2.943 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.48569704 4.98751642 5.44597256 5.83012653 6.12020511 6.31084944
6.4119955 6.44724165 6.44997147 6.45786894 6.50672707 6.62456723
6.82703528 7.11483101 7.4735932 7.87625909 8.28751072 8.66957758
8.98844222 9.21942917 9.35125908 9.38790189 9.34792457 9.26144063
9.16516252 9.09636994 9.08678483 9.15735959 9.31483463 9.55063095
9.84225965 10.15701679 10.45735656 10.70705993 10.87718535 10.95082378
10.92587676 10.81540122 10.64546423 10.45086027 10.2693946 10.13567196
10.0754111 10.10121927 10.21052089 10.38598069 10.59835186 10.81128124
10.98728199 11.09389262]
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.0986701 10.95537979 10.68628459 10.3421188 9.98966661 9.69541101
9.50925641 9.45231034 9.51171553 9.6437979 ]
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 4.961530
x1 0.494829
np.sin(x1) 0.567696
I((x1 - 5) ** 2) -0.019033
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.098670
1 10.955380
2 10.686285
3 10.342119
4 9.989667
5 9.695411
6 9.509256
7 9.452310
8 9.511716
9 9.643798
dtype: float64
