时间序列分析（4）| 滞后算子与季节性模型

本篇内容分为两部分：滞后算子和季节性模型，其中介绍前者既回顾了之前的内容，也为介绍后者做了铺垫。

1 滞后算子

时间序列模型中的自回归项和移动平均项都可以使用滞后算子（Lag Operator）进行方便地表示。如可以被表示为。

• 对于常数而言，它的任意滞后项都还是它本身，即。

• 分配律

滞后算子满足分配律：。

• 结合律

滞后算子满足结合律：

根据以上性质，对于ARMA(, )，有

进而有

上式是ARMA(, )模型的一个特解，它的展开式是一个MA()过程。

滞后算子有如下几个常用的公式：

当时，

当时，

当时，

2 季节性模型

2.1 示例数据

许多时间序列数据都带有一定的季节性。USAccDeaths数据集是R中系统自带的数据，它储存的是美国1973-1978年逐月的意外死亡人口数。

data <- c(USAccDeaths)
plot(USAccDeaths, type = "l")
points(USAccDeaths, pch = 21, col = "red", bg = "white")

从图中可以看出该序列存在明显的周期性。下面再绘制出它的ACF和PACF图象：

acf(data, lag.max = 48)
pacf(data, lag.max = 48)

从ACF图象中可以看出，相差6六个月的月份之间的意外死亡人口数呈现明显负相关（），而相差12个月则呈现明显正相关（），这也表明数据序列的周期性可能与季节因子有关。

由于许多时间序列的周期性都与季节因子相关，因此周期性和季节性两个概念在时间序列分析中常常混用。

2.2 普通季节模型

考虑如下两个模型：

model.1:

model.2:

对于model.1，它的自相关系数的特点是：当为周期的整数倍时，，其余为0，也就是呈现出间断指数递减的态势。假定、，使用arima.sim()函数模拟该序列数据，并绘制序列走势图和ACF图象：

y1 <- arima.sim(model = list(order = c(4,0,0),
                             ar = c(0,0,0,0.5)),
                n = 200)
plot(y1, type = "l")
acf(y1)

对于model.2，它的自相关系数在时存在一个峰值，而在其他处均为0。同样，假定、模拟该序列数据：

y2 <- arima.sim(model = list(order = c(0,0,4),
                             ma = c(0,0,0,0.5)),
                n = 200)
plot(y2, type = "l")
acf(y2)

从数据集USAccDeaths所对应的ACF图象上看，除了在12的整数倍呈现峰值外，在单个周期内自相关系数也呈现递减趋势，而非直接截尾为0，这表明该序列数据在季节性之外还存在非季节性。现实中的数据也往往不会如model.1和model.2所刻画的那么理想。

可以使用以下模型尝试对USAccDeaths数据集进行拟合：

model.3:
model.4:

由于上述模型中有许多滞后值的系数为0，可以使用fixed参数进行设置：

(model.3 <- arima(data, order = c(12,0,0),
                  fixed = c(NA, rep(0,10), NA, NA)))
## arima(x = data, order = c(12, 0, 0), fixed = c(NA, rep(0, 10), NA, NA))
## 
## Coefficients:
##          ar1  ar2  ar3  ar4  ar5  ar6  ar7  ar8  ar9  ar10  ar11    ar12
##       0.3442    0    0    0    0    0    0    0    0     0     0  0.6300
## s.e.  0.0648    0    0    0    0    0    0    0    0     0     0  0.0662
##       intercept
##       9176.4060
## s.e.   685.9749
## 
## sigma^2 estimated as 214051:  log likelihood = -548.73,  aic = 1105.47

(model.4 <- arima(USAccDeaths, order = c(1,0,12),
                  fixed = c(NA, NA, rep(0,10), NA, NA)))
## arima(x = USAccDeaths, order = c(1, 0, 12), fixed = c(NA, NA, rep(0, 10), NA, 
##     NA))
## 
## Coefficients:
##          ar1     ma1  ma2  ma3  ma4  ma5  ma6  ma7  ma8  ma9  ma10  ma11
##       0.7241  0.0028    0    0    0    0    0    0    0    0     0     0
## s.e.  0.0946  0.1432    0    0    0    0    0    0    0    0     0     0
##         ma12  intercept
##       0.7203  8953.6706
## s.e.  0.1365   328.2761
## 
## sigma^2 estimated as 235875:  log likelihood = -552.24,  aic = 1114.47

从AIC和模型简洁性的角度上看，model.3都优于model.4。下面再比较二者的预测值与实际值的差距：

fit.3 <- data - model.3$residuals
fit.4 <- data - model.4$residuals

plot(data, type = "l")
lines(1:72, fit.3, col = "red") # model.3
lines(1:72, fit.4, col = "blue") # model.4

2.3 乘法季节模型

假设示例数据中的非季节性由AR(1)过程控制，该过程使用滞后算子可以书写成；假设季节性由model.1所示的过程控制，其中，该过程使用滞后算子可以书写成。

同理，若示例数据中的非季节性和季节性分别由MA过程控制，那么可以分别写作和。

在上节所示的普通季节模型中，其思想是将控制季节性和非季节性的过程进行相加。如果要考虑二者的相互影响，可以使用乘法季节模型，如：

model.5:

或者，

model.6:

上述两式展开后可以得到ARMA模型的一般形式，并且其阶数高于周期12，但其需要估计的系数并没有增加。

实际运用中并没有必要将上述两式展开，它们可以写成ARIMA(, , )(, , )的形式，其中、、为非季节性过程的参数，、、为非季节性过程的参数，为周期。前面已经提过，对于ARMA模型，、。具体地，model.5可以写作ARIMA(1,0,0)(1,0,0)，model.6可以写作ARIMA(1,0,1)(0,0,1)。

上述形式可以直接使用arima()函数进行估计：

(model.5 <- arima(USAccDeaths, order = c(1,0,0),
                  seasonal = list(order = c(1,0,0), period = 12)))
## arima(x = USAccDeaths, order = c(1, 0, 0), seasonal = list(order = c(1, 0, 0), 
##     period = 12))
## 
## Coefficients:
##          ar1    sar1  intercept
##       0.7579  0.8501  9214.7102
## s.e.  0.0768  0.0498   676.1701
## 
## sigma^2 estimated as 127203:  log likelihood = -533.45,  aic = 1074.89

(model.6 <- arima(USAccDeaths, order = c(1,0,1),
                  seasonal = list(order = c(0,0,1), period = 12)))
## arima(x = USAccDeaths, order = c(1, 0, 1), seasonal = list(order = c(0, 0, 1), 
##     period = 12))
## 
## Coefficients:
##          ar1     ma1    sma1  intercept
##       0.6849  0.0907  0.7060  8947.0658
## s.e.  0.1065  0.1332  0.1287   311.2709
## 
## sigma^2 estimated as 235882:  log likelihood = -552,  aic = 1114.01

从结果上看，model.5要优于model.6。下面对比model.5和普通季节模型中的优胜者model.3之间的预测值差异：

fit.5 <- data - model.5$residuals

plot(data, type = "l")
lines(1:72, fit.3, col = "red", lty = 2) # model.3
lines(1:72, fit.5, col = "red") # model.5

• 上图中，黑色实线为真实数据，红色虚线为使用普通季节模型（model.3）预测的结果，红色实线表示使用乘法季节模型（model.5）预测的结果。可以看出后者的拟合度更高。