原文鏈接：https://www.joinquant.com/view/community/detail/ab2a6ecda285c5415de0e1a43db68914
Statsmodels 是 Python 中一個強大的統計分析包，包含了回歸分析、時間序列分析、假設檢
驗等等的功能。Statsmodels 在計量的簡便性上是遠遠不及 Stata 等軟件的，但它的優點在於可以與 Python 的其他的任務（如 NumPy、Pandas）有效結合，提高工作效率。在本文中，我們重點介紹最回歸分析中最常用的 OLS（ordinary least square）功能。
當你需要在 Python 中進行回歸分析時……
import statsmodels.api as sm！！！
在一切開始之前
上帝導入了 NumPy（大家都叫它囊派？我叫它囊辟），
import numpy as np
便有了時間。
上帝導入了 matplotlib，
import matplotlib.pyplot as plt
便有了空間。
上帝導入了 Statsmodels，
import statsmodels.api as sm
世界開始了。
簡單 OLS 回歸
假設我們有回歸模型
Y=β0+β1X1+⋯+βnXn+ε,
並且有 k 組數據 
。OLS 回歸用於計算回歸系數 βi 的估值 b0,b1,…,bn，使誤差平方
最小化。
statsmodels.OLS 的輸入有 (endog, exog, missing, hasconst) 四個，我們現在只考慮前兩個。第一個輸入 endog 是回歸中的反應變量（也稱因變量），是上面模型中的 y(t), 輸入是一個長度為 k 的 array。第二個輸入 exog 則是回歸變量（也稱自變量）的值，即模型中的x1(t),…,xn(t)。但是要注意，statsmodels.OLS 不會假設回歸模型有常數項，所以我們應該假設模型是
並且在數據中，對於所有 t=1,…,k，設置 x0(t)=1。因此，exog的輸入是一個 k×(n+1) 的 array，其中最左一列的數值全為 1。往往輸入的數據中，沒有專門的數值全為1的一列，Statmodels 有直接解決這個問題的函數：sm.add_constant()。它會在一個 array 左側加上一列 1。（本文中所有輸入 array 的情況也可以使用同等的 list、pd.Series 或 pd.DataFrame。）
確切地說，statsmodels.OLS 是 statsmodels.regression.linear_model 里的一個函數（從這個命名也能看出，statsmodel 有很多很多功能，其中的一項叫回歸）。它的輸出結果是一個 statsmodels.regression.linear_model.OLS，只是一個類，並沒有進行任何運算。在 OLS 的模型之上調用擬合函數 fit()，才進行回歸運算，並且得到 statsmodels.regression.linear_model.RegressionResultsWrapper，它包含了這組數據進行回歸擬合的結果摘要。調用 params 可以查看計算出的回歸系數 b0,b1,…,bn。
簡單的線性回歸
上面的介紹繞了一個大圈圈，現在我們來看一個例子，假設回歸公式是:
我們從最簡單的一元模型開始，虛構一組數據。首先設定數據量，也就是上面的 k 值。
nsample = 100
然后創建一個 array，是上面的 x1 的數據。這里，我們想要 x1 的值從 0 到 10 等差排列。
x = np.linspace(0, 10, nsample)
使用 sm.add_constant() 在 array 上加入一列常項1。
X = sm.add_constant(x)
然后設置模型里的 β0,β1，這里要設置成 1,10。
beta = np.array([1, 10])
然后還要在數據中加上誤差項，所以生成一個長度為k的正態分布樣本。
e = np.random.normal(size=nsample)
由此，我們生成反應項 y(t)。
y = np.dot(X, beta) + e
好嘞，在反應變量和回歸變量上使用 OLS() 函數。
model = sm.OLS(y,X)
然后獲取擬合結果。
results = model.fit()
再調取計算出的回歸系數。
print(results.params)
得到
[ 1.04510666, 9.97239799]
和實際的回歸系數非常接近。
當然，也可以將回歸擬合的摘要全部打印出來。
print(results.summary())
得到
Adj. R-squared：俗稱 R平方。反映了模型的擬合。R平方值的范圍從 0 到 1，其中較高的值通常表示更好的擬合，假設滿足某些條件； const coef：這個系數是你的 Y 截距。這意味着如果 Interest_Rate 和 Unemployment_Rate 系數都為零，則預期輸出（即 Y）將等於 const 系數； Interest_Rate：表示由於利率中一個單位的而變化導致的輸出 Y 的變化（其他一切保持不變）； Unemployment_Rate：表示由於失業率中一個單位的變化導致的產出 Y 的變化（其他一切保持不變）； std err：反映了系數的准確度，它越低，准確度越高； P>|t|：這是 p 值，p 值小於 0.05 被認為是統計學上顯著的； Confidence Interval：這是置信區間，表示我們的系數可能下降的范圍（可能性為 95%）；
中間偏下的 coef 列就是計算出的回歸系數。
我們還可以將擬合結果畫出來。先調用擬合結果的 fittedvalues 得到擬合的 y 值。
y_fitted = results.fittedvalues
然后使用 matplotlib.pyploft 畫圖。首先設定圖軸，圖片大小為 8×6。
fig, ax = plt.subplots(figsize=(8,6))
畫出原數據，圖像為圓點，默認顏色為藍。
ax.plot(x, y, 'o', label='data')
畫出擬合數據，圖像為紅色帶點間斷線。
ax.plot(x, y_fitted, 'r--.',label='OLS')
放置注解。
ax.legend(loc='best')
得到
在大圖中看不清細節，我們在 0 到 2 的區間放大一下，可以見數據和擬合的關系。
加入改變坐標軸區間的指令
ax.axis((-0.05, 2, -1, 25))
高次模型的回歸
假設反應變量 Y 和回歸變量 X 的關系是高次的多項式，即
我們依然可以使用 OLS 進行線性回歸。但前提條件是，我們必須知道 X 在這個關系中的所有次方數；比如，如果這個公式里有一個 
.5項，但我們對此並不知道，那么用線性回歸的方法就不能得到准確的擬合。
雖然 X 和 Y 的關系不是線性的，但是 Y 和 
的關系是高元線性的。也就是說，只要我們把高次項當做其他的自變量，即設 
。那么，對於線性公式
可以進行常規的 OLS 回歸，估測出每一個回歸系數 βi。可以理解為把一元非線性的問題映射到高元，從而變成一個線性關系。
下面以
為例做一次演示。
首先設定數據量，也就是上面的 k 值。
nsample = 100
然后創建一個 array，是上面的 x1 的數據。這里，我們想要 x1 的值從 0 到 10 等差排列。
x = np.linspace(0, 10, nsample)
再創建一個 k×2 的 array，兩列分別為 x1 和 x2。我們需要 x2 為 x1 的平方。
X = np.column_stack((x, x**2))
使用 sm.add_constant() 在 array 上加入一列常項 1。
X = sm.add_constant(X)
然后設置模型里的 β0,β1,β2，我們想設置成 1,0.1,10。
beta = np.array([1, 0.1, 10])
然后還要在數據中加上誤差項，所以生成一個長度為k的正態分布樣本。
e = np.random.normal(size=nsample)
由此，我們生成反應項 y(t)，它與 x1(t) 是二次多項式關系。
y = np.dot(X, beta) + e
在反應變量和回歸變量上使用 OLS() 函數。
model = sm.OLS(y,X)
然后獲取擬合結果。
results = model.fit()
再調取計算出的回歸系數。
print(results.params)
得到
[ 0.95119465, 0.10235581, 9.9998477]
獲取全部摘要
print(results.summary())
得到
擬合結果圖如下
在橫軸的 [0,2] 區間放大，可以看到
啞變量
一般而言，有連續取值的變量叫做連續變量，它們的取值可以是任何的實數，或者是某一區間里的任何實數，比如股價、時間、身高。但有些性質不是連續的，只有有限個取值的可能性，一般是用於分辨類別，比如性別、婚姻情況、股票所屬行業，表達這些變量叫做分類變量。在回歸分析中，我們需要將分類變量轉化為啞變量(dummy variable)。
如果我們想表達一個有 d 種取值的分類變量，那么它所對應的啞變量的取值是一個 d 元組（可以看成一個長度為 d 的向量），其中有一個元素為 1，其他都是 0。元素呈現出 1 的位置就是變量所取的類別。比如說，某個分類變量的取值是 {a,b,c,d}，那么類別 a 對應的啞變量是(1,0,0,0)，b 對應 (0,1,0,0)，c 對應 (0,0,1,0)，d 對應 (0,0,0,1)。這么做的用處是，假如 a、b、c、d 四種情況分別對應四個系數 β0,β1,β2,β3，設 (x0,x1,x2,x3) 是一個取值所對應的啞變量，那么
可以直接得出相應的系數。可以理解為，分類變量的取值本身只是分類，無法構成任何線性關系，但是若映射到高元的 0,1 點上，便可以用線性關系表達，從而進行回歸。
Statsmodels 里有一個函數 categorical() 可以直接把類別 {0,1,…,d-1} 轉換成所對應的元組。確切地說，sm.categorical() 的輸入有 (data, col, dictnames, drop) 四個。其中，data 是一個 k×1 或 k×2 的 array，其中記錄每一個樣本的分類變量取值。drop 是一個 Bool值，意義為是否在輸出中丟掉樣本變量的值。中間兩個輸入可以不用在意。這個函數的輸出是一個k×d 的 array（如果 drop=False，則是k×(d+1)），其中每一行是所對應的樣本的啞變量；這里 d 是 data 中分類變量的類別總數。
我們來舉一個例子。這里假設一個反應變量 Y 對應連續自變量 X 和一個分類變量 Z。常項系數為 10，XX 的系數為 1；Z 有 {a,b,c}三個種類，其中 a 類有系數 1，b 類有系數 3，c 類有系數 8。也就是說，將 Z 轉換為啞變量 (Z1,Z2,Z3)，其中 Zi 取值於 0,1，有線性公式
可以用常規的方法進行 OLS 回歸。
我們按照這個關系生成一組數據來做一次演示。先定義樣本數量為 50。
nsample = 50
設定分類變量的 array。前 20 個樣本分類為 a。
groups = np.zeros(nsample, int)
之后的 20 個樣本分類為 b。
groups[20:40] = 1
最后 10 個是 c 類。
groups[40:] = 2
轉變成啞變量。
dummy = sm.categorical(groups, drop=True)
創建一組連續變量，是 50 個從 0 到 20 遞增的值。
x = np.linspace(0, 20, nsample)
將連續變量和啞變量的 array 合並，並加上一列常項。
-
X = np.column_stack((x, dummy))
-
X = sm.add_constant(X)
定義回歸系數。我們想設定常項系數為 10，唯一的連續變量的系數為 1，並且分類變量的三種分類 a、b、c 的系數分別為 1,3,8。
beta = [10, 1, 1, 3, 8]
再生成一個正態分布的噪音樣本。
e = np.random.normal(size=nsample)
最后，生成反映變量。
y = np.dot(X, beta) + e
得到了虛構數據后，放入 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 同時影響了回歸的結果。初學者往往會誤解這一點，認為啞變量是一個選擇變量：也就是說，上圖中給出的回歸結果，是在只做了一次回歸的情況下完成的，而不是分成３段進行３次回歸。啞變量的取值藏在其他的三個維度中。可以理解成：上圖其實是將高元的回歸結果映射到平面上之后得到的圖。
簡單應用
我們來做一個非常簡單的實際應用。設 x 為上證指數的日收益率，y 為深證成指的日收益率。通過對股票市場的認知，我們認為 x 和 y 有很強的線性關系。因此可以假設模型
並使用 Statsmodels 包進行 OLS 回歸分析。
我們取上證指數和深證成指一年中的收盤價。
-
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
計算兩個指數一年內的日收益率，記載於 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)
將數據轉化為 array 的形式；不要忘記添加常數項。
-
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。最后，畫出數據和擬合線。
-
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')
結語
本篇文章中，我們介紹了 Statsmodels 中很常用 OLS 回歸功能，並展示了一些使用方法。線性回歸的應用場景非常廣泛。在我們量化課堂應用類的內容中，也有相當多的策略內容采用線性回歸的內容。我們會將應用類文章中涉及線性回歸的部分加上鏈接，鏈接到本篇文章中來，形成體系。量化課堂在未來還會介紹 Statsmodel 包其他的一些功能，敬請期待。
到JoinQuant查看代碼並與作者交流討論：【量化課堂】Statsmodels 統計包之 OLS 回歸
更多參考：
https://blog.csdn.net/CoderPai/article/details/83899268
https://www.jianshu.com/p/e45558ccf533