1 時序的平滑化和季節性分解
對時序數據建立復雜模型之前也需要對其進行描述和可視化。在本節中,我們將對時序進行平滑化以探究其總體趨勢,並對其進行分解以觀察時序中是否存在季節性因素。
1.1 通過簡單移動平均進行平滑處理
時序數據集中通常有很顯著的隨機或誤差成分。為了辨明數據中的規律,我們總是希望能夠撇開這些波動,畫出一條平滑曲線。畫出平滑曲線的最簡單辦法是簡單移動平均。比如每個數據點都可用這一點和其前后兩個點的平均值來表示,這就是居中移動平均處理。
時序數據的第一步是畫圖。這里介紹Nile數據集。這一數據集是埃及阿斯旺市在1871年至1970年間所記錄的尼羅河的年度流量。R中有幾個函數都可以做簡單移動平均,包括TTR包中的SMA()函數,zoo包中的rollmean()函數,forecast包中的ma()函數。這里我們用R中自帶的ma()函數來對Nile時序數據進行平滑處理。
如下代碼給出了時序數據的原始數據圖,以及平滑后的圖(對應k=3、7和15),生成的圖像如下圖所示。
opar <- par(no.readonly=TRUE)
par(mfrow=c(2,2)) #把畫布分成2*2四個區域
ylim <- c(min(Nile), max(Nile)) #ylim()分別指定刻度的下限和上限
plot(Nile, main="Raw time series")
plot(ma(Nile, 3), main="Simple Moving Averages (k=3)", ylim=ylim)
plot(ma(Nile, 7), main="Simple Moving Averages (k=7)", ylim=ylim)
plot(ma(Nile, 15), main="Simple Moving Averages (k=15)", ylim=ylim)
par(opar)
結果分析:圖左上角畫出了這一數據集的原始數據信息。從上圖來看,數據總體呈下降趨勢,但不同年份的變動非常大。
從圖像來看,隨着k的增大,圖像變得越來越平滑。因此我們需要找到最能畫出數據中規律的k,避免過平滑或者欠平滑。這里並沒有什么特別的科學理論來指導k的選取,我們只是需要先嘗試多個不同的k,再決定一個最好的k。從本例的圖像來看,尼羅河的流量從1892年到1900年有明顯下降;其他的變動則並不是太好解讀,比如1941年到1961年水量似乎略有上升,但這也可能只是一個隨機波動。
1.2 季節性分解
對於間隔大於1的時序數據(即存在季節性因子),我們需要了解的就不僅僅是總體趨勢了。此時,我們需要通過季節性分解幫助我們探究季節性波動以及總體趨勢。
存在季節性因素的時間序列數據(如月度數據、季度數據等)可以被分解為趨勢因子、季節性因子和隨機因子。趨勢因子(trend component)能捕捉到長期變化;季節性因子(seasonal component)能捕捉到一年內的周期性變化;而隨機(誤差)因子(irregular/error component)則能捕捉到那些不能被趨勢或季節效應解釋的變化。
此時,可以通過相加模型,也可以通過相乘模型來分解數據。在相加模型中,各種因子之和應等於對應的時序值,即:
其中時刻t的觀測值即這一時刻的趨勢值、季節效應以及隨機影響之和。
而相乘模型則將時間序列表示為:(即趨勢項、季節項和隨機影響相乘)
這里通過一個小例子進一步說明相加模型與相乘模型的區別。假設我們有一個時序,記錄了10年來摩托車的月銷量。在可加模型中,11月和12月(聖誕節)的銷量一般會增加500,而1月(一般是銷售淡季)的銷量則會減少200。此時季節性的波動量和當時的銷量無關。在相乘模型中,11月和12月的銷量則會增加20%,1月的銷量減少10%,即季節性的波動量和當時的銷量是成比例的。這也使得在很多時候,相乘模型比相加模型更現實一些。
將時序分解為趨勢項、季節項和隨機項的常用方法是用LOESS光滑做季節性分解。這可以通過R中的stl()函數實現:
stl(ts, s.window=, t.window=)
#其中ts是將要分解的時序,參數s.window控制季節效應變化的速度,t.window控制趨勢項變化的速度。較小的值意味着更快的變化速度。令s.windows="periodic"可使得季節效應在各年間都一樣。這一函數中,參數ts和s.windows是必須提供的。
雖然stl()函數只能處理相加模型,但這也不算一個多嚴重的限制,因為相乘模型總可以通
過對數變換轉換成相加模型:
log(Yt) = log(Trendt * Seasonalt * Irregulart)
= log(Trendt) + log(Seasonalt) + log(Irregulart)
用經過對數變換的序列擬合出的相加模型也總可以再轉化回原始尺度。下面給出一個例子。
R中自帶的AirPassengers序列描述了1949~1960年每個月國際航班的乘客(單位:千)。
序列圖中的第一幅圖。從圖像來看,序列的波動隨着整體水平的增長而增長,即相乘模型更適合這個序列。
序列圖中的第二幅圖是經過對數變換后的序列。這樣序列的波動就穩定了下來,對數變換后的序列就可以用相加模型來擬合了。
用stl()函數做季節性分解,代碼如下:
查看數據集:
plot(AirPassengers) #查看數據集原始圖像
lAirPassengers <- log(AirPassengers)
plot(lAirPassengers, ylab="log(AirPassengers)") #畫出取對數log之后的圖像
fit <- stl(lAirPassengers, s.window="period") #分解時間序列
#用stl()函數做季節性分解,s.window="period"用於季節性分解提取
plot(fit) #畫出季節性分解圖
結果分析:上圖給出了1949~1960年的時序圖、季節效應圖、趨勢圖以及隨機波動項。注意此時將季節效應限定為每年都一樣(即設定s.window="period")。序列的趨勢為單調增長,季節效應表明夏季乘客數量更多(可能因為假期)。每個圖的y軸尺度不同,因此我們通過圖中右側的灰色長條來指示量級,即每個長條代表的量級一樣。
fit$time.series #每個觀測值各分解項的值
# stl()函數返回的對象中有一項是time.series,它包括每個觀測值中的趨勢、季節以及隨機效應的具體組成。此時,直接用fit$time.series則返回對數變換后的時序,
exp(fit$time.series) #對每個觀測值各分解項的值取指數,將結果轉化為原始尺度
結果分析:觀察季節效應可發現,7月的乘客數增長了24%(即乘子為1.24),而11月的乘客數減少了20%(即乘子為0.8)
我們還可以通過兩幅圖來對季節分解進行可視化,即用R中自帶的monthplot()函數和forecast包中的seasonplot()函數來畫圖:
opar <- par(no.readonly=TRUE)
par(mfrow=c(2,1)) #把畫布分成兩個區域
library(forecast)
monthplot(AirPassengers, xlab="", ylab="")
# monthplot()函數繪制一個時間序列的季節性(或其他)子序列,對於每個季節(或其他類別),都會繪制出一個時間序列
seasonplot(AirPassengers, year.labels="TRUE", main="")
# seasonplot()繪制季節圖,這就像一個時間圖,只是數據是根據不同年份的季節來繪制的。
par(opar)
結果分析:這是AirPassengers序列的月度圖(上)和季度圖(下),從兩幅圖中都可以看出總體的增長趨勢以及相似的季節模式。第一幅圖是月度圖,表示的是每個月份組成的子序列(連接所有1月的點、連接所有2月的點,以此類推),以及每個子序列的平均值。從這幅圖來看,每個月的增長趨勢幾乎是一致的。另外,我們還可以看到7月和8月的乘客數量最多。第二幅圖是季節圖(season plot),這幅圖以年份為子序列。從這幅圖中我們也可以觀測到同樣的趨勢性和季節效應。