5.4 高效的预测变量
当我们对时间序列数据做回归时，可以构造一些高效的预测变量。
趋势
很多的时间序列存在趋势。当存在简单的线性趋势时，可以直接使用 \(x_{1,t}=t\) 作为预测变量，
\[
y_{t}= \beta_0+\beta_1t+\varepsilon_t,
\]
其中， \(t=1,\dots,T\)。在tslm()函数中，可以利用trend来指定趋势变量。在5.8节中，我们会详细讨论如何对非线性趋势建模。
虚拟变量
目前，我们讨论的每个预测变量都是数值型变量。但是当某个预测变量为分类变量且只有两个取值时（例如，“是”或“否”）应当怎么处理？例如，当你想要预测日销量时，你想把当天是否为 法定节假日 考虑进来。此时，则需要引入一个预测变量，当天为法定节假日时该变量取值为“是”，否则取值为“否”。
在这种情况下，我们可以通过在多元模型中添加“虚拟变量”来进行处理。当虚拟变量的取值为1时，代表“是”；取值为0时代表“否”。虚拟变量通常也被称为“指示变量”。
虚拟变量也可以用来处理数据中的 离群点 。与简单删除离群点不同的是，虚拟变量剔除了离群点对模型的影响。当该观测值是离群点时，虚拟变量取值为1，在其他观测值处，虚拟变量取值均为0。虚拟变量也可以表示特殊事件是否发生。例如，在预测到达巴西的游客时，我们需要考虑2016年里约热内卢夏季奥运会对游客数量的影响。
如果有两个以上的类别，则可以使用多个虚拟变量（需要注意的是，虚拟变量个数应比类别数少1）对变量进行编码。在tslm()函数中，如果你指定某个分类变量为预测变量，tslm()将自动处理，而不需要再手动添加虚拟变量。
季节性虚拟变量
假设我们想要预测日度数据，且想把星期数（周一、周二等等）作为预测变量。我们可以构造如下的虚拟变量。
| \(d_{1,t}\) | \(d_{2,t}\) | \(d_{3,t}\) | \(d_{4,t}\) | \(d_{5,t}\) | \(d_{6,t}\) | |
|---|---|---|---|---|---|---|
| 周一 | 1 | 0 | 0 | 0 | 0 | 0 |
| 周二 | 0 | 1 | 0 | 0 | 0 | 0 |
| 周三 | 0 | 0 | 1 | 0 | 0 | 0 |
| 周四 | 0 | 0 | 0 | 1 | 0 | 0 |
| 周五 | 0 | 0 | 0 | 0 | 1 | 0 |
| 周六 | 0 | 0 | 0 | 0 | 0 | 1 |
| 周日 | 0 | 0 | 0 | 0 | 0 | 0 |
| 周一 | 1 | 0 | 0 | 0 | 0 | 0 |
| 周二 | 0 | 1 | 0 | 0 | 0 | 0 |
| 周三 | 0 | 0 | 1 | 0 | 0 | 0 |
| 周四 | 0 | 0 | 0 | 1 | 0 | 0 |
| 周五 | 0 | 0 | 0 | 0 | 1 | 0 |
| ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
值得注意的是，编码七个类别只需要六个虚拟变量。当所有虚拟变量都取0时，即可表示第七类（上例中的周日）。
许多初学者会给第七类添加第七个虚拟变量，这会导致模型预测变量之间出现完全共线性，一般被称为“虚拟变量陷阱”，它会导致回归失败。因此如果定性变量有 \(m\) 个类别，只需要引入 \(m-1\) 个虚拟变量。例如对于季度数据，需要引入3个虚拟变量；对于月度数据，需要引入11个虚拟变量；对于日度数据，需要引入6个虚拟变量。
与虚拟变量相关的每个系数的解释是 该类别相对于忽略的类别对模型的影响程度 。在上例中，“周一”的系数 \(d_{1,t}\) 即是与“周日”相比，“周一”对被预测变量的影响。以下是解释估计的虚拟变量系数的示例，其显示了澳大利亚啤酒产量的季度季节性。
在tslm()函数中，如果你设置预测变量为season，tslm()将自动处理，而不需要再手动添加虚拟变量。
示例：澳大利亚啤酒季度产量
图5.14为澳大利亚啤酒的季度产量。
beer2 <- window(ausbeer, start=1992)
autoplot(beer2) + xlab("年份") + ylab("万升")+
theme(text = element_text(family = "STHeiti"))
图 5.14: 澳大利亚啤酒的季度产量。
我们想要预测未来的啤酒产量。因此可以构造一个具有线性趋势和季度虚拟变量的回归模型： \[ y_{t} = \beta_{0} + \beta_{1} t + \beta_{2}d_{2,t} + \beta_3 d_{3,t} + \beta_4 d_{4,t} + \varepsilon_{t}, \] 其中，当 \(t\) 为第 \(i\) 个季度时 \(d_{i,t}=1\) ，否则 \(d_{i,t}=0\) 。模型中没有变量表示第一季度，因为其他季度虚拟变量的系数是这些季度与第一季度之间差异的度量。
fit.beer <- tslm(beer2 ~ trend + season)
summary(fit.beer)
#>
#> Call:
#> tslm(formula = beer2 ~ trend + season)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -42.90 -7.60 -0.46 7.99 21.79
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 441.8004 3.7335 118.33 < 2e-16 ***
#> trend -0.3403 0.0666 -5.11 2.7e-06 ***
#> season2 -34.6597 3.9683 -8.73 9.1e-13 ***
#> season3 -17.8216 4.0225 -4.43 3.4e-05 ***
#> season4 72.7964 4.0230 18.09 < 2e-16 ***
#> ---
#> Signif. codes:
#> 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 12.2 on 69 degrees of freedom
#> Multiple R-squared: 0.924, Adjusted R-squared: 0.92
#> F-statistic: 211 on 4 and 69 DF, p-value: <2e-16需要注意的是，“趋势”和“季节”都不是R工作空间中的实际对象。当采用这种方式指定时，它们由 tslm() 自动创建。
由上述结果可知，每季度平均下降趋势为-0.34兆升。平均而言，第二季度的产量比第一季度低34.7兆升，第三季度的产量比第一季度低17.8兆升，第四季度的产量比第一季度高72.8兆升。
autoplot(beer2, series="真实值") +
autolayer(fitted(fit.beer), series="拟合值") +
xlab("年份") + ylab("万升") +
ggtitle("啤酒的季度产出")+
theme(text = element_text(family = "STHeiti"))+
theme(plot.title = element_text(hjust = 0.5))
图 5.15: 啤酒产出的时间序列图及预测值。
cbind(Data=beer2, Fitted=fitted(fit.beer)) %>%
as.data.frame() %>%
ggplot(aes(x=Data, y=Fitted, colour=as.factor(cycle(beer2)))) +
geom_point() +
ylab("拟合值") + xlab("真实值") +
ggtitle("啤酒的季度产出") +
scale_colour_brewer(palette="Dark2", name="季度") +
geom_abline(intercept=0, slope=1)+
theme(text = element_text(family = "STHeiti"))+
theme(plot.title = element_text(hjust = 0.5))
图 5.16: 实际啤酒产量与预测的啤酒产量的对比。
干预变量
建模时，我们通常需要考虑可能对被预测变量的产生影响的干预因素。例如竞争对手的活动、广告支出、工业行动等等都会对被预测变量产生影响。
当干预因素的影响仅持续一个时期时，我们可以使用“尖峰”变量来描述。尖峰变量的处理方法和处理离群点非常相似，也是构造一个虚拟变量，在干预因素作用期间取值1，在其他地方取0。
干预因素的影响还可能是长期或永久的。如果干预因素导致水平偏移（即序列的值从干预时间点之后突然且永久地改变），那么我们使用“阶梯”变量。阶梯变量在干预产生之前取值为0，从干预产生之后取值为1。
干预因素的另一种长远影响是斜率的变化。此时需采取分段处理，在干预因素产生影响前后斜率是不同的，因此模型是非线性的。我们将会在5.8节详细讨论。
交易日
一个月的交易日数可能会有很大差异，并会对销售数据产生重大影响。为此，可以将每个月的交易日数作为预测变量。
对于月度或季度数据，bizdays()函数可以计算出每个时期内的交易日数。
我们可以引入七个解释变量，每个解释变量定义如下： \[\begin{align*} x_{1} &= \text{当月中周一的数目；}\\ x_{2} &= \text{当月中周二的数目；}\\ & \vdots \\ x_{7} &= \text{当月中周日的数目。} \end{align*}\]
分布滞后
通常情况下，把广告支出作为解释变量会十分有效。但是，广告效应往往会具有滞后性。因此，我们可以使用如下变量： \[\begin{align*} x_{1} &= \text{一个月前的广告支出；}\\ x_{2} &= \text{两个月前的广告支出；}\\ & \vdots \\ x_{m} &= \text{ $m$ 个月前的广告支出。} \end{align*}\]
一般情况下，系数随着滞后阶数的增加而减小。
复活节
复活节与其他大多数的假期不同，因为它不是每年在同一天举行，并且其影响可持续一段时间。在这种情况下，在复活节特定的时间段内，虚拟变量取值为1，在其他时间段内取值为0。
当数据为月度数据时，若复活节在三月份，那么虚拟变量在三月份时取1；同样的，若复活节在四月份时，虚拟变量在四月份时取1。
easter()函数可以计算复活节所代表的虚拟变量。
傅立叶级数
对于季节性虚拟变量，尤其是长季节周期，通常可以采用傅里叶级数。让·巴普蒂斯·约瑟夫·傅里叶是一位出生于18世纪的法国数学家。他表明一定频率的一系列正弦和余弦项可以逼近任何周期函数。我们可以把傅里叶级数用于季节模式中。
当序列的季节周期为 \(m\) 时，其傅里叶级数的前几项为： \[ x_{1,t} = \sin\left(\textstyle\frac{2\pi t}{m}\right), x_{2,t} = \cos\left(\textstyle\frac{2\pi t}{m}\right), x_{3,t} = \sin\left(\textstyle\frac{4\pi t}{m}\right), \] \[ x_{4,t} = \cos\left(\textstyle\frac{4\pi t}{m}\right), x_{5,t} = \sin\left(\textstyle\frac{6\pi t}{m}\right), x_{6,t} = \cos\left(\textstyle\frac{6\pi t}{m}\right), \] 如果数据中存在月度季节性，那么我们使用这些预测变量中的前11个，我们将得到与使用11个虚拟变量完全相同的预测。
当采用傅里叶级数时，尤其当 \(m\) 值很大时，我们通常可以用较少的预测变量得到与采用虚拟变量相同的预测结果。例如，周度数据中 \(m\approx52\) ，因此这对于周度数据非常有效。对于短季节周期数据（例如，季度数据），使用傅里叶级数相比于季节性虚拟变量几乎没有优势。
傅里叶级数可通过 fourier()函数来构造。例如，澳大利亚啤酒产出数据可通过以下方式建模。
fourier.beer <- tslm(beer2 ~ trend + fourier(beer2, K=2))
summary(fourier.beer)
#>
#> Call:
#> tslm(formula = beer2 ~ trend + fourier(beer2, K = 2))
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -42.90 -7.60 -0.46 7.99 21.79
#>
#> Coefficients:
#> Estimate Std. Error t value
#> (Intercept) 446.8792 2.8732 155.53
#> trend -0.3403 0.0666 -5.11
#> fourier(beer2, K = 2)S1-4 8.9108 2.0112 4.43
#> fourier(beer2, K = 2)C1-4 53.7281 2.0112 26.71
#> fourier(beer2, K = 2)C2-4 13.9896 1.4226 9.83
#> Pr(>|t|)
#> (Intercept) < 2e-16 ***
#> trend 2.7e-06 ***
#> fourier(beer2, K = 2)S1-4 3.4e-05 ***
#> fourier(beer2, K = 2)C1-4 < 2e-16 ***
#> fourier(beer2, K = 2)C2-4 9.3e-15 ***
#> ---
#> Signif. codes:
#> 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 12.2 on 69 degrees of freedom
#> Multiple R-squared: 0.924, Adjusted R-squared: 0.92
#> F-statistic: 211 on 4 and 69 DF, p-value: <2e-16fourier()函数的第一个参数是季节性周期的长度 \(m\) ；第二个参数K用来确定包含多少正弦项和余弦项，K的最大值为 \(K=m/2\) 。当K取最大值时，预测结果与使用季节性虚拟变量时获得的预测结果相同。
如果仅使用前两个傅立叶项（ \(x_{1,t}\) 和 \(x_{2,t}\) ），此时季节性模式将遵循简单的正弦波。由于连续的傅里叶项表示前两个傅立叶项的谐波，因此包含傅里叶项的回归模型通常称为谐波回归。