Statsmodels 是 Python 中一個強大的統(tǒng)計分析包冗疮,包含了回歸分析、時間序列分析檩帐、假設(shè)檢
驗等等的功能术幔。Statsmodels 在計量的簡便性上是遠遠不及 Stata 等軟件的,但它的優(yōu)點在于可以與 Python 的其他的任務(wù)(如 NumPy湃密、Pandas)有效結(jié)合诅挑,提高工作效率。在本文中泛源,我們重點介紹最回歸分析中最常用的 OLS(ordinary least square)功能拔妥。
當(dāng)你需要在 Python 中進行回歸分析時……
import statsmodels.api as sm!4锕俊没龙!
在一切開始之前
上帝導(dǎo)入了 NumPy(大家都叫它囊派?我叫它囊辟),
import numpy as np
便有了時間兜畸。
上帝導(dǎo)入了 matplotlib努释,
import matplotlib.pyplot as plt
便有了空間碘梢。
上帝導(dǎo)入了 Statsmodels咬摇,
import statsmodels.api as sm
世界開始了。
簡單 OLS 回歸
假設(shè)我們有回歸模型
Y=β0+β1X1+?+βnXn+ε,
并且有 k 組數(shù)據(jù) [圖片上傳失敗...(image-9b948c-1585191416601)]
煞躬。OLS 回歸用于計算回歸系數(shù) βi 的估值 b0,b1,…,bn肛鹏,使誤差平方
最小化。
statsmodels.OLS 的輸入有 (endog, exog, missing, hasconst) 四個恩沛,我們現(xiàn)在只考慮前兩個在扰。第一個輸入 endog 是回歸中的反應(yīng)變量(也稱因變量),是上面模型中的 y(t), 輸入是一個長度為 k 的 array雷客。第二個輸入 exog 則是回歸變量(也稱自變量)的值芒珠,即模型中的x1(t),…,xn(t)。但是要注意搅裙,statsmodels.OLS 不會假設(shè)回歸模型有常數(shù)項皱卓,所以我們應(yīng)該假設(shè)模型是
并且在數(shù)據(jù)中,對于所有 t=1,…,k部逮,設(shè)置 x0(t)=1娜汁。因此,exog的輸入是一個 k×(n+1) 的 array兄朋,其中最左一列的數(shù)值全為 1掐禁。往往輸入的數(shù)據(jù)中,沒有專門的數(shù)值全為1的一列颅和,Statmodels 有直接解決這個問題的函數(shù):sm.add_constant()傅事。它會在一個 array 左側(cè)加上一列 1。(本文中所有輸入 array 的情況也可以使用同等的 list峡扩、pd.Series 或 pd.DataFrame享完。)
確切地說,statsmodels.OLS 是 statsmodels.regression.linear_model 里的一個函數(shù)(從這個命名也能看出有额,statsmodel 有很多很多功能般又,其中的一項叫回歸)。它的輸出結(jié)果是一個 statsmodels.regression.linear_model.OLS巍佑,只是一個類茴迁,并沒有進行任何運算。在 OLS 的模型之上調(diào)用擬合函數(shù) fit()萤衰,才進行回歸運算堕义,并且得到 statsmodels.regression.linear_model.RegressionResultsWrapper,它包含了這組數(shù)據(jù)進行回歸擬合的結(jié)果摘要。調(diào)用 params 可以查看計算出的回歸系數(shù) b0,b1,…,bn倦卖。
簡單的線性回歸
上面的介紹繞了一個大圈圈洒擦,現(xiàn)在我們來看一個例子,假設(shè)回歸公式是:
我們從最簡單的一元模型開始怕膛,虛構(gòu)一組數(shù)據(jù)熟嫩。首先設(shè)定數(shù)據(jù)量,也就是上面的 k 值褐捻。
nsample = 100
然后創(chuàng)建一個 array掸茅,是上面的 x1 的數(shù)據(jù)。這里柠逞,我們想要 x1 的值從 0 到 10 等差排列昧狮。
x = np.linspace(0, 10, nsample)
使用 sm.add_constant() 在 array 上加入一列常項1。
X = sm.add_constant(x)
然后設(shè)置模型里的 β0,β1板壮,這里要設(shè)置成 1,10逗鸣。
beta = np.array([1, 10])
然后還要在數(shù)據(jù)中加上誤差項,所以生成一個長度為k的正態(tài)分布樣本绰精。
e = np.random.normal(size=nsample)
由此撒璧,我們生成反應(yīng)項 y(t)。
y = np.dot(X, beta) + e
好嘞茬底,在反應(yīng)變量和回歸變量上使用 OLS() 函數(shù)沪悲。
model = sm.OLS(y,X)
然后獲取擬合結(jié)果。
results = model.fit()
再調(diào)取計算出的回歸系數(shù)阱表。
print(results.params)
得到
[ 1.04510666, 9.97239799]
和實際的回歸系數(shù)非常接近殿如。
當(dāng)然,也可以將回歸擬合的摘要全部打印出來最爬。
print(results.summary())
得到
中間偏下的 coef 列就是計算出的回歸系數(shù)涉馁。
我們還可以將擬合結(jié)果畫出來。先調(diào)用擬合結(jié)果的 fittedvalues 得到擬合的 y 值爱致。
y_fitted = results.fittedvalues
然后使用 matplotlib.pyploft 畫圖烤送。首先設(shè)定圖軸,圖片大小為 8×6糠悯。
fig, ax = plt.subplots(figsize=(8,6))
畫出原數(shù)據(jù)帮坚,圖像為圓點,默認顏色為藍互艾。
ax.plot(x, y, 'o', label='data')
畫出擬合數(shù)據(jù)试和,圖像為紅色帶點間斷線。
ax.plot(x, y_fitted, 'r--.',label='OLS')
放置注解纫普。
ax.legend(loc='best')
得到
在大圖中看不清細節(jié)阅悍,我們在 0 到 2 的區(qū)間放大一下,可以見數(shù)據(jù)和擬合的關(guān)系。
加入改變坐標軸區(qū)間的指令
ax.axis((-0.05, 2, -1, 25))
高次模型的回歸
假設(shè)反應(yīng)變量 Y 和回歸變量 X 的關(guān)系是高次的多項式节视,即
我們依然可以使用 OLS 進行線性回歸拳锚。但前提條件是,我們必須知道 X 在這個關(guān)系中的所有次方數(shù)寻行;比如霍掺,如果這個公式里有一個 [圖片上傳失敗...(image-94ece6-1585191416601)]
.5項,但我們對此并不知道寡痰,那么用線性回歸的方法就不能得到準確的擬合抗楔。
雖然 X 和 Y 的關(guān)系不是線性的棋凳,但是 Y 和 [圖片上傳失敗...(image-fd911a-1585191416601)]
的關(guān)系是高元線性的拦坠。也就是說,只要我們把高次項當(dāng)做其他的自變量剩岳,即設(shè) [圖片上傳失敗...(image-af157-1585191416601)]
贞滨。那么,對于線性公式
可以進行常規(guī)的 OLS 回歸拍棕,估測出每一個回歸系數(shù) βi晓铆。可以理解為把一元非線性的問題映射到高元绰播,從而變成一個線性關(guān)系骄噪。
下面以
為例做一次演示。
首先設(shè)定數(shù)據(jù)量蠢箩,也就是上面的 k 值链蕊。
nsample = 100
然后創(chuàng)建一個 array,是上面的 x1 的數(shù)據(jù)谬泌。這里滔韵,我們想要 x1 的值從 0 到 10 等差排列。
x = np.linspace(0, 10, nsample)
再創(chuàng)建一個 k×2 的 array掌实,兩列分別為 x1 和 x2陪蜻。我們需要 x2 為 x1 的平方。
X = np.column_stack((x, x**2))
使用 sm.add_constant() 在 array 上加入一列常項 1贱鼻。
X = sm.add_constant(X)
然后設(shè)置模型里的 β0,β1,β2宴卖,我們想設(shè)置成 1,0.1,10。
beta = np.array([1, 0.1, 10])
然后還要在數(shù)據(jù)中加上誤差項邻悬,所以生成一個長度為k的正態(tài)分布樣本症昏。
e = np.random.normal(size=nsample)
由此,我們生成反應(yīng)項 y(t)拘悦,它與 x1(t) 是二次多項式關(guān)系齿兔。
y = np.dot(X, beta) + e
在反應(yīng)變量和回歸變量上使用 OLS() 函數(shù)。
model = sm.OLS(y,X)
然后獲取擬合結(jié)果。
results = model.fit()
再調(diào)取計算出的回歸系數(shù)分苇。
print(results.params)
得到
[ 0.95119465, 0.10235581, 9.9998477]
獲取全部摘要
print(results.summary())
得到
擬合結(jié)果圖如下
在橫軸的 [0,2] 區(qū)間放大添诉,可以看到
啞變量
一般而言,有連續(xù)取值的變量叫做連續(xù)變量医寿,它們的取值可以是任何的實數(shù)栏赴,或者是某一區(qū)間里的任何實數(shù),比如股價靖秩、時間须眷、身高。但有些性質(zhì)不是連續(xù)的沟突,只有有限個取值的可能性花颗,一般是用于分辨類別,比如性別惠拭、婚姻情況扩劝、股票所屬行業(yè),表達這些變量叫做分類變量职辅。在回歸分析中棒呛,我們需要將分類變量轉(zhuǎn)化為啞變量(dummy variable)。
如果我們想表達一個有 d 種取值的分類變量域携,那么它所對應(yīng)的啞變量的取值是一個 d 元組(可以看成一個長度為 d 的向量)簇秒,其中有一個元素為 1,其他都是 0秀鞭。元素呈現(xiàn)出 1 的位置就是變量所取的類別趋观。比如說,某個分類變量的取值是 {a,b,c,d}气筋,那么類別 a 對應(yīng)的啞變量是(1,0,0,0)拆内,b 對應(yīng) (0,1,0,0),c 對應(yīng) (0,0,1,0)宠默,d 對應(yīng) (0,0,0,1)麸恍。這么做的用處是,假如 a搀矫、b抹沪、c、d 四種情況分別對應(yīng)四個系數(shù) β0,β1,β2,β3瓤球,設(shè) (x0,x1,x2,x3) 是一個取值所對應(yīng)的啞變量融欧,那么
可以直接得出相應(yīng)的系數(shù)∝韵郏可以理解為噪馏,分類變量的取值本身只是分類麦到,無法構(gòu)成任何線性關(guān)系,但是若映射到高元的 0,1 點上欠肾,便可以用線性關(guān)系表達瓶颠,從而進行回歸。
Statsmodels 里有一個函數(shù) categorical() 可以直接把類別 {0,1,…,d-1} 轉(zhuǎn)換成所對應(yīng)的元組刺桃。確切地說粹淋,sm.categorical() 的輸入有 (data, col, dictnames, drop) 四個。其中瑟慈,data 是一個 k×1 或 k×2 的 array桃移,其中記錄每一個樣本的分類變量取值。drop 是一個 Bool值葛碧,意義為是否在輸出中丟掉樣本變量的值借杰。中間兩個輸入可以不用在意。這個函數(shù)的輸出是一個k×d 的 array(如果 drop=False吹埠,則是k×(d+1))第步,其中每一行是所對應(yīng)的樣本的啞變量疮装;這里 d 是 data 中分類變量的類別總數(shù)缘琅。
我們來舉一個例子。這里假設(shè)一個反應(yīng)變量 Y 對應(yīng)連續(xù)自變量 X 和一個分類變量 Z廓推。常項系數(shù)為 10刷袍,XX 的系數(shù)為 1;Z 有 {a,b,c}三個種類樊展,其中 a 類有系數(shù) 1呻纹,b 類有系數(shù) 3,c 類有系數(shù) 8专缠。也就是說雷酪,將 Z 轉(zhuǎn)換為啞變量 (Z1,Z2,Z3),其中 Zi 取值于 0,1涝婉,有線性公式
可以用常規(guī)的方法進行 OLS 回歸哥力。
我們按照這個關(guān)系生成一組數(shù)據(jù)來做一次演示。先定義樣本數(shù)量為 50墩弯。
nsample = 50
設(shè)定分類變量的 array吩跋。前 20 個樣本分類為 a。
groups = np.zeros(nsample, int)
之后的 20 個樣本分類為 b渔工。
groups[20:40] = 1
最后 10 個是 c 類锌钮。
groups[40:] = 2
轉(zhuǎn)變成啞變量。
dummy = sm.categorical(groups, drop=True)
創(chuàng)建一組連續(xù)變量引矩,是 50 個從 0 到 20 遞增的值梁丘。
x = np.linspace(0, 20, nsample)
將連續(xù)變量和啞變量的 array 合并侵浸,并加上一列常項。
X = np.column_stack((x, dummy))
X = sm.add_constant(X)
定義回歸系數(shù)氛谜。我們想設(shè)定常項系數(shù)為 10通惫,唯一的連續(xù)變量的系數(shù)為 1,并且分類變量的三種分類 a混蔼、b履腋、c 的系數(shù)分別為 1,3,8。
beta = [10, 1, 1, 3, 8]
再生成一個正態(tài)分布的噪音樣本惭嚣。
e = np.random.normal(size=nsample)
最后遵湖,生成反映變量。
y = np.dot(X, beta) + e
得到了虛構(gòu)數(shù)據(jù)后晚吞,放入 OLS 模型并進行擬合運算延旧。
result = sm.OLS(y,X).fit()
print(result.summary())
得到
再畫圖出來
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(x, y, 'o', label="data")
ax.plot(x, result.fittedvalues, 'r--.', label="OLS")
ax.legend(loc='best')
得出
這里要指出,啞變量是和其他自變量并行的影響因素槽地,也就是說迁沫,啞變量和原先的 x 同時影響了回歸的結(jié)果。初學(xué)者往往會誤解這一點捌蚊,認為啞變量是一個選擇變量:也就是說集畅,上圖中給出的回歸結(jié)果,是在只做了一次回歸的情況下完成的缅糟,而不是分成3段進行3次回歸挺智。啞變量的取值藏在其他的三個維度中〈盎拢可以理解成:上圖其實是將高元的回歸結(jié)果映射到平面上之后得到的圖赦颇。
簡單應(yīng)用
我們來做一個非常簡單的實際應(yīng)用。設(shè) x 為上證指數(shù)的日收益率赴涵,y 為深證成指的日收益率媒怯。通過對股票市場的認知,我們認為 x 和 y 有很強的線性關(guān)系髓窜。因此可以假設(shè)模型
并使用 Statsmodels 包進行 OLS 回歸分析扇苞。
我們?nèi)∩献C指數(shù)和深證成指一年中的收盤價。
data = get_price(['000001.XSHG', '399001.XSHE'], start_date='2015-01-01', end_date='2016-01-01', frequency='daily', fields=['close'])['close']
x_price = data['000001.XSHG'].values
y_price = data['399001.XSHE'].values
計算兩個指數(shù)一年內(nèi)的日收益率纱烘,記載于 x_pct 和 y_pct 兩個 list 中杨拐。
x_pct, y_pct = [], []
for i in range(1, len(x_price)):
x_pct.append(x_price[i]/x_price[i-1]-1)
for i in range(1, len(y_price)):
y_pct.append(y_price[i]/y_price[i-1]-1)
將數(shù)據(jù)轉(zhuǎn)化為 array 的形式;不要忘記添加常數(shù)項。
x = np.array(x_pct)
X = sm.add_constant(x)
y = np.array(y_pct)
上吧,λu.λv.(sm.OLS(u,v).fit())货裹!全靠你了!
results = sm.OLS(y, X).fit()
print(results.summary())
得到
恩屋吨,y=0.002+0.9991x蜒谤,合情合理,或者干脆直接四舍五入到 y=x至扰。最后鳍徽,畫出數(shù)據(jù)和擬合線。
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(x, y, 'o', label="data")
ax.plot(x, results.fittedvalues, 'r--', label="OLS")
ax.legend(loc='best')
結(jié)語
本篇文章中敢课,我們介紹了 Statsmodels 中很常用 OLS 回歸功能阶祭,并展示了一些使用方法。線性回歸的應(yīng)用場景非常廣泛直秆。在我們量化課堂應(yīng)用類的內(nèi)容中濒募,也有相當(dāng)多的策略內(nèi)容采用線性回歸的內(nèi)容。我們會將應(yīng)用類文章中涉及線性回歸的部分加上鏈接圾结,鏈接到本篇文章中來瑰剃,形成體系。量化課堂在未來還會介紹 Statsmodel 包其他的一些功能筝野,敬請期待晌姚。
轉(zhuǎn):https://zhuanlan.zhihu.com/p/22692029?refer=JoinQuant
到JoinQuant查看代碼并與作者交流討論:【量化課堂】Statsmodels 統(tǒng)計包之 OLS 回歸