時(shí)間序列與R教材

讀入數(shù)據(jù)ts(data,frequency=n,start=x,end=y,names=c(a,b,c,...))n：單位時(shí)間的觀測值數(shù)目，如n=1為年，n=7為周，n=12為月，n=4為季度等x：從第幾個(gè)開始，若為兩個(gè)整數(shù)的向量，如(1986,2)則從1986年的第2個(gè)月開始插入數(shù)據(jù)例：>birthsts<-ts(births,frequency=12,start=c(1946,1))

>ts(births,frequency=12,start=c(1946,1))時(shí)間序列檢驗(yàn)1、自相關(guān)自相關(guān)：

自相關(guān)系數(shù)通過acf()實(shí)現(xiàn)acf(x,lag.max=n,type=c("correlation","covariance","partial"),plot=TRUE)--x

時(shí)間序列--lag.max

最大之后階數(shù)--type

分別為“自相關(guān)系數(shù)”，“自協(xié)方差”，“偏自相關(guān)系數(shù)”--plot

是否要畫出acf圖形，F(xiàn)ALSE為不畫2、偏自相關(guān)偏自相關(guān)：偏自相關(guān)系數(shù)通過pacf()或者acf()中type="partial"需注意的一點(diǎn)是，acf從0開始，pacf從1階開始。3、平穩(wěn)性檢驗(yàn)平穩(wěn)性：序列的均值為常數(shù)，且自協(xié)方差函數(shù)只與時(shí)滯有關(guān)，與時(shí)間點(diǎn)無關(guān)。平穩(wěn)性的目的是使序列不存在周期性和趨勢性。平穩(wěn)化最好的方法是差分，檢驗(yàn)平穩(wěn)化與否最有效的方法是單位根檢驗(yàn)(AugmentedDickey-Fuller:ADF)。差分通過diff()實(shí)現(xiàn)diff(data,differences=d,lag=n)--data

時(shí)間序列--differences

d階差分--lag

周期內(nèi)差分，如12為為今年1月與去年1月的差分(1)時(shí)序圖檢驗(yàn)

平穩(wěn)時(shí)間序列的均值和方差為常數(shù)，可用plot()檢查序列，應(yīng)在一個(gè)常數(shù)值附近隨機(jī)波動(dòng)，且波動(dòng)范圍有界。(2)自相關(guān)圖檢驗(yàn)

平穩(wěn)時(shí)間序列有短期的相關(guān)性，因此，自相關(guān)系數(shù)會(huì)很快遞減為0。(3)單位根檢驗(yàn)fUnitRoots包以下討論幾種時(shí)間序列的平穩(wěn)性(1)確定趨勢的時(shí)間序列

帶隨機(jī)波動(dòng)的線性趨勢序列，對y=0.8x+3加入隨機(jī)擾動(dòng)后形成序列yy>x<-seq(0.1,15,0.02)>y<-0.8*x+3>yy<-ts(jitter(y,amount=0.5))>plot(yy)>acf(yy,lag.max=300)

>pacf(yy,lag.max=300)

對于一個(gè)單增序列來說，acf是非常緩慢下降的序列。(2)周期性時(shí)間序列>x<-seq(0.1,15,0.02)>y<-sin(x)>yy<-ts(jitter(y,amount=0.2))

>plot(yy)>acf(yy,lag.max=300)>pacf(yy,lag.max=300)

acf有偽周期的特點(diǎn)(3)偽周期時(shí)間序列>x<-seq(0.1,20,0.02)

>y<-10*sin(x)/x

>yy<-ts(jitter(y,amount=0.2))

>plot(yy)>acf(yy,lag.max=300)>pacf(yy,lag.max=300)

acf同樣有偽周期的特點(diǎn)(4)隨機(jī)游走Y=ts(rnorm(1000,mean=0,sd=1));

for(iin2:length(Y)){

Y[i]=Y[i]+Y[i-1];

}

plot(Y,

main="隨機(jī)游走",type="b",col="red");

abline(h=0)acf(Y,lag.max=300)pacf(Y,lag.max=300)4、純隨機(jī)檢驗(yàn)（白噪音檢驗(yàn)）

白噪音：時(shí)間序列的每一個(gè)元素都為獨(dú)立同分布，且均值為0。

barlett檢驗(yàn)，統(tǒng)計(jì)量有Q統(tǒng)計(jì)量和LB統(tǒng)計(jì)量，兩者都漸進(jìn)服從卡方分布。H0：該序列為白噪音。p-value>p0則接受原假設(shè)。

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

--x

待檢驗(yàn)序列

--lag

滯后階數(shù)

--type

type="Box-Pierce"為Q統(tǒng)計(jì)量，type="Ljung-Box"為LB統(tǒng)計(jì)量

--fitdf

自由度調(diào)整

值得注意的是，對于滯后階數(shù)的確定，可以寫一循環(huán)遍歷n，看哪個(gè)階數(shù)不通過假設(shè)檢驗(yàn)

以下以白噪音為例，進(jìn)行檢驗(yàn)>D=ts(rnorm(300,0,1))

>par(mfrow=c(3,1))

>plot(D,main="白噪音")

>acf(D,main="acf")

>pacf(D,main="pacf")>LAG=c();LB=c();p=c()

>for(iin1:12){

+Btest=Box.test(D,type="Ljung-Box",lag=i)

+LAG[i]=i+LB[i]=Btest$statistic

+p[i]=Btest$p.value

+result=cbind(LAG,LB,p)

+}

>result

LAG

LB

p

[1,]

10.28839210.5912528

[2,]

20.29009140.8649828

[3,]

31.38220960.7097101

[4,]

41.38341490.8470724

[5,]

51.66484770.8932983

[6,]

63.77557900.7070157

[7,]

74.35137630.7385322

[8,]

84.53751680.8056680

[9,]

95.99163930.7407548

[10,]

106.05279810.8108080

[11,]

117.04457820.7954621

[12,]

127.60513830.8151766可見，各階LB檢驗(yàn)都不能拒絕原假設(shè)，因此可認(rèn)為序列純隨機(jī)。時(shí)間序列的分解

--通過看時(shí)序圖觀察是否其趨勢性和周期性，若季節(jié)性和隨機(jī)部分在整個(gè)時(shí)間周期內(nèi)大致不變，即可采用相加模型1、非季節(jié)性時(shí)間序列分解：非季節(jié)性（趨勢性）+隨機(jī)部分

趨勢性一般用平滑法來實(shí)現(xiàn)，目標(biāo)是消除隨機(jī)性，將曲線轉(zhuǎn)化為平滑的趨勢線

最簡單的為移動(dòng)平均法（包括簡單移動(dòng)平均，加權(quán)移動(dòng)平均，多次移動(dòng)平均），

其中簡單移動(dòng)平均法為：

對x=(x1,x2,x3,x4,x5,...),若計(jì)算跨度n=5，則x6=(x1+...+x5)/5,x7=(x2+...+x6)/5...

若隨機(jī)性較大，則取較大的n，以避免隨機(jī)造成的偏差；若隨機(jī)性較小，則取較小的n，方便跟蹤數(shù)據(jù)趨勢

在R中用TTR程序包的SMA函數(shù)，SMA(data,n=8)king<-scan("/tsdldata/misc/kings.dat",skip=3)kingts<-ts(king)plot(kingts)install.packages("g:/R_cran/Rpkg/TTR_0.22-0.zip")library(TTR)kingstime<-SMA(kingts,n=8)

//周期為8的簡單移動(dòng)平均plot(kingstime)2、季節(jié)性時(shí)間序列分解：非季節(jié)性（趨勢性）+季節(jié)性（周期性）+隨機(jī)部分采用R語言中的decompose()函數(shù)，生成trend+seasonal+random三部分births<-scan("/tsdldata/data/nybirths.dat")birthsts<-ts(births,

frequency=12,start=c(1946,1))plot(birthsts)birthsts_decompose<-decompose(birthsts)plot(birthsts_decompose)時(shí)間序列的預(yù)測1、簡單指數(shù)平滑法

處于恒定水平和無季節(jié)性變動(dòng)的時(shí)間序列，可使用簡單指數(shù)平滑法進(jìn)行短期預(yù)測。其中，恒定水平指的是在某一數(shù)值周圍波動(dòng)。

指數(shù)平滑模型為：F(t)=α×x(t)+(1-α)×F(t-1),

α>0，即用上期的值與上期的預(yù)測值來預(yù)測本期值。

更多時(shí)間序列平滑預(yù)測模型介紹，見《時(shí)間序列平滑預(yù)測法.ppt》例子：1813~1912年倫敦每年降雨量rain<-scan("/tsdldata/hurst/precip1.dat",skip=1)raints<-ts(rain,start=c(1813))plot(raints)--數(shù)據(jù)在25附近隨機(jī)波動(dòng)，且無季節(jié)性和趨勢性，采用簡單指數(shù)平滑法raints_exp<-

HoltWinters(raints,beta=FALSE,gamma=FALSE)Holt-Wintersexponentialsmoothingwithouttrendandwithoutseasonalcomponent.

Call:

HoltWinters(x=raints,beta=FALSE,gamma=FALSE)

Smoothingparameters:

alpha:0.02412151

beta:FALSE

gamma:FALSE

Coefficients:

[,1]

a24.67819--α值為0.024，權(quán)重更多依賴于上一期的預(yù)測值--系數(shù)a為整體的水平基準(zhǔn)值plot(raints_exp)--

raints_exp$fitted和

raints_exp$SSE分別為預(yù)測值（上圖的紅線）及誤差平方和raints_exp<-HoltWinters(raints,beta=FALSE,gamma=FALSE,

l.start=23)--可用l.start設(shè)置初始值水平--以上是對已知樣本的樣本內(nèi)預(yù)測，若想預(yù)測更遠(yuǎn)的值，可基于上邊輸出的結(jié)果，使用forecast.HoltWinters函數(shù)預(yù)測library("forecast")raints_exp_forecast<-forecast.HoltWinters(raints_exp,h=8,

level=c(95))

//h=8表示預(yù)測之后的8個(gè)值，level表示給出95%的預(yù)測區(qū)間（默認(rèn)給出80%和95%）Point

Forecast

Lo95

Hi95

1913

24.6781916.2616933.09470

1914

24.6781916.2592433.09715

1915

24.6781916.2567933.09960

1916

24.6781916.2543433.10204

1917

24.6781916.2519033.10449

1918

24.6781916.2494533.10694

1919

24.6781916.2470133.10938

1920

24.6781916.2445633.11182plot(raints_exp_forecast)--顏色條內(nèi)的為預(yù)測值區(qū)間>raints_exp_forecast$

raints_exp_forecast$method

raints_exp_forecast$model

raints_exp_forecast$level

raints_exp_forecast$mean(預(yù)測值)

raints_exp_forecast$lower

raints_exp_forecast$upper

raints_exp_forecast$x

raints_exp_forecast$xname

raints_exp_forecast$fitted

raints_exp_forecast$residuals--上邊為輸出結(jié)果包含的變量，其中raints_exp_forecast$residuals為樣本內(nèi)預(yù)測的誤差2、霍爾特指數(shù)平滑法holt

有增長或下降趨勢的，但無季節(jié)性因素的時(shí)間序列，可使用該方法短期預(yù)測。

Holt指數(shù)平滑法需要估計(jì)兩個(gè)參數(shù)：水平參數(shù)（水平值變化情況）和斜率參數(shù)（斜率變化情況），表現(xiàn)為α和β，若需預(yù)先設(shè)定，則表現(xiàn)為l.start和b.start例子：1866~1911年女人裙子直徑skirts<-scan("/tsdldata/roberts/skirts.dat",skip=5)skirtsts<-ts(skirts,start=c(1866))plot(skirtsts)--數(shù)據(jù)有上升和下降的趨勢，且無季節(jié)性，考慮采用Holt預(yù)測法skirtsts_holt<-HoltWinters(skirtsts,gamma=FALSE)Holt-Wintersexponentialsmoothingwithtrendandwithoutseasonalcomponent.

Call:

HoltWinters(x=skirtsts,gamma=FALSE)

Smoothingparameters:

alpha:0.8383481

beta:1

gamma:FALSE

Coefficients:

[,1]

a529.308585

b

5.690464plot(skirtsts_holt)library(forecast)skirtsts_holt_forecast<-forecast.HoltWinters(skirtsts_holt,h=10,

level=c(95))

Point

Forecast

Lo95

Hi95

1912

534.9990496.08130

573.9168

1913

540.6895464.71204

616.6670

1914

546.3800422.47258

670.2874

1915

552.0704372.11216

732.0287

1916

557.7609314.86713

800.6547

1917

563.4514251.49103

875.4117

1918

569.1418182.51596

955.7677

1919

574.8323108.348291041.3163

1920

580.5228

29.313621131.7319

1921

586.2132-54.318701226.7452plot(skirtsts_holt_forecast)3、霍爾特-溫特指數(shù)平滑法holt-wintersHolt-Winters模型的基本思想是把具體線性趨勢、季節(jié)變動(dòng)和隨機(jī)變動(dòng)的時(shí)間序列進(jìn)行分解研究，并與指數(shù)平滑法相結(jié)合，分別對長期趨勢(Ut)、趨勢的增量(bt)和季節(jié)變動(dòng)(Ft)作出估計(jì)，然后建立預(yù)測模型，外推預(yù)測值。/view/75e854886529647d272852e0.html

有增長或下降趨勢，且有季節(jié)性變動(dòng)趨勢的時(shí)間序列，可用該方法預(yù)測。

Holt-Winters需要估計(jì)三個(gè)參數(shù)：α、β和γ，分別表示水平變化參數(shù)，斜率變化參數(shù)和季節(jié)部分例子：1987年1月~1993年12月每月銷售情況sales<-scan("/tsdldata/data/fancy.dat")salests<-ts(sales,start=c(1987,1),frequency=12)plot(salests)salests_hw<-HoltWinters(salests)Holt-Wintersexponentialsmoothingwithtrendandadditiveseasonalcomponent.

Call:

HoltWinters(x=salests)

Smoothingparameters:

alpha:0.021882

beta:1

gamma:1

Coefficients:

[,1]

a

31552.6045

b

778.6445

s1

-12994.3112

s2

-13070.8952

s3

-3550.8387

s4

-8961.9894

s5

-11179.7153

s6

-9320.5816

s7

-2480.7650

s8

-712.8248

s9

639.4040

s10

505.2548

s11

15860.4200

s12

73108.0655plot(salests_hw)salests_hw_forecast<-forecast.HoltWinters(salests_hw,h=12,level=c(95))

Point

Forecast

Lo95

Hi95

Jan1994

19336.94

9416.573

29257.30

Feb1994

20039.00

10109.138

29968.86

Mar1994

30337.70

20386.507

40288.89

Apr1994

25705.19

15716.189

35694.20

May1994

24266.11

14218.313

34313.91

Jun1994

26903.89

16772.026

37035.75

Jul1994

34522.35

24277.173

44767.53

Aug1994

37068.94

26677.616

47460.26

Sep1994

39199.81

28626.424

49773.19

Oct1994

39844.30

29050.384

50638.23

Nov1994

55978.11

44923.226

67033.00

Dec1994

114004.40102646.754125362.05plot(salests_hw_forecast)--對原始數(shù)據(jù)取對數(shù)再擬合salests_log<-log(salests)

plot(salests_log)salests_log_hw<-HoltWinters(salests_log)Holt-Wintersexponentialsmoothingwithtrendandadditiveseasonalcomponent.

Call:

HoltWinters(x=salests_log)

Smoothingparameters:

alpha:0.413418

beta:0

gamma:0.9561275

Coefficients:

[,1]

a

10.37661961

b

0.02996319

s1

-0.80952063

s2

-0.60576477

s3

0.01103238

s4

-0.24160551

s5

-0.35933517

s6

-0.18076683

s7

0.07788605

s8

0.10147055

s9

0.09649353

s10

0.05197826

s11

0.41793637

s12

1.18088423salests_log_hw$SSE

[1]2.011491

plot(salests_log_hw)salests_log_hw_forecast<-forecast.HoltWinters(salests_log_hw,h=12,level=c(95))

Point

Forecast

Lo95

Hi95

Jan1994

9.597062

9.267409

9.926715

Feb1994

9.830781

9.47406810.187495

Mar1994

10.47754210.09568010.859403

Apr1994

10.254867

9.84941510.660319

May1994

10.167100

9.73935610.594845

Jun1994

10.375632

9.92670110.824562

Jul1994

10.66424810.19508711.133409

Aug1994

10.71779610.22924111.206350

Sep1994

10.74278210.23557511.249989

Oct1994

10.72823010.20303211.253427

Nov1994

11.12415110.58155911.666743

Dec1994

11.91706211.35761612.476508plot(salests_log_hw_forecast)自回歸移動(dòng)平均模型（ARIMA）指數(shù)平滑法對于預(yù)測來說是非常有幫助的，而且它對時(shí)間序列上面連續(xù)的值之間相關(guān)性沒有要求。在某種情況下，若要考慮數(shù)據(jù)之間的相關(guān)性來創(chuàng)建更好的預(yù)測模型，則會(huì)采用ARIMA模型。自回歸移動(dòng)平均模型（ARIMA）包含一個(gè)確定（explicit）的統(tǒng)計(jì)模型用于處理時(shí)間序列的不規(guī)則部分，它也允許不規(guī)則部分可以自相關(guān)。ARIMA模型的目標(biāo)在于找到合適的p,d,q值1、ARIMA必須是對平穩(wěn)序列做預(yù)測，因此，對于非平穩(wěn)序列，需要對其做差分直到平穩(wěn)。若d階差分才能得到平穩(wěn)序列，則為ARIMA(p,d,q)，其中d為差分的階數(shù)。2、若考慮對平穩(wěn)序列做預(yù)測，用ARMA模型，基本識別過程如下：獲取時(shí)序數(shù)據(jù)|生成模擬數(shù)據(jù)->模型參數(shù)估計(jì)->模型建立->模型結(jié)果評估->模型預(yù)測和應(yīng)用

1)生成模擬數(shù)據(jù)

arima.sim(model,