時間序列分解Decompose

時間序列分解算法和decompose函數(shù)實(shí)現(xiàn)李思亮目錄TOC\o"1-5"\h\z\o"CurrentDocument"時間序列分解算法和 decompose函數(shù)實(shí)現(xiàn) 1\o"CurrentDocument"1數(shù)據(jù)讀入并生成時間序列 2\o"CurrentDocument"2數(shù)據(jù)可視化 4\o"CurrentDocument"3 時間序列分解 7在時間序列分析的過程中，往往需要對時間序列作出初步分析，本文主要采用R語言作為分析平臺，從數(shù)據(jù)的讀入，可視化圖，分解(decompose)為趨勢項(xiàng)，季節(jié)項(xiàng)，隨機(jī)波動等角度對數(shù)據(jù)開展分析的幾個案例。最后對分解算法作出初步描述并探討其預(yù)測預(yù)報中的潛在應(yīng)用。本文的數(shù)據(jù)和部分內(nèi)容主要采用/en/latest中的內(nèi)容，感興趣的讀者可以參考。1數(shù)據(jù)讀入并生成時間序列對于數(shù)據(jù)分析來講，數(shù)據(jù)讀入是一個比較關(guān)鍵的步驟。常用的數(shù)據(jù)讀入函數(shù)有scan，read.table等。下面列舉了幾種常見的數(shù)據(jù)。首先是/tsdldata/misc/kings.dat，中包含了英國國王的壽命從William開始，數(shù)據(jù)來源(HipelandMcleod,1994)。kings<-scan("/tsdldata/misc/kings.dat”,skip=3)Read42itemskings[1]604367505642506568436534473449411335535616436959485986556851334967778167718168707756上述例子中，讀入了連續(xù)42個公國國王的壽命并將其賦給變量‘kings’如果我們希望對讀入數(shù)據(jù)開展分析，下一步就是將其轉(zhuǎn)化為時間序列對象(時間序列類)，R提供了很多函數(shù)用于分析時間序列類數(shù)據(jù)?？梢允褂胻s函數(shù)將變量轉(zhuǎn)化為時間序列類。kingsts<-ts(kings)kingstsTimeSeries:Start=1End=42Frequency=1[1]604367505642506568436534473449411335535616436959485986556851334967778167718168707756對于上述數(shù)據(jù)操作的好處是將數(shù)據(jù)轉(zhuǎn)化為特定的“時間序列類”便于我們使用R中的函數(shù)分析數(shù)據(jù)。有時候我們會按照一定的時間周期來收集數(shù)據(jù)，這個周期可能是季度，月，日，小時，分。在大數(shù)據(jù)時代，有些情況下的數(shù)據(jù)是按照秒來采集收集。這種情況下，我們需要對數(shù)據(jù)的周期或頻率進(jìn)行設(shè)置。這里采用ts函數(shù)中的frequency參數(shù)可以實(shí)現(xiàn)這種功能。比方說，若按1年為一個周期，我們的月度時間序列數(shù)據(jù)應(yīng)為frequency=12，若為季度時間序列數(shù)據(jù)，則可設(shè)置frequency=4。另外，還可以利用start參數(shù)，設(shè)置時間序列的起點(diǎn)，比如若我們一個周期共4個觀測，而第一個數(shù)據(jù)對應(yīng)為1986年的起的第2個觀測，則可使用start=c(1986,2)。有一個紐約市月出生數(shù)量的數(shù)據(jù)集，從1946年1月至1959年12月。數(shù)據(jù)可通過/tsdldata/data/nybirths.dat獲取。我們將其讀入至R中。births<-scan("/tsdldata/data/nybirths.dat")Read168itemsbirthsTS<-ts(births,frequency=12,start=c(1946,1))birthsTSJanFebMarAprMayJunJulAugSepOctNovDec26.66323.59826.93124.74025.80624.36424.47723.90123.17523.22721.67221.87021.43921.08923.70921.66921.75220.76123.47923.82423.10523.11021.75922.07321.93720.03523.59021.67222.22222.12323.95023.50422.23823.14221.05921.57321.54820.00022.42420.61521.76122.87424.10423.74823.26222.90721.51922.02522.60420.89424.67723.67325.32023.58324.67124.45424.12224.25222.08422.99123.28723.04925.07624.03724.43024.66726.45125.61825.01425.11022.96423.98123.79822.27024.77522.64623.98824.73726.27625.81625.21025.19923.16224.70724.36422.64425.56524.06225.43124.63527.00926.60626.26826.46225.24625.18024.65723.30426.98226.19927.21026.12226.70626.87826.15226.37924.71225.68824.99024.23926.72123.47524.76726.21928.36128.59927.91427.78425.69326.88126.21724.21827.91426.97528.52727.13928.98228.16928.05629.13626.29126.98726.58924.84827.54326.89628.87827.39028.06528.14129.04828.48426.63427.73527.13224.92428.96326.58927.93128.00929.22928.75928.40527.94525.91226.61926.07625.28627.66025.95126.39825.56528.86530.00029.26129.01226.99227.897類似的，文件/tsdldata/data/fancy.dat中包含了澳大利亞昆士蘭州的一個海灘小鎮(zhèn)紀(jì)念品商店的月銷量數(shù)據(jù)，時間從1987/01?1993/12(數(shù)據(jù)來源WheelwrightandHyndman,1998)。我們用下述代碼將其讀入至R中。>souvenir<-scan("/tsdldata/data/fancy.dat")Read84items>souvenirTS<-ts(souvenir,frequency=12,start=c(1987,1))>souvenirTSJanFeb MarAprMayJunJulAugSep1987Oct1664.81Nov Dec2397.53 2840.713547.293752.963714.744349.613566.3419885021.822499.816423.48 7600.605198.24 7225.1419756.214806.035900.884951.346179.124752.1519895496.434717.025835.10 12600.085702.63 9957.5828541.725304.786492.436630.807349.628176.6219908573.175921.109690.50 15151.845814.58 12421.2534061.016369.777609.127224.758121.227979.2519918093.064826.648476.70 17914.666470.23 9638.7730114.418821.178722.3710209.4811276.5512552.22 11637.39 13606.89 21822.11 45060.697615.03 9849.69 14558.40 11587.33 9332.56 13082.09 16732.7819888.61 23933.38 25391.35 36024.80 80721.7110243.24 11266.88 21826.84 17357.33 15997.79 18601.53 26155.1528586.52 30505.41 30821.33 46634.38 104660.672數(shù)據(jù)可視化上述章節(jié)中我們讀入并生成了3個時間序列數(shù)據(jù)分別為，kingsts，birthsTS和souvenirTS。數(shù)據(jù)分析的后續(xù)過程為可視化過程，讓有關(guān)人員可以從直觀角度來認(rèn)識數(shù)據(jù)的變化規(guī)律。plot函數(shù)提供了各種類的繪圖功能>plot(kingsts)

上圖中的時間序列看來的話，沒有比較明顯的周期性規(guī)律，直觀判斷上述時間序列類似于一個隨機(jī)波動序列(平穩(wěn)時間序列)，后續(xù)我們會探討如何判斷序列為平穩(wěn)時間序列。我們對紐約市月出生數(shù)量畫圖，結(jié)果為>plot(birthsTS)CDCDy垣q從上述結(jié)果看來，每個月的出生數(shù)量似乎含有一定的季節(jié)變化（周期性），在夏季較高，在冬季較低。另外每個周期內(nèi)的波動幅度（振幅）似乎是固定的，并不隨時間序列的水平（趨勢）變化而變化。隨機(jī)波動隨時間變化也不明顯。我們同樣對澳大利亞昆士蘭州海灘旅游小鎮(zhèn)紀(jì)念品商店的月銷量數(shù)據(jù)畫圖。S+S+心?3+^寸?心9其十導(dǎo)mmqjz8+擋從上圖中可以看出誤差波動和振幅都隨時間水平增加而增加。因此，我們需要對原始時間序列作出對數(shù)變換（取log），以便使用additive模型對數(shù)據(jù)開展分析預(yù)測。19S7 1988 1989 1990 1991 1992 1993 1994Time直觀上看來，上述的對數(shù)化時間序列的季節(jié)變化和隨機(jī)變化隨時間穩(wěn)定在一個范圍內(nèi)，不隨時間序列變化而變化。因此，時間序列對數(shù)化可以用于additive模型。3時間序列分解時間序列分解主要針對非平穩(wěn)時間序列，非平穩(wěn)時間序列包含4種主要趨勢：長期趨勢（Trend），循環(huán)趨勢，季節(jié)趨勢（Seasonal）和隨機(jī)趨勢（Random）?？梢岳胊dditive模型分解長期趨勢項(xiàng)和季節(jié)項(xiàng)。使用R語言的decompose函數(shù)。該函數(shù)可以分解出趨勢項(xiàng)，季節(jié)項(xiàng)和隨機(jī)項(xiàng)并利用additive模型描述時間序列。時間序列分解的主要參考的函數(shù)有decompose函數(shù)和stl函數(shù)。在這里我們先介紹decompose函數(shù)的主要思路。相加分解模型（additive）的主要形式為：xt=Tt+St+et另外還有相乘模型（multiplicative），主要形式為：Xt=Tt*St*et其中，在t時刻，xt為該時刻的觀測值，Tt為該時刻的趨勢值，St為該時刻的季節(jié)（周期）效應(yīng)，et為該時刻的無規(guī)則部分，也就是剩余項(xiàng)。一般而言，理想的分解模型結(jié)果中的剩余項(xiàng)應(yīng)該為一個均值為0

的隨機(jī)變量。我們先利用R的decompose函數(shù)，看看有關(guān)時間序列分解的結(jié)果，再從結(jié)果入手分析其具體數(shù)據(jù)分析處理流程。分解我們之前的birthsTS數(shù)據(jù)，可視化后的結(jié)果為：birthsTSDecompose<-decompose(birthsTS)plot(birthsTSDecompose)Decompositionofadditivetimeseries1946 1948 1950 1952 1954 1956 1958 1960Timedecompose函數(shù)提供了additive和multiplicative2種方法來估算各分解項(xiàng)。本文主要分析additive模型的分解步驟?，F(xiàn)在的已知條件：時間序列，包含了時間序列的頻率(frequency)，起始時間(start)。Stepl：分解出趨勢項(xiàng)。采用中心化移動平均值來計算趨勢項(xiàng)，1(：一號)J…+:+("J+1/f-1)TOC\o"1-5"\h\zT= 2 2 2—,te( ,l— )\o"CurrentDocument"t f 2 2當(dāng)f為奇數(shù)時采用上述計算方法0.5x +x+...+x +0.5xt-(*) t—(f)+i t+(f)—i t+(*)//f f、T= 2 2 2 J,te(—+1,l——)\o"CurrentDocument"t f 2 2當(dāng)f為偶數(shù)時采用上述計算方法，其中Tt為趨勢項(xiàng)，f為時間序列頻率，l為時間序列長度。結(jié)果為長度為l的時間序列，為便于后續(xù)的向量計算，當(dāng)t超出上述下標(biāo)的定義域時，其值為NA，如T1。R中提供filter函數(shù)計算2個序列的線性濾波(卷積)，上述時間序列xt的權(quán)重恰好為一個對稱序列，因此可以通過線性濾波函數(shù)實(shí)現(xiàn)中心化移動平均，見如下代碼：f<-frequency(birthsTS)Ts<-filter(birthsTS,c(0.5,rep(1,f-1),0.5)/f)identical(Ts,birthsTSDecompose$trend)[1]TRUEStep2:計算季節(jié)項(xiàng)。1、 采用將原始時間序列減去趨勢項(xiàng)S=x—T2、 將各個周期內(nèi)相同頻率下的值平均化，得到季節(jié)項(xiàng)figurefigure=YS/f,t(1,f),n =l對/取整，即max(n,nf<l)i=03、 將figure中心化，得到中心化的季節(jié)項(xiàng)Figure，代碼可表述為figure=figure-mean(figure)最終得到的的長度為f的季節(jié)項(xiàng)。假設(shè)時間序列的頻率為f，長度為l,設(shè)xt-Tt為St，則上述算法可以表示為:用表格可