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

01-27

$Y_{t} =level+slopetimes t+irregular_{t}$ 因為上一篇文章老是保存失敗所以分開兩篇……

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

我們還可以用R中自帶的monthplot()函數和forecast包中的seasonplot()函數來畫圖對季節分解進行可視化，代碼和圖形如下：

par(mfrow=c(2,1))nlibrary(forecast)nmonthplot(AirPassengers, xlab="", ylab="")nseasonplot(AirPassengers, year.labels="TRUE", main="")n

上圖第一幅圖是月度圖，表示每個月份的子序列及其平均值。從這幅圖來看，每個月的增長趨勢幾乎是一致的，且7月和8月的乘客數量最多。第二幅圖是季節圖，以年份為子序列，我們可以觀測到同樣的趨勢性和季節效應。

3 指數預測模型

指數模型常用於預測時序未來值，模型相對簡單，但是短期預測效果較好。根據因子不同可能需要選用不同的指數模型：

單指數模型（simple/single exponential model）——不存在趨勢項和季節效應的時間序列；
雙指數模型（double exponential model），又稱Hold指數平滑（Holt exponential smoothing）——同時存在水平項和趨勢項的時間序列；
三指數模型（triple exponential model），又稱Holt-Winters指數平滑（Holt-Winters exponential smoothing）——同時存在水平項、趨勢項以及季節效應的時間序列。

R中自帶的HoltWinters()函數或者forecast包中的ets()函數都可以擬合指數模型，因為ets()函數備選參數更多，更為實用，本節只討論ets()函數，格式如下：

ets(ts, model = "ZZZ")n

其中ts是要分析的時序對象，ZZZ是限定模型的字母選項：

第一個字母代表誤差項，第二個字母代表趨勢項，第三個字母則代表季節項。具體選項如下：

A——相加模型
M——相乘模型
N——無
Z——自動模型

常用模型羅列如下：

ses()、holt()、以及hw()函數都是ets()函數的簡化版本，與相應的模型相匹配。

3.1 單指數平滑

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

假定時序中的觀測值可表示如下：

$Y_{t} =level+irregular_{t}$

則在時間點 $T_{t+1}$ 的預測值可記為：

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

其中 $c_{i} =alpha left( 1-alpha right) ^{i}$ ， $i=0,1,2,...$ ，並且 $0leq alpha leq 1$ ，且 $c_{i}$ 的總和為1。

$alpha$ 參數控制權數的下降速度，越接近1，則近期觀測值的權重越大，反之，則歷史觀測值的權重越大。通常，其值的大小有計算機進行優化選擇。下面舉一個例子。

nhtemp時序中有康涅狄格州紐黑文市從1912年到1971年每一年的平均華氏溫度，時序的折線圖如下：

從圖中可以看出，不存在某種趨勢，也看不出季節性因素，因此可以擬合一個單指數模型。代碼如下：

#導入包n> library(forecast)n#擬合模型n> fit <- ets(nhtemp, model = "ANN")n#查看模型n> fitnETS(A,N,N) nnCall:n ets(y = nhtemp, model = "ANN") nn Smoothing parameters:n alpha = 0.182 nn Initial states:n l = 50.2759 nn sigma: 1.1263nn AIC AICc BIC n265.9298 266.3584 272.2129 n#一步向前預測n> forecast(fit,1)n Point Forecast Lo 80 Hi 80 Lo 95 Hi 95n1972 51.87045 50.42708 53.31382 49.66301 54.0779n#繪製圖形n> plot(forecast(fit, 1), xlab="Year",n+ ylab=expression(paste("Temperature (", degree*F,")",)),n+ main="New Haven Annual Mean Temperature")n得到準確性度量n> accuracy(fit)n ME RMSE MAE MPE MAPE MASEnTraining set 0.1460295 1.126268 0.8951331 0.2418693 1.748922 0.7512497n ACF1nTraining set -0.00653111n> n

從擬合後的模型信息看，α值比較小（α=0.18），說明模型同時考慮了較近和較遠的觀測值。

forecast()函數用於預測時序為了的k步，格式為：

forecast(fit, k)n

本數據集一步向前預測的結果為51.9°F，並且給出了80%及95%置信區間的數據。

forecast包還提供了accuracy()函數，主要度量值如下：

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

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

$Y_{t} =level+slopetimes t+irregular_{t}$

Holt指數平滑主要有兩個參數，包括：

平滑參數 $alpha$ （alpha）控制水平項的指數型下降
$beta$ （beta）控制斜率的指數型下降

以上兩個參數的取值範圍都是[0,1]，參數越大意味著越近的觀測值的權重越大。

Holt-Winters指數光滑可用來擬合有水平項、趨勢項以及季節項的時間序列。此時時序對象的觀測值可表示如下：

$Y_{t} =level+slopetimes t+s_{t}+ irregular_{t}$

Holt-Winters指數光滑有三個參數，包括：

平滑參數 $alpha$ （alpha）控制水平項的指數型下降
$beta$ （beta）控制斜率的指數型下降
光滑參數 $gamma$ （gamma）控制季節項的指數下降

gamma參數的取值範圍同樣是[0,1]，值越大，意味著越近的觀測值的季節效應權重越大。

在此我們用Holt-Winters指數光滑來預測AirPassengers時序中接下來的五個值。代碼和圖形如下：

#導入包n> library(forecast)n#三指數模型擬合n> fit <- ets(log(AirPassengers),model = "AAA")n> fitnETS(A,A,A) nnCall:n ets(y = log(AirPassengers), model = "AAA") nn#光滑參數n Smoothing parameters:n alpha = 0.6534 n beta = 1e-04 n gamma = 1e-04 nn Initial states:n l = 4.8022 n b = 0.01 n s=-0.1047 -0.2186 -0.0761 0.0636 0.2083 0.217n 0.1145 -0.011 -0.0111 0.0196 -0.1111 -0.0905nn sigma: 0.0359nn AIC AICc BIC n-208.3619 -203.5047 -157.8750 n#查看準確性度量n> accuracy(fit)n ME RMSE MAE MPE MAPEnTraining set -0.0006710596 0.03592072 0.02773886 -0.01250262 0.508256n MASE ACF1nTraining set 0.2291672 0.09431354n#未來值預測n> pred <- forecast(fit,5)n#預測值查看n> predn Point Forecast Lo 80 Hi 80 Lo 95 Hi 95nJan 1961 6.103667 6.057633 6.149701 6.033264 6.174070nFeb 1961 6.093102 6.038107 6.148096 6.008995 6.177208nMar 1961 6.233814 6.171126 6.296502 6.137940 6.329688nApr 1961 6.213130 6.143591 6.282668 6.106780 6.319480nMay 1961 6.223273 6.147500 6.299047 6.107388 6.339159n#圖形繪製n> plot(pred, main="Forecast for Air Travel",n+ ylab="Log(AirPassengers)", xlab="Time")n#用原始值預測n> pred$mean <- exp(pred$mean)n> pred$lower <- exp(pred$lower)n> pred$upper <- exp(pred$upper)n> p <- cbind(pred$mean,pred$lower,pred$upper)n> dimnames(p)[[2]] <- c("mean", "Lo 80", "Lo 95", "Hi 80", "Hi 95")n> pn mean Lo 80 Lo 95 Hi 80 Hi 95nJan 1961 447.4958 427.3626 417.0741 468.5774 480.1365nFeb 1961 442.7926 419.0991 407.0741 467.8256 481.6452nMar 1961 509.6958 478.7246 463.0988 542.6706 560.9814nApr 1961 499.2613 465.7230 448.8908 535.2148 555.2839nMay 1961 504.3514 467.5469 449.1637 544.0531 566.3198n> n

3.3 ets()函數和自動預測

ets()函數還有其他功能，比如，擬合有可乘項的指數模型，加入抑制因子以及進行自動預測。

通過過ets(AirPassengers, model="MAM") 函數或hw(AirPassengers,seasonal="multiplicative")函數對原始數據擬合可乘模型；
ets()函數也可以用來擬合抑制項，大部分問題在加入抑制項以後，往往更符合實際情況；
ets()函數自動選取對原始數據擬合優度最高的模型。

以下以Johnson & Johnson數據集給出自動選取最優模型的步驟。

#導入包n> library(forecast)n#自動擬合模型n> fit <- ets(JohnsonJohnson)n#模型查看n> fitnETS(M,A,M) nnCall:n ets(y = JohnsonJohnson) nn Smoothing parameters:n alpha = 0.1481 n beta = 0.0912 n gamma = 0.4908 nn Initial states:n l = 0.6146 n b = 0.005 n s=0.692 1.2644 0.9666 1.077nn sigma: 0.0889nn AIC AICc BIC n166.6964 169.1289 188.5738 n#圖形展示n> plot(forecast(fit), main="Johnson & Johnson Forecasts",n+ ylab="Quarterly Earnings (Dollars)", xlab="Time", flty=2)n> n

這裡軟體自動搜索了一系列模型，並在其中找到最小化擬合標準的模型。

4 ARIMA預測模型

ARIMA模型當中，預測值表示為由最近的真實值和最近的預測誤差組成的線性函數。

首先需要了解以下概念：

滯後階數（lag）
自相關（autocorrelation）
偏自相關（partial autocorrelation）
差分（differencing）
平穩性（stationarity）

首先，滯後階數（lag）即歐美向後追溯的觀測值的數量。0階滯後項（Lag 0）代表沒有移位的時序，一階滯後（Lag 1）代表時序向左移動一位，二階滯後（Lag 2）代表時序向左移動兩位，以此類推。Nile時序不同滯後階數的數據如下：

時序對象可以通過lag(ts,k)函數變成k階滯後。

自相關，顧名思義，就是用來表示時序中各個觀測值之間的相關性。我們將 $AC_{k}$ 記為一些列觀測值（ $Y_{t}$ ）和k時期之前的觀測值（ $Y_{t-k}$ ）之間的相關性。當k取值不同時相關性所構成的圖即自相關函數圖（AutoCorrelation Function plot，ACF圖）。通常，ACF圖可用於為ARIMA模型選擇合適的參數，並評估最終模型的擬合效果。stats程序包中的acf()函數或者forecast包中的Acf()函數可以生成ACF圖。在本書中，我們採用Acf()函數。

偏自相關即當序列 $Y_{t}$ 和 $Y_{t-k}$ 之間的所有值帶來的效應都被移除後，兩個序列間的相關性。對於不同的k也可以畫出偏自相關圖。我們採用Pacf()函數來畫PACF圖。

平穩是指序列的統計性質並不會隨著時間的推移而改變，比如 $Y_{t}$ 的均值和方差都是恆定的。且對於任意滯後階數k，序列的自相關性不改變。

通常，ARIM模型主要用於擬合具有平穩性或者可以被轉換為平穩序列的時間序列，在實際運用中，在擬合前對序列進行變換非常常見。比如前面提到過的對數變換和Box-Cox變換。

從平穩的定義可以看出，具有平穩性的時序肯定不含有趨勢項，對於有趨勢項的非平穩時序可以通過差分來轉換為平穩性序列。

差分就是將時序中的每一個觀測值 $Y_{t}$ 都替換為 $Y_{t-1} -Y_{t}$ 。一次差分可以移除序列中的線性趨勢，二次差分移除二次項趨勢，三次差分移除三次項趨勢。對序列進行兩次以上的差分不太常用。

在R中可以用diff()函數對序列進行差分，格式如下：

diff(ts, differences=d)n

其中d為對序列ts的差分次數，默認d=1。

此外，forecast包中的ndiffs()函數可以找到最有的d值，格式為ndiffs(ts)。

可以通過觀察時序圖直觀判斷一個時序是否平穩，也可以通過R中的tseries包中的adf.test()進行ADF（Augmented Dickey-Fuller）統計檢驗來驗證平穩性假定，如果結果顯著則序列滿足平穩性要求。對於不穩定的序列，如果方差不是常數，則需要對數據進行變換，如果數據存在趨勢項，則需要進行差分。

4.2 ARMA和ARIMA模型

ARMA中的AR表示自回歸（AutoRegressive，AR），MA表示移動平均（Moving Averages，MA）。

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

其中：

$Y_{t}$ 是時序中的任一觀測值
$mu$ 是序列的均值
$beta$ 是權重
$varepsilon _{t}$ 是隨機擾動

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

其中：

$varepsilon$ 是預測的殘差
$theta$ 是權重

將以上兩種模型混合既是ARMA(p,q)模型，綜合表達式為：

此時，序列中的每個觀測值用過去的p個觀測值和q個殘差的線性組合來表示。

ARIMA(p,d,q)模型則意味著時序還被差分了d次，以下結合實例說明建立ARIMA模型的主要步驟。

1. 驗證序列的穩定性

nile序列的各年間的方差是穩定的，因此無需對數據做變換，但是數據中可能存在某種趨勢。

#導入包n> library(forecast)n> library(tseries)nn 『tseries』 version: 0.10-40nn 『tseries』 is a package for time series analysis and computationaln finance.nn See 『library(help="tseries")』 for details.nn#繪製原數據n> plot(Nile)n#確定差分次數n> ndiffs(Nile)n[1] 1n#差分n> dNile <- diff(Nile)n#繪製差分後數據n> plot(dNile)n#進行adf測試n> adf.test(dNile)nntAugmented Dickey-Fuller Testnndata: dNilenDickey-Fuller = -6.5924, Lag order = 4, p-value = 0.01nalternative hypothesis: stationarynnWarning message:nIn adf.test(dNile) : p-value smaller than printed p-valuen> n

原始序列在差分一次後變得平穩，並且通過ADF檢驗。

2. 選擇模型

藉助ACF圖和PACF圖選擇備選模型：

Acf(dNile)nPacf(dNile)n

從前文可知d=1，還需要確定p和q，下表給出了結合ACF和PACF圖選擇p和q的方法。

從上圖可以看出，可以看到在滯後項為一階時有一個比較明顯的自相關，而當滯後階數逐漸增加時，偏相關逐漸減小至零。因此，我們可以考慮ARIMA(0,1,1)模型。

3. 擬合模型

採用arima()函數來擬合ARIMA模型，格式如下：

arima(ts,order=c(q,d,q))n

對Nile序列進行擬合：

#導入包n> library(forecast)n#擬合模型n> fit <- arima(Nile, order = c(0,1,1))n#模型展示n> fitnnCall:narima(x = Nile, order = c(0, 1, 1))nnCoefficients:n ma1n -0.7329ns.e. 0.1143nnsigma^2 estimated as 20600: log likelihood = -632.55, aic = 1269.09n#查看擬合準確性n> accuracy(fit)n ME RMSE MAE MPE MAPE MASE ACF1nTraining set -11.9358 142.8071 112.1752 -3.574702 12.93594 0.841824 0.1153593n> n

4. 模型評價

如果模型合適，模型的殘差應該滿足獨立正態分布。

#繪製Q-Q正態圖n> qqnorm(fit$residuals)n#擬合Q-Q正太線n> qqline(fit$residuals)n#Box.test檢測n> Box.test(fit$residuals, type = "Ljung-Box")nntBox-Ljung testnndata: fit$residualsnX-squared = 1.3711, df = 1, p-value = 0.2416nn> n

從qqnorm()和qqline()的圖形來看，本例的結果還不錯。不知道如何看Q-Q圖的可以複習一下第八章的內容。

從Box.test()函數的結果來看，模型的殘差不顯著，ARIMA模型能夠較好地擬合本數據。

5. 預測

在模型通過評價以後，可以用來做預測。

#預測未來三年的數據n> forecast(fit,3)n Point Forecast Lo 80 Hi 80 Lo 95 Hi 95n1971 798.3673 614.4307 982.3040 517.0605 1079.674n1972 798.3673 607.9845 988.7502 507.2019 1089.533n1973 798.3673 601.7495 994.9851 497.6663 1099.068n

圖中的藍點為預測值的點估計，淺灰色和深灰色區域分別代表80%和95%的置信區間。

4.3 ARIMA的自動預測

auto.arima()函數可以實現最優ARIMA模型的自動選取。代碼如下：

#導入包n> library(forecast)n#自動擬合模型n> fit <- auto.arima(sunspots)n#模型查看n> fitnSeries: sunspots nARIMA(2,1,2) nnCoefficients:n ar1 ar2 ma1 ma2n 1.3467 -0.3963 -1.7710 0.8103ns.e. 0.0303 0.0287 0.0205 0.0194nnsigma^2 estimated as 243.8: log likelihood=-11745.5nAIC=23500.99 AICc=23501.01 BIC=23530.71n#預測n> forecast(fit,3)n Point Forecast Lo 80 Hi 80 Lo 95 Hi 95nJan 1984 40.43784 20.42717 60.44850 9.834167 71.04150nFeb 1984 41.35311 18.26341 64.44281 6.040458 76.66576nMar 1984 39.79670 15.23663 64.35677 2.235319 77.35808n#查看準確性n> accuracy(fit)n ME RMSE MAE MPE MAPE MASE ACF1nTraining set -0.02672716 15.60055 11.02575 NaN Inf 0.4775401 -0.01055012n> n

函數選定ARIMA模型的參數為p=2、d=1和q=2，此時模型的AIC值最小。

5 延伸閱讀

6 小結

預測歷史悠久，從預測天氣到預測大選，預測在各個學科當中都是一個基礎性問題。在本章當中，我們研究了在R中如何生存時間序列、如何判斷時間序列是否存在趨勢或者季節性因素，並探討了兩種常用的預測手段：指數模型和ARIMA模型。

但是需要注意的是，這些方法都假定為了的條件與現在是相似的，而事實去並不是這樣。很多外部因素會影響序列中的趨勢和模式，預測的時間跨度越大，不確定性越大。