R語言與時間序列學習筆記(1)
今天分享的是R語言中時間序列的有關內容。主要有:時間序列的建立,ARMA模型的建立與自相關和偏自相關函式。
一、 時間序列的建立
時間序列的建立函式為:ts().函式的引數列表如下:
ts(data = NA, start = 1, end = numeric(),frequency = 1,
deltat = 1, ts.eps = getOption("ts.eps"), class = , names = )
引數說明:data:這個必須是一個矩陣,或者向量,再或者資料框frame
Frequency:這個是時間觀測頻率數,也就是每個時間單位的資料數目
Start:時間序列開始值,允許第一個個時間單位出現數據缺失
舉例:ts(matrix(c(NA,NA,NA,1:31,NA),byrow=T,5,7),frequency=7,names=c("Sun","Mon ","Tue", "Wen" ,"Thu","Fri"," Sat"))
執行上面的程式碼就可以得到一個日曆:
Sun Mon Tue Wen Thu Fri Sat
NA NA NA 1 2 3 4
5 6 7 8 9 10 11
12 13 14 15 16 17 18
19 20 21 22 23 24 25
26 27 28 29 30 31 NA
在R語言中本身也有不少資料集,比如統計包中的sunspots,你可以通過函式data(sunspots)來呼叫它們。
二、 一些時間序列模型
這裡主要介紹AR,MA,隨機遊走,餘弦曲線趨勢,季節趨勢等
首先介紹一下AR模型:AR模型,即自迴歸(AutoRegressive,AR)模型,數學表示式為: AR :y(t)=a1y(t-1)+...any(t-n)+e(t)
其中,e(t)為均值為0,方差為某值的白噪聲訊號。
那麼產生AR模型的資料,我們就有兩種方法:1、呼叫R中的函式filter(線性濾波器)去產生AR模型;2、根據AR模型的定義自己編寫函式
先說第一種方法:呼叫R中的函式filter(線性濾波器)去產生AR模型
介紹函式filter的用法如下:
filter(x, filter, method = c("convolution", "recursive"),
sides = 2, circular = FALSE, init)
對於AR(2)模型x(t)=x(t-1)--0.9x(t-2)+e(t)
w<-rnorm(550)#我們假定白噪聲的分佈是正態的。
x<-filter(w,filter=c(1,-0.9),"recursive")
#方法:無論是“卷積”或“遞迴”(可以縮寫)。如果使用移動平均選擇“卷積”:如果“遞迴”便是選擇了自迴歸。
再說第二種方法:依據定義自己程式設計產生AR模型,還是以AR(2)模型x(t)=x(t-1)--0.9x(t-2) +e(t)為例,可編寫函式如下:
w<-rnorm(550)
AR<-function(w){
x<-w
x[2]=x[1]+w[1]
for(i in 3:550)
x[i]=x[i-1]-0.9*x[i-2]+w[i]
x
}
呼叫AR(W)即可得到。如果對相同的隨機數,我們可以發現兩個產生的時間序列是一致的。當然對於第二種方法產生的序列需要轉換為時間序列格式,用as.ts()處理。
類似的,我們給出MA,隨機遊走的模擬:
MA模型:
w<-rnorm(500)
v<-filter(w,sides=2,rep(1,3)/3)
隨機遊走:
w<-rnorm(200)
x<-cumsum(w)#累計求和,seeexample:cumsum(1:!0)
wd<-w+0.2
xd<-cumsum(wd)
可以做出相應的圖形:
再說一下季節性模型:
最簡單的季節模型就是一個分段的周期函式。比如說某地區一年的氣溫就是一個季節性模型。利用TSA包裡給出的資料tempdub我們可以發現他就是這樣的模型
給出驗證:
library(TSA)
data(tempdub)
month<-season(tempdub)
model1<-lm(tempdub~month)
summary(model1)
根據R輸出的結果:
Call:
lm(formula = tempdub ~month)
Residuals:
Min 1Q Median 3Q Max
-8.2750 -2.2479 0.1125 1.8896 9.8250
Coefficients:
Estimate Std. Error t valuePr(>|t|)
(Intercept) 16.608 0.987 16.828 < 2e-16 ***
monthFebruary 4.042 1.396 2.896 0.00443 **
monthMarch 15.867 1.396 11.368 < 2e-16 ***
monthApril 29.917 1.396 21.434 < 2e-16 ***
monthMay 41.483 1.396 29.721 < 2e-16 ***
monthJune 50.892 1.396 36.461 < 2e-16 ***
monthJuly 55.108 1.396 39.482 < 2e-16 ***
monthAugust 52.725 1.396 37.775 < 2e-16 ***
monthSeptember 44.417 1.396 31.822 < 2e-16 ***
monthOctober 34.367 1.396 24.622 < 2e-16 ***
monthNovember 20.042 1.396 14.359 < 2e-16 ***
monthDecember 7.033 1.396 5.039 1.51e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’1
Residual standarderror: 3.419 on 132 degrees of freedom
Multiple R-squared:0.9712, Adjusted R-squared: 0.9688
F-statistic: 405.1 on11 and 132 DF, p-value: < 2.2e-16
這裡2月份係數表明了一月份平均氣溫與二月份平均氣溫的差異,以此類推。
在介紹一下一個季節模型:餘弦趨勢μ1=βcos(2pi*f*t+φ)
還是考慮上面氣溫的例子:
驗證:
har<-harmonic(tempdub,1)
model2<-lm(tempdub~har)
summary(model2)
看看結果:
Call:
lm(formula = tempdub ~har)
Residuals:
Min 1Q Median 3Q Max
-11.1580 -2.2756 -0.1457 2.3754 11.2671
Coefficients:
Estimate Std. Error t valuePr(>|t|)
(Intercept) 46.2660 0.3088 149.816 < 2e-16 ***
harcos(2*pi*t)-26.7079 0.4367 -61.154 < 2e-16 ***
harsin(2*pi*t) -2.1697 0.4367 -4.968 1.93e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’1
Residual standarderror: 3.706 on 141 degrees of freedom
Multiple R-squared:0.9639, Adjusted R-squared: 0.9634
F-statistic: 1882 on 2 and 141 DF, p-value: < 2.2e-16
我們可以作圖來看擬合效果:
順便指出季節模型也可以模擬:比如μ1=βcos(2pi*f*t+φ)模型可以模擬如下:
t<-1:500
w<-rnorm(500)
c<-2*cos(2*pi*t/50+0.6*pi+w)
三、 自相關與偏自相關
我們可以根據定義給出自相關係數(ACF)的演算法:
例如資料:
> x<-1:10
>u<-mean(x)
>v<-var(x)
>sum((x[1:9]-u)*(x[2:10]-u))/(9*v) #延遲1
[1] 0.7
>sum((x[1:8]-u)*(x[3:10]-u))/(9*v) #延遲2
[1]0.4121212
>sum((x[1:7]-u)*(x[4:10]-u))/(9*v) #延遲3
[1]0.1484848
在R中也提供了直接計算acf的函式acf(),利用該函式也計算1至3階的acf,結果如下:
>a<-acf(x,3)
> a
Autocorrelationsof series ‘x’, by lag
0 1 2 3
1.0000.700 0.412 0.148
可以看出,是一樣的。
利用acf()可以處理很多階的acf,以太陽黑子數的資料集做例子:
>data(sunspots)
>acf(sunspots) #給出了相應的圖形
>a<-acf(sunspots,6) #為下面做估計做鋪墊,列出前6階的acf
> a
Autocorrelationsof series ‘sunspots’, by lag
0.00000.0833 0.1667 0.2500 0.3333 0.4167 0.5000
1.000 0.922 0.890 0.875 0.864 0.850 0.836
偏自相關:
對於一個平穩AR(p)模型,求出滯後k自相關係數p(k)時,實際上得到並不是x(t)與x(t-k)之間單純的相關關係。因為x(t)同時還會受到中間k-1個隨機變數x(t-1)、x(t-2)、……、x(t-k+1)的影響,而這k-1個隨機變數又都和x(t-k)具有相關關係,所以自相關係數p(k)裡實際摻雜了其他變數對x(t)與x(t-k)的影響。
為了能單純測度x(t-k)對x(t)的影響,引進偏自相關係數的概念。
對於平穩時間序列{x(t)},用數學語言描述就是:
p[(x(t),x(t-k)]|(x(t-1),……,x(t-k+1)={E[(x(t)-Ex(t)][x(t-k)-Ex(t-k)]}/E{[x(t-k)-Ex(t-k)]^2}
這就是滯後k偏自相關係數的定義。
總之,偏自相關就是在試圖解釋在剔除了中間k-1個隨機變數x(t-1)、x(t-2)、……、x(t-k+1)的干擾之後,x(t-k)對x(t)影響的相關程度。
在R語言中,使用函式PACF()可求解
還是使用太陽黑子數的例子:
> b<-pacf(sunspots,6)
> b
Partial autocorrelations of series ‘sunspots’, bylag
0.0833 0.1667 0.2500 0.3333 0.4167 0.5000
0.922 0.272 0.189 0.135 0.064 0.044
最後,我們利用這兩個函式來看看AR(p),MA(q)的自相關函式與偏自相關函式的截尾性與拖尾性。
利用二中所介紹的方法生成AR(2),MA(2)的資料。
AR(2)模型:
w<-rnorm(550)#我們假定白噪聲的分佈是正態的。
x<-filter(w,filter=c(1,-0.9),"recursive")
MA(3)模型:
w<-rnorm(500)
v<-filter(w,sides=2,rep(1,3)/3)
> qq<-pacf(x,5)
Partial autocorrelations of series ‘x’, by lag
1 2 3 4 5
0.532-0.861 -0.082 0.000
可以看出AR(2)模型的偏自相關函式是截尾的(但由於這個是資料,所以出現pacf只能看出趨勢,而不是在2步後直接變為0)
對於MA(3)模型的自相關函式,由於v的第一項與最後一項缺失,不妨擷取v的一部分資料,命名為a,有:
> y<-acf(a,5)
> y
Autocorrelations of series ‘a’, by lag
0 1 2 3 4 5
1.000 0.652 0.397 0.059 0.067 0.035
也可以看出趨勢。
關於給出模型後的引數估計,我們將在下一篇博文中討論。