时间序列绘图
LJJ
2020/2/16
2 时间序列图形
2.1 ts对象
library(pacman)
p_load(tidyverse,fpp2,forecast)y <- ts(c(123,39,78,52,110), start=2012)
y## Time Series:
## Start = 2012
## End = 2016
## Frequency = 1
## [1] 123 39 78 52 110y <- ts(1:50, start=2003, frequency=12)
y## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2003 1 2 3 4 5 6 7 8 9 10 11 12
## 2004 13 14 15 16 17 18 19 20 21 22 23 24
## 2005 25 26 27 28 29 30 31 32 33 34 35 36
## 2006 37 38 39 40 41 42 43 44 45 46 47 48
## 2007 49 502.2 时间图
如果观测频率大于每周一次，可以采用多种方法来处理频率。例如，日观测数据可能具有周季节性（frequency=7）或者具有年度季节性（frequency=365.25）。类似地，一个每分钟观测一次的数据可能具有时季节性（frequency=60），可能是日季节性（frequency=24x60=1440），还可能是周季节性（frequency=24x60x7=10080），甚至可能具有年度周期性（frequency=24x60x365.25=525960）。在我们处理时间序列之前，确定其频率至关重要。
对于时间序列数据而言，我们从最简单的时间图开始。时间图是用将观测值与观测时间点作图，散点之间用直线连接。例如图2.1表示在澳大利亚两个最大的城市之间，Ansett航空公司的每周客流量。
melsyd %>% head()## Time Series:
## Start = c(1987, 26)
## End = c(1987, 31)
## Frequency = 52
## First.Class Business.Class Economy.Class
## 1987.481 1.912 NA 20.167
## 1987.500 1.848 NA 20.161
## 1987.519 1.856 NA 19.993
## 1987.538 2.142 NA 20.986
## 1987.558 2.118 NA 20.497
## 1987.577 2.048 NA 20.770autoplot(melsyd[,"Economy.Class"]) +
ggtitle("墨尔本 - 悉尼经济舱乘客客流量") +
xlab("年份") +
ylab("千")+
theme(plot.title = element_text(hjust = 0.5))
之后我们会频繁使用autoplot()函数，它可以自动画出你传给其第一个参数的内容。在本例中，它将melsyd[,“Economy.Class”]识别为一个时间序列，进而自动生成时间图。
该时间图直观地展现出数据具有的一些特征：
- 由于1989年当地的工业纠纷，当年的客流量为0。
- 在1992年中，由于一部分经济舱被商务舱取代，导致客流量大幅减少。
- 1991年下半年客流量大幅上升。
- 由于假日效应，在每年年初，客流量都会有一定幅度的下降。
- 这是序列存在长期波动，在1987年向上波动，在1988年向下波动，而在1990年和1991年又再次向上波动。
- 在某些时期存在缺失值。
对该数据进行预测建模时，需要考虑上述所有的特征，以便有效预测未来的客流量。
a10 %>% head(20)## Jan Feb Mar Apr May Jun Jul Aug
## 1991 3.526591 3.180891
## 1992 5.088335 2.814520 2.985811 3.204780 3.127578 3.270523 3.737851 3.558776
## 1993 6.192068 3.450857
## Sep Oct Nov Dec
## 1991 3.252221 3.611003 3.565869 4.306371
## 1992 3.777202 3.924490 4.386531 5.810549
## 1993autoplot(a10) +
ggtitle("澳大利亚降糖药物的销量") +
ylab("百万（美元）") +
xlab("年份")+
theme(plot.title = element_text(hjust = 0.5))
显然，图示的时间序列具有明显增长的趋势。同时，在上升趋势中伴随着明显的季节性。在每年年底，由于政府补贴计划，使得降糖药品更便宜，所以人们倾向于在年底囤积药物，从而导致年初的销售额大幅下降。因此，当我们对降糖药物的销量进行预测时，需同时考虑其趋势和季节性因素。
2.3 时间序列模式
我们通常使用例如“趋势”、“季节性”等词语描述时间序列。在深入研究时间序列模式时，应该更精确的定义这些词语。
趋势
当一个时间序列数据长期增长或者长期下降时，表示该序列有趋势。在某些场合，趋势代表着“转换方向”。例如从增长的趋势转换为下降趋势。
季节性
当时间序列中的数据受到季节性因素（例如一年的时间或者一周的时间）的影响时，表示该序列具有季节性。季节性总是一个已知并且固定的频率。由于抗糖尿病药物的成本在年底时会有变化，导致上述抗糖尿药物的月销售额存在季节性。
周期性
当时间序列数据存在不固定频率的上升和下降时，表示该序列有周期性。这些波动经常由经济活动引起，并且与“商业周期”有关。周期波动通常至少持续两年。许多初学者都不能很好的区分季节性和周期，然而这两个概念是完全不同的。当数据的波动是无规律时，表示序列存在周期性；如果波动的频率不变并且与固定长度的时间段有关，表示序列存在季节性。一般而言，周期的长度较长，并且周期的波动幅度也更大。
许多时间序列同时包含趋势、季节性以及周期性。当我们选择预测方法时，首先应该分析时间序列数据所具备的特征，然后再选择合适的预测方法抓取特征。
2.4 季节图
季节图和时间序列图很相似，不同之处是季节图是针对观察数据的“季节性”绘制的。下面的例子是降糖药物的销售情况。
ggseasonplot(a10, year.labels=TRUE, year.labels.left=TRUE) +
xlab("月份")+
ylab("百万（美元）") +
ggtitle("季节图：降糖药物销量")+
theme(plot.title = element_text(hjust = 0.5))
在本例中，在每年一月份降糖药物的销量都会大幅下降。实际上，患者会在12月下旬大量购买降糖药物，但是这部分销量会在一两周后才向政府登记。从上图还可以看出，2008年3月销量大幅下降（其他年份2月份至3月份的销量增加）。2008年6月份销量较少可能是由于销量数据收集不完整导致。
ggseasonplot(a10, polar=TRUE) +
xlab("月份")+
ylab("百万（美元）") +
ggtitle("极坐标季节图：降糖药物销量")+
theme(plot.title = element_text(hjust = 0.5))
2.5 子系列季节图
ggsubseriesplot(a10) +
xlab("月份")+
ylab("百万（每月）") +
ggtitle("子序列季节图：降糖药物销量")+
theme(plot.title = element_text(hjust = 0.5))
图中的水平线表示每月的平均销量。子系列季节图可以清晰的描绘出数据的潜在季节性形态，并且显示了季节性随时间的变化情况。这类图可以很好地查看各时期内数据的变化情况。在本例中，子系列季节图并没有明显地体现数据特性，但是这是观察季节性变化最有用的方式。
2.6 散点图
在此之前，我们所讨论的内容都是单个时间序列的可视化。此外，多个时间序列的可视化也是非常有用的。
我们可以在一张图上绘制两个时间序列的散点图来研究用电量和温度之间的关系。
elecdemand %>% head(10)## Multi-Seasonal Time Series:
## Start: 2014 1
## Seasonal Periods: 48 336 17520
## Data:
## Demand WorkDay Temperature
## [1,] 3.914647 0 18.2
## [2,] 3.672550 0 17.9
## [3,] 3.497539 0 17.6
## [4,] 3.339145 0 16.8
## [5,] 3.204313 0 16.3
## [6,] 3.100043 0 16.6
## [7,] 3.039468 0 16.6
## [8,] 3.012089 0 16.7
## [9,] 3.017270 0 16.2
## [10,] 3.026672 0 16.6elecdemand %>% class()## [1] "msts" "ts"elecdemand %>% as.data.frame() %>%
ggplot(aes(Temperature,Demand)) +
geom_point() +
ylab("用电量 (千兆瓦)") + xlab("温度 (摄氏度)")+
theme(plot.title = element_text(hjust = 0.5))
这个散点图可以很好的帮助我们理解变量之间的相互关系。从图中我们可以看出，当温度很高时，人们会大量的使用空调进行降温，进而导致用电量随之增加；当温度很低时，人们会使用空调取暖，也会使得用电量一定程度上增加。
2.6.1 相关性
需要注意的是，相关系数仅仅衡量了变量之间的线性关系，并且有时会导致错误的结果。在分析变量之间关系时，不仅要看相关系数值，而且要关注生成的图形。
2.6.2 散点图矩阵
当所分析的数据有多个变量时，将每个变量与其他变量进行比较也很有意义。如图2.11所示，表示澳大利亚新南威尔士五个地区的季度游客人数。
visnights %>% head(20)## NSWMetro NSWNthCo NSWSthCo NSWSthIn NSWNthIn QLDMetro QLDCntrl
## 1998 Q1 9.047095 8.565678 5.818029 2.679538 2.977507 12.106052 2.748374
## 1998 Q2 6.962126 7.124468 2.466437 3.010732 3.477703 7.786687 4.040915
## 1998 Q3 6.871963 4.716893 1.928053 3.328869 3.014770 11.380024 5.343964
## 1998 Q4 7.147293 6.269299 2.797556 2.417772 3.757972 9.311460 4.260419
## 1999 Q1 7.956923 9.493901 4.853681 3.224285 3.790760 12.671942 4.186113
## 1999 Q2 6.542243 5.401201 2.759843 2.428489 3.395284 9.582965 4.237806
## 1999 Q3 6.330364 5.542626 2.041827 2.893273 3.626120 11.192969 6.414599
## 1999 Q4 7.509212 6.383474 2.650677 2.814578 3.691498 9.871108 3.710498
## 2000 Q1 7.662491 8.337008 5.367951 2.553072 3.135834 11.709851 3.800144
## 2000 Q2 6.341802 4.923476 3.460795 2.368211 2.824604 9.357777 4.182855
## 2000 Q3 7.827301 5.134491 1.606173 2.934893 3.390675 11.363533 4.173056
## 2000 Q4 9.579562 6.389888 2.563254 2.561272 2.946215 9.923923 4.192215
## 2001 Q1 8.270488 8.259398 4.913846 2.378108 3.117019 11.281072 3.781164
## 2001 Q2 7.240427 5.442669 3.143754 2.881288 2.967460 8.812295 5.061589
## 2001 Q3 6.640490 5.767428 1.776069 2.470918 3.111254 11.639739 5.086359
## 2001 Q4 7.111875 5.611910 2.437988 2.508259 2.778671 11.172069 4.217269
## 2002 Q1 6.827826 9.504843 4.503678 2.539149 3.246055 12.252613 4.405266
## 2002 Q2 6.404992 6.032525 2.738989 2.325263 3.182464 8.556381 3.262611
## 2002 Q3 6.615760 5.824036 2.597074 3.618077 3.550722 11.348967 6.493243
## 2002 Q4 7.226376 6.521765 3.036104 2.556570 3.654710 11.991551 4.357780
## QLDNthCo SAUMetro SAUCoast SAUInner VICMetro VICWstCo VICEstCo
## 1998 Q1 2.137234 2.881372 2.591997 0.8948773 7.490382 2.4420048 3.3819722
## 1998 Q2 2.269596 2.124736 1.375780 0.9792509 5.198178 0.9605047 1.8279401
## 1998 Q3 4.890227 2.284870 1.079542 0.9803289 5.244217 0.7559744 1.3519521
## 1998 Q4 2.621548 1.785889 1.497664 1.5094343 6.274246 1.2716040 1.4934148
## 1999 Q1 2.483203 2.293873 2.247684 0.9635227 9.187422 2.3850583 2.8969289
## 1999 Q2 3.377830 2.197418 1.672802 0.9968803 4.992303 1.3288638 1.5479014
## 1999 Q3 5.577959 2.033513 1.104703 1.0576872 4.746483 0.7588587 0.9140970
## 1999 Q4 4.279066 2.252621 1.503242 0.7705422 4.684532 0.9422794 1.3416072
## 2000 Q1 2.909853 3.199796 2.480801 1.0666867 8.449047 2.2271104 2.6818755
## 2000 Q2 2.415396 1.967285 1.892488 0.7516705 5.391724 1.0137120 1.4295229
## 2000 Q3 5.428682 1.974432 1.196388 1.6079280 4.921812 0.7175117 1.1963615
## 2000 Q4 3.730405 2.379026 1.308082 1.0255240 5.884153 1.1046218 1.8069715
## 2001 Q1 2.111366 2.498923 2.386227 0.9742663 8.894267 2.3544679 2.7279212
## 2001 Q2 2.888621 1.830694 1.564455 1.1099215 5.724432 1.0716824 1.3704442
## 2001 Q3 4.467746 2.094040 1.032725 1.0319447 5.562026 0.7110844 1.2424391
## 2001 Q4 3.841950 2.226933 1.532747 0.9884083 6.807252 0.8922901 1.3825742
## 2002 Q1 2.474732 2.461231 2.211109 0.6846369 7.903871 2.0967922 2.8917124
## 2002 Q2 2.406979 2.421714 1.601205 0.8437067 4.735441 1.2323374 1.8148167
## 2002 Q3 4.786755 2.394183 1.173739 1.2871040 5.815456 0.8790911 0.9761707
## 2002 Q4 3.011368 2.182758 1.664355 1.1469908 6.731291 1.1075884 1.7994384
## VICInner WAUMetro WAUCoast WAUInner OTHMetro OTHNoMet
## 1998 Q1 5.326655 3.075779 3.066555 0.6949954 3.437924 2.073469
## 1998 Q2 4.441119 2.154929 3.334405 0.5576796 2.677081 1.787939
## 1998 Q3 3.815645 2.787286 4.365844 1.0061844 3.793743 2.345021
## 1998 Q4 3.859567 2.752910 4.521996 1.1725514 3.304231 1.943689
## 1999 Q1 4.588755 3.519564 3.579347 0.3981829 3.510819 2.165838
## 1999 Q2 4.070401 3.160430 3.408533 0.5960182 2.871867 1.803940
## 1999 Q3 4.114489 2.707726 3.979314 0.9514562 3.833231 1.612655
## 1999 Q4 3.723363 2.294065 3.364682 0.8324046 3.143406 1.652280
## 2000 Q1 4.958003 3.606988 3.728502 0.5481368 4.126248 2.026803
## 2000 Q2 4.186498 2.812306 3.340795 0.6126219 3.032598 1.716114
## 2000 Q3 3.765293 1.835913 3.755450 0.6568117 4.203969 1.748992
## 2000 Q4 3.827241 3.294711 2.906949 1.0925982 3.426293 1.883513
## 2001 Q1 4.610200 2.733989 2.737205 0.7301631 2.919347 1.867394
## 2001 Q2 3.939684 2.598865 3.633654 0.5994801 3.413370 1.853570
## 2001 Q3 3.927555 2.520268 3.550497 0.8604462 3.461560 1.961257
## 2001 Q4 3.893688 2.853196 3.933320 1.0063039 3.102415 1.411294
## 2002 Q1 5.550154 3.706924 3.815283 1.0006523 3.195784 1.999595
## 2002 Q2 4.100652 2.308935 3.217102 1.0988116 2.816705 1.846532
## 2002 Q3 4.143904 2.381530 3.926514 0.7976717 3.733359 1.749374
## 2002 Q4 4.190439 2.510406 3.696802 0.7382546 3.861258 2.000513autoplot(visnights[,1:5], facets=TRUE) +
ylab("每季度的游客人数（百万）")+
xlab("时间")
theme(plot.title = element_text(hjust = 0.5))## List of 1
## $ plot.title:List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 0.5
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi FALSE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## - attr(*, "class")= chr [1:2] "theme" "gg"
## - attr(*, "complete")= logi FALSE
## - attr(*, "validate")= logi TRUEvisnights[,1:5] %>% as.data.frame() %>% GGally::ggpairs()## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
对于图中的每一块区域，其行变量是垂直轴行的变量，列变量是水平轴的变量。有许多设置可以控制生成的图形的形态。在默认设置中，相关系数在图的右上方显示，散点图在左下方显示，对角线上是密度曲线。
2.7 滞后图
澳大利亚每季度啤酒产量的散点图，横轴表示时间序列的滞后阶数。
beer2 <- window(ausbeer, start=1992)
gglagplot(beer2)
图中不同颜色代表不同季节，每条线都按时间顺序连接。从图中可以看出，滞后四阶和滞后八阶有正相关关系，说明数据具有很强的季节性。二阶滞后图和六阶滞后图显示，第四季度的峰值对应第二季度的最低点。
2.8 自相关
ggAcf(beer2) +
ggtitle('') 
2.8.1 ACF 图中的趋势性和季节性
当数据具有趋势性时，短期滞后的自相关值较大，因为观测点附近的值波动不会很大。时间序列的ACF一般是正值，随着滞后阶数的增加而缓慢下降。
当数据具有季节性时，自相关值在滞后阶数与季节周期相同时（或者在季节周期的倍数）较大。
当数据同时具有趋势和季节性时，我们会观察到组合效应。
aelec <- window(elec, start=1980)
autoplot(aelec) + xlab("年份") + ylab("百万千瓦") +
theme(plot.title = element_text(hjust = 0.5))
ggAcf(aelec, lag=48) +
ggtitle('') 
自相关系数值随着滞后阶数增加而缓慢降低，是因为原时间序列中具有趋势变化，而图中的“圆齿状”形状是来源于原时间序列中的季节性变化。
2.9 白噪声
set.seed(30)
y <- ts(rnorm(50))
autoplot(y) + ggtitle("白噪声") +
theme(plot.title = element_text(hjust = 0.5))
ggAcf(y) +
ggtitle('')
对于白噪声而言，我们期望它的自相关值接近0。但是由于随机扰动的存在，自相关值并不会精确地等于0。对于一个长度为T的白噪声序列而言，我们期望在0.95的置信度下，它的自相关值处于±2/√T之间。我们可以很容易的画出ACF的边界值（图中蓝色虚线）。如果一个序列中有较多的自相关值处于边界之外，那么该序列很可能不是白噪声序列。
在上例中，序列长度T=50，边界为 ±2/√50=±0.28。所有的自相关值均落在边界之内，证明序列是白噪声。
2.10 练习
2.10.1 练习1
使用ggseasonplot()函数和ggsubseriesplot()函数来观察writing、 fancy、 a10和 h02四个序列的季节性特征。
writing %>% ggseasonplot() +
theme(plot.title = element_text(hjust = 0.5))
writing %>% ggsubseriesplot() +
theme(plot.title = element_text(hjust = 0.5))
2.10.2 练习2
使用autoplot()、 ggseasonplot()、 ggsubseriesplot()、 gglagplot()和 ggAcf()这些画图函数，探索序列hsales、 usdeaths、 bricksq、 sunspotarea和 gasoline的特征
月度数据
hsales %>%
autoplot()
hsales %>%
ggseasonplot()
hsales %>%
ggsubseriesplot()
hsales %>%
gglagplot()
hsales %>%
ggAcf()
明显季节性
usdeaths %>%
autoplot()
usdeaths %>%
ggseasonplot()
usdeaths %>%
ggsubseriesplot()
usdeaths %>%
gglagplot()
usdeaths %>%
ggAcf()
明显趋势
bricksq %>% head()## Qtr1 Qtr2 Qtr3 Qtr4
## 1956 189 204 208 197
## 1957 187 214bricksq %>%
autoplot()
bricksq %>%
ggseasonplot()
bricksq %>%
ggsubseriesplot()
bricksq %>%
gglagplot()
bricksq %>%
ggAcf()
sunspotarea %>%
autoplot()
# sunspotarea %>%
# ggseasonplot()
# sunspotarea %>%
# ggsubseriesplot()
sunspotarea %>%
gglagplot()
sunspotarea %>%
ggAcf()
gasoline %>% head()## Time Series:
## Start = 1991.1
## End = 1991.19582477755
## Frequency = 52.1785714285714
## [1] 6.621 6.433 6.582 7.224 6.875 6.947gasoline %>%
autoplot()
gasoline %>%
ggseasonplot() +
coord_flip()
# gasoline %>%
# ggsubseriesplot()
gasoline %>%
gglagplot()
gasoline %>%
ggAcf()