Learn R | 時間序列圖表的可視化

01-28

前言

對於時間序列（Time Series）數據的分析與可視化展示，一直是R語言以及統計學習中重要的一部分，本篇文章將使用ggplot2包和ggfortify包（重點！），對時間序列數據的可視化展示進行詳細的講解。

一、基礎Time Series圖形繪製

1. 前期準備

在R語言中，時間序列數據對應的格式是ts（不同於以往的data.frame）。為此，在進行圖表繪製之前，我們可以使用class()函數進行數據類型的查詢：

# 查詢數據集類型n> class(Nile)n[1] "ts"n

2. Plotting Time Series for ggplot2

之前在學習ggplot2包時，我們使用ggplot()函數也繪製過關於時間序列數據的圖表，但時間值是作為一個變數添加在aes()函數里的，數據集也是普通的『data.frame類型（economics數據集）。而ggplot2是無法直接對ts類型的數據直接進行圖表繪製的，為此，在此之前，我們需要對數據集進行轉化。

# 對ts類型的數據進行df類型轉換n> a <- data.frame(Time=c(time(Nile)),Nile=c(Nile))n> p <- ggplot(a,aes(x=Time,y=Nile))n> p + geom_line(colour = green) + xlab(The Time Series of Date) n+ ylab(The Time Series of Nile)n

經過轉化後，我們可以像對待data.frame的數據一樣，將ts類型的數據進行各種精細化的調整（詳見Learn R | 可視化之ggplot2包（上））

3. Plotting Time Series for ggfortify

ggfortify部分的內容將是本篇文章的精華所在，與ggplot2包相比，ggfortify包對於Time Series圖表的繪製將使用更為簡潔的代碼，並且更容易實現不同需求下的Time Series圖表的繪製。

（1）單變數時間序列

首先，導入包（同時會自動載入ggplot2包）

> library(ggfortify)n

接下來，繪製Time Series基礎圖形（autoplot函數）

> autoplot(Nile,ts.colour = green) + xlab("Date of Nile") n+ ylab("Number of Nile") + ggtitle("The Time Series of Nile")n# ts.linetype可以改變線條形狀n

（2）多變數時間序列

首先，我們載入需要用到的包。

> library(Vars)n

# 需要用到的數據集，包含e、prod、rw和U四個變數n> Canadan e prod rw Un1980 Q1 929.6105 405.3665 386.1361 7.53n1980 Q2 929.8040 404.6398 388.1358 7.70n1980 Q3 930.3184 403.8149 390.5401 7.47n1980 Q4 931.4277 404.2158 393.9638 7.27n1981 Q1 932.6620 405.0467 396.7647 7.37n1981 Q2 933.5509 404.4167 400.0217 7.13n......n

> autoplot(Canada,ts.colour = green) + ggtitle("The Time Series of Canada")n

如果要把四個變數放在同一張表上，則添加（facets = FALSE）選項

> autoplot(Canada,facets = FALSE) + ggtitle("The Time Series of Canada")n

二、Plotting Time Series 進階

1. 使用移動平均進行平滑處理

時序數據集中通常有很顯著的隨機或誤差成分。為了辨明數據中的規律，我們總是希望能夠撇開這些波動，畫出一條平滑曲線。畫出平滑曲線的最簡單辦法是簡單移動平均。比如每個數據點都可用這一點和其前後兩個點的平均值來表示，這就是居中移動平均（centered movingnaverage）。

移動平均的數學表達式為： $S_{t} = ( Y_{t-q} + ... + Y_{t} + ... + Y_{t+q} ) / (2q+1)$

其中， $S_{t}$ 是時間點t的平滑值， $k = 2q + 1$ 是每次用來平均的觀測值的個數，不過，每個時序集中會損失最後的 $(k-1)/2$ 個觀測值

圖形繪製（使用grid.arrange()函數組合圖形）：

> library(forecast)n> a <- autoplot(ma(Nile,3)) n> b <- autoplot(ma(Nile,7)) n> c <- autoplot(ma(Nile,15))n> d <- autoplot(Nile)n> library(gridExtra)n> grid.arrange(d,a,b,c,ncol=2)n

從圖像來看，隨著k的增大，圖像變得越來越平滑。因此我們需要找到最能畫出數據中規律n的k，避免過平滑或者欠平滑。這裡並沒有什麼特別的科學理論來指導k的選取，我們只是需要先嘗試多個不同的k，再決定一個最好的k（關於移動平均的進一步詳細講解，可參考時間序列的相關資料）。

2. 分解時間序列

# stl:Decompose a time series into seasonal, trend and irregular componentsn> autoplot(stl(AirPassengers, s.window = periodic), ts.colour = green)n

3. Plotting with changepoint package

The changepoint package provides a simple approach for identifying shifts in mean and/or variance in a time series.
ggfortify supports cpt object in changepoint package.

> library(changepoint)n> autoplot(cpt.meanvar(Nile),colour=green) n+ ggtitle("The Time Series of Nile")n

4. Plotting with strucchange package

strucchange is an R package for detecting jumps in data.

> library(strucchange)n> autoplot(breakpoints(Nile ~ 1),colour = green)n

四、Time Series的指數預測模型

指數模型是用來預測時序未來值的最常用模型。這類模型相對比較簡單，但是實踐證明它們的短期預測能力較好。

不同的指數模型建模時選取的因子有所不同，主要的指數預測模型包括以下三種：

單指數模型（simple exponential model）：擬合的是只有常數水平項和時間點 $i$ 處隨機項的時間序列，並不考慮時間序列的趨勢項和季節效應。
雙指數模型（double exponential model）：擬合的是有水平項和趨勢項的時序，也叫Holt指數平滑。
三指數模型（triple exponential model)：擬合的是有水平項、趨勢項以及n季節效應的時序。也叫Holt-Winters指數平滑（Holt-Winters exponential smoothing）

1. 單指數平滑（simple exponential model）

單指數平滑根據現有的時序值的加權平均對未來值做短期預測，其中權數選擇的宗旨是使得距離現在越遠的觀測值對平均數的影響越小。

單指數平滑模型假定時序中的觀測值可被表示為： $Y_{t} = level + irregular$

在時間點 $Y_{t+1}$ 的預測值（一步向前預測，1-step ahead forecast）可寫作：

$Y_{t+1} = c_{0}Y_{t} + c_{1} Y_{t-1} + c_{2} Y_{t-2} + ...$

其中， $c_{i} = a (1-a)^{i}$ $a in (0,1)$ 且 $sum_{i=1}^{n}{c_{i} } =1$

一步向前預測可看作當前值和全部歷史值的加權平均。式中α參數控制權數下降的速度，α越接近於1，則近期觀測值的權重越n大；反之，α越接近於0，則歷史觀測值的權重越大。

# 使用forecast包ets函數繪製單指數平滑觀測圖（指定model=ANN)n> AirPassengers_fit <- ets(AirPassengers,model=ANN)n> plot(forecast(AirPassengers_fit,1))n

2. Holt指數平滑與Holt-Winters指數平滑

Holt指數平滑可以對有水平項和趨勢項（斜率）的時序進行擬合。

時刻 $t$ 的觀測值可表示為： $Y_{t} = level + slope times t + irregular_{t}$

Holt-Winters指數光滑可用來擬合有水平項、趨勢項以及季節項的時間序列。

時刻 $t$ 的觀測值可表示為： $Y_{t} = level + slope times t + S_{t} +irregular_{t}$

對nottem數據集進行H-W指數平滑預測（model=AAA)n> fit <- ets(nottem,model=AAA)n> fitnETS(A,A,A) nnCall:n ets(y = nottem, model = "AAA") nn Smoothing parameters:n alpha = 0.0809 n beta = 1e-04 n gamma = 2e-04 nn Initial states:n l = 49.3104 n b = 6e-04 n s=-9.3916 -6.3872 0.4708 7.4714 11.63 12.8874n 9.0503 3.316 -2.6881 -7.0277 -9.9069 -9.4245nn sigma: 2.2578nn AIC AICc BIC n1740.259 1743.016 1799.430 n# alpha：控制水平項的指數型下降n# beta：控制斜率的指數型下降n# gamma：光滑參數控制季節項的指數下降n# 以上三個參數值越大,意味著越近的觀測值的權重越大n# 以上三個參數取值均在[0,1]之間n

# 對接下來五個月的數據進行預測n> nottem_forecast <- forecast(ets(nottem),5)n# 後四列數據分別為80%和95%區間的上下界n> nottem_forecastn Point Forecast Lo 80 Hi 80 Lo 95 Hi 95nJan 1940 39.96893 37.08495 42.85292 35.55826 44.37960nFeb 1940 39.59189 36.70703 42.47675 35.17988 44.00390nMar 1940 42.44175 39.55602 45.32748 38.02840 46.85509nApr 1940 46.72112 43.83452 49.60772 42.30644 51.13580nMay 1940 52.90164 50.01417 55.78911 48.48563 57.31765n> plot(nottem_forecast)n

更新

注1：關於指數預測模型的進一步學習，可以參考時間序列的相關書籍

注2：在時間序列預測模型中，ARIMA預測模型也是一塊非常重要的內容，由於ARIMA模型較為複雜，本篇文章將不涉及這一部分知識，之後會放在時間序列學習的文章里。

五、時間序列其他數據類型的可視化

autoplot函數也可以理解其他的時間序列類型，支持的R包包括以下：

# xts包：xtsn# zoo包：zooregn# tseries包：irtsn# timeSeries包：timeSeriesn

> library(zoo)n> autoplot(as.zooreg(AirPassengers),colour = green)n> library(xts)n> autoplot(as.xts(AirPassengers),colour = green) # 等價於as.zooreg()n

References：

1. Plotting ts objects

2. R語言實戰（第2版）

3. Time Series Graphics

4. RPubs - Plotting Time Series with ggplot2 and ggfortify

5. ggfortify : Extension to ggplot2 to handle some popular packages