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

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

#計算最優(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())

#對于給定范圍內(nèi)的p,q計算擬合得最好的arima模型，這里是對差分好的數(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將年月日的哪個作為索引，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)
    # 對上述函數(shù)求得的值進(jìn)行語義描述
    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ù)分解為趨勢分量，季節(jié)周期和隨機(jī)分量
#對trend進(jìn)行平穩(wěn)定檢驗
testStationarity(trend)

#對序列進(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ù)的自動識別
model = arima_model(trend_diff_2)
model.get_proper_model()
predict_ts = model.properModel.predict()

#還原數(shù)據(jù),因為使用的是乘法模型，將趨勢分量還原之后需要乘以對應(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")

#模型評價
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()