使用ARIMA模型拟合股票資料
ARIMA可拆分為AR , I , MA。
AR模型：
MA模型：
ARMA模型：
ARIMA模型：
ARIMA：差分自回歸移動平均模型
- AR是自回歸，P是自回歸項；
- d 為時間序列稱為平穩時 所做的差分次數
- MA 是移動平均， q 為移動平均項
- 原理： 将非平穩時間序列轉換成平穩時間序列， 然後将因變量僅對它的滞後值（p階）以及随機誤差項的現值和滞後值進行回顧所建立的模型。
ARIMA模組化流程
- 序列平穩化（ 差分确定）
- p和q的階數确定（ACF,PACF）
- 建立模型并預測
- 資訊準則：BIC
模型選擇準則AIC&&BIC
AIC： Akaike information criterion，赤池資訊量。
BIC：Bayesian information criterion，貝葉斯資訊度量
要解決模型選擇的準則問題需要先搞懂似然函數。
百度百科說，似然函數是一種關于統計模型中的參數的函數，表示模型參數中的似然性。
不明覺厲，還是看wikipedia吧
常說的機率是指給定參數後，預測即将發生的事件的可能性。拿硬币這個例子來說，我們已知一枚均勻硬币的正反面機率分别是0.5，要預測抛兩次硬币，硬币都朝上的機率：p(HH;p=0.5)
而似然機率正好與這個過程相反，我們關注的量不再是事件的發生機率，而是已知發生了某些事件，我們希望知道參數應該是多少。
現在我們已經抛了兩次硬币，并且知道了結果是兩次頭朝上，這時候，我希望知道這枚硬币抛出去正面朝上的機率為0.5的機率是多少？正面朝上的機率為0.8的機率是多少？
如果我們希望知道正面朝上機率為0.5的機率，這個東西就叫做似然函數，可以說成是對某一個參數的猜想（p=0.5）的機率，這樣表示成(條件)機率就是
L(pH=0.5|HH) = P(HH|pH=0.5)
為什麼可以寫成這樣？我覺得可以這樣來想：
似然函數本身也是一種機率，我們可以把L(pH=0.5|HH)寫成P(pH=0.5|HH); 而根據貝葉斯公式，P(pH=0.5|HH) = P(pH=0.5,HH)/P(HH)；既然HH是已經發生的事件，理所當然P(HH) = 1,是以：
P(pH=0.5|HH) = P(pH=0.5,HH) = P(HH;pH=0.5).
右邊的這個計算我們很熟悉了，就是已知頭朝上機率為0.5，求抛兩次都是H的機率，即0.5*0.5=0.25。
是以，我們可以safely得到:
L(pH=0.5|HH) = P(HH|pH=0.5) = 0.25.
這個0.25的意思是，在已知抛出兩個正面的情況下，pH = 0.5的機率等于0.25。
再算一下
L(pH=0.6|HH) = P(HH|pH=0.6) = 0.36.
把pH從0~1的取值所得到的似然函數的曲線畫出來得到這樣一張圖：
（來自wikipedia）
可以發現，pH = 1的機率是最大的。
即L(pH = 1|HH) = 1。
那麼最大似然機率的問題也就好了解了。
最大似然機率，就是在已知觀測的資料的前提下，找到使得似然機率最大的參數值。
這就不難了解，在data mining領域，許多求參數的方法最終都歸結為最大化似然機率的問題。
回到這個硬币的例子上來，在觀測到HH的情況下，pH = 1是最合理的（卻未必符合真實情況，因為資料量太少的緣故）。
說得有點多了哈哈，回到正題。
很多參數估計問題均采用似然函數作為目标函數，當訓練資料足夠多時，可以不斷提高模型精度，但是以提高模型複雜度為代價的，同時帶來一個機器學習中非常普遍的問題——過拟合。是以，模型選擇問題在模型複雜度與模型對資料集描述能力（即似然函數）之間尋求最佳平衡。
人們提出許多資訊準則，通過加入模型複雜度的懲罰項來避免過拟合問題，而AIC和BIC就是來解決這個問題的準則。
AIC是衡量統計模型拟合優良性的一種标準，由日本統計學家赤池弘次在1974年提出，它建立在熵的概念上，提供了權衡估計模型複雜度和拟合資料優良性的标準。
AIC = 2k - 2ln(L)
其中k是模型參數個數，L是似然函數。從一組可供選擇的模型中選擇最佳模型時，通常選擇AIC最小的模型。
當兩個模型之間存在較大差異時，差異主要展現在似然函數項，當似然函數差異不顯著時，上式第一項，即模型複雜度則起作用，進而參數個數少的模型是較好的選擇。當模型複雜度提高（k增大）時，似然函數L也會增大，進而使AIC變小，但是k過大時，似然函數增速減緩，導緻AIC增大，模型過于複雜容易造成過拟合現象。目标是選取AIC最小的模型，AIC不僅要提高模型拟合度（極大似然），而且引入了懲罰項，使模型參數盡可能少，有助于降低過拟合的可能性。
BIC貝葉斯資訊準則與AIC相似，用于模型選擇,
BIC = kln(n) - 2ln(L)
其中，k為模型參數個數，n為樣本數量，L為似然函數。kln(n)懲罰項在維數過大且訓練樣本資料相對較少的情況下，可以有效避免出現次元災難現象。
訓練模型時，增加參數數量，也就是增加模型複雜度，會增大似然函數，但是也會導緻過拟合現象，針對該問題，AIC和BIC均引入了與模型參數個數相關的懲罰項，BIC的懲罰項比AIC的大，考慮了樣本數量，樣本數量過多時，可有效防止模型精度過高造成的模型複雜度過高。
好吧，其實AIC和BIC也有點偏題了，回正題哈哈
看一下怎麼用ARIMA咯
首先我們使用ACF和PACF确定ARMA的階數，進而确定ARIMA的模型參數。第二種方式使用bic選擇參數
上代碼，首先導入package
import pandas as pd
import pandas_datareader
import datetime
import matplotlib.pylab as plt
from matplotlib.pylab import style
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
擷取茅台的股票資料來分析
start_date = datetime.datetime(2007, 1, 1)
# 指定股票分析截止日期
end_date = datetime.datetime(2019, 3, 1)
# 股票代碼
stock_code = '600519.SS' # 滬市貴州茅台
stock_df = pandas_datareader.data.DataReader(stock_code, 'yahoo', start_date, end_date)
# 預覽資料
print(stock_df.head())
隻分析收盤價
stock_all = stock_df.Close
plt.plot(stock_all)
plt.title('股票每日收盤價')
plt.show()
資料太多，太遠的資料其實沒啥用
stock_resample = stock_all.resample('W-MON').mean()
stock_train = stock_resample['2014':'2018']
plt.plot(stock_train)
plt.title('周收盤價')
plt.show()
自相關圖
#畫出自相關性圖
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
plot_acf(stock_train,lags=50)
plt.show()
其實從資料圖和自相關圖就可以看到，資料是非平穩的，檢測一下看看
#平穩性檢測
from statsmodels.tsa.stattools import adfuller
print('原始序列的檢驗結果為：',adfuller(stock_train))
#傳回值依次為：adf, pvalue p值，usedlag, nobs, critical values臨界值 , icbest, regresults, resstore
answer：
原始序列的檢驗結果為： (-0.9039522834177104, 0.7866680876785267, 15, 245, {'1%': -3.4573260719088132, '5%': -2.873410402808354, '10%': -2.573095980841316}, 1984.5615330766404)
看到p-value遠遠大于0.05，做一下差分吧
stock_train_diff = stock_train.diff().dropna()
stock_train_diff.plot()
plt.show()
好像是平穩了一些窩。畫出ACF和PACF圖定階吧
plot_acf(stock_train_diff,lags=20) #畫出自相關圖
plt.show()
plot_pacf(stock_train_diff,lags=20) #畫出偏相關圖
plt.show()
通過PACF确定p=4，ACF确定q=3
檢驗結果，p-value達标了，不過差強人意。
print(u'差分序列的ADF 檢驗結果為： ', adfuller(stock_train_diff)) #平穩性檢驗
#一階差分後的序列的時序圖在均值附近比較平穩的波動， 自相關性有很強的短期相關性， 機關根檢驗 p值小于 0.05
是以說一階差分後的序列是平穩序列
#對一階差分後的序列做白噪聲檢驗
from statsmodels.stats.diagnostic import acorr_ljungbox
print(u'差分序列的白噪聲檢驗結果：',acorr_ljungbox(stock_train_diff, lags= 1)) #傳回統計量和 p 值
# 差分序列的白噪聲檢驗結果：(array([13.84221767]), array([0.00019882])) p值為第二項，
最後拟合模型并預測，把預測結果和訓練資料結合concat
model = ARIMA(stock_train,(4,1,3)).fit()
pre = model.forecast(len(stock_resample['2019':]))[0]
pre_result = pd.Series(pre,index = stock_resample['2019':].index)
result_1 = pd.concat([stock_train,pre_result],axis = 0)
第二種方法是使用bic對模型進行定階
#對模型進行定階
from statsmodels.tsa.arima_model import ARIMA
pmax = 5 #一般階數不超過 length /10
qmax = 5
bic_matrix = []
for p in range(pmax +1):
temp= []
for q in range(qmax+1):
try:
temp.append(ARIMA(stock_train_diff, (p, 1, q)).fit().bic)
except:
temp.append(None)
bic_matrix.append(temp)
bic_matrix = pd.DataFrame(bic_matrix) #将其轉換成Dataframe 資料結構
p,q = bic_matrix.stack().idxmin() #先使用stack 展平， 然後使用 idxmin 找出最小值的位置
print(u'BIC 最小的p值 和 q 值：%s,%s' %(p,q)) # BIC 最小的p值 和 q 值：0,2
同樣地，拟合模型并預測
model2 = ARIMA(stock_train,(p,1,q)).fit()
pre2 = model2.forecast(9)[0]
pre_result2 = pd.Series(pre2,index = stock_resample['2019':].index)
result_2 = pd.concat([stock_train,pre_result],axis = 0)
最後，畫圖，看看兩種方法選擇的階數拟合的模型預測效果如何
plt.plot(result_1,c = 'r')
plt.plot(stock_resample['2019':],c = 'b')
plt.show()
plt.plot(result_2,c = 'r')
plt.plot(stock_resample['2019':],c = 'b')
plt.show()
說實話，好像基本一樣，可能是本身模型預測效果就不好，這裡也沒做很多的資料處理，具體問題，再想想吧哈哈