摘要：對(duì)時(shí)間序列模型進(jìn)行優(yōu)化首先將時(shí)序數(shù)據(jù)分解為趨勢(shì)分量，季節(jié)周期分量和隨機(jī)分量對(duì)趨勢(shì)分量使用模型進(jìn)行擬合季節(jié)周期分量則使用歷史同期分量隨機(jī)分量則是使用歷史同類的平均值進(jìn)行預(yù)測(cè)使用面向?qū)ο蟮姆绞剑瑯?gòu)造模型的類，自動(dòng)選取最優(yōu)的模型參數(shù)定義的類計(jì)算最

**對(duì)時(shí)間序列模型進(jìn)行優(yōu)化
1.首先將時(shí)序數(shù)據(jù)分解為趨勢(shì)分量，季節(jié)周期分量和隨機(jī)分量
2.對(duì)趨勢(shì)分量使用ARIMA模型進(jìn)行擬合
3.季節(jié)周期分量則使用歷史同期分量
4.隨機(jī)分量則是使用歷史同類的平均值進(jìn)行預(yù)測(cè)
5.使用面向?qū)ο蟮姆绞剑瑯?gòu)造模型的類，自動(dòng)選取最優(yōu)的模型參數(shù)**

import numpy as np
import pandas as pd
from datetime import datetime
import matplotlib.pylab  as plt
from statsmodels.tsa.stattools import adfuller
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from statsmodels.graphics.tsaplots import plot_acf,plot_pacf 
import sys
from dateutil.relativedelta import relativedelta
from  copy import deepcopy
from statsmodels.tsa.arima_model import ARMA
import warnings
warnings.filterwarnings("ignore")```

class arima_model:

def __init__(self,ts,maxLag = 9):
    self.data_ts = ts
    self.resid_ts = None
    self.predict_ts = None
    self.forecast_ts = None
    self.maxLag = maxLag
    self.p = maxLag
    self.q = maxLag
    self.properModel = None
    self.bic = sys.maxsize

#計(jì)算最優(yōu)的ARIMA模型，將相關(guān)結(jié)果賦給相應(yīng)的屬性
def get_proper_model(self):
    self._proper_model()
    self.predict_ts = deepcopy(self.properModel.predict())
    self.resid_ts = deepcopy(self.properModel.resid)
    self.forecast_ts = deepcopy(self.properModel.forecast())

#對(duì)于給定范圍內(nèi)的p,q計(jì)算擬合得最好的arima模型，這里是對(duì)差分好的數(shù)據(jù)進(jìn)行擬合，故差分恒為0
def _proper_model(self):
    for p in np.arange(self.maxLag):
        for q in np.arange(self.maxLag):
            model = ARMA(self.data_ts, order = (p,q))
            try:
                results_ARMA = model.fit(disp = -1, method = "css")
            except:
                continue
            bic = results_ARMA.bic
            
            if  bic < self.bic:
                self.p = p
                self.q = q
                self.properModel = results_ARMA
                self.bic = bic
                self.resid_ts = deepcopy(self.properModel.resid)
                self.predict_ts = self.properModel.predict()

#參數(shù)確定模型
def certain_model(self,p,q):
    model = ARMA(self.data_ts,order = (p,q))
    try:
        self.properModel = model.fit(disp = -1,method = "css")
        self.p = p
        self.q = q
        self.bic = self.properModel.bic
        self.predict_ts = self.properModel.predict()
        self.resid_ts = deepcopy(self.properModel.resid)
        self.forecast_ts = self.properModel.forecast()
    except:
        print ("You can not fit the model with this parameter p,q")```
        

dateparse = lambda dates:pd.datetime.strptime(dates,"%Y-%m")
#paese_dates指定日期在哪列 index_dates將年月日的哪個(gè)作為索引，date_parser將字符串轉(zhuǎn)為日期
f = open("D:福建AirPassengers.csv")
data = pd.read_csv(f, parse_dates=["Month"],index_col="Month",date_parser=dateparse)
ts = data["#Passengers"]

def draw_ts(timeSeries,title):
    f = plt.figure(facecolor = "white")
    timeSeries.plot(color = "blue")
    plt.title(title)
    plt.show()

def seasonal_decompose(ts):
    from statsmodels.tsa.seasonal import seasonal_decompose
    decomposition = seasonal_decompose(ts, model = "multiplicative")
    trend = decomposition.trend
    seasonal = decomposition.seasonal
    residual = decomposition.resid
    draw_ts(ts,"origin")
    draw_ts(trend,"trend")
    draw_ts(seasonal,"seasonal")
    draw_ts(residual,"residual")
    return trend,seasonal,residual
def testStationarity(ts):
    dftest = adfuller(ts)
    # 對(duì)上述函數(shù)求得的值進(jìn)行語(yǔ)義描述
    dfoutput = pd.Series(dftest[0:4], index=["Test Statistic","p-value","#Lags Used","Number of Observations Used"])
    for key,value in dftest[4].items():
        dfoutput["Critical Value (%s)"%key] = value
#     print ("dfoutput",dfoutput)        
    return dfoutput

ts_log = np.log(ts)
trend,seasonal,residual = seasonal_decompose(ts_log)
seasonal_arr = seasonal
residual = residual.dropna()
residual_mean = np.mean(residual.values)
trend = trend.dropna()

代碼運(yùn)行如下：

#將原始數(shù)據(jù)分解為趨勢(shì)分量，季節(jié)周期和隨機(jī)分量
#對(duì)trend進(jìn)行平穩(wěn)定檢驗(yàn)
testStationarity(trend)

#對(duì)序列進(jìn)行平穩(wěn)定處理
trend_diff_1 = trend.diff(1)
trend_diff_1 = trend_diff_1.dropna()
draw_ts(trend_diff_1,"trend_diff_1")
testStationarity(trend_diff_1)
trend_diff_2 = trend_diff_1.diff(1)
trend_diff_2 = trend_diff_2.dropna()
draw_ts(trend_diff_2,"trend_diff_2")
testStationarity(trend_diff_2)

#使用模型擬合趨勢(shì)分量
#使用模型參數(shù)的自動(dòng)識(shí)別
model = arima_model(trend_diff_2)
model.get_proper_model()
predict_ts = model.properModel.predict()

#還原數(shù)據(jù),因?yàn)槭褂玫氖浅朔Ｐ?，將趨?shì)分量還原之后需要乘以對(duì)應(yīng)的季節(jié)周期分量和隨機(jī)分量
diff_shift_ts = trend_diff_1.shift(1)
diff_recover_1 = predict_ts.add(diff_shift_ts)
rol_shift_ts = trend.shift(1)
diff_recover = diff_recover_1.add(rol_shift_ts)
recover = diff_recover["1950-1":"1960-6"] * seasonal_arr["1950-1":"1960-6"] * residual_mean
log_recover = np.exp(recover)
draw_ts(log_recover,"log_recover")

#模型評(píng)價(jià)
ts_quantum = ts["1950-1":"1960-6"]
plt.figure(facecolor = "white")
log_recover.plot(color = "blue",label = "Predict")
ts_quantum.plot(color = "red", label = "Original")
plt.legend(loc = "best")
plt.title("RMSE %.4f" % np.sqrt(sum((ts_quantum - log_recover) ** 2) / ts_quantum.size))
plt.show()