時間序列

函數	程序包	用途
ts()	stats	生成時序對象
plot()	graphics	畫出時間序列的折線圖
start()	stats	返回時間序列的開始時間
end()	stats	返回時間序列的結束時間
frequency()	stats	返回時間序列中時間點的個數
window()	stats	對時序對象取子集
ma()	forecast	擬合一個簡單的移動平均模型
stl()	stats	用LOESS光滑將時序分解為季節項、趨勢項和隨機項
monthplot()	stats	畫出時序中的季節項
seasonplot()	forecast	生成季節圖
HoltWinters()	stats	擬合指數平滑模型
forecast()	forecast	預測時序的未來值
accuracy()	forecast	返回時序的擬合優度度量
ets()	forecast	擬合指數平滑模型，同時也可以自動選取最優模型
lag()	stats	返回取過指定滯后項后的時序
Acf()	forecast	估計自相關函數
Pacf()	forecast	估計偏自相關函數
diff()	base	返回取過滯后項和（或）差分后的序列
ndiffs()	forecast	找到最優差分次數以移除序列中的趨勢項
adf.test()	tseries	對序列做ADF檢驗以判斷其是否平穩
arima()	stats	擬合ARIMA模型
Box.test()	stats	進行Ljung-Box檢驗以判斷模型的殘差是否獨立
bds.test()	tseries	進行BDS檢驗以判斷序列中的隨機變量是否服從獨立同分布
auto.arima()	forecast	自動選擇ARIMA模型

---

生成時序對象

一個數值型向量或數據框中的一列可通過ts()函數存儲為時序對象。

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

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

我們可以通過window()函數生成原始數據的一個子時序。

subset <- window(data, start=, end=)

---

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

時序數據集中通常由很顯著的隨機或誤差成分。為了辨明數據中的規律，我們總是希望能夠撇開這些波動，畫出一條平滑曲線。畫出平滑曲線的最簡單的辦法是簡單移動平均。

比如每個數據點都可用這一點和其前后兩個點的平均值來表示，這就是居中移動平均(centered moving average)。

\[S_t=\frac{Y_{t-q}+\cdots+Y_t+\cdots+Y_{t+q}}{2q+1} \]

其中，\(S_t\) 是時間點 \(t\) 的平滑值，\(k=2q+1\) 是每次用來平均的觀測值的個數。

居中移動平均法的代價是，每個時序集中我們會損失最后的 \(\frac{k-1}{2}\) 個觀測值。

R 中有幾個函數都可以做簡單移動平均，包括 TTR 包中的SMA()函數，zoo 包中的rollmean()函數，forecast 包中的ma()函數。

使用自帶的尼羅河數據集(Nile)為例子。

library(forecast)
opar <- par(no.readonly=TRUE)
par(mfrow=c(2,2))
ylim <- c(min(Nile), max(Nile))
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)

【注】報錯的話可能需要先安裝 forecast 包。

install.packages('forecast')

然后我們就可以得到下圖。

可以發現隨着 \(k\) 的增大，圖像變得越來越平滑。因此我們需要找到最能畫出數據中規律的 \(k\)，避免過平滑或者欠平滑，這里並沒有什么特別的科學理論來指導選取，只能嘗試多個不同的 \(k\) 再決定。

---

季節性分解

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

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

而我們可以通過相加模型或相乘模型來分解數據。

相加模型 \(Y_t=Trend_t+Seasonal_t+Irregular_t\)

相乘模型 \(Y_t=Trend_t\times Seasonal_t \times Irregular_t\)

很多時候，相乘模型比相加模型更現實一些。

將時序分解為趨勢項、季節項和隨機項的常用方法是用 LOESS光滑 做季節性分解。

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

其中，ts是將要分解的時序，參數s.window控制季節效應變化的速度，t.window控制趨勢項變化的速度。較小的值意味着更快的變化速度。令s.window="periodic"可使得季節效應在各年間都一樣。這一函數中，參數ts和s.window是必須提供的。

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

\[log(Y_t)=log(Trend_t \times Seasonal_t \times Irregular_t)=log(Trend_t)+log(Seasonal_t)+log(Irregular_t) \]

以自帶的 AirPassengers 數據集為例子。

> plot(AirPassengers)
> lAirPassengers <- log(AirPassengers)
> plot(lAirPassengers, ylab="log(AirPassengers)")
> fit <- stl(lAirPassengers, s.window="period")
> plot(fit)
> fit$time.series
            seasonal    trend     remainder
Jan 1949 -0.09164042 4.829389 -0.0192493585
Feb 1949 -0.11402828 4.830368  0.0543447685
...output omitted ...
> exp(fit$time.series)
          seasonal    trend remainder
Jan 1949 0.9124332 125.1344 0.9809347
Feb 1949 0.8922327 125.2571 1.0558486
Mar 1949 1.0159924 125.3798 1.0362293
Apr 1949 0.9860703 125.6345 1.0412930
May 1949 0.9850875 125.8897 0.9757094
...output omitted ...

得到的圖如下。

上圖中，序列的趨勢為單調增長，季節效應表明夏季乘客數量更多。每個圖的 y 軸尺度不同，因此通過圖中右側的灰色長條來指示量級，即每個長條代表的量級一樣。

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

通過R中自帶的monthplot()函數和forecast包中的seasonplot()函數可以對季節分解進行可視化。

par(mfrow=c(2,1))
library(forecast)
monthplot(AirPassengers, xlab="", ylab="")
seasonplot(AirPassengers, year.labels="True", main="")

得到下圖。

第一幅圖是月度圖，表示的是每個月份組成的子序列（連接所有 1 月的點、連接所有 2 月的點，以此類推），以及每個子序列的平均值。

第二幅圖是季節圖，這幅圖以年份為子序列。

指數預測模型

指數模型是用來預測時序未來值的最常用模型。它們的短期預測能力較好。不同指數模型建模時選用的因子可能不同。

單指數模型(simple/single exponential model)擬合的是只有常數水平項和時間點\(i\)處隨機項的時間序列，這時認為時間序列不存在趨勢項和季節效應。

雙指數模型(double exponential model; 也叫Holt指數平滑， Holt exponential smoothing)擬合的是有水平項和趨勢項的時序。

三指數模型(triple exponential model; 也叫Holt-Winters指數平滑， Holt-Winters exponential smoothing)擬合的是有水平項、趨勢項以及季節效應的時序。

而使用R中自帶的HoltWinters()函數或者forecast包中的ets()函數可以擬合指數模型。但是ets()函數的備選參數更多，因此更實用。

ets(ts, model="zzz")

其中ts是要分析的時序，限定模型的字母有三個。第一個字母代表誤差項，第二個字母代表趨勢項，第三個字母則代表季節項。可選的字母包括：相加模型(A)、相乘模型(M)、無(N)、自動選擇(Z)。

類型	參數	函數
單指數	水平項	ets(ts, model="ANN")
		ses(ts)
雙指數	水平項、斜率	ets(ts, model="AAN")
		holt(ts)
三指數	水平項、斜率、季節項	ets(ts, model="AAA")
		hw(ts)

ses()、holt()和hw()都是ets()函數的便捷包裝(convenience wrapper)，函數中有事先默認設定的參數值。

單指數平滑

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

單指數平滑模型假定時序中的觀測值可被表示為：

\[Y_t=level+irregular_r \]

在時間點\(Y_{t+1}\)的預測值（一步向前預測, 1-step ahead forecast）可寫作

\[Y_{t+1}=c_0 Y_t+c_1 Y_{t-1}+c_2 Y_{t-2}+\cdots \]

其中\(c_i=\alpha (1-\alpha)^i, i=0,1,2,\cdots\)並且\(0\leq \alpha \leq 1\)。

權數\(c_i\)的總和為1，則一步向前預測可看作當前值和全部歷史值的加權平均。

式中\(\alpha\)參數控制權數下降的速度，\(\alpha\) 越接近於 1，則近期觀測值的權重越大；反之，\(\alpha\) 越接近於 0，則歷史觀測值的權重越大。

為最優化某種擬合標准，\(\alpha\) 的實際值一般由計算機選擇，常見的擬合標准是真實值和預測值之間的殘差平方和。

以自帶的康涅狄格州紐黑文地區的年平均氣溫數據集 nhtemp 為例子。

> library(forecast)
> fit <- ets(nhtemp,model ="ANN")
> fit
ETS(A,N,N) 

Call:
 ets(y = nhtemp, model = "ANN") 

  Smoothing parameters:
    alpha = 0.1819 

  Initial states:
    l = 50.2762 

  sigma:  1.1455

     AIC     AICc      BIC 
265.9298 266.3584 272.2129 
> forecast(fit,1)
     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
1972       51.87031 50.40226 53.33835 49.62512 54.11549
> plot(forecast(fit,1), xlab="Year",ylab=expression(paste("Temperature(", degree*F,")",)),main="New Haven Annual Mean Temperature")
> accuracy(fit)
                    ME     RMSE       MAE       MPE   MAPE      MASE
Training set 0.1460657 1.126268 0.8951225 0.2419373 1.7489 0.7512408
                     ACF1
Training set -0.006441923

可以看到 \(\alpha\) 值比較小(\(\alpha\)=0.18)說明預測時同時考慮了離現在較近和較遠的觀測值，這樣的 \(\alpha\) 值可以最優化模型在給定數據集上的擬合效果。

其中forecast()函數用於預測時序未來的 \(k\) 步，其形式為 \(forecast(fit,k)\) 。得到的預測值是 51.87°F，其 95% 的置信區間為 49.63°F 到 54.12°F，還給出了 80% 的置信區間。

accuracy()函數展示了時序預測中最主流的幾個准確性度量。下表中 \(e_t\) 代表第 \(t\) 個觀測值的誤差項（隨機項），即\((Y_t- \hat{Y_i})\)。

度量標准	簡寫	定義
平均誤差	ME	\(mean(e_t)\)
平均殘差平方和的平方根	RMSE	\(sqrt(mean(e_t^2))\)
平均絕對誤差	MAE	$mean(
平均百分比誤差	MPE	\(mean(100 \times e_t / Y_t)\)
平均絕對百分誤差	MAPE	$mean(
平均絕對標准化誤差	MASE	$mean(

其中\(q_t=e_t/(1/(T-1)\times sum(|y_t-y_{t-1}|))\)，T 是觀測值的個數，對 \(t=2\) 到 \(t=T\) 求累加。

一般來說，平均誤差和平均百分比誤差用處不大，因為正向和負向的誤差會抵消掉。平均絕對百分誤差給出了誤差在真實值中的占比，它沒有單位，因此可以用於比較不同時序間的預測准確性；但它同時假定測量尺度中存在一個真實為零的點。RMSE 是相對最有名、最常用的。

Holt指數平滑和Holt-Winters指數平滑

Holt指數平滑可以對有水平項和趨勢項（斜率）的時序進行擬合。時刻 t 的觀測值可表示為：

\[Y_t=level+slope \times t +irregular_t \]

平滑參數 \(\alpha\)(alpha) 控制水平項的指數型下降，beta 控制斜率的指數型下降。同樣，兩個參數的有效范圍都是 [0,1]，參數取值越大意味着越近的觀測值的權重越大。

Holt-Winters指數光滑可用來擬合有水平項、趨勢項以及季節項的時間序列。此時，模型可表示為：

\[Y_t=level+slope \times t +s_t+irregular_t \]

其中\(s_t\)代表時刻\(t\)的季節效應。除 alpha 和 beta 參數外，gamma 光滑參數控制季節項的指數下降。gamma 參數的取值范圍同樣是 [0,1]，gamma 值越大，意味着越近的觀測值的季節效應權重越大。

仍然使用之前的航空乘客數據集作為例子。

> library(forecast)
> fit <- ets(log(AirPassengers), model="AAA")
> fit
ETS(A,A,A) 

Call:
 ets(y = log(AirPassengers), model = "AAA") 

  Smoothing parameters:
    alpha = 0.6975 
    beta  = 0.0031 
    gamma = 1e-04 

  Initial states:
    l = 4.7925 
    b = 0.0111 
    s = -0.1045 -0.2206 -0.0787 0.0562 0.2049 0.2149
           0.1146 -0.0081 -0.0059 0.0225 -0.1113 -0.0841

  sigma:  0.0383

      AIC      AICc       BIC 
-207.1694 -202.3123 -156.6826 
> accuracy(fit)
                       ME       RMSE        MAE         MPE      MAPE      MASE
Training set -0.001830684 0.03606976 0.02770885 -0.03435608 0.5079142 0.2289192
                   ACF1
Training set 0.05590461
> pred <- forecast(fit,5)
> pred
         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
Jan 1961       6.109335 6.060306 6.158365 6.034351 6.184319
Feb 1961       6.092542 6.032679 6.152405 6.000989 6.184094
Mar 1961       6.236626 6.167535 6.305718 6.130960 6.342292
Apr 1961       6.218531 6.141239 6.295823 6.100323 6.336738
May 1961       6.226734 6.141971 6.311498 6.097100 6.356369
> plot(pred, main="Forecast for Air Travel", ylab="Log(AirPassengers)", xlab="Time")
> pred$mean <- exp(pred$mean)
> pred$lower <- exp(pred$lower)
> pred$upper <- exp(pred$upper)
> p <- cbind(pred$mean, pred$lower, pred$upper)
> dimnames(p)[[2]] <- c("mean", "Lo 80", "Lo 95", "Hi 80", "Hi 95")
> p
             mean    Lo 80    Lo 95    Hi 80    Hi 95
Jan 1961 450.0395 428.5065 417.5279 472.6544 485.0826
Feb 1961 442.5448 416.8301 403.8280 469.8459 484.9735
Mar 1961 511.1312 477.0088 459.8775 547.6945 568.0971
Apr 1961 501.9652 464.6289 446.0019 542.3017 564.9506
May 1961 506.1001 464.9691 444.5667 550.8694 576.1504

可以看到三個光滑參數。趨勢項為 0.0031，它小意味着近期觀測值的斜率不需要更新。

接着預測了接下來五個月的乘客量。圖如下。

此時的預測基於對數變換后的數值，因此要通過冪變換得到預測的乘客量。pred$mean包含了點估計值，其他兩個包含了 80% 和 95% 置信區間的下界以及上屆。cbind()函數用於整合所有結果。

ets()函數和自動預測

ets()函數還可以用來擬合有可乘項的指數模型，加入抑制因子(dampening component)，以及進行自動預測。

時序預測一般假定序列的長期趨勢是一直向上的，而一個抑制項則使得趨勢項在一段時間內靠近一條水平漸近線。在很多問題中，一個有抑制項的模型往往更符合實際情況。

也可以直接使用ets()函數自動選取對原始數據擬合優度最高的模型。

ARIMA預測模型

ARIMA 模型主要用於擬合具有平穩性（或可以被轉換為平穩序列）的時間序列。在平穩的時序中，序列的統計性質並不會隨着時間的推移而改變。

一般來說，擬合 ARIMA 模型前都需要變換序列的值以保證方差為常數。之前用到的對數變換就是一種常用的變換方法，另外常見的還有 Box-Cox 變換。

由於一般假定平穩性時序有常數均值，這樣的序列中肯定不含有趨勢項。非平穩的時序可以通過差分來轉換為平穩性序列。具體來說，差分就是將時序中的每一個觀測值\(Y_t\)都替換為\(Y_{t-1}-Y_t\)。注意，對序列的一次差分可以移除序列中的線性趨勢，二次差分移除二次項趨勢，三次差分移除三次項趨勢。在實際操作中，對序列進行兩次以上的差分通常都是不必要的。

我們可通過diff()函數對序列進行差分，即diff(ts,differences=d)，其中 d 即對序列ts的差分次數，默認值為\(d=1\)。forecast 包中的ndiffs()函數可以幫助我們找到最優的 \(d\) 值，語句為ndiffs(ts)。

平穩性一般可以通過時序圖直觀判斷。如果方差不是常數，則需要對數據做變換；如果數據中存在趨勢項，則需要對其進行差分。

也可以通過ADF(Augmented Dickey-Fuller)統計檢驗來驗證平穩性假定。R中的 tseries 包的adf.test()可以用來做 ADF 檢驗，語句為adf.test(ts)。如果結果顯著，則認為序列滿足平穩性。

相關概念

滯后階數(lag）即我們向后追溯的觀測值的數量。

時序可以通過lag(ts,k)函數變成\(k\)階滯后，其中ts指代目標序列，\(k\)代表滯后項階數。

自相關(autocorrelation)度量時序中各個觀測值之間的相關性。\(AC_k\)即一系列觀測值(\(Y_t\))和\(k\)時期之前的觀測值(\(Y_{t-k}\))之間的相關性。這樣，\(AC_1\)就是一階滯后序列和0階滯后序列間的相關性……這些相關性(\(AC_1,AC_2,\cdots,AC_k\))構成的圖即自相關函數圖（AutoCorrelation Function plot, ACF圖）。ACF 圖可用於為ARIMA模型選擇合適的參數，並評估最終模型的擬合效果。

stats程序包中的acf()函數或者 forecast 包中的Acf()函數可以生成ACF圖。

偏自相關(partial autocorrelation)即當序列\(Y_t\)和\(Y_{t-k}\)之間的所有值(\(Y_{t-1},Y_{t-2},\cdots,Y_{t-k+1}\))帶來的效應都被移除后，兩個序列間的相關性。我們也可以對不同的\(k\)值畫出偏自相關圖。

stats 程序包中的pacf()函數和 forecast 包中的Pacf()函數都可以用來畫 PACF 圖。PACF 圖也可以用來找到ARIMA 模型中最適宜的參數。

ARMA和ARIMA模型

在一個\(p\)階自回歸模型中，序列中的每一個值都可以用它之前\(p\)個值的線性組合來表示：

\[AR(p):Y_t=\mu+\beta_1 Y_{t-1}+\beta_2 Y_{t-2}+\cdots+\beta_p Y_{t-p}+\epsilon_t \]

其中\(Y_t\)是時序中的任一觀測值，\(\mu\)是序列的均值，\(\beta\)是權重，\(\epsilon_t\)是隨機擾動。

在一個\(q\)階移動平均模型中，時序中的每個值都可以用之前的\(q\)個殘差的線性組合來表示：

\[MA(q):Y_t=\mu-\theta_1 \epsilon_{t-1}-\theta_2 \epsilon_{t-2}-\cdots- \theta_q \epsilon_{t-q}+\epsilon_t \]

其中\(\epsilon\)是預測的殘差，\(\theta\)是權重。這里的移動平均和之前說的簡單移動平均不是一個概念。

這兩種方法的混合即\(ARMA(p,q)\)模型，即：

\[Y_t=\mu+\beta_1 Y_{t-1}+\beta_2 Y_{t-2}+\cdots+\beta_p Y_{t-p}-\theta_1 \epsilon_{t-1}-\theta_2 \epsilon_{t-2}-\cdots- \theta_q \epsilon_{t-q}+\epsilon_t \]

---

\(ARIMA(p,d,q)\)模型意味着時序被差分了\(d\)次，且序列中的每個觀測值都是用過去的\(p\)個觀測值和\(q\)個殘差的線性組合表示的。

建立ARIMA模型的步驟包括：

• 確保時序是平穩的；
• 找到一個（或幾個）合理的模型（即選定可能的\(p\)值和\(q\)值）；
• 擬合模型；
• 從統計假設和預測准確性等角度評估模型；
• 預測。

---

選擇ARIMA模型的方法：

模型	ACF	PACF
ARIMA(p,d,0)	逐漸減小到零	在p階后減小到零
ARIMA(0,d,q)	q階后減小到零	逐漸減小到零
ARIMA(p,d,q)	逐漸減小到零	逐漸減小到零

---

擬合模型

可以用arima()函數擬合一個ARIMA模型，其表達式為arima(ts,order=c(p,d,q))。

函數可以返回移動平均項的系數以及模型的AIC值。可以比較AIC值來得到最合理的模型，比較的准則是AIC值越小越好。另外使用accuracy()函數來度量准確性也可以幫助選擇模型。

---

模型評價

一般來說，一個模型如果合適，那模型的殘差應該滿足均值為0的正態分布，並且對於任意的滯后階數，殘差自相關系數都應該為零。也就是說，模型的殘差應該滿足獨立正態分布（即殘差間沒有關聯）。

可以使用下面的代碼進行檢驗。

> qqnorm(fit$residuals)
> qqline(fit$residuals)
> Box.test(fit$residuals, type="Ljung-Box")

    Box-Ljung test

data:  fit$residuals
X-squared = 1.3711, df = 1, p-value = 0.2416

qqnorm()和qqline()函數輸出的是 Q-Q 圖。如果數據滿足正態分布，則數據中的點會落在圖中的線上。

Box.test()函數可以檢驗殘差的自相關系數是否都為零。如果模型的殘差沒有通過顯著性檢驗，那么我們可以認為殘差的自相關系數為零。ARIMA模型能較好地擬合本數據。

---

預測

如果模型殘差不滿足正態性假設或零自相關系數假設，則需要調整模型、增加參數或改變差分次數。

使用 forecast 包中的forecast()函數可以實現對接下來幾年的預測。

使用plot()函數則可以畫出預測圖。圖中黑色的點為預測點的點估計，淺灰色和深灰色區域分別代表 80% 和95% 的置信區間。

ARIMA的自動預測

可以使用 forecast 包中的auto.arima()實現最優ARIMA模型的自動選取。

【參考】

[1] 《R語言實戰》（第2版）