1. 程式人生 > >季節性ARIMA模型【R語言】

季節性ARIMA模型【R語言】

季節性的ARIMA模型可以預測含有季節性,趨勢性的時間序列。他的形式如下 在這裡插入圖片描述

這裡m是每一季節的週期值。季節項與非季節項的模型非常相近。但是季節項中包含了季節週期性。例如對於ARIMA(1,1,1)(1,1,1)4模型能夠寫成: 在這裡插入圖片描述 ACF與PACF 對於AR與MA模型的季節項,我們將會在ACF和PACF的lags上看到差異。例如,ARIMA(0,0,0)(0,0,1)12模型,我們將會在ACF的lag12處看到一個spike(長釘,即非零顯著的形象描述),而在其他地方看不到spike。PACF將會在具有周期性的位置出現指數衰減,即lags 12,24,36…類似的,一個ARIMA(0,0,0)(1,0,0)12模型則會顯示出在ACF圖的週期性位置顯示出指數衰減,而在PACF圖中的lag12處看到一個spike。我以中國人民大學統計學院易丹輝教授在時間序列課堂上給出的資料book19為例進行分析,這裡用到的軟體是R,課上易丹輝教授用的是EViews6.0。 R程式: data=read.csv("~/Downloads/Book19.csv") #讀取本地csv檔案 data=ts(data[,2],frequency=12,start=c(1990,1)) #[,2]是隻讀取第2列資料,start是在資料上增加起始時間,frequency是增加週期,ts為轉換成ts格式。 library(forecast) #載入程式包forecast,如果提示之前沒有安裝,則用install.package(“forecast”)進行安裝,再載入forecast包。 tsdisplay(data,xlab=“year”,ylab=“Retail index”) #顯示自相關係數ACF,偏自相關係數PACF,與資料原始影象。

在這裡插入圖片描述

data_d1=diff(data,1) 從上圖可以看出自相關係數呈現拖尾性,即序列Yt與之前的時刻有強烈的相關性,而且相關項的範圍是比較大的,說明序列本身有趨勢,為了消除這種趨勢我們取每一項與前一項的差。 tsdisplay(data_d1,xlab=“year”,ylab=“Retail index”) 在這裡插入圖片描述 可以看出來,趨勢基本消除了。序列較平穩。而PACF在季節性週期lags12處有spike,因此進一步做季節差分。 data_d1D1=diff(data_diff1,12) tsdisplay(data_d1D1,xlab=“year”,ylab=“Retail index”) 在這裡插入圖片描述

所以對於非季節項,我們只做了一階非季節差分,故d=1,從PACF看出p=2或3,從ACF看出q=1 對於季節項,我們也只做了一階季節差分,故D=1,再看在lag12,24處看是否有spike,從PACF看出P=1,從ACF看出Q=1 從而選擇模型一共有:Arima(2,1,1)(1,1,1)[12]和Arima(3,1,1)(1,1,1)[12]

fit1 <- Arima(data, order=c(2,1,1),seasonal=c(1,1,1)) fit1 #顯示AICc=1110.09 tsdisplay(residuals(fit1)) #顯示殘差 在這裡插入圖片描述

fit2 <- Arima(data, order=c(3,1,1),seasonal=c(1,1,1)) fit2 #顯示 AICc=1112.3 tsdisplay(residuals(fit1)) #顯示殘差 在這裡插入圖片描述 plot(forecast(fit1,h=12)) #選擇AICc值小的,輸出我們最終的預測結果~ 在這裡插入圖片描述 還有一個更加傻瓜的做法就是使用auto.arima函式自動定p,d,q,P,D,Q等引數。 arima1<-auto.arima(data,trace=T) 如下圖: 在這裡插入圖片描述

with drift代表有趨勢,所以最終模型可以加上d=1去除趨勢。 最終的模型可以是: Arima(2,1,0)(1,1,0)[12] 其AICc=1107.86

plot(forecast(Arima(data,order=c(2,1,0),seasonal=c(1,1,0)),h=12)) 在這裡插入圖片描述

這個結果是各個模型中AICc最小的。