1、時(shí)間序列分解算法和decompose函數(shù)實(shí)現(xiàn)李思亮 55531469目錄時(shí)間序列分解算法和decompose函數(shù)實(shí)現(xiàn)11 數(shù)據(jù)讀入并生成時(shí)間序列12 數(shù)據(jù)可視化33 時(shí)間序列分解6在時(shí)間序列分析的過程中，往往需要對(duì)時(shí)間序列作出初步分析，本文主要采用R語言作為分析平臺(tái)，從數(shù)據(jù)的讀入，可視化圖，分解（decompose）為趨勢項(xiàng)，季節(jié)項(xiàng)，隨機(jī)波動(dòng)等角度對(duì)數(shù)據(jù)開展分析的幾個(gè)案例。最后對(duì)分解算法作出初步描述并探討其預(yù)測預(yù)報(bào)中的潛在應(yīng)用。本文的數(shù)據(jù)和部分內(nèi)容主要采用/en/latest/中的內(nèi)

2、容，感興趣的讀者可以參考。1 數(shù)據(jù)讀入并生成時(shí)間序列對(duì)于數(shù)據(jù)分析來講，數(shù)據(jù)讀入是一個(gè)比較關(guān)鍵的步驟。常用的數(shù)據(jù)讀入函數(shù)有scan，read.table等。下面列舉了幾種常見的數(shù)據(jù)。首先是 and Mcleod, 1994）。> kings <- scan("Read 42 items> kings 1 60 43 67 50 56 42 50 65 68 43 65 34 47 34 49 41 13 35 53 56 16 43 69 59 48 59 86 55 68 51 33 49 67 77 81 67 71 81 68 70 77 56上述例子中，讀入了

3、連續(xù)42個(gè)公國國王的壽命并將其賦給變量kings如果我們希望對(duì)讀入數(shù)據(jù)開展分析，下一步就是將其轉(zhuǎn)化為時(shí)間序列對(duì)象（時(shí)間序列類），R提供了很多函數(shù)用于分析時(shí)間序列類數(shù)據(jù)。可以使用ts函數(shù)將變量轉(zhuǎn)化為時(shí)間序列類。> kingsts <- ts(kings)> kingstsTime Series:Start = 1 End = 42 Frequency = 1 1 60 43 67 50 56 42 50 65 68 43 65 34 47 34 49 41 13 35 53 56 16 43 69 59 48 59 86 55 68 51 33 49 67 77 81 67 7

4、1 81 68 70 77 56對(duì)于上述數(shù)據(jù)操作的好處是將數(shù)據(jù)轉(zhuǎn)化為特定的“時(shí)間序列類”便于我們使用R中的函數(shù)分析數(shù)據(jù)。有時(shí)候我們會(huì)按照一定的時(shí)間周期來收集數(shù)據(jù)，這個(gè)周期可能是季度，月，日，小時(shí)，分。在大數(shù)據(jù)時(shí)代，有些情況下的數(shù)據(jù)是按照秒來采集收集。這種情況下，我們需要對(duì)數(shù)據(jù)的周期或頻率進(jìn)行設(shè)置。這里采用ts函數(shù)中的frequency參數(shù)可以實(shí)現(xiàn)這種功能。比方說，若按1年為一個(gè)周期，我們的月度時(shí)間序列數(shù)據(jù)應(yīng)為frequency=12，若為季度時(shí)間序列數(shù)據(jù)，則可設(shè)置frequency=4。另外，還可以利用start參數(shù)，設(shè)置時(shí)間序列的起點(diǎn)，比如若我們一個(gè)周期共4個(gè)觀測，而第一個(gè)數(shù)據(jù)對(duì)應(yīng)為1986

5、年的起的第2個(gè)觀測，則可使用start=c(1986,2)。有一個(gè)紐約市月出生數(shù)量的數(shù)據(jù)集，從1946年1月至1959年12月。數(shù)據(jù)可通過 獲取。我們將其讀入至R中。> births <- scan("Read 168 items> birthsTS<-ts(births,frequency=12,start=c(1946,1)> birthsTS Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec1946 26.663 23.598 26.931 24.740 25.806 24.364 24.477 23.

6、901 23.175 23.227 21.672 21.8701947 21.439 21.089 23.709 21.669 21.752 20.761 23.479 23.824 23.105 23.110 21.759 22.0731948 21.937 20.035 23.590 21.672 22.222 22.123 23.950 23.504 22.238 23.142 21.059 21.5731949 21.548 20.000 22.424 20.615 21.761 22.874 24.104 23.748 23.262 22.907 21.519 22.0251950

7、22.604 20.894 24.677 23.673 25.320 23.583 24.671 24.454 24.122 24.252 22.084 22.9911951 23.287 23.049 25.076 24.037 24.430 24.667 26.451 25.618 25.014 25.110 22.964 23.9811952 23.798 22.270 24.775 22.646 23.988 24.737 26.276 25.816 25.210 25.199 23.162 24.7071953 24.364 22.644 25.565 24.062 25.431 2

8、4.635 27.009 26.606 26.268 26.462 25.246 25.1801954 24.657 23.304 26.982 26.199 27.210 26.122 26.706 26.878 26.152 26.379 24.712 25.6881955 24.990 24.239 26.721 23.475 24.767 26.219 28.361 28.599 27.914 27.784 25.693 26.8811956 26.217 24.218 27.914 26.975 28.527 27.139 28.982 28.169 28.056 29.136 26

9、.291 26.9871957 26.589 24.848 27.543 26.896 28.878 27.390 28.065 28.141 29.048 28.484 26.634 27.7351958 27.132 24.924 28.963 26.589 27.931 28.009 29.229 28.759 28.405 27.945 25.912 26.6191959 26.076 25.286 27.660 25.951 26.398 25.565 28.865 30.000 29.261 29.012 26.992 27.897類似的，文件 中包含了澳大利亞昆士蘭州的一個(gè)海灘小

10、鎮(zhèn)紀(jì)念品商店的月銷量數(shù)據(jù)，時(shí)間從1987/011993/12(數(shù)據(jù)來源Wheelwright and Hyndman, 1998)。我們用下述代碼將其讀入至R中。> souvenir <- scan("Read 84 items> souvenirTS<-ts(souvenir,frequency=12,start=c(1987,1)> souvenirTS Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec1987 1664.81 2397.53 2840.71 3547.29 3752.96 3714.74

11、 4349.61 3566.34 5021.82 6423.48 7600.60 19756.211988 2499.81 5198.24 7225.14 4806.03 5900.88 4951.34 6179.12 4752.15 5496.43 5835.10 12600.08 28541.721989 4717.02 5702.63 9957.58 5304.78 6492.43 6630.80 7349.62 8176.62 8573.17 9690.50 15151.84 34061.011990 5921.10 5814.58 12421.25 6369.77 7609.12 7

12、224.75 8121.22 7979.25 8093.06 8476.70 17914.66 30114.411991 4826.64 6470.23 9638.77 8821.17 8722.37 10209.48 11276.55 12552.22 11637.39 13606.89 21822.11 45060.691992 7615.03 9849.69 14558.40 11587.33 9332.56 13082.09 16732.78 19888.61 23933.38 25391.35 36024.80 80721.711993 10243.24 11266.88 21826

13、.84 17357.33 15997.79 18601.53 26155.15 28586.52 30505.41 30821.33 46634.38 104660.672 數(shù)據(jù)可視化上述章節(jié)中我們讀入并生成了3個(gè)時(shí)間序列數(shù)據(jù)分別為，kingsts，birthsTS和souvenirTS。數(shù)據(jù)分析的后續(xù)過程為可視化過程，讓有關(guān)人員可以從直觀角度來認(rèn)識(shí)數(shù)據(jù)的變化規(guī)律。plot函數(shù)提供了各種類的繪圖功能> plot(kingsts)上圖中的時(shí)間序列看來的話，沒有比較明顯的周期性規(guī)律，直觀判斷上述時(shí)間序列類似于一個(gè)隨機(jī)波動(dòng)序列（平穩(wěn)時(shí)間序列），后續(xù)我們會(huì)探討如何判斷序列為平穩(wěn)時(shí)間序列。我們對(duì)紐

14、約市月出生數(shù)量畫圖，結(jié)果為> plot(birthsTS)從上述結(jié)果看來，每個(gè)月的出生數(shù)量似乎含有一定的季節(jié)變化（周期性），在夏季較高，在冬季較低。另外每個(gè)周期內(nèi)的波動(dòng)幅度（振幅）似乎是固定的，并不隨時(shí)間序列的水平（趨勢）變化而變化。隨機(jī)波動(dòng)隨時(shí)間變化也不明顯。我們同樣對(duì)澳大利亞昆士蘭州海灘旅游小鎮(zhèn)紀(jì)念品商店的月銷量數(shù)據(jù)畫圖。從上圖中可以看出誤差波動(dòng)和振幅都隨時(shí)間水平增加而增加。因此，我們需要對(duì)原始時(shí)間序列作出對(duì)數(shù)變換（取log），以便使用additive模型對(duì)數(shù)據(jù)開展分析預(yù)測。直觀上看來，上述的對(duì)數(shù)化時(shí)間序列的季節(jié)變化和隨機(jī)變化隨時(shí)間穩(wěn)定在一個(gè)范圍內(nèi)，不隨時(shí)間序列變化而變化。因此，時(shí)間

15、序列對(duì)數(shù)化可以用于additive模型。3 時(shí)間序列分解時(shí)間序列分解主要針對(duì)非平穩(wěn)時(shí)間序列，非平穩(wěn)時(shí)間序列包含4種主要趨勢：長期趨勢（Trend），循環(huán)趨勢，季節(jié)趨勢（Seasonal）和隨機(jī)趨勢（Random）。可以利用additive模型分解長期趨勢項(xiàng)和季節(jié)項(xiàng)。使用R語言的decompose函數(shù)。該函數(shù)可以分解出趨勢項(xiàng)，季節(jié)項(xiàng)和隨機(jī)項(xiàng)并利用additive模型描述時(shí)間序列。時(shí)間序列分解的主要參考的函數(shù)有decompose函數(shù)和stl函數(shù)。在這里我們先介紹decompose函數(shù)的主要思路。相加分解模型（additive）的主要形式為：xt=Tt+St+et另外還有相乘模型（multiplic

16、ative），主要形式為：xt=Tt*St*et其中，在t時(shí)刻，xt為該時(shí)刻的觀測值，Tt為該時(shí)刻的趨勢值，St為該時(shí)刻的季節(jié)（周期）效應(yīng)，et為該時(shí)刻的無規(guī)則部分，也就是剩余項(xiàng)。一般而言，理想的分解模型結(jié)果中的剩余項(xiàng)應(yīng)該為一個(gè)均值為0的隨機(jī)變量。我們先利用R的decompose函數(shù)，看看有關(guān)時(shí)間序列分解的結(jié)果，再從結(jié)果入手分析其具體數(shù)據(jù)分析處理流程。分解我們之前的birthsTS數(shù)據(jù)，可視化后的結(jié)果為：> birthsTSDecompose<-decompose(birthsTS)> plot(birthsTSDecompose)decompose函數(shù)提供了additive

17、和multiplicative 2種方法來估算各分解項(xiàng)。本文主要分析additive模型的分解步驟。現(xiàn)在的已知條件：時(shí)間序列，包含了時(shí)間序列的頻率（frequency），起始時(shí)間（start）。Step1：分解出趨勢項(xiàng)。采用中心化移動(dòng)平均值來計(jì)算趨勢項(xiàng)，當(dāng)f為奇數(shù)時(shí)采用上述計(jì)算方法當(dāng)f為偶數(shù)時(shí)采用上述計(jì)算方法，其中Tt為趨勢項(xiàng)，f為時(shí)間序列頻率，l為時(shí)間序列長度。結(jié)果為長度為l的時(shí)間序列，為便于后續(xù)的向量計(jì)算，當(dāng)t超出上述下標(biāo)的定義域時(shí)，其值為NA，如T1。R中提供filter函數(shù)計(jì)算2個(gè)序列的線性濾波（卷積），上述時(shí)間序列xt的權(quán)重恰好為一個(gè)對(duì)稱序列，因此可以通過線性濾波函數(shù)實(shí)現(xiàn)中心化移動(dòng)平

18、均，見如下代碼：> f<-frequency(birthsTS)> Ts<-filter(birthsTS,c(0.5,rep(1,f-1),0.5)/f)> identical(Ts,birthsTSDecompose$trend)1 TRUE Step2：計(jì)算季節(jié)項(xiàng)。1、 采用將原始時(shí)間序列減去趨勢項(xiàng)2、 將各個(gè)周期內(nèi)相同頻率下的值平均化，得到季節(jié)項(xiàng)figure3、 將figure中心化，得到中心化的季節(jié)項(xiàng)Figure，代碼可表述為figure=figure-mean(figure)最終得到的的長度為f的季節(jié)項(xiàng)。假設(shè)時(shí)間序列的頻率為f，長度為l,設(shè)xt-Tt為St，則上述算法可以表示為:用表格可以表示為-S1S2SfSf+1S2f+1Snf+1-Sl figure1fugure2figuref注:-為NA數(shù)據(jù)，其中St序列本身或已包含了NA數(shù)據(jù)。Step3：計(jì)算周期性季節(jié)項(xiàng)Seasonal，周期性季節(jié)項(xiàng)的計(jì)算則為將figure重復(fù)至長度為l的序列。Step4：計(jì)算隨機(jī)項(xiàng)Random，隨機(jī)項(xiàng)的計(jì)算方法為4 總結(jié)和展望根據(jù)上述的分析可以看出R內(nèi)置的decompose函數(shù)的核心為趨勢項(xiàng)（step1）和季節(jié)項(xiàng)（step2）。關(guān)于趨勢