时序分析会用到的函数

函数	程序包	用途
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模型

---

生成时序对象

一个数值型向量或数据框中的一列可通过ts()函数存储为时序对象。

myseries <- ts(data, start=, end=, frequency=)

其中，myseries是所生成的时序对象，data是原始的包含观测值的数值型向量，start参数和end参数（可选）给出时序的起始时间和终止时间，frequency参数（可选）为每个单位所包含的观测值数量（如frequency=1对应年度数据，frequency=12对应月度数据，frequency=4对应季度数据）。

我们可以通过window()函数生成原始数据的一个子时序。

subset <- window(data, start=, end=)

---

通过简单移动平均进行平滑处理

时序数据集中通常由很显著的随机或误差成分。为了辨明数据中的规律，我们总是希望能够撇开这些波动，画出一条平滑曲线。画出平滑曲线的最简单的办法是简单移动平均。

比如每个数据点都可用这一点和其前后两个点的平均值来表示，这就是居中移动平均(centered moving average)。

\[S_t=\frac{Y_{t-q}+\cdots+Y_t+\cdots+Y_{t+q}}{2q+1}\]

其中，\(S_t\)是时间点\(t\)的平滑值，\(k=2q+1\)是每次用来平均的观测值的个数。

居中移动平均法的代价是，每个时序集中我们会损失最后的\(\frac{k-1}{2}\)个观测值。

R中有几个函数都可以做简单移动平均，包括TTR包中的SMA()函数，zoo包中的rollmean()函数，forecast包中的ma()函数。

使用自带的尼罗河数据集(Nile)为例子。

library(forecast)
opar <- par(no.readonly=TRUE)
par(mfrow=c(2,2))
ylim <- c(min(Nile), max(Nile))
plot(Nile,main = "Raw time series")
plot(ma(Nile,3), main = "Simple Moving Averages (k=3)", ylim=ylim)
plot(ma(Nile,7), main = "Simple Moving Averages (k=7)", ylim=ylim)
plot(ma(Nile,15), main = "Simple Moving Averages (k=15)", ylim=ylim)
par(opar)

【注】报错的话可能需要先安装forecast包。

install.packages('forecast')

然后我们就可以得到下图。

可以发现随着\(k\)的增大，图像变得越来越平滑。因此我们需要找到最能画出数据中规律的\(k\)，避免过平滑或者欠平滑，这里并没有什么特别的科学理论来指导选取，只能尝试多个不同的\(k\)再决定。

---

季节性分解

存在季节性因素的时间序列数据（如月度数据、季度数据等）可以被分解为趋势因子、季节性因子和随机因子。

趋势因子(trend component)能捕捉到长期变化；季节性因子（seasonal component)能捕捉到一年内的周期性变化；而随机（误差）因子(irregular/error component)则能捕捉到那些不能被趋势或季节效应解释的变化。

而我们可以通过相加模型或相乘模型来分解数据。

相加模型 \(Y_t=Trend_t+Seasonal_t+Irregular_t\)

相乘模型 \(Y_t=Trend_t\times Seasonal_t \times Irregular_t\)

很多时候，相乘模型比相加模型更现实一些。

将时序分解为趋势项、季节项和随机项的常用方法是用LOESS光滑做季节性分解。

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

其中，ts是将要分解的时序，参数s.window控制季节效应变化的速度，t.window控制趋势项变化的速度。较小的值意味着更快的变化速度。令s.window="periodic"可使得季节效应在各年间都一样。这一函数中，参数ts和s.window是必须提供的。

虽然stl()函数只能处理相加模型，但相乘模型总可以通过对数变换转换成相加模型：

\[log(Y_t)=log(Trend_t \times Seasonal_t \times Irregular_t)=log(Trend_t)+log(Seasonal_t)+log(Irregular_t)\]

以自带的AirPassengers数据集为例子。

> plot(AirPassengers)
> lAirPassengers <- log(AirPassengers)
> plot(lAirPassengers, ylab="log(AirPassengers)")
> fit <- stl(lAirPassengers, s.window="period")
> plot(fit)
> fit$time.seriesseasonal    trend     remainder
Jan 1949 -0.09164042 4.829389 -0.0192493585
Feb 1949 -0.11402828 4.830368  0.0543447685
...output omitted ...
> exp(fit$time.series)seasonal    trend remainder
Jan 1949 0.9124332 125.1344 0.9809347
Feb 1949 0.8922327 125.2571 1.0558486
Mar 1949 1.0159924 125.3798 1.0362293
Apr 1949 0.9860703 125.6345 1.0412930
May 1949 0.9850875 125.8897 0.9757094
...output omitted ...

得到的图如下。

上图中，序列的趋势为单调增长，季节效应表明夏季乘客数量更多。每个图的y轴尺度不同，因此通过图中右侧的灰色长条来指示量级，即每个长条代表的量级一样。

stl()函数返回的对象中有一项是time.series，它包括每个观测值中的趋势、季节以及随机效应的具体组成。此时，直接用fit$time.series则返回对数变换后的时序，而通过exp(fit$time.series)可将结果转化为原始尺度。

通过R中自带的monthplot()函数和forecast包中的seasonplot()函数可以对季节分解进行可视化。

par(mfrow=c(2,1))
library(forecast)
monthplot(AirPassengers, xlab="", ylab="")
seasonplot(AirPassengers, year.labels="True", main="")

得到下图。

第一幅图是月度图，表示的是每个月份组成的子序列（连接所有1月的点、连接所有2月的点，以此类推），以及每个子序列的平均值。

第二幅图是季节图，这幅图以年份为子序列。

指数预测模型

指数模型是用来预测时序未来值的最常用模型。它们的短期预测能力较好。不同指数模型建模时选用的因子可能不同。

单指数模型(simple/single exponential model)拟合的是只有常数水平项和时间点\(i\)处随机项的时间序列，这时认为时间序列不存在趋势项和季节效应。

双指数模型(double exponential model; 也叫Holt指数平滑， Holt exponential smoothing)拟合的是有水平项和趋势项的时序。

三指数模型(triple exponential model; 也叫Holt-Winters指数平滑， Holt-Winters exponential smoothing)拟合的是有水平项、趋势项以及季节效应的时序。

而使用R中自带的HoltWinters()函数或者forecast包中的ets()函数可以拟合指数模型。但是ets()函数的备选参数更多，因此更实用。

ets(ts, model="zzz")

其中ts是要分析的时序，限定模型的字母有三个。第一个字母代表误差项，第二个字母代表趋势项，第三个字母则代表季节项。可选的字母包括：相加模型(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()函数的便捷包装(convenience wrapper)，函数中有事先默认设定的参数值。

单指数平滑

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

单指数平滑模型假定时序中的观测值可被表示为：

\[Y_t=level+irregular_r\]

在时间点\(Y_{t+1}\)的预测值（一步向前预测, 1-step ahead forecast）可写作

\[Y_{t+1}=c_0 Y_t+c_1 Y_{t-1}+c_2 Y_{t-2}+\cdots\]

其中\(c_i=\alpha (1-\alpha)^i, i=0,1,2,\cdots\)并且\(0\leq \alpha \leq 1\)。

权数\(c_i\)的总和为1，则一步向前预测可看作当前值和全部历史值的加权平均。

式中\(\alpha\)参数控制权数下降的速度，\(\alpha\)越接近于1，则近期观测值的权重越大；反之，\(\alpha\)越接近于0，则历史观测值的权重越大。

为最优化某种拟合标准，\(\alpha\)的实际值一般由计算机选择，常见的拟合标准是真实值和预测值之间的残差平方和。

以自带的康涅狄格州纽黑文地区的年平均气温数据集nhtemp为例子。

> library(forecast)
> fit <- ets(nhtemp,model ="ANN")
> fit
ETS(A,N,N) Call:ets(y = nhtemp, model = "ANN") Smoothing parameters:alpha = 0.1819 Initial states:l = 50.2762 sigma:  1.1455AIC     AICc      BIC
265.9298 266.3584 272.2129
> forecast(fit,1)Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
1972       51.87031 50.40226 53.33835 49.62512 54.11549
> plot(forecast(fit,1), xlab="Year",ylab=expression(paste("Temperature(", degree*F,")",)),main="New Haven Annual Mean Temperature")
> accuracy(fit)ME     RMSE       MAE       MPE   MAPE      MASE
Training set 0.1460657 1.126268 0.8951225 0.2419373 1.7489 0.7512408ACF1
Training set -0.006441923

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

其中forecast()函数用于预测时序未来的\(k\)步，其形式为\(forecast(fit,k)\)。得到的预测值是51.87°F，其95%的置信区间为49.63°F到54.12°F，还给出了80%的置信区间。

accuracy()函数展示了时序预测中最主流的几个准确性度量。下表中\(e_t\)代表第\(t\)个观测值的误差项（随机项），即\((Y_t- \hat{Y_i})\)。

度量标准	简写	定义
平均误差	ME	\(mean(e_t)\)
平均残差平方和的平方根	RMSE	\(sqrt(mean(e_t^2))\)
平均绝对误差	MAE	\(mean(|e_t|)\)
平均百分比误差	MPE	\(mean(100 \times e_t / Y_t)\)
平均绝对百分误差	MAPE	\(mean(|100 \times e_t / Y_t|)\)
平均绝对标准化误差	MASE	\(mean(|q_t|)\)

其中\(q_t=e_t/(1/(T-1)\times sum(|y_t-y_{t-1}|))\)，T是观测值的个数，对\(t=2\)到\(t=T\)求累加。

一般来说，平均误差和平均百分比误差用处不大，因为正向和负向的误差会抵消掉。平均绝对百分误差给出了误差在真实值中的占比，它没有单位，因此可以用于比较不同时序间的预测准确性；但它同时假定测量尺度中存在一个真实为零的点。RMSE是相对最有名、最常用的。

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

Holt指数平滑可以对有水平项和趋势项（斜率）的时序进行拟合。时刻t的观测值可表示为：

\[Y_t=level+slope \times t +irregular_t\]

平滑参数\(\alpha\)(alpha)控制水平项的指数型下降，beta控制斜率的指数型下降。同样，两个参数的有效范围都是[0,1]，参数取值越大意味着越近的观测值的权重越大。

Holt-Winters指数光滑可用来拟合有水平项、趋势项以及季节项的时间序列。此时，模型可表示为：

\[Y_t=level+slope \times t +s_t+irregular_t\]

其中\(s_t\)代表时刻\(t\)的季节效应。除alpha和beta参数外，gamma光滑参数控制季节项的指数下降。gamma参数的取值范围同样是[0,1]，gamma值越大，意味着越近的观测值的季节效应权重越大。

仍然使用之前的航空乘客数据集作为例子。

> library(forecast)
> fit <- ets(log(AirPassengers), model="AAA")
> fit
ETS(A,A,A) Call:ets(y = log(AirPassengers), model = "AAA") Smoothing parameters:alpha = 0.6975 beta  = 0.0031 gamma = 1e-04 Initial states:l = 4.7925 b = 0.0111 s = -0.1045 -0.2206 -0.0787 0.0562 0.2049 0.21490.1146 -0.0081 -0.0059 0.0225 -0.1113 -0.0841sigma:  0.0383AIC      AICc       BIC
-207.1694 -202.3123 -156.6826
> accuracy(fit)ME       RMSE        MAE         MPE      MAPE      MASE
Training set -0.001830684 0.03606976 0.02770885 -0.03435608 0.5079142 0.2289192ACF1
Training set 0.05590461
> pred <- forecast(fit,5)
> predPoint Forecast    Lo 80    Hi 80    Lo 95    Hi 95
Jan 1961       6.109335 6.060306 6.158365 6.034351 6.184319
Feb 1961       6.092542 6.032679 6.152405 6.000989 6.184094
Mar 1961       6.236626 6.167535 6.305718 6.130960 6.342292
Apr 1961       6.218531 6.141239 6.295823 6.100323 6.336738
May 1961       6.226734 6.141971 6.311498 6.097100 6.356369
> plot(pred, main="Forecast for Air Travel", ylab="Log(AirPassengers)", xlab="Time")
> 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", "Lo 95", "Hi 80", "Hi 95")
> pmean    Lo 80    Lo 95    Hi 80    Hi 95
Jan 1961 450.0395 428.5065 417.5279 472.6544 485.0826
Feb 1961 442.5448 416.8301 403.8280 469.8459 484.9735
Mar 1961 511.1312 477.0088 459.8775 547.6945 568.0971
Apr 1961 501.9652 464.6289 446.0019 542.3017 564.9506
May 1961 506.1001 464.9691 444.5667 550.8694 576.1504

可以看到三个光滑参数。趋势项为0.0031，它小意味着近期观测值的斜率不需要更新。

接着预测了接下来五个月的乘客量。图如下。

此时的预测基于对数变换后的数值，因此要通过幂变换得到预测的乘客量。pred$mean包含了点估计值，其他两个包含了80%和95%置信区间的下界以及上届。cbind()函数用于整合所有结果。

ets()函数和自动预测

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

时序预测一般假定序列的长期趋势是一直向上的，而一个抑制项则使得趋势项在一段时间内靠近一条水平渐近线。在很多问题中，一个有抑制项的模型往往更符合实际情况。

也可以直接使用ets()函数自动选取对原始数据拟合优度最高的模型。

ARIMA预测模型

ARIMA模型主要用于拟合具有平稳性（或可以被转换为平稳序列）的时间序列。在平稳的时序中，序列的统计性质并不会随着时间的推移而改变。

一般来说，拟合ARIMA模型前都需要变换序列的值以保证方差为常数。之前用到的对数变换就是一种常用的变换方法，另外常见的还有Box-Cox变换。

由于一般假定平稳性时序有常数均值，这样的序列中肯定不含有趋势项。非平稳的时序可以通过差分来转换为平稳性序列。具体来说，差分就是将时序中的每一个观测值\(Y_t\)都替换为\(Y_{t-1}-Y_t\)。注意，对序列的一次差分可以移除序列中的线性趋势，二次差分移除二次项趋势，三次差分移除三次项趋势。在实际操作中，对序列进行两次以上的差分通常都是不必要的。

我们可通过diff()函数对序列进行差分，即diff(ts,differences=d)，其中d即对序列ts的差分次数，默认值为\(d=1\)。forecast包中的ndiffs()函数可以帮助我们找到最优的\(d\)值，语句为ndiffs(ts)。

平稳性一般可以通过时序图直观判断。如果方差不是常数，则需要对数据做变换；如果数据中存在趋势项，则需要对其进行差分。

也可以通过ADF(Augmented Dickey-Fuller)统计检验来验证平稳性假定。R中的tseries包的adf.test()可以用来做ADF检验，语句为adf.test(ts)。如果结果显著，则认为序列满足平稳性。

相关概念

滞后阶数(lag）即我们向后追溯的观测值的数量。

时序可以通过lag(ts,k)函数变成\(k\)阶滞后，其中ts指代目标序列，\(k\)代表滞后项阶数。

自相关(autocorrelation)度量时序中各个观测值之间的相关性。\(AC_k\)即一系列观测值(\(Y_t\))和\(k\)时期之前的观测值(\(Y_{t-k}\))之间的相关性。这样，\(AC_1\)就是一阶滞后序列和0阶滞后序列间的相关性……这些相关性(\(AC_1,AC_2,\cdots,AC_k\))构成的图即自相关函数图（AutoCorrelation Function plot, ACF图）。ACF图可用于为ARIMA模型选择合适的参数，并评估最终模型的拟合效果。

stats程序包中的acf()函数或者forecast包中的Acf()函数可以生成ACF图。

偏自相关(partial autocorrelation)即当序列\(Y_t\)和\(Y_{t-k}\)之间的所有值(\(Y_{t-1},Y_{t-2},\cdots,Y_{t-k+1}\))带来的效应都被移除后，两个序列间的相关性。我们也可以对不同的\(k\)值画出偏自相关图。

stats程序包中的pacf()函数和forecast包中的Pacf()函数都可以用来画PACF图。PACF图也可以用来找到ARIMA模型中最适宜的参数。

ARMA和ARIMA模型

在一个\(p\)阶自回归模型中，序列中的每一个值都可以用它之前\(p\)个值的线性组合来表示：

\[AR(p):Y_t=\mu+\beta_1 Y_{t-1}+\beta_2 Y_{t-2}+\cdots+\beta_p Y_{t-p}+\epsilon_t\]

其中\(Y_t\)是时序中的任一观测值，\(\mu\)是序列的均值，\(\beta\)是权重，\(\epsilon_t\)是随机扰动。

在一个\(q\)阶移动平均模型中，时序中的每个值都可以用之前的\(q\)个残差的线性组合来表示：

\[MA(q):Y_t=\mu-\theta_1 \epsilon_{t-1}-\theta_2 \epsilon_{t-2}-\cdots- \theta_q \epsilon_{t-q}+\epsilon_t\]

其中\(\epsilon\)是预测的残差，\(\theta\)是权重。这里的移动平均和之前说的简单移动平均不是一个概念。

这两种方法的混合即\(ARMA(p,q)\)模型，即：

\[Y_t=\mu+\beta_1 Y_{t-1}+\beta_2 Y_{t-2}+\cdots+\beta_p Y_{t-p}-\theta_1 \epsilon_{t-1}-\theta_2 \epsilon_{t-2}-\cdots- \theta_q \epsilon_{t-q}+\epsilon_t\]

---

\(ARIMA(p,d,q)\)模型意味着时序被差分了\(d\)次，且序列中的每个观测值都是用过去的\(p\)个观测值和\(q\)个残差的线性组合表示的。

建立ARIMA模型的步骤包括：

• 确保时序是平稳的；
• 找到一个（或几个）合理的模型（即选定可能的\(p\)值和\(q\)值）；
• 拟合模型；
• 从统计假设和预测准确性等角度评估模型；
• 预测。

---

选择ARIMA模型的方法：

模型	ACF	PACF
ARIMA(p,d,0)	逐渐减小到零	在p阶后减小到零
ARIMA(0,d,q)	q阶后减小到零	逐渐减小到零
ARIMA(p,d,q)	逐渐减小到零	逐渐减小到零

---

拟合模型

可以用arima()函数拟合一个ARIMA模型，其表达式为arima(ts,order=c(p,d,q))。

函数可以返回移动平均项的系数以及模型的AIC值。可以比较AIC值来得到最合理的模型，比较的准则是AIC值越小越好。另外使用accuracy()函数来度量准确性也可以帮助选择模型。

---

模型评价

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

可以使用下面的代码进行检验。

> qqnorm(fit$residuals)
> qqline(fit$residuals)
> Box.test(fit$residuals, type="Ljung-Box")Box-Ljung testdata:  fit$residuals
X-squared = 1.3711, df = 1, p-value = 0.2416

qqnorm()和qqline()函数输出的是Q-Q图。如果数据满足正态分布，则数据中的点会落在图中的线上。

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

---

预测

如果模型残差不满足正态性假设或零自相关系数假设，则需要调整模型、增加参数或改变差分次数。

使用forecast包中的forecast()函数可以实现对接下来几年的预测。

使用plot()函数则可以画出预测图。图中黑色的点为预测点的点估计，浅灰色和深灰色区域分别代表80%和95%的置信区间。

ARIMA的自动预测

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

【参考】

[1] 《R语言实战》（第2版）