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.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.986
Model: OLS Adj. R-squared: 0.985
Method: Least Squares F-statistic: 1073.
Date: Sun, 17 Jan 2021 Prob (F-statistic): 1.45e-42
Time: 14:02:15 Log-Likelihood: 2.8926
No. Observations: 50 AIC: 2.215
Df Residuals: 46 BIC: 9.863
Df Model: 3
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 4.8931 0.081 60.296 0.000 4.730 5.056
x1 0.5257 0.013 42.001 0.000 0.500 0.551
x2 0.4256 0.049 8.650 0.000 0.327 0.525
x3 -0.0217 0.001 -19.711 0.000 -0.024 -0.019
==============================================================================
Omnibus: 1.103 Durbin-Watson: 2.517
Prob(Omnibus): 0.576 Jarque-Bera (JB): 0.434
Skew: -0.140 Prob(JB): 0.805
Kurtosis: 3.361 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.35159891 4.81989064 5.25321007 5.62836208 5.93052254 6.15567389
6.31126519 6.4149883 6.4918711 6.57016553 6.67670587 6.8324999
7.04927743 7.32756312 7.65659024 8.01606967 8.3795233 8.71863433
9.00789963 9.22881985 9.37293877 9.4432327 9.45362142 9.42668121
9.38993567 9.37133362 9.39465721 9.47561461 9.61925914 9.81915908
10.05845433 10.31262655 10.55352736 10.75400314 10.89235622 10.95590975
10.94309052 10.86368739 10.73724302 10.58984306 10.44983039 10.34314839
10.2890786 10.29707274 10.3651999 10.48046369 10.6209376 10.75936742
10.867649 10.92144575]
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.88966221 10.74317584 10.49867726 10.19420234 9.87981962 9.60537185
9.40827271 9.30434664 9.28395448 9.31435345]
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.893105
x1 0.525667
np.sin(x1) 0.425605
I((x1 - 5) ** 2) -0.021660
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.889662
1 10.743176
2 10.498677
3 10.194202
4 9.879820
5 9.605372
6 9.408273
7 9.304347
8 9.283954
9 9.314353
dtype: float64
