forked from statsmodels/statsmodels
-
Notifications
You must be signed in to change notification settings - Fork 0
/
tut_ols.py
107 lines (85 loc) · 2.69 KB
/
tut_ols.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
'''OLS Example with joint significance tests
'''
import numpy as np
#from scipy import stats
import statsmodels.api as sm
import matplotlib
#..matplotlib.use('Qt4Agg')#, warn=True) #for Spyder
import matplotlib.pyplot as plt
from statsmodels.sandbox.regression.predstd import wls_prediction_std
#fix a seed for these examples
np.random.seed(9876789)
# OLS non-linear curve but linear in parameters
# ---------------------------------------------
nsample = 50
sig = 0.5
x1 = np.linspace(0, 20, nsample)
X = np.c_[x1, np.sin(x1), (x1-5)**2, np.ones(nsample)]
beta = [0.5, 0.5, -0.02, 5.]
y_true = np.dot(X, beta)
y = y_true + sig * np.random.normal(size=nsample)
plt.figure()
plt.plot(x1, y, 'o', x1, y_true, 'b-')
res = sm.OLS(y, X).fit()
print res.params
print res.bse
print res.predict()
prstd, iv_l, iv_u = wls_prediction_std(res)
plt.plot(x1, res.fittedvalues, 'r--.')
plt.plot(x1, iv_u, 'r--')
plt.plot(x1, iv_l, 'r--')
#@savefig tut_ols_0.png
plt.title('blue: true, red: OLS')
print res.summary()
#OLS with dummy variables
#------------------------
sig = 1.
#suppose observations from 3 groups
xg = np.zeros(nsample, int)
xg[20:40] = 1
xg[40:] = 2
print xg
dummy = (xg[:,None] == np.unique(xg)).astype(float)
#use group 0 as benchmark
X = np.c_[x1, dummy[:,1:], np.ones(nsample)]
beta = [1., 3, -3, 10]
y_true = np.dot(X, beta)
y = y_true + sig * np.random.normal(size=nsample)
res2 = sm.OLS(y, X).fit()
print res2.params
print res2.bse
print res.predict()
prstd, iv_l, iv_u = wls_prediction_std(res2)
plt.figure()
plt.plot(x1, y, 'o', x1, y_true, 'b-')
plt.figure()
plt.plot(x1, y, 'o', x1, y_true, 'b-')
plt.plot(x1, res2.fittedvalues, 'r--.')
plt.plot(x1, iv_u, 'r--')
plt.plot(x1, iv_l, 'r--')
#@savefig tut_ols_1.png
plt.title('blue: true, red: OLS')
print res.summary()
R = [[0, 1, 0, 0],
[0, 0, 1, 0]]
# F test joint hypothesis R * beta = 0
# i.e. coefficient on both dummy variables equal zero
print res2.f_test(R)
# strongly rejected Null of identical constant in 3 groups
#..<F test: F=124.19050615860911, p=2.87411973729e-019, df_denom=46, df_num=2>
help(res2.f_test)
# t test for Null hypothesis effects of 2nd and 3rd group add to zero
R = [0, 1, -1, 0]
print res2.t_test(R)
# don't reject Null at 5% confidence level (note one sided p-value)
#..<T test: effect=1.0363792917100714, sd=0.52675137730463362, t=1.9674923243925513, p=0.027586676754860262, df_denom=46>
# OLS with small group effects
beta = [1., 0.3, -0.0, 10]
y_true = np.dot(X, beta)
y = y_true + sig * np.random.normal(size=nsample)
res3 = sm.OLS(y, X).fit()
print res3.f_test(R)
# don't reject Null of identical constant in 3 groups
#..<F test: F=1.9715385826285652, p=0.15083366806, df_denom=46, df_num=2>
#..plt.draw()
#..plt.show()