Python Statsmodels 統(tǒng)計(jì)包之 OLS 回歸JoinQuant 1 年前 Statsmodels 是 Python 中一個(gè)強(qiáng)大的統(tǒng)計(jì)分析包，包含了回歸分析、時(shí)間序列分析、假設(shè)檢 驗(yàn)等等的功能。Statsmodels 在計(jì)量的簡(jiǎn)便性上是遠(yuǎn)遠(yuǎn)不及 Stata 等軟件的，但它的優(yōu)點(diǎn)在于可以與 Python 的其他的任務(wù)（如 NumPy、Pandas）有效結(jié)合，提高工作效率。在本文中，我們重點(diǎn)介紹最回歸分析中最常用的 OLS（ordinary least square）功能。 當(dāng)你需要在 Python 中進(jìn)行回歸分析時(shí)…… import statsmodels.api as sm！?。?/p> 在一切開(kāi)始之前 上帝導(dǎo)入了 NumPy（大家都叫它囊派？我叫它囊辟）， import numpy as np 便有了時(shí)間。 上帝導(dǎo)入了 matplotlib， import matplotlib.pyplot as plt 便有了空間。 上帝導(dǎo)入了 Statsmodels， import statsmodels.api as sm 世界開(kāi)始了。 簡(jiǎn)單 OLS 回歸 假設(shè)我們有回歸模型 Y=β0+β1X1+?+βnXn+ε, 并且有 k 組數(shù)據(jù) 。OLS 回歸用于計(jì)算回歸系數(shù) βi 的估值 b0,b1,…,bn，使誤差平方 最小化。 statsmodels.OLS 的輸入有 (endog, exog, missing, hasconst) 四個(gè)，我們現(xiàn)在只考慮前兩個(gè)。第一個(gè)輸入 endog 是回歸中的反應(yīng)變量（也稱(chēng)因變量），是上面模型中的 y(t), 輸入是一個(gè)長(zhǎng)度為 k 的 array。第二個(gè)輸入 exog 則是回歸變量（也稱(chēng)自變量）的值，即模型中的x1(t),…,xn(t)。但是要注意，statsmodels.OLS 不會(huì)假設(shè)回歸模型有常數(shù)項(xiàng)，所以我們應(yīng)該假設(shè)模型是 并且在數(shù)據(jù)中，對(duì)于所有 t=1,…,k，設(shè)置 x0(t)=1。因此，exog的輸入是一個(gè) k×(n+1) 的 array，其中最左一列的數(shù)值全為 1。往往輸入的數(shù)據(jù)中，沒(méi)有專(zhuān)門(mén)的數(shù)值全為1的一列，Statmodels 有直接解決這個(gè)問(wèn)題的函數(shù)：sm.add_constant()。它會(huì)在一個(gè) array 左側(cè)加上一列 1。（本文中所有輸入 array 的情況也可以使用同等的 list、pd.Series 或 pd.DataFrame。） 確切地說(shuō)，statsmodels.OLS 是 statsmodels.regression.linear_model 里的一個(gè)函數(shù)（從這個(gè)命名也能看出，statsmodel 有很多很多功能，其中的一項(xiàng)叫回歸）。它的輸出結(jié)果是一個(gè) statsmodels.regression.linear_model.OLS，只是一個(gè)類(lèi)，并沒(méi)有進(jìn)行任何運(yùn)算。在 OLS 的模型之上調(diào)用擬合函數(shù) fit()，才進(jìn)行回歸運(yùn)算，并且得到 statsmodels.regression.linear_model.RegressionResultsWrapper，它包含了這組數(shù)據(jù)進(jìn)行回歸擬合的結(jié)果摘要。調(diào)用 params 可以查看計(jì)算出的回歸系數(shù) b0,b1,…,bn。 簡(jiǎn)單的線(xiàn)性回歸 上面的介紹繞了一個(gè)大圈圈，現(xiàn)在我們來(lái)看一個(gè)例子，假設(shè)回歸公式是: 我們從最簡(jiǎn)單的一元模型開(kāi)始，虛構(gòu)一組數(shù)據(jù)。首先設(shè)定數(shù)據(jù)量，也就是上面的 k 值。 nsample = 100 然后創(chuàng)建一個(gè) array，是上面的 x1 的數(shù)據(jù)。這里，我們想要 x1 的值從 0 到 10 等差排列。 x = np.linspace(0, 10, nsample) 使用 sm.add_constant() 在 array 上加入一列常項(xiàng)1。 X = sm.add_constant(x) 然后設(shè)置模型里的 β0,β1，這里要設(shè)置成 1,10。 beta = np.array([1, 10]) 然后還要在數(shù)據(jù)中加上誤差項(xiàng)，所以生成一個(gè)長(zhǎng)度為k的正態(tài)分布樣本。 e = np.random.normal(size=nsample) 由此，我們生成反應(yīng)項(xiàng) y(t)。 y = np.dot(X, beta) + e 好嘞，在反應(yīng)變量和回歸變量上使用 OLS() 函數(shù)。 model = sm.OLS(y,X) 然后獲取擬合結(jié)果。 results = model.fit() 再調(diào)取計(jì)算出的回歸系數(shù)。 print(results.params) 得到 [ 1.04510666, 9.97239799] 和實(shí)際的回歸系數(shù)非常接近。 當(dāng)然，也可以將回歸擬合的摘要全部打印出來(lái)。 print(results.summary()) 得到 中間偏下的 coef 列就是計(jì)算出的回歸系數(shù)。 我們還可以將擬合結(jié)果畫(huà)出來(lái)。先調(diào)用擬合結(jié)果的 fittedvalues 得到擬合的 y 值。 y_fitted = results.fittedvalues 然后使用 matplotlib.pyploft 畫(huà)圖。首先設(shè)定圖軸，圖片大小為 8×6。 fig, ax = plt.subplots(figsize=(8,6)) 畫(huà)出原數(shù)據(jù)，圖像為圓點(diǎn)，默認(rèn)顏色為藍(lán)。 ax.plot(x, y, 'o', label='data') 畫(huà)出擬合數(shù)據(jù)，圖像為紅色帶點(diǎn)間斷線(xiàn)。 ax.plot(x, y_fitted, 'r--.',label='OLS') 放置注解。 ax.legend(loc='best') 得到 在大圖中看不清細(xì)節(jié)，我們?cè)?0 到 2 的區(qū)間放大一下，可以見(jiàn)數(shù)據(jù)和擬合的關(guān)系。 加入改變坐標(biāo)軸區(qū)間的指令 ax.axis((-0.05, 2, -1, 25)) 高次模型的回歸 假設(shè)反應(yīng)變量 Y 和回歸變量 X 的關(guān)系是高次的多項(xiàng)式，即 我們依然可以使用 OLS 進(jìn)行線(xiàn)性回歸。但前提條件是，我們必須知道 X 在這個(gè)關(guān)系中的所有次方數(shù)；比如，如果這個(gè)公式里有一個(gè) （x）^2.5項(xiàng)，但我們對(duì)此并不知道，那么用線(xiàn)性回歸的方法就不能得到準(zhǔn)確的擬合。 雖然 X 和 Y 的關(guān)系不是線(xiàn)性的，但是 Y 和
的關(guān)系是高元線(xiàn)性的。也就是說(shuō)，只要我們把高次項(xiàng)當(dāng)做其他的自變量，即設(shè)
那么，對(duì)于線(xiàn)性公式 可以進(jìn)行常規(guī)的 OLS 回歸，估測(cè)出每一個(gè)回歸系數(shù) βi?？梢岳斫鉃榘岩辉蔷€(xiàn)性的問(wèn)題映射到高元，從而變成一個(gè)線(xiàn)性關(guān)系。 下面以
為例做一次演示。 首先設(shè)定數(shù)據(jù)量，也就是上面的 k 值。 nsample = 100 然后創(chuàng)建一個(gè) array，是上面的 x1 的數(shù)據(jù)。這里，我們想要 x1 的值從 0 到 10 等差排列。 x = np.linspace(0, 10, nsample) 再創(chuàng)建一個(gè) k×2 的 array，兩列分別為 x1 和 x2。我們需要 x2 為 x1 的平方。 X = np.column_stack((x, x**2)) 使用 sm.add_constant() 在 array 上加入一列常項(xiàng) 1。 X = sm.add_constant(X) 然后設(shè)置模型里的 β0,β1,β2，我們想設(shè)置成 1,0.1,10。 beta = np.array([1, 0.1, 10]) 然后還要在數(shù)據(jù)中加上誤差項(xiàng)，所以生成一個(gè)長(zhǎng)度為k的正態(tài)分布樣本。 e = np.random.normal(size=nsample) 由此，我們生成反應(yīng)項(xiàng) y(t)，它與 x1(t) 是二次多項(xiàng)式關(guān)系。 y = np.dot(X, beta) + e 在反應(yīng)變量和回歸變量上使用 OLS() 函數(shù)。 model = sm.OLS(y,X) 然后獲取擬合結(jié)果。 results = model.fit() 再調(diào)取計(jì)算出的回歸系數(shù)。 print(results.params) 得到 [ 0.95119465, 0.10235581, 9.9998477] 獲取全部摘要 print(results.summary()) 得到
擬合結(jié)果圖如下
在橫軸的 [0,2] 區(qū)間放大，可以看到
啞變量 一般而言，有連續(xù)取值的變量叫做連續(xù)變量，它們的取值可以是任何的實(shí)數(shù)，或者是某一區(qū)間里的任何實(shí)數(shù)，比如股價(jià)、時(shí)間、身高。但有些性質(zhì)不是連續(xù)的，只有有限個(gè)取值的可能性，一般是用于分辨類(lèi)別，比如性別、婚姻情況、股票所屬行業(yè)，表達(dá)這些變量叫做分類(lèi)變量。在回歸分析中，我們需要將分類(lèi)變量轉(zhuǎn)化為啞變量(dummy variable)。 如果我們想表達(dá)一個(gè)有 d 種取值的分類(lèi)變量，那么它所對(duì)應(yīng)的啞變量的取值是一個(gè) d 元組（可以看成一個(gè)長(zhǎng)度為 d 的向量），其中有一個(gè)元素為 1，其他都是 0。元素呈現(xiàn)出 1 的位置就是變量所取的類(lèi)別。比如說(shuō)，某個(gè)分類(lèi)變量的取值是 {a,b,c,d}，那么類(lèi)別 a 對(duì)應(yīng)的啞變量是(1,0,0,0)，b 對(duì)應(yīng) (0,1,0,0)，c 對(duì)應(yīng) (0,0,1,0)，d 對(duì)應(yīng) (0,0,0,1)。這么做的用處是，假如 a、b、c、d 四種情況分別對(duì)應(yīng)四個(gè)系數(shù) β0,β1,β2,β3，設(shè) (x0,x1,x2,x3) 是一個(gè)取值所對(duì)應(yīng)的啞變量，那么
可以直接得出相應(yīng)的系數(shù)?？梢岳斫鉃?，分類(lèi)變量的取值本身只是分類(lèi)，無(wú)法構(gòu)成任何線(xiàn)性關(guān)系，但是若映射到高元的 0,1 點(diǎn)上，便可以用線(xiàn)性關(guān)系表達(dá)，從而進(jìn)行回歸。 Statsmodels 里有一個(gè)函數(shù) categorical() 可以直接把類(lèi)別 {0,1,…,d-1} 轉(zhuǎn)換成所對(duì)應(yīng)的元組。確切地說(shuō)，sm.categorical() 的輸入有 (data, col, dictnames, drop) 四個(gè)。其中，data 是一個(gè) k×1 或 k×2 的 array，其中記錄每一個(gè)樣本的分類(lèi)變量取值。drop 是一個(gè) Bool值，意義為是否在輸出中丟掉樣本變量的值。中間兩個(gè)輸入可以不用在意。這個(gè)函數(shù)的輸出是一個(gè)k×d 的 array（如果 drop=False，則是k×(d+1)），其中每一行是所對(duì)應(yīng)的樣本的啞變量；這里 d 是 data 中分類(lèi)變量的類(lèi)別總數(shù)。 我們來(lái)舉一個(gè)例子。這里假設(shè)一個(gè)反應(yīng)變量 Y 對(duì)應(yīng)連續(xù)自變量 X 和一個(gè)分類(lèi)變量 Z。常項(xiàng)系數(shù)為 10，XX 的系數(shù)為 1；Z 有 {a,b,c}三個(gè)種類(lèi)，其中 a 類(lèi)有系數(shù) 1，b 類(lèi)有系數(shù) 3，c 類(lèi)有系數(shù) 8。也就是說(shuō)，將 Z 轉(zhuǎn)換為啞變量 (Z1,Z2,Z3)，其中 Zi 取值于 0,1，有線(xiàn)性公式
可以用常規(guī)的方法進(jìn)行 OLS 回歸。 我們按照這個(gè)關(guān)系生成一組數(shù)據(jù)來(lái)做一次演示。先定義樣本數(shù)量為 50。 nsample = 50 設(shè)定分類(lèi)變量的 array。前 20 個(gè)樣本分類(lèi)為 a。 groups = np.zeros(nsample, int) 之后的 20 個(gè)樣本分類(lèi)為 b。 groups[20:40] = 1 最后 10 個(gè)是 c 類(lèi)。 groups[40:] = 2 轉(zhuǎn)變成啞變量。 dummy = sm.categorical(groups, drop=True) 創(chuàng)建一組連續(xù)變量，是 50 個(gè)從 0 到 20 遞增的值。 x = np.linspace(0, 20, nsample) 將連續(xù)變量和啞變量的 array 合并，并加上一列常項(xiàng)。 X = np.column_stack((x, dummy))X = sm.add_constant(X) 定義回歸系數(shù)。我們想設(shè)定常項(xiàng)系數(shù)為 10，唯一的連續(xù)變量的系數(shù)為 1，并且分類(lèi)變量的三種分類(lèi) a、b、c 的系數(shù)分別為 1,3,8。 beta = [10, 1, 1, 3, 8] 再生成一個(gè)正態(tài)分布的噪音樣本。 e = np.random.normal(size=nsample) 最后，生成反映變量。 y = np.dot(X, beta) + e 得到了虛構(gòu)數(shù)據(jù)后，放入 OLS 模型并進(jìn)行擬合運(yùn)算。 result = sm.OLS(y,X).fit()print(result.summary()) 得到
再畫(huà)圖出來(lái) 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') 得出
這里要指出，啞變量是和其他自變量并行的影響因素，也就是說(shuō)，啞變量和原先的 x 同時(shí)影響了回歸的結(jié)果。初學(xué)者往往會(huì)誤解這一點(diǎn)，認(rèn)為啞變量是一個(gè)選擇變量：也就是說(shuō)，上圖中給出的回歸結(jié)果，是在只做了一次回歸的情況下完成的，而不是分成３段進(jìn)行３次回歸。啞變量的取值藏在其他的三個(gè)維度中。可以理解成：上圖其實(shí)是將高元的回歸結(jié)果映射到平面上之后得到的圖。 簡(jiǎn)單應(yīng)用 我們來(lái)做一個(gè)非常簡(jiǎn)單的實(shí)際應(yīng)用。設(shè) x 為上證指數(shù)的日收益率，y 為深證成指的日收益率。通過(guò)對(duì)股票市場(chǎng)的認(rèn)知，我們認(rèn)為 x 和 y 有很強(qiáng)的線(xiàn)性關(guān)系。因此可以假設(shè)模型
并使用 Statsmodels 包進(jìn)行 OLS 回歸分析。 我們?nèi)∩献C指數(shù)和深證成指一年中的收盤(pán)價(jià)。 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'].valuesy_price = data['399001.XSHE'].values 計(jì)算兩個(gè)指數(shù)一年內(nèi)的日收益率，記載于 x_pct 和 y_pct 兩個(gè) 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ù)項(xiàng)。 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。最后，畫(huà)出數(shù)據(jù)和擬合線(xiàn)。 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é)語(yǔ) 本篇文章中，我們介紹了 Statsmodels 中很常用 OLS 回歸功能，并展示了一些使用方法。線(xiàn)性回歸的應(yīng)用場(chǎng)景非常廣泛。在我們量化課堂應(yīng)用類(lèi)的內(nèi)容中，也有相當(dāng)多的策略?xún)?nèi)容采用線(xiàn)性回歸的內(nèi)容。我們會(huì)將應(yīng)用類(lèi)文章中涉及線(xiàn)性回歸的部分加上鏈接，鏈接到本篇文章中來(lái)，形成體系。量化課堂在未來(lái)還會(huì)介紹 Statsmodel 包其他的一些功能，敬請(qǐng)期待。 |
|
|
來(lái)自： 習(xí)悟齋 > 《編程開(kāi)發(fā)》
