---
# System prepended metadata

title: 時間序列分析與預測

---

# 時間序列分析與預測

Author : Xin

## 前言

關於時間序列分析與預測，強烈推薦有興趣的人閱讀下列書籍：
Hyndman, R.J., & Athanasopoulos, G. (2021) Forecasting: principles and practice, 3rd edition, OTexts: Melbourne, Australia. [(英文版連結)](https://otexts.com/fpp3/index.html) [(簡中版連結)](https://otexts.com/fpp3cn/)
這本書不僅是免費公開的，而且可以在網路上直接閱讀。書中有大量現實數據例子，同時每個章節皆有教學影片和程式碼。

## 時間序列分析

對下列章節，我們使用Air Passengers這個資料集作為範例，可以在[(資料連結)](https://www.kaggle.com/datasets/rakannimer/air-passengers/data)中下載到檔案。
首先載入需要用到的Python套件和設定：
```python=
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import warnings
from statsmodels.stats.diagnostic import acorr_ljungbox
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.seasonal import STL
from sklearn.metrics import mean_squared_error
from pmdarima.arima import auto_arima

warnings.filterwarnings("ignore")
```

### 時間序列組成

通常會使用下列術語來描述時間序列：
- ***趨勢(Trend)***
    這個很容易理解，基本上就是數據變化的方向。例如由於暖化因素，地球溫度長期上升。
    趨勢並不一定是線性的，在某些時候也可能由增長趨勢轉變為遞減趨勢。
- ***季節性(Seasonal)***
    當時間序列受到某個==固定==頻率的因素影響時，則此序列具季節性。例如一年四季的溫度變化就是季節性的一個例子。
- ***週期性(Cyclic)***
    當時間序列受到某個==非固定==頻率的因素影響時，則此序列具週期性。例如景氣循環就是週期性的一個例子。

總結來說，當時間序列波動是無規律且隨機的，則其具有週期性。如果波動頻率具有規律，則其具有季節性。對於如何分辨季節性和週期性之間的不同，也有以下簡易方法：
- 對於季節性，我們可以預測出何時數據會達到高點和低點。比如說每日平均溫度在夏天時達到高點、冬天時達到低點。
- 對於週期性，則無法判斷。舉例來說，對景氣循環週期，無從預知何時景氣好、何時景氣差。

接下來先來看看數據長甚麼樣子，首先載入資料並調整設定。
```python=+
data = pd.read_csv(**填入檔案位置**)
data['Date'] = pd.to_datetime(data['Month'])
data = data.drop(columns = 'Month')
data = data.set_index('Date')
data = data.rename(columns = {'#Passengers':'Passengers'})
```
資料為每月的乘客數量統計，接下來將數據圖像化
```python=+
data.plot(title = "Numbers of Passengers")
```

|![output_3_1](https://hackmd.io/_uploads/SkATPKsR1e.png)|
|:--:|
|從圖中可看出，此時間序列具有增加的趨勢和每12個月循環的季節性|

### 平穩時間序列

由於時間序列的特性，自然的會想使用過去的歷史數據來預測未來的數據。因此當我們在預測時，會希望序列的統計性質不隨著時間改變，這樣就代表未來的數據會表現出和過去相同的規律。
對於這樣的時間序列，我們把它叫做平穩(Stationary)時間序列。通常來說，會使用下列的方式來定義平穩時間序列：
- 時間序列$~y_t~$是平穩的條件為：
    1. 期望值$~\mathbb{E}[y_t]~$不隨著時間改變。
    2. 自共變異數(Autocovariance)不隨著時間改變。
    對於$~y_t~$和$~y_{t-k}~$，通常會將時間間隔$~k~$稱作為 lag。lag 為$~k~$的自共變異數定義為
    $$\operatorname{Cov}(y_t , y_{t-k})$$ 上述式子理論上會跟計算的時間點$~t~$有關，平穩的條件要求其只與$~k~$有關而不隨著時間$~t~$改變。若將$~k=0~$帶入，則得出變異數$$\sigma^2 = \sigma_t^2 = \operatorname{Cov}(y_t , y_{t})$$ 不隨時間改變。

時間序列的自相關函數(Autocorrelation Function 或簡稱為ACF)的定義為
$$\operatorname{ACF}(k) = \frac{\operatorname{Cov}(y_t , y_{t-k})}{\operatorname{Var}[y_t]} $$ 由上述平穩定義可知，平穩時間序列的ACF只與$~k~$有關不隨時間而改變。$~k = 0~$時ACF永遠為1。

對於時間序列，有以下事實：
1. 有趨勢的時間序列非平穩：
    一個簡單的判別方法便是檢查時間序列的移動平均(Moving Average)。若其隨著時間而變動，則與平穩定義第一條相違背。
    ```python=+
    plt.figure()
    line1, = plt.plot(data.rolling(12).mean(), label = 'Moving Average')             
    line2, = plt.plot(data.rolling(12).std(), label = 'Moving Standard Deviation')
    plt.legend(handles = [line1, line2], loc='upper left')
    plt.show()
    ```
    |![output_4_0](https://hackmd.io/_uploads/rJ50DKoAkx.png)|
    |:--:|
    |從圖中可看出，移動平均隨著時間而改變|
1. 季節性的時間序列非平穩：
    若序列是平穩的，則其自相關函數會隨著lag的增加而迅速遞減至零。但是季節性序列的自相關函數會在lag和週期倍數相等時變大，因此與平穩定義第二條相違背。之後會用例子展示這點。
1. 若要確定數據是否平穩，可以同時使用ADF test和KPSS test來檢定，這些檢定可以用來檢驗數據是否需要差分。[(參考範例)](https://www.statsmodels.org/dev/examples/notebooks/generated/stationarity_detrending_adf_kpss.html)

### ACF和PACF

我們也可以從序列的ACF圖中看出一些端倪。
- 當數據具有趨勢時，對於較小的lag，ACF值通常為正且較大，之後隨著lag增加而緩慢遞減。
- 當數據具有季節性時，ACF值在lag為季節性的週期的倍數時會變大。
- 當數據兼具趨勢和季節性時，則兩種特徵會混合出現。

可以使用下列程式碼畫出ACF圖
```python=+
plot_acf(data , lags = 25); print()
```
|![output_7_1](https://hackmd.io/_uploads/HkdJuFoAJe.png)|
|:--:|
|從圖中可看出，此序列具有趨勢和季節性的混合特徵。ACF值在lag為12的倍數時變大|

- ACF反映出了不同$~k~$之下，$~y_t~$和$~y_{t-k}~$之間的關係。假設我們知道$~y_t~$和$~y_{t-1}~$相關，則$~y_t~$和$~y_{t-2}~$也相關。但這可能是因為$~y_{t}~,~y_{t-2}~$都與$~y_{t-1}~$相關，而不是因為$~y_{t-2}~$本身帶有資訊可以預測$~y_{t}~$。
- 為了解決此問題，提出了偏自相關(Partial Autocorrelations)的概念。它的目的是在移除 lag 1, 2, ..., $k-1~$的影響下，衡量$~y_t~$和$~y_{t-k}~$的相關性。PACF圖在之後決定ARIMA模型階數時至關重要。我們在之後會看到這點。

### 白噪音(White Noise)

基本上白噪音代表了不可控的隨機因素。白噪音$~\varepsilon_t~$的性質如下：
1. 對於任意$~t~,\mathbb{E}[\varepsilon_t] = 0$ 
1. 對於任意$~t~,\operatorname{Var}[\varepsilon_t] = \sigma^2 < \infty$
1. 對於任意$~t~$和$~k~,\operatorname{Cov}(\varepsilon_t , \varepsilon_{t-k}) = 0$

根據上述性質，可知白噪音是平穩序列。理想中，我們希望將數據分解為可控制的部分以及白噪音：
$$y_t = \text{Our Model} ~+~ \varepsilon_t $$ 當要預測$~y_t~$時，可計算$~y_t~$的期望值：
$$\mathbb{E}[y_t] = \mathbb{E}[\text{Our Model}]$$ 若是白噪音服從常態分布，則在95%信心水準下，$~y_t~$會落在$~\mathbb{E}[y_t]\pm 2\sigma~$的範圍內。
當建立好模型後，需要確認殘差值 $$\varepsilon_t = y_t - \text{Our Model}$$ 是不是白噪音。此時可以利用Ljung-Box檢定來測試。如果殘差不是白噪音，那麼就必須重新選定一個新的模型。

### 去除趨勢及季節性

如果希望時間序列變平穩，就必須想辦法去除趨勢及季節性。
- 去除趨勢：有一個簡單的方法能夠去除趨勢，那就是對資料做差分。假設數據為$$y_1,y_2,\cdots,y_n$$ 那麼差分後的數據變為
$$y_2 - y_1 ~,~ y_3 - y_2 ~,~ y_4 - y_3 ~,~ \cdots ~,~ y_n - y_{n-1}$$ 也就是說，我們將$~y_t~$的值替換成了$~y_t - y_{t-1}~$。注意到由於第一個資料$~y_1~$並無前一項，差分後的數據會少掉一筆資料，因此剩下$~n-1~$筆數據。
考慮一個簡單的例子，假設數據符合方程式$~y_t = t~$，則此時間序列顯然有線性增長的趨勢。因為每個$~y_t~$在差分後都變為
$$y_t - y_{t-1} = t - (t - 1) = 1$$ 所以差分後數據的趨勢消失了。
    可以使用以下程式碼來對數據進行差分：
    ```python=+
    diff_data = data.diff()
    diff_data.plot()
    ```
    |![output_5_1](https://hackmd.io/_uploads/BkwgdtiCJg.png)|
    |:--:|
    |從圖中可看出，經過差分後數據的趨勢消失了|

- 去除季節性：與去除趨勢時相似，也可以利用差分將季節性去掉。以我們的資料為例，季節性的週期為12個月。那麼對$~y_t~$就做季節性差分$~y_t - y_{t-12}~$，也就是說，將$~y_t~$和去年同一月分的數據做差分。這樣得出的結果都會是同一季節之間互相比較。如同之前一樣，由於第一個季節並無前一季，因此差分後的數據會少掉第一季的資料，因此剩下$~n-12~$筆數據。
可以使用以下程式碼來消除季節性：
    ```python=+
    lag = 12
    lagged_data = data.shift(lag)
    deseason_by_lag = data - lagged_data
    deseason_by_lag.plot()
    ```
    |![output_6_1](https://hackmd.io/_uploads/SkZ-uKoA1g.png)|
    |:--:|
    |從圖中可看出，經過季節性差分後數據的季節性消失了|

### 時間序列分解

對於時間序列分解，通常會將其分解為$~S_t~ , ~T_t~ , ~R_t~$的函數。其中$~S_t~$代表季節項、$~T_t~$代表趨勢和週期項、$~R_t~$代表殘差項。對於具體的分解方法，有以下兩種模型：
1. 加法模型：
    加法模型假設序列是由多種成分相加而得，具體來說，
    $$
    y_t = S_t + T_t + R_t
    $$ 當季節性波動的幅度或趨勢週期項波動的幅度不隨時間序列的值變化而變化時，則比較適合使用加法模型。
1. 乘法模型：
    乘法模型假設序列是由多種成分相乘而得，具體來說，
    $$
    y_t = S_t \times T_t \times R_t
    $$ 當季節性波動的幅度或趨勢週期項波動的幅度與時間序列的值成比例時，則比較適合使用乘法模型。

以我們的資料來說，由於季節性波動的峰值高度和整體寬度隨著時間序列的值增加而增加，因此比較適合使用乘法模型。
|![output_3_1](https://hackmd.io/_uploads/SkATPKsR1e.png)|
|:--:|

對於乘法模型，可以將數據經過對數變換轉為加法模型：
$$y_t = S_t \times T_t \times R_t\Longleftrightarrow\log y_t = \log S_t + \log T_t + \log R_t$$ 可以使用套件的STL分解法幫我們自動分解完成，下面我們使用了對數轉換將模型換為加法模型。
```python=+
decomposition = STL(np.log(data["Passengers"]), period = 12).fit()
figure = decomposition.plot()
```
|![output_8_0](https://hackmd.io/_uploads/rkhfuKoAyg.png)|
|:--:|
|AirPassengers數據的STL分解|

我們也可以使用STL分解的結果來消除趨勢和季節性，例如將原本的數據減掉STL中的Season項就得到了去除季節性的資料。


## 時間序列模型

對於預測方法，我們介紹以下幾個經典的模型。

### 自回歸模型 (AutoRegressive models)

顧名思義，就是對變量自身做回歸。一個$~p~$階的自回歸模型，稱為**AR($p$)**，可以表示如下
$$y_t = c + \phi_1y_{t-1} + \phi_2y_{t-2} + \cdots +\phi_py_{t-p} + \varepsilon_t$$ 其中$~\varepsilon_t~$是白噪音。這就相當於用前$~p~$個  lag 項對$~y_t~$做回歸模型。

一般來說，AR($p$)模型有以下特徵：
- ACF圖中呈現指數型遞減，或類似於正弦的波動。
- PACF圖中在lag$~p~$時有顯著突起，但lag更大時就不再有。

使用下列程式碼來模擬AR(1)模型
$$y_t = 0.8y_{t-1} + \varepsilon_t $$
```python=+
from statsmodels.tsa.arima_process import ArmaProcess

AR_1 = ArmaProcess([1 , -0.8] , [1])
simulated_AR_1 = AR_1.generate_sample(nsample=1000)
plt.plot(simulated_AR_1)
plt.show()
```

|![output_11_0](https://hackmd.io/_uploads/ryhm_KsRkg.png)|
|:--:|
|AR(1)模型的模擬|

```python=+
plot_acf(simulated_AR_1, lags = 10)
plt.show()
plot_pacf(simulated_AR_1, lags = 10)
plt.show()
```

|![output_12_0](https://hackmd.io/_uploads/HJrEdYjAkg.png)|![output_12_1](https://hackmd.io/_uploads/rk04OKoRJx.png)|
|:--:|:--:|
|ACF圖|PACF圖|

從圖中可觀察出
- ACF圖中呈現出指數型遞減。
- PACF圖中在lag 1 時有顯著突起，但lag更大時就不再有。

### 移動平均模型 (Moving Average models)

不同於使用變量的歷史值回歸，移動平均模型有點像是使用歷史預測誤差來建立回歸模型。一個$~q~$階的移動平均模型，稱為**MA($q$)**，可以表示如下
$$y_t = c + \varepsilon_t + \theta_1\varepsilon_{t-1} + \theta_2\varepsilon_{t-2} + \cdots +\theta_p\varepsilon_{t-p}$$ 其中$~\varepsilon_t~$是白噪音。

- 任何AR($p$)模型其實都可以用MA($\infty$)表示：
    例如，AR(1)模型
    $$y_t = \phi_1y_{t-1} + \varepsilon_t$$ 假定$~|\phi_1| < 1~$，則其可以表示為
    \begin{align}
    y_t &= \phi_1y_{t-1} + \varepsilon_t = \phi_1(\phi_1y_{t-2}+\varepsilon_{t-1}) + \varepsilon_t\\
    &= \phi_1^2y_{t-2} + \varepsilon_t + \phi_1\varepsilon_{t-1}\\
    &= \phi_1^2(\phi_1y_{t-3} + \varepsilon_{t-2}) + \varepsilon_t + \phi_1\varepsilon_{t-1}\\
    &= \cdots\\
    &= \varepsilon_t + \phi_1\varepsilon_{t-1} + \phi_1^2\varepsilon_{t-2} + \phi_1^3\varepsilon_{t-3} + \cdots
    \end{align}


一般來說，MA($q$)模型有以下特徵：
- ACF圖中在lag$~q~$時有顯著突起，但lag更大時就不再有。
- PACF圖中呈現指數型遞減，或類似於正弦的波動。

使用下列程式碼來模擬MA(1)模型
$$y_t = \varepsilon_t + 0.8\varepsilon_{t-1}$$
```python=+
MA_1 = ArmaProcess([1] , [1 , 0.8])
simulated_MA_1 = MA_1.generate_sample(nsample=1000)
plt.plot(simulated_MA_1)
plt.show()
```

|![output_13_0](https://hackmd.io/_uploads/Sy2SdYsCJg.png)|
|:--:|
|MA(1)模型的模擬|

```python=+
plot_acf(simulated_MA_1, lags = 10)
plt.show()
plot_pacf(simulated_MA_1, lags = 10)
plt.show()
```

|![output_14_0](https://hackmd.io/_uploads/BJwLdKoCJg.png)|![output_14_1](https://hackmd.io/_uploads/ry50_Fo0yx.png)|
|:--:|:--:|
|ACF圖|PACF圖|

從圖中可觀察出
- ACF圖中在lag 1 時有顯著突起，但lag更大時不再有。
- PACF圖中呈現出類似於正弦的波動。

### ARIMA (AutoRegressive Integrated Moving Average)

ARIMA基本上就是將AR模型和MA模型結合起來，只不過多了一個差分的步驟。假設$~y_t^\prime~$代表經過差分$~d~$次之後的數據，那麼ARIMA($p$, $d$, $q$) 模型可以表示為
$$y_t^\prime = c + \phi_1y_{t-1}^\prime + \phi_2y_{t-2}^\prime + \cdots +\phi_py_{t-p}^\prime + 
\theta_1\varepsilon_{t-1} + \theta_2\varepsilon_{t-2} + \cdots +\theta_p\varepsilon_{t-p}+\varepsilon_t$$ 許多模型都是ARIMA模型的特殊情況：
- ARIMA($p$, 0, 0) 就是AR($p$)。
- ARIMA(0, 0, $q$) 就是MA($q$)。
- ARIMA(0, 1, 0) 為隨機漫步(Random Walk)，也就是說
$$y_t - y_{t-1} = c + y_t^\prime = \varepsilon_t \Longrightarrow y_t =c +  y_{t-1} + \varepsilon_t$$ 此時我們對$~y_t~$的預測只能是前一時刻的值加上常數
$$\mathbb{E}[y_t\mid y_1,y_2,...,y_{t-1}] = c + y_{t-1} + \mathbb{E}[\varepsilon_t] = c+ y_{t-1}$$ 股市價格的時間序列常會出現此類模型，因此它是無法被預測的。此即《漫步華爾街》（A Random Walk Down Wall Street）中提到的***隨機漫步假說***。

為了方便，常會將ARIMA($p$, $d$, $q$)表示為
$$\underbrace{(1-\phi_1B -\cdots - \phi_pB^p)}_{\text{AR}(p)}\underbrace{(1-B)^d}_{d~次差分} y_t = c + \underbrace{(1 + \theta_1B + \cdots +\theta_qB^q)\varepsilon_t}_{\text{MA}(q)}$$ 其中$~B~$代表延遲算子，意即$~By_t = y_{t-1}$。類似的，$$B^2y_t = B(By_t) = By_{t-1} = y_{t-2}$$以此類推，可知 $B^ky_t =y_{t-k}$。

- 對於($p$, $d$, $q$)的選擇並不是一件容易的事，首先必須進行差分直到數據變得平穩，之後使用ACF和PACF來輔助選擇。
- 幸好，已經有一些演算法可以幫助我們自動挑選ARIMA模型。
- 在開頭介紹的書中很好的介紹了選擇ARIMA模型的步驟([連結](https://otexts.com/fpp3cn/arima-r-cn.html))

### SARIMA(Seasonal ARIMA)

ARIMA 模型也可以用於季節性的數據。ARIMA$(p, d, q)(P, D, Q)_{m}$ 模型多引入了季節項，其中$~m~$代表一季的資料數，以我們的資料為例就是$~m = 12$。具體公式有點複雜，這裡只舉例說明。

- 例如，ARIMA$(1,1,1)(1,1,1)_{12}$可表示為
$$\underbrace{(1-\phi_1B)}_{\text{AR}(1)}\underbrace{(1-\Phi_1B^{12})}_{\text{Seasonal AR}}\underbrace{(1-B^{12})}_{\text{季節性差分}}\underbrace{(1-B)}_{\text{差分}}y_t = c + \underbrace{(1+\theta_1B)}_{\text{MA}(1)}\underbrace{(1+\Theta_1B^{12})}_{\text{Seasonal MA}}\varepsilon_t $$

AR 或 MA 模型中的季節性部分也可以從ACF和PACF中觀察到。
- 例如，ARIMA$(0,0,0)(0,0,1)_{12}$ 會顯示出:
    1. ACF在 lag 12 上有顯著突起，但沒有其他的顯著突起。
    1. PACF在 lag 為12的倍數時出現指數性遞減。
- 同樣的，ARIMA$(0,0,0)(1,0,0)_{12}$ 會顯示出：
    1. ACF在 lag 為12的倍數時出現指數性遞減。
    1. PACF在 lag 12 上有顯著突起，但沒有其他的顯著突起。

## 實際預測

首先將數據分為訓練集和測試集，其中訓練集由1949年1月至1958年12月的數據組成，而測試集由1959年1月至1960年12月的數據組成。接著使用 log 函數將乘法模型變為加法模型。
```python=+
tr_start , tr_end = '1949-01-01','1958-12-01'
te_start , te_end = '1959-01-01','1960-12-01'
train = np.log(data[tr_start:tr_end]["Passengers"])
test = np.log(data[te_start:te_end]["Passengers"])
```

### SARIMA

首先，可以直接利用`auto_arima`來自動選擇SARIMA模型。
```python=+
model_autoARIMA = auto_arima(train ,
                      trace = True,
                      m = 12,
                      error_action = 'ignore',  
                      suppress_warnings = True, 
                      approximation = False,
                      stepwise = False)
```

### STL + ARIMA
也可利用STL分解去除數據的季節性，之後再選擇ARIMA模型。
```python=+
decomposition = STL(train , period = 12).fit()
deseason_bySTL = train - decomposition.seasonal
```

對數據進行一次差分後，可以利用ADF和KPSS檢定以確認數據是否平穩
```python=+
from statsmodels.tsa.stattools import kpss
from statsmodels.tsa.stattools import adfuller

def kpss_test(timeseries):
    print("Results of KPSS Test:")
    kpsstest = kpss(timeseries, 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)

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)
    
diff_data = deseason_bySTL.diff().dropna()
adf_test(diff_data)
kpss_test(diff_data)
```
由ADF和KPSS之檢定結果得知，在顯著水平為0.05的條件下，差分一次後的數據是平穩的。因此我們可以只進行一次差分( $d = 1$ )。所以只需從 ARIMA($p$ , 1 , $q$) 中挑選適當的模型。
接著使用`auto_arima`自動挑選ARIMA模型
```python=+
model_autoARIMA = auto_arima(deseason_bySTL ,
                             trace = True,
                             d = 1,
                             error_action = 'ignore',  
                             suppress_warnings = True, 
                             pproximation=False,
                             stepwise = False)
```

`auto_arima` 函數得出了ARIMA(0,1,1) 模型。為了確保殘差值是白噪音，對殘差做Ljung-Box檢定。
```python=+
residual = model_autoARIMA.resid()
results = acorr_ljungbox(residual)
```
檢驗$~p~$值非常接近 1，說明沒有顯著證據能推翻殘差是白噪音這個假說。
假設季節性不變，那麼就可以利用過去同一季的資料來當作下一季的季節性值。以下程式碼就是在做這件事：
```python=+
from datetime import timedelta
from datetime import datetime
from dateutil.relativedelta import relativedelta

m = 12
data1 = {
    "date": pd.date_range(start = te_start, end = te_end, freq = "MS"),
    "season" : None
}

df1 = pd.DataFrame(data1)
df1.set_index("date", inplace=True)

data2 = decomposition.seasonal
df2 = pd.DataFrame(data2)

result_merge = pd.concat([df2, df1])

test_start, test_end = datetime.strptime(te_start, "%Y-%m-%d"), datetime.strptime(te_end, "%Y-%m-%d")
diff = relativedelta(test_end, test_start)
diff_months = diff.months + diff.years * 12

for d in range(1, diff_months + 2):
    k = m - d + m * np.floor((d - 1) / m)
    result_merge.loc[test_start + relativedelta(months = d - 1)] = result_merge.loc[test_start - relativedelta(months = k + 1)]
```
之後，將模型預測出的結果和季節性相加就得出了預測值。
```python=+
df = pd.DataFrame(model_autoARIMA.predict(test.shape[0]))
df.columns = ["season"]
pred = result_merge[test_start : test_end].add(df)
```
得出結果後，再將預測值和測試集的數據做比較。
```python=+
error = np.sqrt(mean_squared_error(np.exp(test), np.exp(pred)))
print('Test RMSE: %.4f' % error)
pd.DataFrame({'test': np.exp(test), 'pred': np.exp(pred.squeeze())}).plot();plt.show()
```

|![output_26_1](https://hackmd.io/_uploads/HJWddKjA1l.png)|
|:--:|
|預測值和測試集數據之比較 (Test RMSE: 20.2322)|

將預測值和原本的資料集合併顯示。
```python=+
plt.figure()
line1, = plt.plot(data, label = 'Original Data')             
line2, = plt.plot(np.exp(pred), label = 'Forecast')
plt.legend(handles = [line1, line2], loc='upper left')
plt.show()
```

|![output_27_0](https://hackmd.io/_uploads/HyFdOKsCJg.png)|
|:--:|
|預測值和原始數據|

發現預測值有正確地顯示出季節性，且擬合的還不錯。
前面寫了這麼多，其實也有內建函數 `STLForecast` 可以處理上述所做的事。
```python=+
from statsmodels.tsa.forecasting.stl import STLForecast
from statsmodels.tsa.arima.model import ARIMA

stlf = STLForecast(train, ARIMA, model_kwargs = dict(order = (0, 1, 1), trend = "t") , period = 12)
stlf_res = stlf.fit()
forecast = stlf_res.forecast(24)

plt.figure()
line1, = plt.plot(data, label = 'Original Data')             
line2, = plt.plot(np.exp(forecast), label = 'Forecast with STLForecast')
plt.legend(handles = [line1, line2], loc='upper left')
plt.show()  
```

|![output_28_0](https://hackmd.io/_uploads/SkGtOFo01l.png)|
|:--:|
|預測值和原始數據|