《R語言實戰》第四部分第十五章-時間序列學習筆記（I）

01-28

之前的大部分章節內容，我們所探討的都是橫截面數據（cross-sectional data），所謂橫截面數據就是在一個給定的時間點測量變數值。與之相對應的是縱向數據（longitudinal data），所謂縱向數據就是隨著時間變化反覆測量的變數值。

本章我們將重點研究在給定的一段時間內有規律記錄的觀測數據。對於這樣的數據我們規定如下：

T表示時間序列中觀測值的個數；
$Y_{i} left( i = 1, 2, 3...t,...T right)$ 表示時間序列；
$Y_{t}$ 為 $Y$ 在時間點 $t$ 的值。

下圖為兩個完全不同時間序列的實例：

對於時間序列的研究主要包括兩個基本問題：

數據描述：即這段時間內都發生了什麼
預測未來：即接下來將會發生什麼

對時序數據的研究在現實世界中有著廣泛的應用：

經濟學家嘗試通過時序分析理解並預測金融市場；
城市規劃者基於時序數據預測未來的交通需求；
氣候學家通過時序數據預測全球氣候變化；
公司需要時序分析來預測產品的需求及未來銷量；
醫療保健人員需要根據時序數據研究疾病傳播範圍及某區域內可能出現的病例數；
地震學家通過時序數據預測地震。

描述時序數據和預測未來值得方法很多，R語言更有很多其他軟體所不具備的精細時序分析工具，本章將對一些最常用的時序數據描述和預測方法以及所對應的R函數進行介紹，下表是本章中按出現順序給出了這些函數。

不同類型的時序數據集特點各不相同，適用不同的模型和分析方法，下表給出了本章中分析的幾個時序數據集，這些數據集在R中都可以直接找到。

1 在R中生成時序對象

在分析時間序列之前，我們需要將分析對象轉成時間序列對象（time-series object），可以將其理解為R中一種特殊的數據結構，它包含觀測值、起始時間、終止時間以及周期（如月、季度或年）等。我們分析時間序列所依賴的分析、建模和繪圖方法只對時間序列對象才起作用。比如一個數值型向量或者數據框中一列（數值型）可通過以下代碼轉換為時序對象：

myseries <- ts(data, start=, end=, frequency=)n

其中：

myseries是所生成的時序對象；
data是原始的包含觀測值的數值型向量；
start參數和end參數（可選）給出時序的起始時間和終止時間；
frequency為每個單位時間所包含的觀測值數量（如frequency=1對應年度數據，frequency=12對應月度數據，frequency=4對應季度數據）。

下列代碼給出了一個示例：

> sales <- c(18, 33, 41, 7, 34, 35, 24 ,25, 21, 25, 20, 22, 31, 40, 29, 25, 21, 22, 54, 31, 25, 26, 35)n#生成時序對象n> tsales <- ts(sales, start = c(2003,1),frequency = 12)n> tsalesn Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Decn2003 18 33 41 7 34 35 24 25 21 25 20 22n2004 31 40 29 25 21 22 54 31 25 26 35 n#繪製時序對象圖形n> plot(tsales)n#獲得這個對象的信息n> start(tsales)n[1] 2003 1n> end(tsales)n[1] 2004 11n> frequency(tsales)n[1] 12n#對時序對象取子集n> tsales.subset <- window(tsales, start= c(2003,5), end = c(2004, 6))n> tsales.subsetn Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Decn2003 34 35 24 25 21 25 20 22n2004 31 40 29 25 21 22 n> n

2 時序的平滑化和季節性分解

在建模之前對時序數據進行描述和可視化等同於對橫截面數據集分析和建模之前的描述性統計和畫圖一樣。本節主要對時序進行平滑化，並對其進行分解以探求是否存在季節性因素。

2.1 通過簡單移動平均進行平滑處理

本節以Nile數據集為例，該數據集1871年至1970年間埃及阿斯旺市所記錄的尼羅河的年度流量，首先對其畫圖，代碼如下：

> install.packages("forecast")n> library(forecast)n> opar <- par(no.readonly = TRUE)n> par(mfrow = c(2,2))n> ylim <- c(min(Nile), max(Nile))n> plot(Nile, main = "Raw time series")n> plot(ma(Nile, 3), main="Simple Moving Averages (k=3)", ylim=ylim)n> plot(ma(Nile, 7), main="Simple Moving Averages (k=7)", ylim=ylim)n> plot(ma(Nile, 15), main="Simple Moving Averages (k=15)", ylim=ylim)n> par(opar)n> n

左上圖要原始數據集的圖形，數據總體呈下降趨勢，但是不同年份變化非常大。為了找出其中的規律，撇開時序數據集當中的隨機或誤差成分，我們通常採用簡單移動平均畫出一條平滑曲線。其中居中移動平均（centered moving average）就是將每一個數據點用這一點和其前後兩個點的平均值來表示，數學表達式如下：

其中 $S_{t}$ 是時間點 $t$ 的平滑值， $k=2q+1$ 是每次用來平均的觀測值的個數，通常設為一個奇數。顯然，居中移動平均法是有代價的，會損失最後的 $left( k-1 right) /2$ 個觀測值。

R中TTR包中的SMA()函數、zoo包中的rollmean()函數以及forecast包中的ma()函數都可以做簡單移動平均，上例中用R自帶的ma()函數對Nile時序數據進行平滑處理，上圖右上及第二列的平滑後的圖形，對應k=3、7和15。不難發現，隨著k的增大，圖像變得越來越平滑。為了避免過平滑或者欠平滑，需要選擇一個合適的k。對於k的選擇，沒有特別好的理論或者方法來進行指導，只有本著對數據的理解進行不斷的嘗試多試幾個不同的數值。

2.2 季節性分解

存在季節性因素（周期性）的時間序列數據（如月度數據、季度數據等）可以被分解為趨勢因子、季節性因子和隨機因子。其中：

趨勢因子（trend component）能捕捉到長期變化；
季節性因子（seasonalcomponent）能捕捉到一年內的周期性變化；
隨機（誤差）因子（irregular/error component）則能捕捉到那些不能被趨勢或季節效應解釋的變化。

可以通過相加模型或者相乘模型來分解數據，在相加模型當中：

在相乘模型當中：

下圖給出了對應的實例：

圖（a）是一個隨機因子；
圖（b）是一個向上的趨勢因子疊加一個隨機因子；
圖（c）是一個季節性因子疊加一個隨機因子；
圖（d）是季節性因子、隨機因子及向上趨勢因子的疊加，適用相加模型；
圖（e）是三個因子相乘模型的疊加，波動與趨勢成正比。

在R中常用LOESS光滑做季節性分解，一般通過R中的stl()函數實現，格式如下：

stl(ts, s.window=, t.window=)n

其中，ts是將要分解的時序，s.window控制季節效益變化的速度，t.window控制趨勢項變化的速度。如果令s.window="periodic"可使得季節效應在各年間都一樣。

stl()函數只能處理相加模型，不過通過相乘模型可以通過對數變換轉換成相加模型：

用經過對數變換的序列擬合出的相加模型也總可以再轉化回原始尺度。

以R中自帶的AirPassengers時序數據集為例，該數據集描述了1949~1960年每個月國際航班的乘客（單位：千）。首先進行繪圖，原序列圖形及對數變換後的圖形如下：

從原序列圖形上看，序列的波動隨著整體水平的增長而增長，適合於相乘模型。經過對數變換後的下圖序列的波動就平穩下來了，變換的序列就可以用相加模型來擬合。以下用相加模型對變換後的時序數據進行處理分析，代碼和圖形如下：

#畫出原時間序列n> plot(AirPassengers)n#對數變換n> lAirPassengers <- log(AirPassengers)n#畫出對數變換後的時間序列n> plot(lAirPassengers, ylab="log(AirPassengers)")n#分解時間序列n> fit <- stl(lAirPassengers, s.window="period")n#畫出分解序列的結果n> plot(fit)n> fit$time.seriesn seasonal trend remaindernJan 1949 -0.09164042 4.829389 -0.0192493585nFeb 1949 -0.11402828 4.830368 0.0543447685nMar 1949 0.01586585 4.831348 0.0355884457nApr 1949 -0.01402759 4.833377 0.0404632511nMay 1949 -0.01502478 4.835406 -0.0245905300nJun 1949 0.10978976 4.838166 -0.0426814256nJul 1949 0.21640041 4.840927 -0.0601151688nAug 1949 0.20960587 4.843469 -0.0558624690n... output omitted ...n> exp(fit$time.series)n seasonal trend remaindernJan 1949 0.9124332 125.1344 0.9809347nFeb 1949 0.8922327 125.2571 1.0558486nMar 1949 1.0159924 125.3798 1.0362293nApr 1949 0.9860703 125.6345 1.0412930nMay 1949 0.9850875 125.8897 0.9757094nJun 1949 1.1160434 126.2377 0.9582166nJul 1949 1.2415994 126.5866 0.9416561nAug 1949 1.2331919 126.9088 0.9456692n... output omitted ...n

上圖為1949~1960年的時序圖、季節效應圖、趨勢圖以及隨機波動項。由於我們設定了s.window="period"，所以趨勢為單調增長，且每年都一樣。

此外，stl()函數返回的對象中有一項是time.series，它包含了每個觀測值中的趨勢、季節以及隨機效應的具體組成，fit$time.series返回對數變換後的時序，exp(fit$time.series)可將結果轉化為原始尺度。