R语言中的时间序列分析·

1 数据集说明

AirPassengers 1949~1960年每月乘坐飞机的乘客数

JohnsonJohnson Johnson&Johnson每股季度收入

nhtemp 康涅狄格州纽黑文地区从1912年至1971年每年的平均气温

Nile 尼罗河的流量

sunspots 1749年~1983年月平均太阳黑子数

2 相关包

xts、forecast、tseries、directlabels

3 在R中生成时序对象

首先将分析对象转成时间序列对象（time-series object）。

xts包提供的xts类支持时间间隔规律的和时间间隔不规律的时间序列，还拥有很多用于操作时序数据的函数。

另外：基础R附带的ts用于储存时间间隔规律的耽搁时间序列，mts’用于储存时间间隔均匀的多个时间序列；zoo包提供的类可以储存时间间隔不规律的时间序列，xts包提供zoo类的超集，包含更多函数；tsbox包提供的函数可以将数据框转化为任意一种时间序列格式，也可以将一种时间序列格式转化为另一种时间序列格式。

创建xts时间序列：

library(xts)myseries<-xts(data,index)

# data是观测值的数值型向量；index表示观测值观测时间的日期向量。

3.1 生成时间序列对象

# 从2018年1月以来的两年的月度销售数据library(xts)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)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)date<-seq(from=as.Date("2018/1/1"),to=as.Date("2019/12/1"),by="month")sales.xts<-xts(sales,date)sales.xts

# 结果

xts格式的时间序列对象可使用方括号[]来设定子集：

# 返回2018年以来的所有数据sales.xts["2018"]

# 结果

# 返回2018年3月到2019年5月的所有数据sales.xts["2018-3/2019-5"]

函数apply用于对时间序列对象的每个不同时间段执行一个函数，在将时间序列聚集合成更大的时间段时很实用。

newseries<-apply.period(x,FUN,…)

# period可以是daily、weekly、monthly、quarterly或yearly；x是一个xts时间序列对象，FUN是要应用的函数，…是传递给FUN的参数。

# 返回8个季度销售总量的时间序列quarterlies<-apply.quarterly(sales.xts,sum)quarterlies

# 结果

# sum可以替换为mean、median、min、max或其他任何返回单个值的函数

3.2 绘制时间序列图

forecast包中的autoplot()可将时序数据绘制为ggplot2图形

library(ggplot2)library(forecast)autoplot(sales.xts)

# 结果

图 3-1 销量数据的时序图，autoplot()提供的默认格式

autoplot(sales.xts)+geom_line(color="blue")+scale_x_date(date_breaks="1 months",date_labels="%b %y")+labs(x="",y="Sales",title='Customized Time Series Plot')+theme_bw()+theme(axis.text.x=element_text(angle=90,vjust=0.5,hjust=1),panel.grid.minor.x=element_blank())

# 选项date_breaks指定刻度间的距离，值可以是“1 day”“2 weeks”“5 years”等；“%b%y”指定了月份（3个字母）和年份（2位数），以及两者之间的空格；最后，选择了黑白主题，x轴标签旋转了90度，取消了垂直的小网格线。

# 结果

图 3-2 销量数据时序图，定义了颜色、更美观的标签和更清晰的主题元素

4 时序的平滑化和季节项分解

对时序数据建立复杂模型之前需要对其进行描述和可视化。我们对时序进行平滑化以探究其总体趋势，并对其进行分解以观察时序中是否存在季节项。

4.1 通过简单移动平均进行平滑处理（Nile数据集）

处理时序数据的第一步是画图。使用数据集Nile，是埃及阿斯旺市在1871年至1970年间所记录的尼罗河的年度流量。由图，数据总体呈下降趋势，但不同年份的变动非常大。（图 4-1 ）

时序数据集中通常有很显著的随机或误差成分。为了辨明数据中的规律，我们希望能够撇开这些波动，画出一条平滑曲线。画出平滑曲线的最简单方法是简单移动平均。比如，每个数据点都可用这一点和其前后两个点的平均值来表示，这就是居中移动平均。

R中可以做简单移动平均的函数：TTR包中的SMA()函数；zoo包中的rollmean()函数；forecast包中的ma()函数。

# 时序数据的原始数据图、平滑后的图library(ggplot2)library(forecast)theme_set(theme_bw())ylim<-c(min(Nile),max(Nile))autoplot(Nile)+ggtitle(“Raw time series”)+scale_y_continuous(limits=ylim)autoplot(ma(Nile,3))+ggtitle(“Simple Moving Average (k=3)”)+scale_y_continuous(limits=ylim)autoplot(ma(Nile,7))+ggtitle(“simple Moving Average (k=7)”)+scale_y+continuous(limits=ylim)autoplot(ma(Nile,15))+ggtitle(“simple Moving Average (k=15)”)+scale_y+continuous(limits=ylim)

# 结果

图 4-1 1871~1970年阿斯旺水站观测到的尼罗河的年流量

图 4-2 平滑水平k=3

图 4-3 平滑水平k=7

图 4-4 平滑水平k=15

# 随着k的增大，图像变得越来越平滑。

# 尼罗河的·流量从1892年到1900年有明显下降；其他的变动不是太好解读，比如1941年到1961年流量似乎略有上升，但这也可能只是一个随机波动。

我们需要找到最能画出数据中规律的k，避免过平滑或者欠平滑。没有什么特别的科学理论来指导k的选取，只需要多次尝试不同的k，再决定一个最好的k。

对于间隔大于1的时序数据（存在季节想），需要了解的不仅是总体趋势。此时需要通过季节项分解帮助我们探究季节性波动以及总体趋势。

4.2 季节项分解

存在季节项的时序数据（如月度数据、季度数据等）可以被分解为趋势项、季节项和随机项。趋势项（trend component）能捕捉到长期变化；季节项（seasonal component）能捕捉到一年内的周期性变化；随机（误差）项（irregular/error component）能捕捉到那些不能被趋势项或季节项解释的变化。

可以通过相加模型或相乘模型来分解数据。在相加模型中，各项之和等于对应的时序值，时刻t的观测值等于这一时刻的趋势项、季节项以及随机项之和；而在相乘模型中，时刻t的观测值等于趋势项、季节项和随机项相乘。

下面举例说明相加模型与相乘模型的区别：假设有一个时序，记录了10年来摩托车的月销量。在具有季节项的相加模型中，11月和12月的销量一般会增加500（圣诞节销售高峰），而1月（一般是销售淡季）的销量则会减少200.此时，季节性的波动量和当时的销量无关。在具有季节项的相乘模型中，11月和12月的销量会增加20%，1月的销量会减少10%，即季节项的波动量和当时的销量是成比例的，这与相加模型不一样。这也使得在很多时候，相乘模型比相加模型更切合实际。

将时序分解为趋势项、季节项和随机项的常用方法是用loess平滑做季节项分解，通过stl()函数实现：

stl(ts, s.window=, t.window=)

# ts是将要分解的时序；s.window控制季节项变化的速度；t.window控制趋势项变化的速度。设置s.windows=”periodic”可使得季节项宰割年检都一样。参数ts和s.windows必须提供。

# stl()函数只能处理相加模型，但是相乘模型可以通过对数变换转换成相加模型。

例 （使用数据集AirPassengers）

该时序数据集描述了1949年~1960年每个月国际航班的乘客数（单位：千）。

library(ggplot2)library(forecast)autoplot(AirPassengers)

# 结果

图 4-5 AirPassengers时序图

# 序列的波动随着整体水平的增长而增长，相乘模型更适合这个序列。

# 对数变换AirPassengers1<-log(AirPassengers)autoplot(AirPassengers1,ylab="log(AirPassengers)")

# 结果

# 分解时间序列## 季节项分解，将季节项限定为每年都一样fit<-stl(AirPassengers1,s.window=”period”)autoplot(fit)

# 结果

图 时序图、季节项分解图、趋势项分解图、随机项分解图

# 时序的趋势为单调增长，季节项表明（可能为假期）夏季乘客更多。

# 每个图的y轴尺度不同，通过图中右侧的灰色长条来指示量级，即每个长条代表的量级一样

# 返回对数变换后的时序fit$time.series

# 结果（部分）

# stl()函数返回的对象中有一项是time.series，包括每个观测值中的各分解项——趋势项、季节项以及随机项的值。

# 将结果转化为原始尺度exp(fit$time$series)

# 结果（部分）

# 观察季节项，7月的乘客数增长了24%（乘子为1.24），11月的乘客数减少了20%（乘子为0.8）

例 通过forecast包提供的其它工具对季节项分解进行可视化（月度图和季节图）

library(forecast)library(ggplot2)library(directlabels)# 月度图ggmonthplot(AirPassengers)+labs(title=”Month plot: AirPassengers”,x=””,y=”Passengers (thousands)”)

# 结果

图 AirPassengers时序的月度图。月度图显示了按月划分的子序列（从1949年到1960年，每年的所有1月的点连接起来，所有2月的点连接起来，以此类推），以及每个字序列的平均值。每个月的增长趋势几乎一致，大多数乘客在7月和8月乘坐飞机。

# 季节图p<-ggseasonplot(AirPassengers)+geom_point()+labs(title=”Seasonal plot: AirPassengers”,x=””,y=”Passengers (thousands)”)direct.label(p)

# 结果

图 季节图显示每年的子序列。我们可以从图中看到类似的规律，包括每年乘客的增长趋势和相同的季节性规模。默认情况下，ggplot2包将为年份变量创建图例。使用directlabels包可以将年份标签直接放置在图中时序的每条线旁边。

5 指数预测模型

指数模型是用来预测时序未来值的最常用模型。这类模型相对比较简单，实践证明它们的短期预测能力良好。不同指数模型建模时选用的因子可能不同。比如，单指数模型（simple/single exponential model）拟合的是只有水平项和时间点i处随机项的时间序列，这时认为时间序列不存在趋势项和季节项；双指数模型（double exponential model；也叫Holt指数平滑，Holt exponential smoothing）拟合的是有水平项和趋势项的时序；三指数模型（triple exponential model；也叫Holt-Winters指数平滑，Holt-Winters exponential smoothing）拟合的是有水平项、趋势项以及季节项的时序。

Forecast包中的ets()函数可以拟合指数模型：

ets(ts,model=”zzz”)

# ts是要分析的时序，限定模型的字母有3个，第1个字母代表随机项，第2个字母代表趋势项，第3个字母代表季节项。可选的字母包括：A（相加模型）、M（相乘模型）、N（无）、z（自动选择）。

用于拟合单指数、双指数和三指数预测模型的函数：

类型	参数	函数
单指数	水平项	ets(ts,model=”ANN”) ses(ts)
双指数	水平项、趋势项	Ets(ts,model=”AAN”) holt(ts)
三指数	水平项、趋势项、季节项	ets(ts,model=”AAA”) hw(ts)

函数ses()、holt()、hw()都是函数ets()的便捷包装，函数中有事先默认设定的参数值。

5.1 单指数平滑

单指数平滑根据现有的时序值的加权平均对未来值做短期预测，权数选择的宗旨是使得距离现在越远的观测值对平均数的影响越小。

一步向前预测可看做当前值和全部历史值的加权平均。

例 （使用数据集nhtemp）

nhtemp时序中有康涅狄格州纽黑文市从1912年至1971年每年的平均气温。

# 拟合模型library(forecast)fit<-ets(nhtemp,model=”ANN”)fit

# A表示可加误差，NN表示时序中不存在趋势项和季节项。

# 结果

# α值比较小（α=0.18），说明预测时同时考虑了离现在较近和较远的观测值，这样的α值可以最优化模型在给定数据集上的拟合效果。

# 一步向前预测forecast(fit,1)

# 结果

# 一步向前预测的结果是51.87031℉，其95%的置信区间为49.62512℉到54.11549℉。

# 输出准确度度量accuracy(fit)

# 结果

# 预测准确度度量

一般来说，平均误差和平均百分比误差用处不大，因为正向和负向的误差会抵消掉

RMSE给出了平均误差平方和的平方根，即1.126℉；平均绝对百分比误差给出了误差在真实值中的占比，它没有单位，可以用于比较不同时序间的预测准确度，但它同时假定测量尺度中存在一个真实为零的点（比如每天的乘客数），但温度中并没有一个真实为零的点，因此这里不能使用这个统计量；平均绝对标准化误差是最新的一种准确度评估指标，通常用于比较不同尺度的时序间的预测准确度。这几种预测准确度评估指标中，并不存在某种最优评估指标，不过相对来说，RMSE最有名、最常用。

5.2 Holt指数平滑和Holt-Winters指数平滑

Holt指数平滑可以对有水平项和趋势项（斜率）的时序进行拟合。

Ets函数中的平滑参数alpha控制水平项的指数型下降，beta控制趋势项的指数型下降。两个参数的有效范围都是[0,1]，参数取值越大意味着越近的观测值的权重越大。

Holt-Winters指数平滑可用来拟合含有水平项、趋势项以及季节项的时间序列。此时，模型可表示为：Yt=level+slope*t+st+irregulart

st代表时刻t的季节项；ets()函数的参数gamma控制季节项的指数下降，取值范围是[0,1]，取值越大意味着越近的观测值的季节项权重越大。

例 用指数模型预测未来的乘客数（Holt-Winters指数平滑）

# 有水平项、趋势项以及季节项的指数平滑模型## 平滑参数llibrary(forecast)fit<-ets(log(AirPassengers),model=”AAA”)fit

# 结果

# 给出了3个平滑参数，即水平项0.6795、趋势项0.0031、季节项1e-04。趋势项的值很小，意味着近期观测值的趋势不需要更新

accuracy(fit)

# 结果

## 未来值预测pred<-forecast(fit,5)pred

# 结果

autoplot(pred)+labs(title=”Forecast for Air Travel”,y=”Log(AirPassengers)”,x=”Time”)

# 结果

基于Holt-Winters指数平滑模型的5年航线乘客数预测（对数变换后，单位：千）

## 用原始尺度预测pred$mean<-exp(pred$mean)pred$lower<-exp(pred$lower)pred$upper<-exp(pred$upper)p<-cbind(pred$mean,pred$lower,pred$upper)dimnames(p)[[2]]<-c("mean","Lo 80","Hi 80","Lo95","Hi 95")p

# 结果

5.3 ets()函数和自动预测

ets()函数还可以用来拟合有可乘项的指数模型，加入抑制项（dampening component），以及进行自动检测。

通过ets(AirPassengers,model=”MAM”)函数对原始数据拟合可乘模型。此时，仍假定趋势项可加，但季节项和误差项可乘。当采用可乘模型时，准确度统计量和预测值都基于原始尺度（以千为单位的乘客数），这也是它的一个明显优势。

ets()函数也可以用来拟合抑制项。时序预测一般假定序列的长期趋势是一直向上的（如房地产市场），而一个抑制项则使得趋势项在一段时间内靠近一条水平渐近线。在很多问题中，一个有抑制项的模型往往更符合实际预测。

最后，可以通过ets()函数最懂选取对原始数据拟合优度最高的模型。

例 自动选取最优拟合模型，使用ets()函数进行自动指数预测（使用数据集JohnsonJohnson）

library(forecast)fit<-ets(JohnsonJohnson)fit

# 结果

# 没有指定模型，软件自动搜索了一系列模型，并在其中找到最小化拟合标准（默认为对数似然）的模型。所选模型同时有可乘趋势项、季节项和误差项（随机项）。

# 下8个季度的预测autoplot(forecast(fit))+labs(x=”Time”,y=”Quarterly Earnings (Dollars)”,title=”Johnson and Johnson Forecasts”)

# 结果

带趋势项和季节项的可乘指数平滑模型，其中预测值由虚线表示，80%和95%置信区间分别由淡蓝色和深蓝色表示

6 ARIMA预测模型

在ARIMA预测模型中，预测值表示为由最近的真实值和最近的预测误差（残差）组成的线性函数。

6.1 概念介绍

时序的滞后阶数：向后追溯的观测值的数量

表 Nile时序的不同滞后阶数

0阶滞后项（Lag 0）代表没有移位的时序，一阶滞后（Lag 1）代表时序向左移动一位，二阶滞后（Lag 2）代表时序向左移动两位，以此类推。

Lag(ts,k)

# 把时序变成k阶滞后，ts指代目标序列，k为滞后项阶数。

自相关度量时序中各个观测值之间的相关性。相关性构成的图即自相关函数图（autocorrelation function plot，ACF图）。ACF图可用于为ARIMA模型选择合适的参数，并评估最终模型的拟合效果。

foecast包中的Acf()函数可以生成ACF图：

acf(ts)

# ts是原始时序

偏自相关图（PACF）使用forecast包中的Pacf()函数绘制：

Pacf(ts)

PACF图也可以用来找到ARIMA模型中最适宜的参数

ARIMA模型主要用于拟合具有平稳性（或可以被转换为平稳序列）的时间序列。在一个平稳的时序中，序列的统计性质并不会随着时间的推移而改变；另外，对任意滞后阶数k，序列的自相关性不改变。

一般来说，拟合ARIMA模型前需要变换序列的值以保证方差为常数，比如对数变换、Box-Cox变换等。

一般假定平稳性时序有常数均值，这样的序列中肯定不含有趋势项。非平稳的时序可以通过差分来转换为平稳性序列。对序列的一次差分可以移除序列中的线性趋势，二次差分移除二次项趋势，三次差分移除三次项趋势。在实际操作中，对序列进行两次以上的差分通常都是不必要的。

通过函数diff()对序列进行差分：

diff(ts,differences=d)

# d是对序列ts的差分次数，默认值为d=1

# forecast包中的ndiffs()函数可以帮我们找到最优的d值：ndiffs(ts)

平稳性一般可以通过时序图直观判断。如果方差不是常数，需要对数据做变换；如果数据中存在趋势项，需要对其进行差分。也可以通过ADF（Augmented Dickey-Fuller）统计检验来验证平稳性假定。tseries包中的adf.test()函数可以用来做ADF检验：adf.test(ts)。如果结果显著，则认为序列满足平稳性。

总之，通过ACF图和PACF图来为ARIMA模型选定参数。平稳性是ARIMA模型中的一个重要假设，可以通过数据变换和差分使得序列满足平稳性嘉定。

下面我们就可以拟合出有自回归（autoregressive，AR）项、移动平均（moving averages，MA）项或者两者都有（ARMA）的模型了，最后，检验有ARMA项的ARIMA模型，并对其进行差分以保证平稳性。

6.2 ARMA和ARIMA模型

ARIMA(p,d,q)模型：时序被差分了d次，且序列中的每个观测值都是用过去的p个观测值和q个残差的线性组合表示的。预测是“无误差的”或完整的。

建立ARIMA模型的步骤：

（1）确保时序是平稳的；

（2）找到一个（或几个）合理的模型（即选定可能的p值和q值）；

（3）拟合模型；

（4）从统计假设和预测准确度等角度评估模型的拟合效果；

（5）预测’

例 对Nile时序拟合ARIMA模型

（1）验证序列的平稳性

画出序列的折线图并判别其平稳性

library(forecast)library(tseries)autoplot(Nile)

# 结果

1871年~1970年在阿斯旺地区测量的尼罗河的年流量时序图

ndiffs(Nile)

# 结果

dNile<-diff(Nile)autoplot(dNile)

# 结果

被差分一次后的时序图，差分后原始图中下降的趋势被移除了

adf.test(dNile)

# 结果

# 对差分后的时序做ADF检验，结果显示时序此时是平稳的，可以继续进行下一步

（2）选择模型

通过ACF图和PACF图来选择备选模型

autoplot(Acf(dNile))autoplot(Pacf(dNile))

# 结果

# 需要为ARIMA模型指定参数p、d、q，从前文可得d=1

选择ARIMA模型的方法

在ACF图中，在滞后项为一阶时有一个比较明显的自相关；当滞后阶数逐渐增加时，偏相关逐渐减小至零。因此，可以考虑ARIMA(0,1,1)模型

（3）拟合模型

使用arima()函数拟合ARIMA模型：

arima(ts,order=c(q,d,q))# 拟合ARIMA模型library(forecast)fit<-arima(Nile,order=c(0,1,1))fit

# 结果

accuracy(fit)

# 结果

# 指定了d=1，函数将对时序做一阶差分，因此直接将模型应用于原始时序即可。函数返回了移动平均项的系数（-0.73）以及模型的AIC值。如果还有其他备选模型，则可以通过比较AIC值来得到最合理的模型，比较的准则是AIC值越小越好。对百分比误差的绝对值做平均的结果是13%。

（4）模型评估

一般来说，一个模型如果合适，那么模型的残差应该满足均值为0的正态分布，并且对于任意的滞后阶数，残差自相关系数都应该为零。也就是说，模型的残差应该满足独立正态分布（残差间没有关联）。

# 模型拟合评估library(ggplot2)## 提取残差df<-data.frame(resid=as.numeric(fit$residuals))## 创建Q-Q图ggplot(df,aes(sample=resid))+stat_qq()+stat_qq_line()+labs(title=”Normal Q-Q Plot”)

# 结果

判断序列残差是否满足正态性假设的正态Q-Q图

# 计算值如果服从正态分布，则数据中的点会落在图中的线上，显然本例的结果还不错## 检验所有滞后阶数的自相关系数是否为零Box.test(fit$residuals,type=”Ljung-Box”)

# 结果

# Box.test()函数可以检验残差的自相关系数是否都为零。本例中，模型的残差没有通过显著性检验，可以认为残差的自相关系数为零。ARIMA模型能够较好地拟合本数据。

（5）预测

如果模型残差不满足正态性假设或自相关系数为零的假设，则需要调整模型、增加参数或改变差分次数。选定模型后，就可以用它来做预测。

# 用ARIMA模型做预测（forecast()函数）forecast(fit,3)

# 结果

autoplot(forecast(fit,3))+labs(x=”Year”,y=”Annual Flow”)

# 结果

用ARIMA(0,1,1)模型对Nile序列做接下来三年的预测。

# 图中的黑线代表点估计，浅蓝色和深蓝色区域分别代表80%置信区间和95%的置信区间

6.3 ARIMA模型的自动预测

使用forecast包中的auto.arima()函数实现最优ARIMA模型的自动选取

# ARIMA模型的自动预测（使用sunspots数据集）library(forecast)fit<-auto.arima(sunspots)fit

# 结果

# 函数选定ARIMA模型的参数为p=2、d=1、q=2，与其他模型相比，这种情况下得到的模型的AIC值最小。

forecast(fit,3)

# 结果

accuracy(fit)

# 结果

# 由于序列中存在值为零的观测值，MPE和MAPE这两个准确度度量都失效了（这也是这两个统计量的一个缺陷）