今天分享的是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)

> qq

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

也可以看出趨勢。

關於給出模型後的引數估計，我們將在下一篇博文中討論。