《R语言实战 (第2版)》试读：15.1　在R中生成时序对象

本章内容  生成时间序列   分解时间序列   建立预测模型   预测未来值  全球变暖的速度有多快？十年后会产生什么影响？除了9.6节中的重复测量方差分析外，前面各章节探讨的都是横截面（cross-sectional）数据。在横截面数据集中，我们是在一个给定的时间点测量变量值。与之相反，纵向（longitudinal）数据则是随着时间的变化反复测量变量值。若持续跟踪某一现象，可能会获得很多了解。 本章，我们将研究在给定的一段时间内有规律地记录的观测值。对于这样的观测值，我们可以将其整合成形如Y1, Y2, Y3, …, Yt, …, YT的时间序列，其中Yt为Y在时间点t的值，T是时间序列中观测值的个数。 图15-1中是两个完全不同的时间序列（简称时序）。左边的序列为1960～1980年每股Johnson & Johnson的季度收入（单位：美元），数据集中共有84个观测值，即观测值依次对应21年间的每个季度。右边的序列为1749～1983年瑞士联邦观测台和东京天文观测台所观测到的月均相对太阳黑子数。太阳黑子时序更长一些，共有2820个观测值，对应235年间的每个月。 对时序数据的研究包括两个基本问题：对数据的描述（这段时间内发生了什么）以及预测（接下来将会发生什么）。对于Johnson & Johnson数据，我们可能有如下疑问。  Johnson & Johnson股价在这段时间内有变化吗？   数据会受到季度影响吗？股价是不是存在某种固定的季度变化？   我们可以预测未来的股价吗？如果可以的话，准确率有多高？  对于太阳黑子数据，我们可能有如下疑问。  哪个统计模型可以更好地描述太阳黑子的活动？   是不是有些模型可以更好地拟合数据？   在一个给定时间内，太阳黑子的数目是否是可预测的？在多大程度上可预测？  正确预测股价的能力关系到我能否早点退休去一个热带岛屿，而预测太阳黑子活动的能力则关系到我能否在这个热带岛屿上保持手机通信畅通。 图15-1　(a) 1960～1980年Johnson & Johnson每股季度收入（美元）的时序图； (b) 1749～1983年月均相对太阳黑子数的时序图 对时间序列数据未来值进行预测是基本的人类活动，对时序数据的研究在现实世界中也有着广泛的应用。经济学家尝试通过时序分析理解并预测金融市场；城市规划者基于时序数据预测未来的交通需求；气候学家通过时序数据预测全球气候变化；公司需要时序分析来预测产品的需求及未来销量；医疗保健人员需要根据时序数据研究疾病传播范围及某区域内可能出现的病例数；地震学家通过时序数据预测地震。在这些研究中，对于历史时间序列的分析都是必不可少的。由于不同类型的时间序列数据可能需要不同的方法，本章将研究多个不同的时序数据集。 描述时序数据和预测未来值的方法有很多，而R软件具备很多其他软件都不具备的精细时序分析工具。本章将介绍一些最常用的时序数据描述和预测方法以及对应的R函数。表15-1按在本章中的出现顺序给出了这些函数。 表15-1　时序分析会用到的函数 函　　数 程　序　包 用　　途 ts() stats 生成时序对象 plot() graphics 画出时间序列的折线图 start() stats 返回时间序列的开始时间 end() stats 返回时间序列的结束时间 frequency() stats 返回时间序列中时间点的个数 window() stats 对时序对象取子集 ma() forecast 拟合一个简单的移动平均模型 stl() stats 用LOESS光滑将时序分解为季节项、趋势项和随机项 monthplot() stats 画出时序中的季节项 seasonplot() forecast 生成季节图 （续） 函　　数 程　序　包 用　　途 HoltWinters() stats 拟合指数平滑模型 forecast() forecast 预测时序的未来值 accuracy() forecast 返回时序的拟合优度度量 ets() forecast 拟合指数平滑模型，同时也可以自动选取最优模型 lag() stats 返回取过指定滞后项后的时序 Acf() forecast 估计自相关函数 Pacf() forecast 估计偏自相关函数 diff() base 返回取过滞后项和（或）差分后的序列 ndiffs() forecast 找到最优差分次数以移除序列中的趋势项 adf.test() tseries 对序列做ADF检验以判断其是否平稳 arima() stats 拟合ARIMA模型 Box.test() stats 进行Ljung-Box检验以判断模型的残差是否独立 bds.test() tseries 进行BDS检验以判断序列中的随机变量是否服从独立同分布 auto.arima() forecast 自动选择ARIMA模型 表15-2给出了我们将分析的几个时序数据集，这些数据集在R中都可以找到，它们各有特点，适用的模型也各不相同。 表15-2　本章用到的数据集 时间序列 描　　述 AirPassengers 1949～1960年每个月乘坐飞机的乘客数 JohnsonJohnson 每股Johnson & Johnson股份每季度的收入 nhtemp 康涅狄格州纽黑文地区从1912年至1971年每年的平均气温 Nile 尼罗河的流量 sunspots 1749～1983年每月太阳黑子的数量 本章首先介绍生成、操作时序数据的方法，对它们进行描述并画图，将它们分解成水平、趋势、季节性和随机（误差）等四个不同部分。在此基础上，我们采用不同的统计模型对其进行预测。将要介绍的方法包括基于加权平均的指数模型，以及基于附近数据点和预测误差间关联的自回归积分移动平均（ARIMA）模型。我们还将介绍模型拟合和预测准确性的评价指标。最后，本章将给出关于时间序列的更多参考书目，以便读者继续学习。 15.1　在R中生成时序对象 在R中分析时间序列的前提是我们将分析对象转成时间序列对象（time-series object），即R中一种包括观测值、起始时间、终止时间以及周期（如月、季度或年）的结构。只有将数据转成时间序列对象后，我们才能用各种时序方法对其进行分析、建模和绘图。 一个数值型向量或数据框中的一列可通过ts()函数存储为时序对象： myseries <- ts(data, start=, end=, frequency=) 其中myseries是所生成的时序对象，data是原始的包含观测值的数值型向量，start参数和end参数（可选）给出时序的起始时间和终止时间，frequency为每个单位时间所包含的观测值数量（如frequency=1对应年度数据，frequency=12对应月度数据，frequency=4对应季度数据）。 代码清单15-1给出了一个示例。数据中包含从2003年1月开始，两年内的月度销量数据。 代码清单15-1　生成时序对象 > sales <- c(18, 33, 41, 7, 34, 35, 24, 25, 24, 21, 25, 20, 22, 31, 40, 29, 25, 21, 22, 54, 31, 25, 26, 35) > tsales <- ts(sales, start=c(2003, 1), frequency=12) > tsales Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec 2003 18 33 41 7 34 35 24 25 24 21 25 20 2004 22 31 40 29 25 21 22 54 31 25 26 35 > plot(tsales) > start(tsales) [1] 2003 1 > end(tsales) [1] 2004 12 > frequency(tsales) [1] 12 > tsales.subset <- window(tsales, start=c(2003, 5), end=c(2004, 6)) > tsales.subset Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec 2003 34 35 24 25 24 21 25 20 2004 22 31 40 29 25 21 在代码清单中，ts()函数被用于生成时序对象➊。生成后，我们可以将它以图像的形式显示在屏幕上，图15-2为时序数据对应的图像。我们可以用第3章中介绍的方法将图像变得更精炼，比如通过plot(tsales, type="0", pch=19)可以将线条变为连接起来的实心点。生成时序对象后，我们可以通过start()、end()、frequency()函数查看其性质➋，也可以通过window()函数生成原始数据的一个子时序➌。 图15-2　由代码清单15-1生成的销量数据的时序图。时间轴中的小数点用来表示一年中的某个位置，比如2003.5表示2003年7月1日（即2003年全年的一半）