McK & Note

＃／Advanced Econometrics II

發表於   |   分類於   |  

Cover

既上一篇 Advanced Econometrics I 介紹了線性迴歸模型後，冬季學期的下半堂課則是針對時間序列資料的分析方法。技術上和上半堂課老實說沒有太多區別，但針對週期性的分析和處理是這門課最有趣的特點。

冬季學期的安排

我們在冬季學期上了三堂課，其中一門 Collaborative Industry Project 其實是做兩個 projects，所以實際上學習只有這兩門課程：

  • Advanced Econometrics II（時間序列）
  • Data Analytics and Mining（SAS 建模）

誠如之前在〈這半年發生的事〉所提到的，這兩門課都不是太偏理論的課程。這一部分是因為當時我們全班都在忙 project，所以教授們的教學方向也更偏操作和應用，或許還很佛心地避免考太過艱深的理論。不過時間序列是個滿有趣的資料型態，我也能就我所學到的皮毛談一談學習心得。如果對時間序列或 ARIMA 模型有點認識的讀者，恐怕會覺得我的分享太淺了，不過如果不介意童言童語的話，還請不吝指正。

老師和教學方法

這門課的老師是 Professor Garland Durham，教學方法是板書、教材和 R 並用，教材主要是之前提過的 Time Series Analysis and Its Applications: With R Examples，由於時間關係，我們只上了前三章多一點。R programming 方面，我們用的套件是作者 Professor Stoffer 寫的 astsa。為了上這門課，我也上了一輪 Datacamp 上的相關課程，其中以 ARIMA Modeling with R 和 Manipulating Time Series Data in R with xts & zoo 最為實用。其他課程其實也不錯，只適合這門課教的內容重疊較少。換句話說，學了以上兩門課程，就和我們所學有八成左右的重疊，剩下的兩成大多是理論的證明和推導，有興趣的讀者也能去讀教材。

事後想想，其實這門課的最終目標，就是讓我們了解如何利用 ARIMA 模型建模，並理解背後原理。前者雖然不難，實務上也大多是靠自動函數解決，但要理解後者還是得花不少時間，這和使用 lm 函數是一樣的道理。我覺得學到後來，印象最深的是時間序列資料的特性，以及處理上要留意的謬誤。

資料產生過程和週期性

時間序列，顧名思義就是和時間有關的資料，更精確一點說，是根據時間排列的資料，所以諸如每日氣溫、銷量，每季度的 CPI、M1，或每年的 GDP 都是時間序列資料。雖然呈現上，時間序列資料和一般資料一樣，也是由行和列構成，只是其中一行是時間，但兩者之間在資料產生過程（data generating process，DGP）上有根本差異。

由於時間序列是根據資料產生順序所排列的一組資料，每一筆資料間有分先後，而且在先後之間或許還有週期性（periodicity，或季節性 seasonality）。如果不考慮先後這點，回想一下線性迴歸中資料的處理，就會發現不論資料產生的時間，觀測值都是根據大小來排列，比方說排成一條斜線；但如果考慮到先後，七月的資料，就得排在二月之後。而且不同時間內，資料可能有週期性，例如七月的氣溫持續比二月高。考慮到這兩個基本特性，資料處理和建模的方法就需要改變。

反過來說，如果我們用一般的方法處理時間序列資料，可能會遺漏一些重要的特性。比方說，如果用一般的線性迴歸法解決季節性的問題，將年和月分別當作連續和不連續的變因建模的話，這樣雖然能區分每個月份的特性，也能得出不同年份的整體趨勢，但這麼做又忽略了順序這項特質，特別是「前者是否對後者有影響」的相關性。畢竟今年八月的氣溫，比起去年八月，可能和今年七月的氣溫更有關聯，上述方法無法從月份之間平均值的區別掌握這層關聯性，也無法從每年的平均值捕捉到兩個月之間的關係。

所以透過另一套方法分析時間序列資料就顯得很重要。某種程度上，我喜歡把時間序列分析想成是一套 data manupulation 的方法。就建模技術上，時間序列資料和一般資料的區別並不大，就設計和假設上所需考慮的面向也很相近，但考慮到上面提到的特性，我們需要一套能正確分析關係的流程，時間序列分析也就應運而生。

訊號和雜訊

如果將上面的敘述用更精準的數學式表達的話，可以先將時間序列資料的組成分為兩者：訊號和雜訊（signal and noise），兩者的區別和線性迴歸提過的一樣，前者是分析的目標，屬於真正的資料產生過程，後者是可忽略、也就是長期平均值為零的雜訊。不過和一般迴歸不同的是，在一般迴歸裡看待兩者的觀點，給我的感覺比較接近「在模型下區分訊號和噪音」，但在時間序列裡考慮到週期性和相關性，觀點比較接近「歸納幾種特定的 DGP 後，用這些關係設計建模方法」。這算是一個滿微妙的區別。

我記得上學期在一般迴歸裡，一開始我認識訊號和雜音的差異很簡單：教授在黑板上點幾個點，在中間畫一條線，說這條線就是模型，點和線之間的距離是雜訊，並接著談最小平方法（OLS）等等；但在這門課裡，我們是先從訊號和雜訊如何構成 AR、MA、RW 等數列講起，才接著說明各種數列的特性，最後沿用之前學過的建模方法。雖然兩種說法對訊號和雜訊的看法完全相同，但後者接近 bottom-up 的做法，或許也說明了為什麼我覺得時間序列的處理，某部分很像融合不同考量的 data manipulation。

在以下的模型表達式中，$x_{t}$ 代表訊號，$w_{t}$ 代表雜訊，兩者共有的 $t$ 則表示時間順序。這邊先提幾個我自己歸納的資料生成規則：

  • 順序：$t$ 決定資料順序，例如一個序列 $X$ 包含從 $x_{1}$ 到 $x_{t}$ 等資料。
  • 關聯：每一筆資料可能包含前一筆資料中不同比例和深度的訊號和雜訊，例如 $x_{t} = \alpha_{1} \cdot x_{t-1} + \alpha_{2} \cdot x_{t-3}$ 或 $x_{t} = \beta_{1} \cdot w_{t-2} + \beta_{2} \cdot w_{t-4}$ 或兩者混合。
  • 平穩：每種數列都需要考慮其週期或平穩性，另外 $w$ 可以被簡單當作白噪音。

這些只是一些簡單的規則，更嚴謹的定義和術語可以參考 Wikipedia 上的時間序列。如果讀者還記得 Advanced Econometrics I 中提到的五個假設，用 Mean Independence 和 Homoskedasticity 其實很好理解時間序列的假設，也能從而推敲何謂 de-trending 和 GARCH 模型所解決的問題。

San Jose

幾種不同的數列

理解了 $x_{t}$ 和 $w_{t}$ 的簡單規則後，以下談由兩者組成的各種數列。理解這些數列的目的，是了解不同資料產生過程所帶來的各種特性，這些會直接影響之後的建模和預測。

White Noise 白噪音

資料產生過程：$w \sim N\left( 0,\sigma_{w}^{2}\right)$

這裡假設 $w$ 符合平均值為零和狹義的平穩性（stationarity），所以我們上課時直接將 $w$ 當作符合 $N(0,\sigma^{2}_{w})$ 分佈的數列。根據 Wikipedia 上的介紹，白噪聲在功率、頻寬等方面有一些有趣的假設，相信這和統計上的理論和方法設計很有關係，但我所上過的幾門課程並沒有詳細說明這點。不過現實生活中有很多時間序列資料，如金融數據，並不呈正態分佈，有著肥尾（fat-tailed）或變異數隨時間變化等特徵，所以白噪音也不全面適用。

Random Walk 隨機漫步

資料產生過程：$x_{t} = c + x_{t-1} + w_{t}$

隨機漫步模型（RW）的特點是累加 $w$，數列本身的期望值不為零，需要透過各項相減還原回 $w$。當 $c$ 不為零時，模型也被稱為 Random Walk with Drift。我上過的課程裡並沒有特別強調 RW，不過透過相減追求平穩性這項特質，在實務上是消除 Drift 滿常見的處理方法。另外 RW 這個詞，相信讓不少讀者聯想到《漫步華爾街》這本書，以及背後的有效市場假說（EMH）。理所當然地，在隨機漫步假說的簡易模型中也能發現 $w$ 的存在，不過這就超出我這個偽金融學生的守備範圍了。

Autoregressive 自迴歸

資料產生過程：$x_{t} = c + \sum _{i=1}^{p}\varphi _{i}x_{t-i} + w_{t} \\= c+ (\varphi _{1} x_{t-1}+\varphi _{2}x_{t-2}\ldots + \varphi _{p}x_{t-p}) + w_{t} $

自迴歸模型（AR）的特點是每一項都帶有過去項的影響，這項設計讓它可以包涵「第 $t$ 項數受過去訊號影響」這項假設。特別的是，由於 AR 模型中每一項都能往以前追溯，針對 AR 模型的平均數和變異數分析，可以先利用迭代法還原回最低項，再計算平均數、確認平穩性等等。

Moving Average 移動平均

資料產生過程：$x_{t} = c + \sum _{i=0}^{p}\theta _{i}w_{t-i} \\= c+ (\theta _{1}w_{t-1}+\theta _{2}w_{t-2}\ldots + \theta _{p}w_{t-p}) + w_{t} $

移動平均模型（MA）和 AR 相似的地方在於包含過去項的關聯，不過 MA 包含的是雜訊而非訊號，有點像 RW 的變形。由於過去項都是雜訊，MA 模型的平均數和變異數分析比較簡單，圖形上比起 AR 看起來也相對穩定。從這邊也能看出 AR 和 MA 的組成不同，兩者可以互補成 ARMA 模型，用於包含過去訊號和雜訊對往後的影響。

ARMA 和省略的內容

資料產生過程：$x_{t} = c + \sum _{i=1}^{p}\varphi _{i}x_{t-i} + \sum _{j=0}^{q}\theta _{j}w_{t-j}$

這就是時間序列分析中使用的 ARMA 模型了，為了區分兩者的範圍，這邊的 $w$ 項數有所更動，不需要和 $x$ 一樣。為了簡化表達，針對這個只包含 $w$、$x$ 和 $c$ 的數學式，我們可以用 $B$ 表達數項中的滯後關係，把它寫成：

$$ (1-\varphi _{1}B-\varphi _{2}B^2- \ldots -\varphi _{p}B^{p}) \cdot x_{t} = c + (1+\theta _{1}B+\theta _{2}B^2+ \ldots +\theta _{q}B^{q}) \cdot w_{t} $$

其中 $ B^{q} \cdot w_{t} = w_{t-q}$，左邊的連加變連減是因為移項。$B$ 可以讀作 Backward，有時也會用 $L$ 代替，稱為 Lag Operator。這只是一種用來簡化 $\sum$ 的寫法，也能讓滯後關係更直觀。我們還可以進一步簡化成下式：

$$ \varphi (B) \cdot x_{t} = c + \theta (B) \cdot w_{t} $$

利用這個模型，時間序列分析能包涵上述特性所帶來的影響。簡化到這一步，相信不難看出我們如何從界定訊號 $x$ 和雜訊 $w$，到考慮到前後影響，加入滯後等設計。不過，雖然這很接近 ARIMA 建模中所使用的模型了，和上一門課一樣，我省略了很多重要的概念和證明。雖然理解這些模型，已經足以使用等一下會提到的函數，以及初步解讀結果，但如果想更精準分析和預測時間序列資料，這些是值得學習的概念：

  • Wold’s Theorem：ARMA 模型適用於時間序列分析的基礎，大意是任何平穩的時間序列資料都能表達為 AR 和 MA 的合體。
  • ACF／PACF：分析數列中前後關聯的方法，用於分辨數列的平穩性。由於時間序列的擬合過程是一段不斷拆解，直到成果接近 $w$ 的過程，利用 ACF 或 PACF 分辨 AR、MA、ARMA 或 WN 很重要。
  • Causality／Invertibility：利用滯後因子分析 AR 或 MA 數列的方法，判斷平穩性和是否收斂。某些情況下 AR 和 MA 可透過 matching coefficient 方法互換。

以上是三個我覺得比較重要的概念，想學得更完整的讀者，可以參考 Wharton 的 Stat910: Time Series Analysis，上面有很詳盡的講義和 R code。

ARIMA 建模

時間序列資料的建模並不難，只要掌握不同模型的特徵就能入門。如前所述，ARIMA 建模本質上就是迴歸分析，只是需要考慮到和一般資料不同的架構。以下先從如何利用 arima.sim() 產生數列，再從另一端談如何建模。

產生資料

了解了不同模型之後，這裡談一下如何產生不同數列。產生時間序列資料的方法有三種：

  1. 利用迴圈產生數列
  2. 利用 filter() 函數
  3. 利用 arima.sim() 函數

第一個方法很簡單，利用 rnorm() 等函數產生必要的雜訊，並定義 $x_{0}$ 等基本項後，就能用 for 迴圈生成向量了。第二種老實說我也不太會用⋯⋯我試著讀過說明文件，但當時沒有很認真想弄懂。這邊想談的主要是第三種，arima.sim() 的使用方法。

在 arima.sim() 模型裡可以定義的變數如下：

  • n：項數
  • model：由 list(order = (p,d,q)) 組成，後三項分別是 AR、差分、MA 的滯後項係數
  • 其他的變數包括雜訊產生方法（例如 rnorm）等等

其中由 (p,d,q) 所產生的模型，就稱為 ARIMA(p,d,q)。前面沒提到的差分項 $d$，讀者可以參考 Wikipedia 上的 ARIMA 模型。ARIMA 的（數階）差分即為 ARMA 模型。由 arima.sim() 產生的數列為 ts 資料型態，想多了解如何處理 ts 的讀者可以參考前面推薦的 Manipulating Time Series Data in R with xts & zoo。

以下是一些 arima.sim() 的範例：

$$ AR(2): x_{t} = 1.5 x_{t−1} − 0.75 x_{t−2} + w_{t} $$
1
2
3
ar_2 <- arima.sim(model = list(order = (2,0,0),
                               ar = c(1.5, -0.75)),
                  n = 1000)
$$ MA(1): x_{t} = w_{t} - 0.8 w_{t-1} $$
1
2
3
ma_1 <- arima.sim(model = list(order = (0,0,1),
                               ma = -0.8),
                  n = 1000)
$$ ARMA(1,1): x_{t} = 0.9 x_{t-1} + w_{t} - 0.4 w_{t-1} $$
1
2
3
4
arma_1_1 <- arima.sim(model = list(order = (1,0,1),
                                   ar = 0.9,
                                   ma = -0.4),
                      n = 1000)
$$ ARIMA(1,1,0): x_{t} - x_{t-1} = y_{t} = 0.9 y_{t-1} + w_{t} $$
1
2
3
arima_1_1_0 <- arima.sim(model = list(order = (1,1,0),
                                      ar = 0.9,
                         n = 1000)

arima 函數和 SARIMA

了解了如何利用 arima.sim() 產生時間序列資料後，就能利用 arima() 分析。arima() 的使用方法和 arima.sim() 大同小異，唯一區別是 arima() 是用來分析資料，後者是用來產生資料。所以在 arima() 裡需要指定的變量包括：

  • y：ts 類別的時間序列資料
  • order：由 list(p,q,r) 組成，即三者的滯後深度
  • seasonal：由 list(P,Q,R) 組成，即整個 ARMA 模型的季節性
  • xreg：外在變量
  • 其他選項包括是否包括 mean、drift 等等

其中 seasonal 是前面沒提到的季節性。延續前面提到的 ARMA 模型，Seasonal ARIMA（SARIMA）模型的型態如下：

加入了 $\Phi$ 和 $\Theta$ 之後，SARIMA 可以包涵長期（$ B^{s} $）的週期性，例如每隔十二個月的變化可以表示為：

可以看出 $\theta$ 決定了前項和後項間 $w$ 的關係，但 $\Phi$ 決定了前後間隔十二項 $x$ 的關係。雖然這個表達式看起來有點複雜，剛入門的讀者可以把這想成是包含季節性的設計。總而言之，在 arima() 函數裡，我們需要指定最多七項（p、d、q、P、D、Q、period）變量，以供 arima() 測定係數。但問題來了，在測試這些數值時，怎樣才算好的模型呢？

模型選擇

我們所學的簡易工作流程如下：

  1. 先用 plot() 或 ts.plot() 觀察散點或線條圖
  2. 如果圖形看起來不太像白噪音，例如斜率不為零或平均不為零，就利用 log()、差分函數 diff()、線性迴歸 lm() 等方式 de-trend 數據
  3. 利用 acf()、pacf() 和 acf2() 等模型判定資料產生過程
  4. 帶入不同的變量，觀察 AIC、BIC 等數值，直到找出最佳模型

可以看出前三項是我前面沒提到，但也滿重要的步驟。我前面有稍微提到「時間序列的擬合過程是一段不斷拆解，直到成果接近 $w$ 的過程」，所以這就是那段讓我感覺像 data manipulation 的步驟，也是需要掌握理論的原因。老實說，第四個步驟就只是 grid-search 和 performance comparison，所以除了一組一組係數慢慢試，也能寫個 for 迴圈批次測試。所以更進階的用法涉及兩個函數：sarima() 跟 auto.arima()。

sarima() 是 astsa 套件中的函數，和 arima() 的功能基本上差不多，不過 sarima() 能提供更詳盡的報告，像是殘差的 ACF、Q-Q plot、Ljung–Box test 等等，也就多了不少可比較的面向。當然，要理解這些報告的意義為何，就和理解 lm() 的報告一樣，最好先熟悉理論。

另一個更強大的函數是 forecast 套件裡的 auto.arima()，顧名思義，這個函數直接包辦了第四步驟的自動化，省下不少一組一組係數測試的時間。不過有時 auto.arima() 返回的結果，不一定是最理想的，所以在更嚴謹的建模中，通常會用 auto.arima() 判斷 AR、MA 係數等建模方向後，再用 sarima() 確認最佳模型。這當然一部分是因為模型間能比較的面向不只一個，所以有時也不存在絕對正確的模型。這也是建模的樂趣（？）之一。

San Diego Zoo

重視時間序列的特性，否則⋯⋯

以上就是我們所學的內容，差不多就是了解 auto.arima() 的運作原理，以及學會一些處理不同狀況、改進模型的知識。除此之外，我真心覺得整堂課學下來最重要的一件事，是如果資料和時間相關，請謹慎處理，不要直接當作一般（截面）資料。讓我講得更清楚一點：

我：如果資料產生的過程和時間相關，別忽略時間序列分析，也別直接當一般資料處理。
某：所以我不能把時間當一個 variable 丟進 lm() 裡就好嗎？
我：不行，因為這樣會忽略許多重要的性質。
某：就算我先用 as.Date() 改變資料型態？
我：還是不行，這跟改不改 Date 沒什麼關聯。平面的 lm() 無法捕捉上述幾種模型的特質。
某：那如果直接用 lm() 發現結果正合我意⋯⋯
我：不！你不能用錯誤的分析方法得出偏頗的結論！
某：可是我完全不會用 arima() 也不了解時間序列分析，怎麼辦？
我：⋯⋯人生很長，你可以慢慢學。
某：我同意。我現在就來去嗑一點 lm() 模型。

會想特別強調這件很基本的事，是因為在之前提到的 Iowa Case Competition 裡，確實有參賽隊伍將時間序列資料當一般資料建模。他們將時間作為一項 variable，和其他的外在因素丟進 lm() 分析以後，得出了「氣溫和銷量呈正相關，而且很顯著」的結論。雖然氣溫的確可能和某些商品的銷量相關，但我們分析的是 3M 的銷售資料⋯⋯於事主管當場就問那支隊伍「氣候暖化那年，我們的銷售額怎麼不增反減？」

聰明的讀者應該一看就知道問題所在。氣溫之所以會和銷量相關，是因為兩者在年復一年的狀況下，都存在季節性等循環，所以他們共同符合「某個時段會上升或下降」的特徵，但這不代表撇除季節性之後，兩者還是有所關聯。也就是說，氣溫和銷量的相似，在於它們有類似的結構，這就像公寓鄰居間雖然房型差不多，但不能因此就說每家的生活很相似一樣。所以，容我再說一次，時間序列分析某種程度上很像 data manipulation，透過假設分離出不同的比較面向後，得出合理的擬合。如果想包含外在因素，請善用 arima() 模型中的 xreg 變量。

眺望遠方的應用

最後，我想說這些內容只算入門。上述幾種數列，和現實中的資料還是有些差異，像是雜訊的分佈不一定常態，這使得擬合的精度不是很高，預測上更是困難重重，誤差大到不實用。所以現今金融市場的應用，又比這篇文章的介紹深入許多，有些是在時間序列的架構下，改善對不同假設的支援，像是前面提到的 GARCH 模型；有些則是結合 Machine Learning 技術，強化模型的 robustness，或是用 RNN 等模型處理時間序列資料的特性。觀察和理解不同技術如何改進時間序列分析，真的是滿有趣的學習經驗，也讓我對仍在和金融工程奮鬥的朋友們充滿敬意。

　
走走拍拍：朋友家好動的 Zoey、San Jose、San Diego Zoo