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

2、容，感興趣的讀者可以參考。1 數據讀入并生成時間序列對于數據分析來講，數據讀入是一個比較關鍵的步驟。常用的數據讀入函數有scan，read.table等。下面列舉了幾種常見的數據。首先是 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、連續42個公國國王的壽命并將其賦給變量kings如果我們希望對讀入數據開展分析，下一步就是將其轉化為時間序列對象（時間序列類），R提供了很多函數用于分析時間序列類數據。可以使用ts函數將變量轉化為時間序列類。> 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對于上述數據操作的好處是將數據轉化為特定的“時間序列類”便于我們使用R中的函數分析數據。有時候我們會按照一定的時間周期來收集數據，這個周期可能是季度，月，日，小時，分。在大數據時代，有些情況下的數據是按照秒來采集收集。這種情況下，我們需要對數據的周期或頻率進行設置。這里采用ts函數中的frequency參數可以實現這種功能。比方說，若按1年為一個周期，我們的月度時間序列數據應為frequency=12，若為季度時間序列數據，則可設置frequency=4。另外，還可以利用start參數，設置時間序列的起點，比如若我們一個周期共4個觀測，而第一個數據對應為1986

5、年的起的第2個觀測，則可使用start=c(1986,2)。有一個紐約市月出生數量的數據集，從1946年1月至1959年12月。數據可通過 獲取。我們將其讀入至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類似的，文件 中包含了澳大利亞昆士蘭州的一個海灘小

10、鎮紀念品商店的月銷量數據，時間從1987/011993/12(數據來源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 數據可視化上述章節中我們讀入并生成了3個時間序列數據分別為，kingsts，birthsTS和souvenirTS。數據分析的后續過程為可視化過程，讓有關人員可以從直觀角度來認識數據的變化規律。plot函數提供了各種類的繪圖功能> plot(kingsts)上圖中的時間序列看來的話，沒有比較明顯的周期性規律，直觀判斷上述時間序列類似于一個隨機波動序列（平穩時間序列），后續我們會探討如何判斷序列為平穩時間序列。我們對紐

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

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

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

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

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