創新狀態空間模型

國

立 政 治 大 學

‧

N a tio na

l C h engchi U ni ve rs it y

第四節 創新狀態空間模型

Hyndman et al. [6] 在 其 書 中， 介 紹 指 數 平 滑 法 如 何 表 示 成 狀 態 空 間 模 型， 記 為 (E,T,S)，其中 E 為誤差項 (Error)，T 為趨勢 (Trend)，S 為季節性 (Seasonal)。針對誤差 項可以分成加法誤差項的模型以及具有乘法誤差項的模型。舉例來說 ETS(A, N, N) 代 表加法誤差且不含趨勢及季節的指數平滑模型，也就是具有加法誤差項的簡單指數平 滑模型。

Hyndman et al. [6] 將指數平滑法寫成狀態空間模型的型式，公式如下:

觀測方程式: y_t = H^Txt+_εt (3.4.1) 狀態方程式: x_t = Ax_t−1+_gεt−1 (3.4.2) 在觀測方程式中，εt 為誤差項，代表無法被觀測到的部分，也稱為 “innovations”，通常 假設其服從變異數為σ²的常態分配。狀態向量 x_t+₁= (t, bt, st, s_t−1, ..., s_t−m+₁)^{T 包含了} 趨勢及季節。於狀態方程式中，Ax_t−1描述過去的狀態 x_t−1對當前狀態 x_t 的影響。向 量 g 決定了 innovations 對狀態的附加影響程度，稱為 “persistence vector”。H 與 g 的維 度皆為 k x 1, A 為 k × k 矩陣。

Hyndman et al. [6] 將 (3.4.1) (3.4.2) 這種模型稱為 “創新狀態空間模型 (innovations state space model)”。該模型中，觀測方程式與狀態方程式中的誤差項只差了一個乘項g

，不同於 (3.3.1) 與 (3.3.2) 中 εt 與 w_t−1有互相獨立的性質。而由於兩誤差項有線性關 係，以往卡爾曼濾波的更新步驟，在此並不能適用。以下舉兩個例子說明指數平滑法 如何改寫成狀態空間模型。

‧

加法誤差項 Holt-Winters 模型: ETS(A,A,A)

透過改寫 (3.1.1) 式至 (3.1.4) 式的 Holt-Winters 法，可以得到具有加法誤差項、加法 趨勢項及加法季節項的 Holt-Winters 模型，稱為 ETS(A,A,A)，模型公式如下:

觀測方程式: y_t = ℓ_t−1+b_t−1+s_t−m+_εt (3.4.3)

‧

None (N) Additive (A)

None (N)

‧

‧ 國

立 政 治 大 學

‧

N a tio na

l C h engchi U ni ve rs it y

第 四 章

創新狀態空間模型的估計

本章主要討論 (3.4.1)(3.4.2) 式這種狀態空間模型的狀態更新，以下第一節中，我們 先介紹 Hyndman et al. [6] 如何估計參數。第二節提出類似卡爾曼濾波的估計方法，討論 系統狀態的更新。過往使用的卡爾曼濾波更新步驟，其推導過程有利用到觀測與狀態 方程式中的誤差項互相獨立的性質，該性質在創新狀態空間模型不復存在。因此，狀 態更新的式子將有所不同。

第一節 平滑係數選取

關於平滑係數 (_{α, β, γ}) 初始值的選取有幾種方法，第一種是參考過去的經驗，第 二種是試誤法。針對創新狀態空間模型 (3.4.1)(3.4.2) 有第三種方法，利用 MLE 來求得 MSE 最小的初始值。Hyndman et al. [6] 給出一般性的指數平滑模型:

y_t =w(x_t−1) +r(x_t−1)_εt, x_t = f(x_t−1) +g(x_t−1)_εt−1, 有了指數平滑模型便能求得概似函數 (Likelihood function):

L^∗(_{θ, x0}) = nlog(

∑

n t=₁

ε²_t) +2

∑

n t=₁

log|^r(x_t−1)|

透過最小化L^{∗可以求得平滑係數}θ= (_{α, β, γ})^T與初始狀態 x₀= (ℓ₀, b₀, s₀, s₋₁, ..., s_−m+₁)^T。 這邊我們令 w(x_t−1) = H_t^Tx_t−1、r(x_t−1) = 1 、 f(x_t−1) = Ax_t−1，可得到指數平滑模

‧

‧

‧ 國

立 政 治 大 學

‧

N a tio na

l C h engchi U ni ve rs it y

Proof.

P_t⁻ =E[(x_t− ˆx⁻_t )(x_t− ˆx⁻_t )^T|y_t−1]

=E[(Ax_t−1+_ωt−1−A ˆx_t−1)(Ax_t−1+_ωt−1−A ˆx_t−1)^T|y_t−1]

=E[ [A(x_t−1− ^ˆxt−1) +ω_t−1][A(x_t−1− ^ˆxt−1) +ω_t−1]^T|^yt−1]

=AP_t−1A^T+Q

+ E[A(x_t−1− ˆx_t−1)_ωt^T₋₁|y_t−1] +E[ [A(x_t−1− ˆx_t−1)_ω^T_t−1]^T |y_t−1]

=AP_t−1A^T+Q+_αA[(1−^zt−1)ˆx_t−1ˆp_t−1−^Pt−1H_t−1]_2×1[1 0] + _αA([(1−z_t−1)ˆx_t−1ˆp_t−1−P_t−1H_t−1]_2×1[1 0])^T (by Lemma 4)

在 P_t⁻ 部分比原先的卡爾曼濾波多出了兩個期望值 E[A(x_t−1− ˆx_t−1)_ω^T_t−1|y_t−1]與 E[ [A(x_t−1− ˆx_t−1)_ω^T_t−1]^T |y_t−1]，互為轉置矩陣，因此只要計算其中一個即可。為 了讓推導過程簡潔明瞭，會把 E[(x_t−1− ˆx_t−1)_ω^T_t−1|y_t−1] 分成三個 Lemma 來推導，

E[(x_t−1− ˆx_t−1)_ωt^T₋₁|y_t−1]推導部份則為 Lemma 4。

Lemma 1:

E[x_t−1w^T_t−1|^yt−1] =_αyt−1ˆx_t−1−αE[(H_t^T₋₁x_t−1)x_t−1|^yt−1]

Proof.

E[x_t−1w^T_t−1|^yt−1] =E[ [ℓ_t−1

p_t−1 ]

[_αεt−1 0]|^yt−1]

=_αyt−1ˆx_t−1−αE[(H_t^T₋₁x_t−1)x_t−1|y_t−1]

Lemma 2:

ˆx_t−1E[w^T_t−1|^yt−1] =_{α ˆxt−1}y_t−1−α ˆx_t−1[1 1]ˆx_t−1

‧

‧ 國

立 政 治 大 學

‧

N a tio na

l C h engchi U ni ve rs it y

Proof.

E[(x_t−1− ˆx_t−1)_ω^T_t−1|y_t−1] =E[x_t−1w^T_t−1|y_t−1]− ˆx_t−1E[w^T_t−1|y_t−1]

=_αyt−1ˆx_t−1−αE[(H_t^T₋₁x_t−1)x_t−1|y_t−1] (by Lemma 1 and 2)

− (α ˆx_t−1y_t−1−α ˆx_t−1[1 1]ˆx_t−1)

=_{α ˆxt−1}[1 1]ˆx_t−1−αE[(H_t^T₋₁x_t−1)x_t−1|y_t−1]

=_{α ˆxt−1}[1 1]ˆx_t−1−α(ˆx_t−1ˆx_t^T₋₁+P_t−1)H_t−1 (by Lemma 3)

=_{α ˆxt−1}[1 1]ˆx_t−1−α( [

ˆx_t−1ℓˆ_t−1 ˆx_t−1ˆp_t−1 ]

H_t−1)

− αPt−1H_t−1

=_{α ˆxt−1}^ˆℓ_t−1+_{α ˆxt−1}ˆp_t−1−α ˆx_t−1ℓˆ_t−1+_{α ˆxt−1}ˆp_t−1z_t−1

− αPt−1H_t−1

=_{α ˆxt−1}ˆp_t−1−α ˆxt−1ˆp_t−1z_t−1−αPt−1H_t−1

=_α[(1−^zt−1)ˆp_t−1ˆx_t−1−^Pt−1H_t−1]

將 Lemma 4 帶回 P_t⁻ 完成推導。

‧ 國

立 政 治 大 學

‧

N a tio na

l C h engchi U ni ve rs it y

第三節 模擬

本 節 以 模 擬 實 驗 來 討 論 上 一 節 的 理 論 推 導 之 正 確 性。 模 擬 資 料 採 用 R 語 言

“smooth” [11] 套件中的 “sim.es” 函數生成 50 筆 ETS(A, N, A) 資料，週期設為 5 ，初始 狀態設為(ℓ₀, b₀, s₀) = (0.5, 0.5, 0.5)，設定誤差項服從常態分配，初始值設為 0，再使 用 (4.2.3)-(4.2.4) 進行估計。模擬結果顯示，(4.2.3)-(4.2.4) 能適當進行狀態更新。

圖 4-1 模擬結果

‧ 國

立 政 治 大 學

‧

N a tio na

l C h engchi U ni ve rs it y

第 五 章 實證研究

第一節 資料介紹與分析步驟

本論文使用美國某精品服飾店資料，圖 5-1 紀錄 2007 年 1 月 1 日到 2007 年 12 月 30 日每小時的來店人數。因為 4 月 8 日復活節、11 月 22 日感恩節及 12 月 25 日聖誕節三 日未營業，本研究沿用施佩吟 [15] 的資料處理做法，針對放假當日以及放假隔日，用 前五週同一時間資料的平均值進行差補。因為營業時間非 24 小時，且剛開店以及關店 前的人潮不穩定，因此只採用每日的 12 點至 18 點共 7 小時，在雙季節指數平滑模型上 設定的雙季節變為小循環為 7，大循環為 49(7 小時 *7 天)，也就是說經過 7 個時間點就 會出現小更新，經過 49 個時間點會出現大更新。

圖 5-1 來店人數

‧ 國

立 政 治 大 學

‧

N a tio na

l C h engchi U ni ve rs it y

施佩吟 [15] 從每日每時的來店人數資料中，擷取每日同一時段的資料進行分析與預 測。例如時間點為 12 時的資料，這樣的時間序列共有 364 筆資料，透過這樣針對不同 時間點的資料進行切割，所得到的時間序列仍有每週內的週期性，但無每日內的週期 性。本論文不依時間點將資料切割，如此可以保有每日內的週期性資訊，優點請見下 一節。

本研究使用 ETS(A,N,N)、ETS(A,A,N)、ETS(A,A,A)、DSHW(7,49) ((3.4.7)-(3.4.11)) 與包含一個外生變數 day-of-week(DOW) 的 ETS(A,N,N) 模型 ((4.2.1)-(4.2.2)) 進行預測，

並對 Holt-Winters 模型與 DSHW(7, 49) (式 (3.4.7)-式 (3.4.11)) 的殘差做 ACF 圖。

分析工具使用統計程式語言 R ，預測模型採用 “forecast” [5] 套件中的 “dshw” 函數與

“ets” 函數，以向前一期預測 (one-step ahead forecast) 來進行，從 t=1274 開始預測，得 到 t=1275 到 t=2548 的預測值，再用 “forecast” 套件中的 “accuracy” 函數進行評估，選 擇 RMSE、MAE、MAPE 作為預測指標。

‧ 國

立 政 治 大 學

‧

N a tio na

l C h engchi U ni ve rs it y

第二節 模型分析

先針對上一節提到的切割問題作討論。表 5-1(a) 使用未經切割的資料，以向前一期 預測的方式所得到的結果，再將預測結果進行切割。表 5-1(b) 為切割時間序列後的結 果。比較雙季節指數平滑模型 (ds) 的 MAPE 、RMSE 與 MAE ，原始結果明顯優於切割 後結果，Holt-Winters 模型 (hw) 則是部分時段稍微優於切割後結果。

由於每個時段的時間序列都是以日為單位，不用考慮一日內的週期性。從表 5-1(b) 切割後結果可以看出，雙季節指數平滑模型的效果比 Holt-Winters 模型差，在表 5-1(a) 原始結果中，雙季節指數平滑模型的效果明顯優於 Holt-Winters 模型。

表 5-1(a) 原始結果

hwMAE dsMAE hwRMSE dsRMSE hwMAPE dsMAPE hr12 40.7785 12.8239 53.6892 18.3876 75.5364 23.0775 hr13 17.0110 16.5195 22.4340 22.7794 19.3047 18.4896 hr14 21.2681 18.1779 28.4254 24.6962 21.2290 17.6267 hr15 20.9838 17.3916 27.9896 23.9497 18.6574 14.1633 hr16 19.3669 16.9492 26.6916 22.8884 15.7448 12.9095 hr17 19.7582 17.0975 26.9127 24.3946 14.2363 12.6561 hr18 23.8890 18.7278 30.8040 26.1964 17.4071 14.2650

表 5-1(b) 切割後結果

hwMAE dsMAE hwRMSE dsRMSE hwMAPE dsMAPE hr12 12.6006 13.7750 17.9290 18.5899 21.5247 24.5132 hr13 18.4042 24.8268 25.9697 34.4116 19.6669 27.0415 hr14 21.9140 25.2312 31.8695 35.5935 19.4884 23.7716 hr15 23.9387 29.9521 34.9250 41.6261 18.7113 25.0338 hr16 25.3270 32.2892 37.0469 44.0801 18.4891 23.9149 hr17 27.0031 31.2908 39.6480 44.3243 18.2696 21.8242

‧ 國

立 政 治 大 學

‧

N a tio na

l C h engchi U ni ve rs it y

從表 5-2 可以看出，ETS(A,N,N) 跟 ETS(A,A,N) 的預測結果差距不大，因為 ETS(A,N,N) 與 ETS(A,A,N) 兩種模型的差別在於有沒有趨勢影響，剛好這筆來店人數資料沒有明顯 的趨勢，可以預期會有這樣的預測結果。ETS(A,N,N)+DOW 中的 DOW 代表外生變數

“day-of-week” (星期一至星期五設為 0，六日設為 1)，這個模型我們使用上一節所推導 的時間更新 (time update) 與測量更新 (measurement update) 來進行迭代更新。他比前兩 個模型多了一個外生變數 DOW，在預測上有略為改進。以上這三個模型都沒有考慮季 節性。

而 ETS(A,A,A) 考慮了季節性，預測結果明顯改善，再看到雙季節指數平滑模型的 三個預測指標又明顯下降，可見該筆資料的確具有雙季節影響，適合使用雙季節指數 平滑模型來預測。

表 5-2 實驗結果

RMSE MAE MAPE ETS(A,N,N) 47.03177 31.91049 38.11577 ETS(A,A,N) 47.04344 31.90937 38.12456 ETS(A,N,N)+DOW 46.66465 31.59696 37.26808 ETS(A,A,A) 32.43201 23.29364 26.01653 DSHW(_{7, 49}) _{23.43904 16.81247 16.16967}

在圖5-2 雙季節指數平滑模型部分，除了少部分出現高估之外，在季節循環上有不 錯的預測結果。在 Holt-Winters 模型預測上，出現不少低估的情況。在一開始虛線包住 的區段為平日週三，卻出現假日的人潮，經查詢之後發現該日為 7 月 4 日美國獨立紀念 日。

看到圖5-3，黑色為真實來店人數，在 12 月 30 日附近，因為接近跨年的關係，店家 的人潮快速上升。Holt-Winters 模型在低點經常出現低估的情況，還出現小於 0 的預測 值。雖然雙季節指數平滑模型在高點低點都有高估低估的狀況，但與真實值相距不遠。

‧ 國

立 政 治 大 學

‧

N a tio na

l C h engchi U ni ve rs it y

圖 5-2 Holt-Winters vs. DSHW 1-250

圖 5-3 Holt-Winters vs. DSHW 900-1274

‧ 國

立 政 治 大 學

‧

N a tio na

l C h engchi U ni ve rs it y

圖 5-4 縮小區間範圍至 1200-1274 ，後期出現震幅放大的情形，Holt-Winters 模型在 這段時間的走勢中有晚一期的現象，人數下降時還保持上升的趨勢，雙季節指數平滑 模型在這方面表現穩定，幾乎能抓到上升下降的時機。

圖 5-4 Holt-Winters vs. DSHW 1200-1274

‧ 國

立 政 治 大 學

‧

N a tio na

l C h engchi U ni ve rs it y

圖 5-5 為兩模型的殘差圖，先看到 Holt-Winters 殘差圖，每隔一段時間殘差會異常 變大，在這裡推測可能是因為第二個季節因素導致預測結果不理想，在頭尾有較大 的變動，稍後會與雙季節指數平滑模型一起討論。再看到雙季節指數平滑模型殘差 圖，其他時間的殘差明顯比 Holt-Winters 模型小且穩定，頭尾部分一樣有不規則變動。

Holt-Winters 模型與雙季節指數平滑模型兩種方法的殘差在t ∈ (^{300, 800})之間有較小變 異，但在頭 (1, 250)、尾 (900, 1274) 部分，兩者的殘差就有明顯的起伏。在雙季節指數 平滑模型殘差圖第 237 的殘差突然增加且為正值，表示在該點的預測出現高估的情況，

推測該商店可能發生什麼事。

圖 5-5 殘差圖

‧ 國

立 政 治 大 學

‧

N a tio na

l C h engchi U ni ve rs it y

從圖 5-6 Holt-Winters 殘差的 ACF 圖可以看出預測結果還有一個季節影響，週期約 為 49，符合上面討論的長週期長度。DSHW 殘差的 ACF 則落在標準內，但仍能看出具 有些微週期性，或許是平日假日的季節影響。

圖 5-6 殘差 ACF

‧ 國

立 政 治 大 學

‧

N a tio na

l C h engchi U ni ve rs it y

第 六 章

結論與建議

在三種傳統的指數平滑模型比較中，Holt-Winters 模型明顯優於簡單指數平滑模型 與 Holt’s linear 模型。接著與 Taylor 的雙季節指數平滑模型作比較，雙季節指數平滑更 優於 Holt-Winters 模型，對於年末突然增長的時段依舊保有不錯的預測結果。

本論文將卡爾曼濾波器的推導邏輯，應用到觀測與狀態方程式中，誤差項具有線性 關係的創新狀態空間模型，推導得到一組新的更新方程式。因為創新狀態空間模型是 指數平滑法的一種模型表示法，目前文獻中，對於創新狀態的更新步驟，基本上仍沿 用指數平滑法的操作步驟，在處理外生變數的情況較為棘手。有了類似卡爾曼濾波器 的更新方法後，處理外生變數就容易許多。未來只要擴展狀態向量的維度，根據指數 平滑法的公式做調整，將相對應的參數放入矩陣 H_t ，就能套用到其他的指數平滑法。

若能套用至雙季節指數平滑模型，加上國定假日或競爭商家的資訊作為外生變數，預 期能得到更好的預測結果。

由於研究資料的時間長度只有一年，無法比較歷年的季節變化，另外從圖 6-1 可以 看出該筆資料可能存在兩個以上季節循環。設定星期一至星期日為一個區段，圖 6-1 中 以虛線作為分段點，可以看出每個區段包含七個尖點，前四個尖點約為 100 人，也就 是星期一至星期四的人數約為 100 人，第五至第七尖點明顯上升，因此一週七天又能 再分成兩個小季節循環。星期六的人數明顯高於其他天，或許能將星期六單獨視為一 個季節循環。

‧ 國

立 政 治 大 學

‧

N a tio na

l C h engchi U ni ve rs it y

圖 6-1 多季節循環

‧ 國

立 政 治 大 學

‧

N a tio na

l C h engchi U ni ve rs it y

參考文獻

[1] S. Tom Au, Guang-Qin Ma, and Shu-Ngai Yeung. Automatic forecasting of double sea-sonal time series with applications on mobility network traffic prediction. 2011.

[2] Robert Goodell Brown. Statistical forecasting for inventory control. McGraw-Hill, New York, 1959.

[3] Phillip G Gould, Anne B Koehler, J Keith Ord, Ralph D Snyder, Rob J Hyndman, and Farshid Vahid-Araghi. Forecasting time series with multiple seasonal patterns. European Journal of Operational Research, 2008.

[4] Charles C Holt. Forecasting seasonals and trends by exponentially weighted moving av-erages. Carnegie Institute of Technology, Graduate school of Industrial Administration, 1957.

[5] Rob J Hyndman and Yeasmin Khandakar. Automatic time series forecasting: the forecast package for R. Journal of Statistical Software, 26(3):1–22, 2008.

[6] Rob J. Hyndman, Anne B. Koehler, J. Keith Ord, and Ralph D. Snyder. Forecasting with Exponential Smoothing:The State Space Approach. Springer, 2008.

[7] Mohamed A. Ismail, Alyaa R. Zahran, and Eman M. Abd El-Metaal. Forecasting hourly electricity demand in egypt using double seasonal autoregressive integrated moving aver-age model. 2015.

[8] R. E. Kalman. A new approach to linear filtering and prediction problems. Journal of

在文檔中 應用雙季節指數平滑模型於來店人數預測之研究 - 政大學術集成 (頁 16-0)