緒論

國

立 政 治 大 學

‧

N a tio na

l C h engchi U ni ve rs it y

數於此模型中。然而，對於單一誤差的狀態空間模型，目前文獻上並無討論如何以類 似卡爾曼濾波器的方式，來進行狀態的更新。普遍看到的預測步驟，仍是原本的指數 平滑法。

本研究的主要貢獻有兩點。首先是關於指數平滑法在來店人數預測的比較。因該筆 資料具有雙季節性，使用 Taylor [12] 的雙季節指數平滑法，在預測上明顯優於單季節指 數平滑法 (Holt-Winters)。我們也與施佩吟 [15] 擷取每日特定時段的作法進行比較，發 現擷取後的時間序列，因為沒有考慮每日內的週期性，在預測表現上確實比雙季節指 數平滑法差。第二點貢獻是針對單一誤差的狀態空間模型，本研究試著提出類似卡爾 曼濾波器的迭代更新。這裡考慮的是不具任何趨勢或季節的簡單指數平滑法加上一個 外生變數的模型，我們推導出時間更新 (time update) 與測量更新 (measurement update)，

並以模擬實驗評估推導結果的正確性。

本論文分成以下流程: 第一章緒論介紹研究動機與目的；第二章是文獻回顧；第三 章研究方法介紹兩種指數平滑法、卡爾曼濾波以及指數平滑的狀態空間模型；第四章 說明如何使用卡爾曼濾波來對指數平滑法對應的狀態空間模型做狀態更新，並與傳統 的指數平滑模型做比較，最後以模擬的方式來驗證推導的正確性；第五章實證研究介 紹本研究採用的資料，利用該筆資料進行預測與評估；第六章結論與建議針對研究結 果提出結論。

‧ 國

立 政 治 大 學

‧

N a tio na

l C h engchi U ni ve rs it y

第 二 章 文獻回顧

指數平滑是一種預測的方法，不少有效的預測方法都能看見指數平滑的身影。這些 方法都具有將過去的觀測值做加權組合的特性。在設定權重的部分，新的觀測值會比 舊的觀測值來的多。因此“指數平滑”可以反映隨著觀測時間增加，權重呈指數下降 的事實。最早的指數平滑法可追溯到 1944 年，由美國海軍的分析師 Brown [2] 所發現。

同一時期，Holt [4] 也在美國海軍研究辦公室 (ONR) 發展指數平滑法，但 Holt 沒有公開 發表該方法，只記錄在 ONR 備忘錄，直到 2004 年才發表。Holt 關於加法性和乘法性 的季節指數平滑則是由他的學生 Winters [14] 提出。因此具有季節性的 Holt’s linear 法也 被稱為 Holt-Winters 法。

Holt-Winters 法是針對單一週期時間序列的預測方法，這個方法在許多實際應用上 有很好的預測表現。然而，隨著科技進步，我們能夠記錄大量且多樣的資訊，單位尺 度也更小。如果以小時資料來看，每 24 小時會循環一次，又因為 7 日為一週，經過 168 小時又會出現一個循環，對電力需求來說，冬季夏季的使用量大不同，一年 8760 小時又能分出一個循環。對於這種雙週期，甚至多週期的資料，Holt-Winters 法有其 限制。Taylor [12] 為了解決短期電力預測的問題，在 Holt-Winters 法中加入第二個季節 項，提出雙季節指數平滑法，並與季節性 ARIMA 模型作比較。

雙季節指數平滑法在預測短期電力需求上有不錯的表現，如 Ismail, Zahran & El-Metaal [7]，使用雙季節整合移動平均自我迴歸模型 (DSARIMA) 與雙季節指數平滑模 型，預測埃及的電力需求並比較兩種方法的效果。Souza, Barros & Miranda [10] 一樣使

‧ 國

立 政 治 大 學

‧

N a tio na

l C h engchi U ni ve rs it y

用雙季節指數平滑法預測電力需求，並考慮假日以及溫度的影響。另一個例子是 Au, Ma & Yeung [1]，預測基地台每小時的手機網路流量，建立 DSARIMA 的自動化預測流 程，與雙季節指數平滑模型進行比較。

‧ 國

立 政 治 大 學

‧

N a tio na

l C h engchi U ni ve rs it y

第 三 章 研究方法

第一節 Holt-Winters 法

Holt-Winters 法同時考慮了三個因子，分別為水平 (Level)、成長 (Growth) 與季節 (Seasonal)，能夠更有效的預測具有季節性的時間序列資料。

本論文採用具有加法趨勢項及加法季節項的 Holt-Winters 法，定義如下:

(Level) ℓ_t =_α(y_t−s_t−m) + (1−α)(ℓ_t−1+b_t−1) (3.1.1) (Growth) b_t =_β(ℓ_t− ℓt−1) + (1−β)b_t−1 (3.1.2) (Seasonal) s_t =_γ(y_t− ℓt−1−^bt−1) + (1−γ)s_t−m (3.1.3) (Forecast) ˆy_t+_h|t = ℓ_t+hb_t+s_t−m+h⁺_m (3.1.4)

其中α 、β 、γ 分別為水平影響、趨勢影響、季節影響的平滑參數，(3.1.3) 式中的

m 代表一個季節週期的長度，平滑參數反映了水平、趨勢、季節項對於資料的適應速 度。 ˆy_t+h|t 代表向前 h 期預測，其中 h⁺_m = [(h−¹)mod m] +1

‧ 國

立 政 治 大 學

‧

N a tio na

l C h engchi U ni ve rs it y

第二節 雙季節指數平滑法

現實生活中的時間序列可能包含不同長度的季節週期，例如一年的用電量，在一天 24 小時的週期中，上下班時間有著固定的週期，一週七天中又有假日與非假日的週期。

若將時間拉長至兩年以上，又會有冬季夏季的週期，當資料存在著多種不同的週期時，

Holt-Winters 法在這種情況下很難發揮最大功用。Taylor [12] 為了解決短期電力預測的 問題，在 Holt-Winters 法中加入第二個季節項，提出可以處理雙季節性的模型。

不同於 Holt-Winters 法，雙季節指數平滑法有兩個季節週期，m₁為短週期，m₂為 長週期，則公式可改寫成:

ℓ_t = (1−α)(ℓ_t−1+b_t−1) +_α (y_t− (^s(_t−¹⁾_m1 +s⁽_t−²⁾_m

2) ) (3.2.1) b_t =b_t−1+_β(y_t− ℓt−1−b_t−1− (s⁽_t−¹⁾_m1 +s⁽_t−²⁾_m2) ) (3.2.2) s⁽_t¹⁾ = (1−γ_d1)s⁽_t¹₋⁾_m

1+_γd1(y_t− ℓ_t−1−b_t−1−s⁽_t−²⁾_m2) (3.2.3) s⁽_t²⁾ = (1−γ_d2)s⁽_t²₋⁾_m

2+γ_d2(y_t− ℓt−1−^bt−1−^s(_t−¹⁾_m1) (3.2.4) ˆy_t = ℓ_t−1+b_t−1+s⁽_t−¹⁾_m1+s⁽_t−²⁾_m2 (3.2.5) (3.2.1) (3.2.2) 中的 α 、β 為水平項及趨勢項的平滑參數，(3.2.3)(3.2.4) 中的 γd₁ 及 γd₂ 分別為短週期季節項 s⁽_t¹⁾ 及長週期季節項 s⁽_t²⁾ 為水平項及趨勢項的平滑參數。

在長周期中有 m₂ 個季節項，每 m₂ 單位時間更新一次。在較短的周期內有 m₁個 季節項，每 m₁ 個單位時間更新一次。在這模型中並沒有要求 m₂與 m₁有倍數關係。

但是，如果 k =m₂/m₁ ，則較長周期內有 k 個較短的周期。因此，對於單位時間為小 時的資料，七天 (168 小時) 為一個長週期，在經過 168 小時之後，會有 168 個季節項 (s⁽_t+¹⁾₁, s⁽_t¹₊⁾₂, ..., s⁽_t+¹⁾₁₆₈) 會更新一次，一天 (24 小時) 則為短週期，經過 24 小時就會有 24 個季節項(s⁽_t+²⁾₁, s⁽_t+²⁾₂, ..., s⁽_t+²⁾₂₄)^{會更新一次。}

在參數估計的部分，對於長週期的 168 個季節項，每一項的平滑參數皆為γd2，對 於短週期的 24 個季節項，每一項的平滑參數皆為γ_d1。

‧ 國

立 政 治 大 學

‧

N a tio na

l C h engchi U ni ve rs it y

第三節 卡爾曼濾波

Kalman [8] 發表著名的論文，描述了離散數據線性濾波問題的遞歸解決方案。自那 時以來，由於數值計算技術的進步，卡爾曼濾波 (Kalman filter) 一直是廣泛研究和應用 的主題，特別是在自主或輔助導航領域。卡爾曼濾波是一組數學方程，它提供了一種 有效的計算（遞歸）方法來估計過程的狀態，以最小化平方誤差的均值。該濾波器在幾 個方面非常強大：它支持對過去，現在和甚至未來狀態的估計，並且即使未知建模系 統的確切性質，也可以這樣做。

卡爾曼濾波是一種高效率的遞歸濾波器，它能夠從一系列包含雜訊的觀測值中，估 計動態系統的狀態。卡爾曼濾波會根據每個觀測值在不同時間下的值，考慮各時間下 的聯合分布，再對未知變數進行估計，因此會比只以單一觀測值為基礎的估計方式要 準。

卡爾曼濾波的一個典型實例是針對一組有限的、包含噪聲且根據物體位置的觀察序 列（可能有偏差），預測出物體的位置、坐標及速度。在很多工程應用（如雷達、電腦 視覺）中都可以找到它的身影。

首先定義卡爾曼濾波的狀態空間模型:

y_t = Hx_t +εt (3.3.1)

x_t =Ax_t−1+Bu_t−1+w_t−1 (3.3.2) 隨機變數 w 與ε 分別為程序雜訊 (process noise) 以及量測雜訊 (measurement noise)

p(w) ∼^N(0, Q) (3.3.3) p(_ε)∼ N(0, R) (3.3.4) 在 這 邊 假 設 w 與 ε 互相獨立且服從常態分配，Q 為程序雜訊變異量 (process noise covariance)，R 為量測雜訊變異量 (measurement noise covariance)。H 為觀測模型，將狀 態空間映射到測量空間，在現實中，Q 、R 、H 會隨著更新而變化，此處都先假設為定 值。

首先定義 ˆx⁻_t ∈ ^Rn 是已知第 t 步之前的狀態估計下，第 t 步的先驗 ( a priori ) 狀態估 計。而 ˆx_t ∈ ^Rn 為已知觀測值 y_t 下，第 t 步的後驗 ( a posteriori ) 狀態估計。接著定義

‧

time update )，第二是測量更新 ( measurement update )。時間更新負責估計當前狀態以及 誤差共變異數，以獲得下一時間的先驗估計，也就是更新 ˆx⁻_t 與 P_t⁻，測量更新則是負

時間更新方程式 (time update equations)

ˆx⁻_t =A ˆx_t−1+Bu_t−1 (3.3.9) P_t⁻ =AP_t−1A^T+Q (3.3.10) 由於卡爾曼濾波是依賴於前一步的估計值(ˆx_t−1, P_t−1)，因此在初始階段需自定一組初 始值(ˆx₀, P₀)

‧ 國

立 政 治 大 學

‧

N a tio na

l C h engchi U ni ve rs it y

測量更新方程式 (measurement update equations) K_t = ^P

t−H^T

HP_t⁻H^T +R (3.3.11) ˆx_t =ˆx⁻_t +K_t(y_t−H ˆx⁻_t ) (3.3.12) P_t =(I−^KtH)P_t⁻ (3.3.13)

卡爾曼增益是測量更新首先計算的項目，有了卡爾曼增益之後，便能求得最佳化的 後驗估計。測量更新方程式計算完之後，將後驗估計值(ˆx_k, P_k)^{帶入時間更新方程式，}

求得下一期的先驗估計，形成一個循環。

圖 3-1 卡爾曼濾波

此範例引用自 Greg Welch and Gary Bishop [13]

綠色直線為真實值，黑點為加上白噪音後的觀測值，藍線為卡爾曼濾波模擬的結果，

可以發現其收斂速度非常快

‧ 國

立 政 治 大 學

‧

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]

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

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