时间序列分析：ts/mts基础数据类型

时间序列是以时间为索引的数据。时间索引为等间距递增的时间序列称为规则时间序列，其他为不规则时间序列。由于R语言核心组件的构建坚持以通用性为原则，因此R的基础包中只定义了规则时间序列，而不规则时间序列的定义和处理则留给开发人员做扩展。

本文先介绍规则时间序列。

1 ts数据类型

R语言中基本的时间序列对象为ts，在stats基本包中定义，由同名构造函数ts()产生。使用前先看看它的参数：

args(ts)

## function (data = NA, start = 1, end = numeric(), frequency = 1, 
##     deltat = 1, ts.eps = getOption("ts.eps"), class = if (nseries > 
##         1) c("mts", "ts", "matrix") else "ts", names = if (!is.null(dimnames(data))) colnames(data) else paste("Series", 
##         seq(nseries))) 
## NULL

参数中最重要的是data、frequency、start和end。

• data 即时间序列中的观测值，可以是向量或矩阵，或者转换为向量或矩阵的数据类型（如data.frame）。默认为NA。
• frequency 是一个时间周期中的间隔频率。如果设为12，那么时间序列将自动识别为12个自然月；如果设置为4，识别为4个季度。ts另外一个参数deltat=1/frequency，两者只能设置其中一个。
• start和end 开始时间和结束时间，长度为1或2。如果长度为2则第二个值设定的是周期中的具体值，从1开始，不大于frequency。

dt <- sample(1:100, 24)
ts(dt, freq = 12, start = 2010)

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2010  28  80  85  29  96  45  30  69  70  93  60  42
## 2011  63  16  82  81  27  57  10  78  97  51  36  13

ts(dt, freq = 4, start = 2010)

##      Qtr1 Qtr2 Qtr3 Qtr4
## 2010   28   80   85   29
## 2011   96   45   30   69
## 2012   70   93   60   42
## 2013   63   16   82   81
## 2014   27   57   10   78
## 2015   97   51   36   13

ts(dt, freq = 1, start = 2010)

## Time Series:
## Start = 2010 
## End = 2033 
## Frequency = 1 
##  [1] 28 80 85 29 96 45 30 69 70 93 60 42 63 16 82 81 27 57 10 78 97 51 36
## [24] 13

ts(dt, freq = 7, start = 1)

## Time Series:
## Start = c(1, 1) 
## End = c(4, 3) 
## Frequency = 7 
##  [1] 28 80 85 29 96 45 30 69 70 93 60 42 63 16 82 81 27 57 10 78 97 51 36
## [24] 13

ts(dt, freq = 12, start = c(2010, 3))

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2010          28  80  85  29  96  45  30  69  70  93
## 2011  60  42  63  16  82  81  27  57  10  78  97  51
## 2012  36  13

frequency=1/4/12时ts显示形式比较合理，但其他周期数据不甚理想。这是因为ts的本意就是用来存储规则时间序列数据，而年度、季度和月度数据是最常用的。以7天为周期的统计数据相对少见，要想获得一个日历表要费点周折：

dd <- weekdays(as.Date("2016-10-17") + 0:6)
nn <- which(weekdays(as.Date("2016-10-01")) == dd)
print(ts(1:30, freq = 7, start = c(1, nn)), calendar = T)

##   p1 p2 p3 p4 p5 p6 p7
## 1                 1  2
## 2  3  4  5  6  7  8  9
## 3 10 11 12 13 14 15 16
## 4 17 18 19 20 21 22 23
## 5 24 25 26 27 28 29 30

data如果是向量，时间序列的class属性为“ts”；如果是matrix，则对象还具有mts（具有相同时间索引的多列数据）和matrix类属性。mts对象用于记录一个时间点上有多个观测指标的数据：

data(presidents, package = "datasets")
data(EuStockMarkets, package = "datasets")
head(presidents)

## [1] NA 87 82 75 63 50

head(EuStockMarkets)

##          DAX    SMI    CAC   FTSE
## [1,] 1628.75 1678.1 1772.8 2443.6
## [2,] 1613.63 1688.5 1750.5 2460.2
## [3,] 1606.51 1678.6 1718.0 2448.2
## [4,] 1621.04 1684.1 1708.1 2470.4
## [5,] 1618.16 1686.6 1723.1 2484.7
## [6,] 1610.61 1671.6 1714.3 2466.8

class(presidents)

## [1] "ts"

class(EuStockMarkets)

## [1] "mts"    "ts"     "matrix"

函数as.ts()和is.ts()用于ts数据转换和判断。

2 ts辅助函数

2.1 获取时间序列属性（tsp）

每个ts对象都有特定的的start、end和frequency，这三者合称时间序列属性，即tsp（time series properties），可以用tsp()函数获取。R还提供了helper函数用于获取ts的相应属性。

(ts1 <- ts(dt, start = 2015, frequency = 12))

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2015  28  80  85  29  96  45  30  69  70  93  60  42
## 2016  63  16  82  81  27  57  10  78  97  51  36  13

tsp(ts1)

## [1] 2015.000 2016.917   12.000

start(ts1)

## [1] 2015    1

end(ts1)

## [1] 2016   12

frequency(ts1)

## [1] 12

deltat(ts1)

## [1] 0.08333333

2.2 获取时间索引

ts对象的时间索引可以使用time()函数提取，返回值为向量：

(tt <- time(ts1))

##           Jan      Feb      Mar      Apr      May      Jun      Jul
## 2015 2015.000 2015.083 2015.167 2015.250 2015.333 2015.417 2015.500
## 2016 2016.000 2016.083 2016.167 2016.250 2016.333 2016.417 2016.500
##           Aug      Sep      Oct      Nov      Dec
## 2015 2015.583 2015.667 2015.750 2015.833 2015.917
## 2016 2016.583 2016.667 2016.750 2016.833 2016.917

ts的时间索引并不是我们想要的整数或因子形式，cycle()函数可能更有用：

(ck <- cycle(ts1))

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2015   1   2   3   4   5   6   7   8   9  10  11  12
## 2016   1   2   3   4   5   6   7   8   9  10  11  12

注意time和cycle函数获得的也是ts类对象，其tsp和原ts对象完全相同：

class(tt)

## [1] "ts"

tsp(tt)

## [1] 2015.000 2016.917   12.000

tsp(tt) == tsp(ts1)

## [1] TRUE TRUE TRUE

class(ck)

## [1] "ts"

tsp(ck)

## [1] 2015.000 2016.917   12.000

tsp(ck) == tsp(ts1)

## [1] TRUE TRUE TRUE

3 ts/mts对象操作

3.1 截取时间窗口

首先想到的是subset函数，它可以用于取数据子集，很多情况下它获得的数据类型和原对象一样，如data.frame类数据：

str(airquality)

## 'data.frame':	153 obs. of  6 variables:
##  $ Ozone  : int  41 36 12 18 NA 28 23 19 8 NA ...
##  $ Solar.R: int  190 118 149 313 NA NA 299 99 19 194 ...
##  $ Wind   : num  7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
##  $ Temp   : int  67 72 74 62 56 66 65 59 61 69 ...
##  $ Month  : int  5 5 5 5 5 5 5 5 5 5 ...
##  $ Day    : int  1 2 3 4 5 6 7 8 9 10 ...

str(subset(airquality, !is.na(Ozone)))

## 'data.frame':	116 obs. of  6 variables:
##  $ Ozone  : int  41 36 12 18 28 23 19 8 7 16 ...
##  $ Solar.R: int  190 118 149 313 NA 299 99 19 NA 256 ...
##  $ Wind   : num  7.4 8 12.6 11.5 14.9 8.6 13.8 20.1 6.9 9.7 ...
##  $ Temp   : int  67 72 74 62 66 65 59 61 74 69 ...
##  $ Month  : int  5 5 5 5 5 5 5 5 5 5 ...
##  $ Day    : int  1 2 3 4 6 7 8 9 11 12 ...

但subset函数用于ts对象时获得的仅仅是一个普通向量：

str(presidents)

##  Time-Series [1:120] from 1945 to 1975: NA 87 82 75 63 50 43 32 35 60 ...

tsp(presidents)

## [1] 1945.00 1974.75    4.00

tt <- time(presidents)
(pp <- subset(presidents, subset = tt > 1960 & tt < 1965))

##  [1] 62 61 57 72 83 71 78 79 71 62 74 76 64 62 57 80 73 69 69

str(pp)

##  num [1:19] 62 61 57 72 83 71 78 79 71 62 ...

ts对象子集的提取有一个专门函数window()，即窗口函数，其产生和ts函数类似：

window(x, start = NULL, end = NULL, frequency = NULL, deltat = NULL, extend = FALSE,
    ...)

常用的参数是start和end。如设置frequency/deltat，要注意真正取到的是哪些值，而且有可能会出错：

window(presidents, start = 1960, end = c(1969, 4))

##      Qtr1 Qtr2 Qtr3 Qtr4
## 1960   71   62   61   57
## 1961   72   83   71   78
## 1962   79   71   62   74
## 1963   76   64   62   57
## 1964   80   73   69   69
## 1965   71   64   69   62
## 1966   63   46   56   44
## 1967   44   52   38   46
## 1968   36   49   35   44
## 1969   59   65   65   56

window(presidents, start = 1960, end = c(1969, 4), frequency = 1)

## Time Series:
## Start = 1960 
## End = 1969 
## Frequency = 1 
##  [1] 71 72 79 76 80 71 63 44 36 59

window(presidents, start = 1960, end = c(1969, 4), frequency = 2)

## Time Series:
## Start = c(1960, 1) 
## End = c(1969, 2) 
## Frequency = 2 
##  [1] 71 61 72 71 79 62 76 62 80 69 71 69 63 56 44 38 36 35 59 65

window(presidents, start = 1960, end = c(1969, 4), frequency = 3)

##      Qtr1 Qtr2 Qtr3 Qtr4
## 1960   71   62   61   57
## 1961   72   83   71   78
## 1962   79   71   62   74
## 1963   76   64   62   57
## 1964   80   73   69   69
## 1965   71   64   69   62
## 1966   63   46   56   44
## 1967   44   52   38   46
## 1968   36   49   35   44
## 1969   59   65   65   56

3.2 合并数据

ts对象y可以用cbind函数能有效合并，前提是frequency必需一样。rbind不行。

ts1 <- ts(sample(1:100, 12), start = 1, freq = 12)
ts2 <- ts(rnorm(12), start = 1, freq = 12)
ts3 <- ts(runif(12), start = c(1, 4), freq = 12)
ts4 <- ts(rnorm(12), start = 2, freq = 4)
cbind(ts1, ts2)

##       ts1         ts2
## Jan 1  20 -1.38615633
## Feb 1  83 -0.44157612
## Mar 1  55 -0.49842737
## Apr 1  72 -0.88409143
## May 1  63  1.25211918
## Jun 1  41  1.13397750
## Jul 1  69  0.68428102
## Aug 1   5  0.22880012
## Sep 1  17 -0.23620775
## Oct 1  98  1.08904161
## Nov 1  75  2.09957979
## Dec 1  43  0.03089808

cbind(ts1, ts3)

##       ts1        ts3
## Jan 1  20         NA
## Feb 1  83         NA
## Mar 1  55         NA
## Apr 1  72 0.35403318
## May 1  63 0.93236546
## Jun 1  41 0.38389533
## Jul 1  69 0.93882046
## Aug 1   5 0.20153206
## Sep 1  17 0.27942223
## Oct 1  98 0.43420658
## Nov 1  75 0.24039645
## Dec 1  43 0.70481722
## Jan 2  NA 0.02078979
## Feb 2  NA 0.47112608
## Mar 2  NA 0.77757363

cbind(ts1, ts4)

## Error in .cbind.ts(list(...), .makeNamesTs(...), dframe = FALSE, union = TRUE): not all series have the same frequency

rbind(ts1, ts3)

##           [,1]       [,2]       [,3]       [,4]       [,5]       [,6]
## ts1 20.0000000 83.0000000 55.0000000 72.0000000 63.0000000 41.0000000
## ts3  0.3540332  0.9323655  0.3838953  0.9388205  0.2015321  0.2794222
##           [,7]      [,8]       [,9]       [,10]      [,11]      [,12]
## ts1 69.0000000 5.0000000 17.0000000 98.00000000 75.0000000 43.0000000
## ts3  0.4342066 0.2403965  0.7048172  0.02078979  0.4711261  0.7775736

3.3 前/后推移：lag函数

ts1

##   Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1  20  83  55  72  63  41  69   5  17  98  75  43

lag(ts1, k = 4)

##   Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 0                                  20  83  55  72
## 1  63  41  69   5  17  98  75  43

lag(ts1, k = -4)

##   Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1                  20  83  55  72  63  41  69   5
## 2  17  98  75  43

3.4 ts同比差值：diff函数

diff(ts1, lag = 1)

##   Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1  63 -28  17  -9 -22  28 -64  12  81 -23 -32

ts1 - lag(ts1, -1)

##   Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1  63 -28  17  -9 -22  28 -64  12  81 -23 -32

diff(ts1, lag = -4)

## Error in diff.ts(ts1, lag = -4): bad value for 'lag' or 'differences'

ts1 - lag(ts1, 4)

##   Jan Feb Mar Apr May Jun Jul Aug
## 1 -43  42 -14  67  46 -57  -6 -38

diff只能是当前值减去历史值，即lag>0；用lag函数更灵活些。

4 ts/mts对象作图

R为ts类型对象的plot函数设置了特殊的参数：

args(plot.ts)

## function (x, y = NULL, plot.type = c("multiple", "single"), xy.labels, 
##     xy.lines, panel = lines, nc, yax.flip = FALSE, mar.multi = c(0, 
##         5.1, 0, if (yax.flip) 5.1 else 2.1), oma.multi = c(6, 
##         0, 5, 0), axes = TRUE, ...) 
## NULL

plot.type参数相对有用：

plot(EuStockMarkets)
plot(EuStockMarkets, plot = "single", col = 1:4)
legend("topleft", colnames(EuStockMarkets), lty = 1, col = 1:4, text.col = 1:4,
    box.col = NA, inset = 0.05)

---

作者: ZGUANG@LZU

Created: 2016-10-25 二 11:12

Emacs 24.4.1 (Org mode 8.2.10)