文章目錄

• 一、準備工作
• • 1、數據準備
  • 2、基本概念
• 二、數據處理
• • 1、模式識別
  • 2、參數估計
  • 3、診斷性檢驗
  • • 1 殘差序列
    • 2 Ljung-Box 檢驗
  • 4、預測

一、準備工作

1、數據準備

所使用的數據是TSA包中的co2數據，如果沒有這個包的話，可以先裝一下

install.packages("TSA") # 安裝包 TSA

會有讓你選鏡像的過程，隨便選就行了。下載好之后，導入并查看數據

library(TSA) data(co2) win.graph(width = 4.875,height = 3,pointsize = 8) plot(co2,ylab='CO2') #繪制原始數據

可以看到，原始數據明顯有一個向上的趨勢和一個周期趨勢。

2、基本概念

赤池信息準則（Akaike’s(1973) Information Criterion, AIC）是建立在熵的概念基礎上，可以權衡所估計模型的復雜度和此模型擬合數據的優良性。
$AIC=?2log(極大似然估計值)+2k$
其中，如果模型包含截距或常數項，那么k=p+q+1；否則k=p+q。AIC越小越好。

Ljung-Box檢驗即LB檢驗、隨機性檢驗，用來檢驗m階滯后范圍內序列的自相關性是否顯著，或序列是否為白噪聲（或者統計量服從自由度為m的卡方分布）。若是白噪聲數據，則該數據沒有價值提取，即不用繼續分析了。

二、數據處理

拿到一個序列之后，首先判斷它是不是平穩時間序列，如果是就進行模式識別；如果不是就扣除趨勢項將其變成一個平穩時間序列。接著做模式識別、參數估計、模型診斷和預測。

ps: 這是從老師課件上找的流程圖，個人感覺模式識別部分，不應包含參數d，因為含有d的一般是ARIMA(p,d,q)模型，它是非平穩模型，而上一步已將非平穩時間序列轉成平穩時間序列了，d應該是在上一步確認的。也有可能是，扣除趨勢項和模式識別的界限根本不可能分那么細，但是又要用流程圖表示出來，所以才這么寫的。

1、模式識別

一般來講，模式識別就是判別出ARIMA(p,d,q)中的各階數p,d,q。模式識別常用的方法有：acf, pacf, eacf

首先來看它的自相關函數

acf(as.vector(co2),lag.max = 36) #自相關函數

季節自相關關系十分顯著：在滯后12，24，36，……上表現出很強的相關性。

plot(diff(co2),ylab='1st Diff. of CO2',xlab='Year') #一次差分，消除整體上升趨勢

可以看到，經過一次差分后，序列中的整體上升趨勢已經消除。再來看其樣本自相關函數

acf(as.vector(diff(co2)),lag.max = 36) #一次差分后的自相關函數

一次差分后，序列中仍存在強烈的季節性；應用季節差分法應該可以得到更為簡約的模型。

plot(diff(diff(co2),lag=12),xlab='Year',ylab='1st & seasonal Diff.')

繪制其自相關函數

acf(as.vector(diff(diff(co2),lag=12)),lag.max=36,ci.type='ma')

可以看到，經過一次差分和季節差分后的時間序列已經消除了季節性的大部分影響。根據樣本自相關函數可以看到，除了在滯后1和12上具有自相關性外，經一次和季節差分后的序列幾乎不再具有自相關性，所建模型只需要在滯后1和12上具有自相關性即可。

綜上，考慮構建乘法季節 $ARIMA(0,1,1)\times (0,1,1{)}_{12}$模型。

2、參數估計

模型建立后，需要估計模型的參數。乘法季節ARIMA模型只是一般ARIMA模型的特例。

m1.co2=arima(co2,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12)) print(m1.co2) ------------------------------------------- Coefficients:ma1 sma1-0.5792 -0.8206 s.e. 0.0791 0.1137sigma^2 estimated as 0.5446: log likelihood = -139.54, aic = 283.08

上面第一行代碼便得到了參數的極大似然估計值，參數估值的標準差為0.5446，對數似然值為-139.54，AIC=283.08。模型的參數估值均為高度顯著，進而將對該模型加以檢驗。

ps:為什么能根據這些指標值說明參數估值高度顯著，標準是多少？

3、診斷性檢驗

1 殘差序列

首先觀察殘差的時間序列圖

plot(window(rstandard(m1.co2),start=c(1995,2)),ylab='Standardized Resi.',type='o'); abline(h=0)

除了序列中間存在某些異常行為外，殘差圖中并沒有表明模型有任何主要的不規則性。然后對殘差的樣本自相關函數進一步觀察

acf(as.vector(window(rstandard(m1.co2),start=c(1995,2))),lag.max=36)

統計上顯著的相關系數位于滯后22，其值僅為-0.194，相關性非常小，且滯后22上的依賴關系難以給出合理的解釋。除了滯后22處的邊緣顯著以外，該模型似乎已捕捉到了序列中依賴關系的本質。

注：在acf前打個print即可輸出滯后各階的自相關函數的值。

2 Ljung-Box 檢驗

下面進行Ljung-Box檢驗：

win.graph(width=3,height =3,pointsize = 8) hist(window(rstandard(m1.co2),start=c(1995,2)),xlab='Standardized Resi.',ylab='Frequency')

直方圖的形狀像鐘形，但并不標準。對模型進行Ljung-Box檢驗，給出自由度為22的x=25.59, p=0.27，表明該模型已捕獲時間序列中的依賴關系。

why?

R 語言進行 Ljung-Box 檢驗的函數如下：

Box.test(x, lag = 1, type = c("Box-Pierce", "Ljung-Box"), fitdf = 0)

• x: 一個時間序列，殘差檢驗時，一般是殘差
• lag: 基于自相關因子得出的lag值
• type: Ljung-Box 檢驗就設置為 Ljung-Box
• fitdf: 如果x是一系列殘差，則需要減去自由度。

調用函數后，我們關心的就是p值，如果p > 0.05，則說明是白噪聲序列，是純隨機性序列。否則數據不是白噪聲，具有研究價值。

示例如下：

x <- rnorm (100) Box.test(x, lag = 5) Box.test(x, lag = 10, type = "Ljung") a=Box.test(resid(m1.xpole),type="Ljung",lag=20,fitdf=11)

接著繪制分位數-分位數圖（qq圖）

win.graph(width=5,height =5,pointsize = 8) qqnorm(window(rstandard(m1.co2),start=c(1995,2))) abline(c(0,0),c(1,1),col='red')

QQ 圖的上尾部，再次出現了一個異常值。但是，Shapiro-Wilk正態性檢驗給出的統計量W=0.982，進而得到p=0.11，且在任何顯著水平上正態性都未被拒絕。

作為對模型的進一步檢驗，考慮用$ARIMA(0,1,2)\times (0,1,1{)}_{12}$模型進行過度擬合。

m2.co2=arima(co2,order=c(0,1,2),seasonal=list(order=c(0,1,1),period=12)) print(m1.co2) print(m2.co2) -------------------------------- arima(x = co2, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12)) Coefficients:ma1 sma1-0.5792 -0.8206 s.e. 0.0791 0.1137 sigma^2 estimated as 0.5446: log likelihood = -139.54, aic = 283.08 -------------------------------- arima(x = co2, order = c(0, 1, 2), seasonal = list(order = c(0, 1, 1), period = 12))Coefficients:ma1 ma2 sma1-0.5714 -0.0165 -0.8274 s.e. 0.0897 0.0948 0.1224sigma^2 estimated as 0.5427: log likelihood = -139.52, aic = 285.05

可以看到，${\theta }_1$和$\theta$的估計變化很小（考慮標準差的大小時）。新參數${\theta }_2$的估值在統計上與零無異。AIC 已經增加的情況下，sigma^2和對數似然值均無顯著變化。所以使用$ARIMA(0,1,2)\times (0,1,1{)}_{12}$模型是過度擬合，根據“奧卡姆剃刀原理”，如無必要，勿增實體。所以使用$ARIMA(0,1,1)\times (0,1,1{)}_{12}$模型即可。

4、預測

前置時間設為2年，進行預測2年并繪制預測值與預測極限。

win.graph(width = 4.875,height = 3,pointsize = 8) plot(m1.co2,n1=c(2003,1),nahead=24,xlab='Year',type='o',ylab='CO2 Levels')

前置時間設為1年，進行預測4年并繪制預測值與預測極限。

win.graph(width = 4.875,height = 3,pointsize = 8) plot(m1.co2,n1=c(2004,1),n.ahead=48,xlab='Year',type='b',ylab='CO2 Levels')

ps: 上面參數的含義，n1=c(2004,1)代表從2004年1月開始（實際數據到2005年12月結束），n.ahead=48代表預測48個值（一年12個值，所以是4年）

總結

以上是生活随笔為你收集整理的R 语言做时间序列分析的实例（模式识别、拟合、检验、预测）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。