本文主要分析時間序列的趨勢和季節性，分解時間序列，實現預測模型。

時間序列預測是基于時間數據進行預測的任務。它包括建立模型來進行觀測，并在諸如天氣、工程、經濟、金融或商業預測等應用中推動未來的決策。

本文主要介紹時間序列預測并描述任何時間序列的兩種主要模式(趨勢和季節性)。并基于這些模式對時間序列進行分解。最后使用一個被稱為Holt-Winters季節方法的預測模型，來預測有趨勢和/或季節成分的時間序列數據。

為了涵蓋所有這些內容，我們將使用一個時間序列數據集，包括1981年至1991年期間墨爾本(澳大利亞)的溫度。

導入庫和數據

首先，導入運行代碼所需的下列庫。除了最典型的庫之外，該代碼還基于statsmomodels庫提供的函數，該庫提供了用于估計許多不同統計模型的類和函數，如統計測試和預測模型。

from?datetime?import?datetime,?timedelta
import?matplotlib.pyplot?as?plt
import?numpy?as?np
import?pandas?as?pd
import?statsmodels.api?as?sm
from?statsmodels.tsa.stattools?import?adfuller,?kpss
from?statsmodels.tsa.api?import?ExponentialSmoothing
%matplotlib?inline

下面是導入數據的代碼。數據由兩列組成，一列是日期，另一列是1981年至1991年間墨爾本(澳大利亞)的溫度。

#?date
numdays?=?365*10?+?2
base?=?'2010-01-01'
base?=?datetime.strptime(base,?'%Y-%m-%d')
date_list?=?[base?+?timedelta(days=x)?for?x?in?range(numdays)]
date_list?=?np.array(date_list)
print(len(date_list),?date_list[0],?date_list[-1])

#?temp
x?=?np.linspace(-np.pi,?np.pi,?365)
temp_year?=?(np.sin(x)?+?1.0)?*?15
x?=?np.linspace(-np.pi,?np.pi,?366)
temp_leap_year?=?(np.sin(x)?+?1.0)
temp_s?=?[]
for?i?in?range(2010,?2020):
????if?i?==?2010:
????????temp_s?=?temp_year?+?np.random.rand(365)?*?20
????elif?i?in?[2012,?2016]:
????????temp_s?=?np.concatenate((temp_s,?temp_leap_year?*?15?+?np.random.rand(366)?*?20?+?i?%?2010))
????else:
????????temp_s?=?np.concatenate((temp_s,?temp_year?+?np.random.rand(365)?*?20?+?i?%?2010))
print(len(temp_s))

#?df
data?=?np.concatenate((date_list.reshape(-1,?1),?temp_s.reshape(-1,?1)),?axis=1)
df_orig?=?pd.DataFrame(data,?columns=['date',?'temp'])
df_orig['date']?=??pd.to_datetime(df_orig['date'],?format='%Y-%m-%d')
df?=?df_orig.set_index('date')
df.sort_index(inplace=True)
df

?

?

可視化數據集

在我們開始分析時間序列的模式之前，讓我們將每個垂直虛線對應于一年開始的數據可視化。

ax?=?df_orig.plot(x='date',?y='temp',?figsize=(12,6))
xcoords?=?['2010-01-01',?'2011-01-01',?'2012-01-01',?'2013-01-01',?'2014-01-01',
???????????'2015-01-01',?'2016-01-01',?'2017-01-01',?'2018-01-01',?'2019-01-01']
for?xc?in?xcoords:
????plt.axvline(x=xc,?color='black',?linestyle='--')
ax.set_ylabel('temperature')

圖1 溫度時間序列

在進入下一節之前，讓我們花點時間看一下數據。這些數據似乎有一個季節性的變化，冬季溫度上升，夏季溫度下降(南半球)。而且氣溫似乎不會隨著時間的推移而增加，因為無論哪一年的平均氣溫都是相同的。

時間序列模式

時間序列預測模型使用數學方程(s)在一系列歷史數據中找到模式。然后使用這些方程將數據[中的歷史時間模式投射到未來。

有四種類型的時間序列模式:

趨勢:數據的長期增減。趨勢可以是任何函數，如線性或指數，并可以隨時間改變方向。

季節性:以固定的頻率(一天中的小時、星期、月、年等)在系列中重復的周期。季節模式存在一個固定的已知周期

周期性:當數據漲跌時發生，但沒有固定的頻率和持續時間，例如由經濟狀況引起。

噪音:系列中的隨機變化。

大多數時間序列數據將包含一個或多個模式，但可能不是全部。這里有一些例子，我們可以識別這些時間序列模式:

維基百科年度受眾(左圖):在這張圖中，我們可以確定一個增加的趨勢，受眾每年線性增加。

美國用電量季節性圖(中圖):每條線對應的是一年，因此我們可以觀察到每年的用電量重復出現的季節性。

ibex35的日收盤價(右圖):這個時間序列隨著時間的推移有一個增加的趨勢，以及一個周期性的模式，因為有一些時期ibex35由于經濟原因而下降。

圖2 從左至右，維基百科的年度受眾，美國用電量的季節性圖，IBEX 35日關閉次數

如果我們假設對這些模式進行加法分解，我們可以這樣寫:

Y[t] = t [t] + S[t] + e[t]

其中Y[t]為數據，t [t]為趨勢周期分量，S[t]為季節分量，e[t]為噪聲，t為時間周期。

另一方面，乘法分解可以寫成:

Y[t] = t [t] *S[t] *e[t]

當季節波動不隨時間序列水平變化時，加法分解是最合適的方法。相反，當季節成分的變化與時間序列水平成正比時，則采用乘法分解更為合適。

分解數據

平穩時間序列被定義為其不依賴于觀察到該序列的時間。因此具有趨勢或季節性的時間序列不是平穩的，而白噪聲序列是平穩的。從數學意義上講，如果一個時間序列的均值和方差不變，且協方差與時間無關，那么這個時間序列就是平穩的。有不同的例子來比較平穩和非平穩時間序列。一般來說，平穩時間序列不會有長期可預測的模式。

為什么穩定性很重要呢?

平穩性已經成為時間序列分析中許多實踐和工具的常見假設。其中包括趨勢估計、預測和因果推斷等。因此，在許多情況下，需要確定數據是否是由固定過程生成的，并將其轉換為具有該過程生成的樣本的屬性。

如何檢驗時間序列的平穩性呢?

我們可以用兩種方法來檢驗。一方面，我們可以通過檢查時間序列的均值和方差來手動檢查。另一方面，我們可以使用測試函數來評估平穩性。

有些情況可能會讓人感到困惑。例如一個沒有趨勢和季節性但具有周期行為的時間序列是平穩的，因為周期的長度不是固定的。

查看趨勢

為了分析時間序列的趨勢，我們首先使用有30天窗口的滾動均值方法分析隨時間推移的平均值。

def?analyze_stationarity(timeseries,?title):
????fig,?ax?=?plt.subplots(2,?1,?figsize=(16,?8))

????rolmean?=?pd.Series(timeseries).rolling(window=30).mean()?
????rolstd?=?pd.Series(timeseries).rolling(window=30).std()
????ax[0].plot(timeseries,?label=?title)
????ax[0].plot(rolmean,?label='rolling?mean');
????ax[0].plot(rolstd,?label='rolling?std?(x10)');
????ax[0].set_title('30-day?window')
????ax[0].legend()
?
????rolmean?=?pd.Series(timeseries).rolling(window=365).mean()?
????rolstd?=?pd.Series(timeseries).rolling(window=365).std()
????ax[1].plot(timeseries,?label=?title)
????ax[1].plot(rolmean,?label='rolling?mean');
????ax[1].plot(rolstd,?label='rolling?std?(x10)');
????ax[1].set_title('365-day?window')
?
pd.options.display.float_format?=?'{:.8f}'.format
analyze_stationarity(df['temp'],?'raw?data')
????ax[1].legend()

圖3 滾動均值和標準差

在上圖中，我們可以看到使用30天窗口時滾動均值是如何隨時間波動的，這是由數據的季節性模式引起的。此外，當使用365天窗口時，滾動平均值隨時間增加，表明隨時間略有增加的趨勢。

這也可以通過一些測試來評估，如Dickey-Fuller (ADF)和Kwiatkowski, Phillips, Schmidt和Shin (KPSS):

ADF檢驗的結果(p值低于0.05)表明，存在的原假設可以在95%的置信水平上被拒絕。因此，如果p值低于0.05，則時間序列是平穩的。

def?ADF_test(timeseries):
????print("Results?of?Dickey-Fuller?Test:")
????dftest?=?adfuller(timeseries,?autolag="AIC")
????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)

ADF_test(df)

Results of Dickey-Fuller Test:
Test Statistic                  -3.69171446
p-value                          0.00423122
Lags Used                       30.00000000
Number of Observations Used   3621.00000000
Critical Value (1%)             -3.43215722
Critical Value (5%)             -2.86233853
Critical Value (10%)            -2.56719507
dtype: float64

KPSS檢驗的結果(p值高于0.05)表明，在95%的置信水平下，不能拒絕的零假設。因此如果p值低于0.05，則時間序列不是平穩的。

def?KPSS_test(timeseries):
????print("Results?of?KPSS?Test:")
????kpsstest?=?kpss(timeseries.dropna(),?regression="c",?nlags="auto")????
????kpss_output?=?pd.Series(
????????kpsstest[0:3],?index=["Test?Statistic",?"p-value",?"Lags?Used"]
????)
????for?key,?value?in?kpsstest[3].items():
????????kpss_output["Critical?Value?(%s)"?%?key]?=?value
????print(kpss_output)

KPSS_test(df)

Results of KPSS Test:
Test Statistic           1.04843270
p-value                  0.01000000
Lags Used               37.00000000
Critical Value (10%)     0.34700000
Critical Value (5%)      0.46300000
Critical Value (2.5%)    0.57400000
Critical Value (1%)      0.73900000
dtype: float64

雖然這些檢驗似乎是為了檢查數據的平穩性，但這些檢驗對于分析時間序列的趨勢而不是季節性是有用的。

統計結果還顯示了時間序列的平穩性的影響。雖然兩個檢驗的零假設是相反的。ADF檢驗表明時間序列是平穩的(p值> 0.05)，而KPSS檢驗表明時間序列不是平穩的(p值> 0.05)。但這個數據集創建時帶有輕微的趨勢，因此結果表明，KPSS測試對于分析這個數據集更準確。

為了減少數據集的趨勢，我們可以使用以下方法消除趨勢:

df_detrend?=?(df?-?df.rolling(window=365).mean())?/?df.rolling(window=365).std()

analyze_stationarity(df_detrend['temp'].dropna(),?'detrended?data')
ADF_test(df_detrend.dropna())

圖4 對時間序列去趨勢后的滾動均值和標準差

檢查季節性

正如在之前從滑動窗口中觀察到的，在我們的時間序列中有一個季節模式。因此應該采用差分方法來去除時間序列中潛在的季節或周期模式。由于樣本數據集具有12個月的季節性，我使用了365個滯后差值:

df_365lag?=??df?-?df.shift(365)

analyze_stationarity(df_365lag['temp'].dropna(),?'12?lag?differenced?data')
ADF_test(df_365lag.dropna())

圖5 對時間序列進行差分后的滾動均值和標準差

現在，滑動均值和標準差隨著時間的推移或多或少保持不變，所以我們有一個平穩的時間序列。

上面方法合在一起的代碼如下：

df_365lag_detrend?=??df_detrend?-?df_detrend.shift(365)

analyze_stationarity(df_365lag_detrend['temp'].dropna(),?'12?lag?differenced?de-trended?data')
ADF_test(df_365lag_detrend.dropna())

圖6 對時間序列進行去趨勢和差分后的滾動均值和標準差

分解模式

基于上述模式的分解可以通過' statmodels '包中的一個有用的Python函數seasonal_decomposition來實現:

def?seasonal_decompose?(df):
????decomposition?=?sm.tsa.seasonal_decompose(df,?model='additive',?freq=365)
?
????trend?=?decomposition.trend
????seasonal?=?decomposition.seasonal
????residual?=?decomposition.resid
?
????fig?=?decomposition.plot()
????fig.set_size_inches(14,?7)
????plt.show()
?
????return?trend,?seasonal,?residual
?
seasonal_decompose(df)

圖7 時間序列分解

在看了分解圖的四個部分后，可以說，在我們的時間序列中有很強的年度季節性成分，以及隨時間推移的增加趨勢模式。

時序建模

時間序列數據的適當模型將取決于數據的特定特征，例如，數據集是否具有總體趨勢或季節性。請務必選擇最適數據的模型。

我們可以使用下面這些模型：

Autoregression (AR)

Moving Average (MA)

Autoregressive Moving Average (ARMA)

Autoregressive Integrated Moving Average (ARIMA)

Seasonal Autoregressive Integrated Moving-Average (SARIMA)

Seasonal Autoregressive Integrated Moving-Average with Exogenous Regressors (SARIMAX)

Vector Autoregression (VAR)

Vector Autoregression Moving-Average (VARMA)

Vector Autoregression Moving-Average with Exogenous Regressors (VARMAX)

Simple Exponential Smoothing (SES)

Holt Winter’s Exponential Smoothing (HWES)

由于我們的數據中存在季節性，因此選擇HWES，因為它適用于具有趨勢和/或季節成分的時間序列數據。

這種方法使用指數平滑來編碼大量的過去的值，并使用它們來預測現在和未來的“典型”值。指數平滑指的是使用指數加權移動平均(EWMA)“平滑”一個時間序列。

在實現它之前，讓我們創建訓練和測試數據集:

y?=?df['temp'].astype(float)
y_to_train?=?y[:'2017-12-31']
y_to_val?=?y['2018-01-01':]
predict_date?=?len(y)?-?len(y[:'2017-12-31'])

下面是使用均方根誤差(RMSE)作為評估模型誤差的度量的實現。

def?holt_win_sea(y,?y_to_train,?y_to_test,?seasonal_period,?predict_date):
?   # 完整數據和代碼，在公眾號「機器學習研習院」后臺回復【墨爾本】獲取
????fit1?=?ExponentialSmoothing(y_to_train,?seasonal_periods=seasonal_period,?trend='add',?seasonal='add').fit(use_boxcox=True)
????fcast1?=?fit1.forecast(predict_date).rename('Additive')
????mse1?=?((fcast1?-?y_to_test.values)?**?2).mean()
????print('The?Root?Mean?Squared?Error?of?additive?trend,?additive?seasonal?of?'+?
??????????'period?season_length={}?and?a?Box-Cox?transformation?{}'.format(seasonal_period,round(np.sqrt(mse1),?2)))

????y.plot(marker='o',?color='black',?legend=True,?figsize=(10,?5))
????fit1.fittedvalues.plot(style='--',?color='red',?label='train')
????fcast1.plot(style='--',?color='green',?label='test')
????plt.ylabel('temp')
????plt.title('Additive?trend?and?seasonal')
????plt.legend()
????plt.show()
?
holt_win_sea(y,?y_to_train,?y_to_val,?365,?predict_date)

The?Root?Mean?Squared?Error?of?additive?trend,?additive?seasonal?of?period?season_length=365?and?a?Box-Cox?transformation?6.27

圖8 Holt-Winters指數平滑法的結果

從圖中我們可以觀察到模型是如何捕捉時間序列的季節性和趨勢的，在異常值的預測上則存在一些誤差。

總結

在本文中，我們通過一個基于溫度數據集的實際示例來介紹趨勢和季節性。除了檢查趨勢和季節性之外，我們還看到了如何降低它，以及如何創建一個基本模型，利用這些模式來推斷未來幾天的溫度。

了解主要的時間序列模式和學習如何實現時間序列預測模型是至關重要的，因為它們有許多應用。

編輯：黃飛

?