時間序列分析

 

 

名稱

說明

ARIMA模型

  • 一般的ARIMA($p,d,q$)模型如下:

    \begin{displaymath}(1-\phi_1B-\cdots-\phi_pB^p)(1-B)^dZ_t=
\theta_0+(1-\theta_1B...
...-\theta_qB^q)a_t\mbox{\raisebox{-1.2mm}{\scriptsize {$\circ$}}}\end{displaymath}
  • 在模型中, 

            $\phi_i, i=1,\cdots,p$為AR的參數, 位階(order)為$p$

            $\theta_j, j=1,\cdots,q$ 為MA的參數, 位階為$q$

            $d$決定模型須經過幾次的差分動作。

            $Z_t$則為在時間$t$得到的觀察值。

  $a_t$稱為白噪音(white noise), 一般取$a_t$滿足$E(a_t)=0$, $Var(a_t)$為一常數。

模型的建立

  • 在給定一時間序列資料, 我們通常使用下列步驟來決定一可能的模型。

先繪製時間序列圖, 並且選擇適當的轉換方式。對任何的時間序列分析, 首先採取的步驟皆是先畫出資料的圖形。透過圖形的檢驗, 可提供關於資料的一些資訊, 像是趨勢、週期、離群值或者變異等。對時間序列資料有時需要經過變異數平穩的轉換或者差分。常用的轉換如對資料取log或開平方。

計算並檢驗原始資料的ACF(autocorrelation function)及 PACF (partial autocorrelation function)圖形由資料的ACF及PACF圖形通常可推估$p$, $q$$d$值。

          $d>0$, 檢定所決定的趨勢項$\theta_0$

 

要利用R來分析時間序列資料時, 需先下載名為“ ts”的資料庫。


加拿大山貓數量資料分析

 

 

名稱

說明

加拿大山貓

(lynx)

  • 在R的資料庫中加拿大山貓(lynx)數量的數據資料名為“ lynx”。

  • 此筆資料記錄加拿大西北部馬更些(Mackenzie)流域附近每年的山貓數量, 資料所收集的時間是西元1821至1934年。

 

  • ACF:acf(x, lag.max = NULL, type = c("correlation", "covariance", "partial"), plot = TRUE)

  • PACF:pacf(x, lag.max, plot, na.action, ...)

  • AR模型:ar(x, aic = TRUE, order.max = NULL, method=c("yule-walker", "burg", "ols", "mle", "yw"))

            參數說明

  •  x:代表數據名稱。

  • lag.max:設定圖形的最大lag數。

  • type:設定acf圖形計算的型式。內定值為“correlation

  • aic:是否利用AIC值來設定AR模型位階的選擇。內定值為“TRUE”,若改為“FALSE”,則模型的選擇依照設定的最大位階,及指令中的“order.max”來計算AR模型的參數估計值。

  • order.max:模型配適的最大位階數。內定值為10*log10(N),N代表觀測值的數目。

  • method:利用何種方法來配適模型。內定值為"yule-walker"

利用“ plot”指令先畫出加拿大山貓數量的時間序列圖, 如圖1。

       圖 2為加拿大山貓數量的ACF及PACF圖。

  

圖 1 加拿大山貓數量的時間序列圖

  

     

圖 10.2 加拿大山貓數量的ACF及PACF圖

 

      範例程式碼:  

> data(lynx)
> plot(lynx,main="The yearly number of lynx")
> acf(lynx)
> pacf(lynx)

由圖 3的常態Q-Q圖來看, 資料並不是很符合常態分佈, 因此對原始資料作對數轉換。轉換後的資料由圖 4的常態Q-Q圖可明顯看出較原始資料更符合常態分佈。對經過對數轉換後的數據, 利用“plot”、“ acf”及“ pacf”指令分別畫出時間序列圖、ACF及PACF圖, 如圖 5及圖 6。

圖 3 加拿大山貓數量的常態Q-Q圖

圖 4 經對數轉換後, 加拿大山貓數量的常態Q-Q圖

圖 5 經對數轉換後, 加拿大山貓數量的時間序列圖

  

     

圖  6 經對數轉換後, 加拿大山貓數量的ACF及PACF圖

 

      範例程式碼:  

> qqnorm(lynx,pch=20,main="Normal Q-Q plot of lynx")
> qqline(lynx)
> tlynx<-log(lynx)
> qqnorm(tlynx,pch=20,main="Normal Q-Q plot of log(lynx)")
> qqline(tlynx)
> plot(tlynx,main="The natural logarithms of yearly number 

+ of lynx")
> acf(tlynx)

> pacf(tlynx)

 

發現經過轉換後之數據仍不平穩, 由圖 5的序列圖及圖10.6的ACF圖形可粗略 看出其具有週期性, 且週期大約在9至10之間。因為ACF呈現指數下降, PACF在 lag=2或11之後大致都落在水平線之間。經由這些圖形的判斷, 我們決定配適一個AR($p$)的模型。利用R中的指令`` ar'' 可得圖10.7的結果, R自動配適出AR(11)的模型, 並給出每一係數值。

圖  7 經對數轉換後, 加拿大山貓數量的模型配適結果

 

      範例程式碼:  

> mlynx<-ar(tlynx)
> mlynx

 

 


太陽黑子數據資料分析

 

 

名稱

說明

太陽黑子

(sunspots)

  • 在R的資料庫中太陽黑子的數據資料名“sunspots”。

  • 此筆資料記錄從西元1749到1983年間太陽黑子之月及年平均數, 我們僅針對太陽黑子的年平均數作討論, 並設此筆資料為$Z_t$$t$表時間。

 

  • 差分:diff(x, lag=1, differences=1, ...)

  • ARIMA模型:
    arima(x, order = c(0, 0, 0), seasonal = list(order = c(0, 0, 0), 
    period = NA))

    參數說明

    •  x:代表數據名稱。

    • order :分別代表ARIMA(p,d,q)之位階值。

    • seasonal:代表模型是否有週期性,是否須經過差分的轉換。

     

先畫出其時間序列圖, 如圖 8。利用R的指令“ acf”及“ pacf”可分別畫出太陽黑子的ACF及PACF圖形, 如圖 9。

圖 8 太陽黑子的時間序列圖


圖 9 太陽黑子的ACF及PACF圖
 

      範例程式碼:  

> data(sunspot)
> plot(sunspot.year,pch=20)

> acf(sunspot.year)
> pacf(sunspot.year) 
     

 

由圖 8的時間序列圖觀察到變異數可能不為一 常數值, 所以考慮將此筆資料做開方 轉換。分別畫原始資料$Z_t$$\sqrt{Z_t}$的常態Q-Q圖, 如圖 10。 經轉換後之太陽黑子資料的時間序列圖如圖 11。對轉換後之數據畫出ACF及PACF圖, 參考圖 12。


圖 10 太陽黑子原始資料$Z_t$$\sqrt{Z_t}$的常態Q-Q圖

圖 11 經開方轉換後, 太陽黑子時間序列圖

圖 12 經開方轉換後, 太陽黑子的ACF及PACF圖
 

      範例程式碼:  

> sqsunspot<-sqrt(sunspot.year)

> par(mfrow=c(1,2))
> qqnorm(sunspot.year,pch=20,main=expression(paste("

+ The normal Q-Q plot",Z[t])))
> qqline(sunspot.year)
> qqnorm(sqsunspot,pch=20,main=expression(paste("

+ The normal Q-Q plot",sqrt(Z[t]))))
> qqline(sqsunspot)

> par(mfrow=c(1,1))
> plot(sqsunspot,main=expression(paste("The time plot of ",

+  sqrt(Z[t]))))
> par(mfrow=c(1,2))
> acf(sqsunspot)
> pacf(sqsunspot)

 



許多學者已試過AR(2)模型和AR(9)模型, 我們同樣利用此二模型來配適轉換後之數據。利用“ar”指令可分別得到AR(2)和AR(9)參數估計值, 如圖 13及圖 14。


圖 13 經開方轉換後, 太陽黑子AR(2)模型配適結果



 
圖 14 經開方轉換後, 太陽黑子AR(9)模型配適結果
 

      範例程式碼:  

> ar(sqsunspot,order=c(2))

> ar(sqsunspot,order=c(9))

 

由圖 12之ACF可約略看出資料具有週期性, 在一些文獻探討中, 提出週期可能為11。故先將資料去除週期性, 利用“diff”指令對資料作差分, 差分後可畫出時間序列圖、ACF及PACF圖, 如圖 15及圖 16。 最後考慮ARIMA(2,0,0)$\times$(1,1,0)$_{11}$模型, 利用“ arima”指令來配適太陽黑子資料 , 結果如圖 17。

圖 15 經差分後, 太陽黑子時間序列圖


圖 16 經差分後, 太陽黑子的ACF及PACF圖

 圖 17 太陽黑子的模型配適結果

 

      範例程式碼:  

> dsunspot<-diff(sqsunspot,lag=11)

> plot(dsunspot,pch=20,main=expression(paste("

+ The time plot of ",(1-B^11)*sqrt(Z[t]))))

> par(mfrow=c(1,2))
> acf(dsunspot)
> pacf(dsunspot)

> arima(sqsunspot,order=c(2,0,0),

+ seasonal=list(order=c(1,1,0),period=11))