Python Statsmodels 統(tǒng)計包之 OLS 回歸JoinQuant 1 年前 Statsmodels 是 Python 中一個強(qiáng)大的統(tǒng)計分析包，包含了回歸分析、時間序列分析、假設(shè)檢 驗(yàn)等等的功能。Statsmodels 在計量的簡便性上是遠(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)行回歸分析時…… import statsmodels.api as sm?。?！ 在一切開始之前 上帝導(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ù) 。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ù)項(xiàng)，所以我們應(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 有很多很多功能，其中的一項(xiàng)叫回歸）。它的輸出結(jié)果是一個 statsmodels.regression.linear_model.OLS，只是一個類，并沒有進(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 可以查看計算出的回歸系數(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 上加入一列常項(xiàng)1。 X = sm.add_constant(x) 然后設(shè)置模型里的 β0,β1，這里要設(shè)置成 1,10。 beta = np.array([1, 10]) 然后還要在數(shù)據(jù)中加上誤差項(xià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)取計算出的回歸系數(shù)。 print(results.params) 得到 [ 1.04510666, 9.97239799] 和實(shí)際的回歸系數(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ù)，圖像為圓點(diǎn)，默認(rèn)顏色為藍(lán)。 ax.plot(x, y, 'o', label='data') 畫出擬合數(shù)據(jù)，圖像為紅色帶點(diǎn)間斷線。 ax.plot(x, y_fitted, 'r--.',label='OLS') 放置注解。 ax.legend(loc='best') 得到 在大圖中看不清細(xì)節(jié)，我們在 0 到 2 的區(qū)間放大一下，可以見數(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)行線性回歸。但前提條件是，我們必須知道 X 在這個關(guān)系中的所有次方數(shù)；比如，如果這個公式里有一個 （x）^2.5項(xiàng)，但我們對此并不知道，那么用線性回歸的方法就不能得到準(zhǔn)確的擬合。 雖然 X 和 Y 的關(guān)系不是線性的，但是 Y 和
的關(guān)系是高元線性的。也就是說，只要我們把高次項(xiàng)當(dāng)做其他的自變量，即設(shè)
那么，對于線性公式 可以進(jìn)行常規(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 上加入一列常項(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)，所以生成一個長度為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)取計算出的回歸系數(shù)。 print(results.params) 得到 [ 0.95119465, 0.10235581, 9.9998477] 獲取全部摘要 print(results.summary()) 得到
擬合結(jié)果圖如下
在橫軸的 [0,2] 區(qū)間放大，可以看到
啞變量 一般而言，有連續(xù)取值的變量叫做連續(xù)變量，它們的取值可以是任何的實(shí)數(shù)，或者是某一區(qū)間里的任何實(shí)數(shù)，比如股價、時間、身高。但有些性質(zhì)不是連續(xù)的，只有有限個取值的可能性，一般是用于分辨類別，比如性別、婚姻情況、股票所屬行業(yè)，表達(dá)這些變量叫做分類變量。在回歸分析中，我們需要將分類變量轉(zhuǎn)化為啞變量(dummy variable)。 如果我們想表達(dá)一個有 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 點(diǎn)上，便可以用線性關(guān)系表達(dá)，從而進(jì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。常項(xiàng)系數(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ī)的方法進(jìn)行 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 合并，并加上一列常項(xiàng)。 X = np.column_stack((x, dummy))X = sm.add_constant(X) 定義回歸系數(shù)。我們想設(shè)定常項(xiàng)系數(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 模型并進(jìn)行擬合運(yùn)算。 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é)者往往會誤解這一點(diǎn)，認(rèn)為啞變量是一個選擇變量：也就是說，上圖中給出的回歸結(jié)果，是在只做了一次回歸的情況下完成的，而不是分成３段進(jìn)行３次回歸。啞變量的取值藏在其他的三個維度中。可以理解成：上圖其實(shí)是將高元的回歸結(jié)果映射到平面上之后得到的圖。 簡單應(yīng)用 我們來做一個非常簡單的實(shí)際應(yīng)用。設(shè) x 為上證指數(shù)的日收益率，y 為深證成指的日收益率。通過對股票市場的認(rèn)知，我們認(rèn)為 x 和 y 有很強(qiáng)的線性關(guān)系。因此可以假設(shè)模型
并使用 Statsmodels 包進(jìn)行 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'].valuesy_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ù)項(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。最后，畫出數(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 包其他的一些功能，敬請期待。 |
|
|
