Holt 两参数指数平滑适用于对含有线性趋势的序列进行修匀。他的基本思想是假定序列有一个比较固定的线性趋势——每期都递增或递减 \(r\)，那么第 \(t\) 期的估计值就应该等于第 \(t-1\) 期的观察值加上每期固定的趋势变动值，即
\[ \hat{x}_{t} = x_{t-1} + r \]
但是由于随机因素的影响，每期的递增或递减值不会恒定为 \(r\)，它会随着时间上下波动，所以趋势序列其实是一个随机序列 \(\{r_{t}\}\)，因而
\[ \hat{x}_{t} = x_{t-1} + r_{t-1} \]
考虑用第t期的观察值和第 \(t\) 期的估计值加权平均数作为第 \(t\) 期的修匀值：
\[ \tilde{x}_{t}=\alpha x_{t}+(1-\alpha)\hat{x}_{t} \]
\[ \Rightarrow \tilde{x}_{t}=\alpha x_{t}+(1-\alpha)(x_{t-1}+r_{t-1}), \quad 0< \alpha <1 \]
因为趋势序列 \(\{r_{t}\}\) 也是一个随机序列，为了使修匀序列 \(\tilde{x}\) 更平滑，我们对 \(\{r_{t}\}\) 也进行一次修匀处理：
\[ r_{t}=\beta(\tilde{x}_{t}-\tilde{x}_{t-1})+(1-\beta)r_{t-1}, \quad 0<\beta<1 \]
这就是 Holt 两参数指数平滑的构造思想，它的平滑公式为：
\[ \left\{\begin{matrix} \tilde{x}_{t}=\alpha x_{t}+(1-\alpha)(x_{t-1}+r_{t-1}) \\ r_{t}=\beta(\tilde{x}_{t}-\tilde{x}_{t-1})+(1-\beta)r_{t-1} \end{matrix}\right. \]
式中，\(\alpha\)，\(\beta\) 为两个平滑系数，也称平滑参数，它们满足 \(0<\alpha,\beta<1\)。
此外，还需要确定两个序列的初始值：
平滑序列的初始值 \(\tilde{x}_{0}\)。最简单的方法是指定 \(\tilde{x}_{t}=x_{1}\)；
趋势序列的初始值 \(r_{0}\)。最简单的方法是指定一个区间长度n，用这段区间的平均趋势作为趋势初始值：
\[ r_{0}=\frac{x_{n+1}-x_{1}}{n} \]
假定最后一期的修匀值为 \(\tilde{x}_{T}\)，那么使用 Holt 两参数指数平滑方法，向前 \(l\) 期的预测值为：\(\tilde{x}_{T}(l)=\tilde{x}_{T} + l \cdot r_{T}\)。
数据：1990-2020年安徽省老年人口数
来源：安徽统计年鉴
# 读取数据
data <- read.csv("data/安徽省老年人口数.csv")
x <- ts(data$老年人口数, start = 1990) # 转换为时间序列数据（本质是一个原子向量）
head(data)
## 时间 老年人口数
## 1 1990 306.2601
## 2 1991 340.7778
## 3 1992 352.8546
## 4 1993 368.1400
## 5 1994 376.8960
## 6 1995 398.6179# 画出时序图
plot(
x,
pch = "*", type = "o", lwd = 2, col = "blue",
main = "", xlab = "时间(年)", ylab = "老年人口数（万人）"
)
# 进行Holt两参数平滑
x_fit <- HoltWinters(x, gamma=FALSE)
x_fit # 查看拟合结果
## Holt-Winters exponential smoothing with trend and without seasonal component.
##
## Call:
## HoltWinters(x = x, gamma = FALSE)
##
## Smoothing parameters:
## alpha: 0.9727709
## beta : 0.06830705
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 915.84161
## b 24.45372
x_fit$fitted # 查看拟合值
## Time Series:
## Start = 1992
## End = 2020
## Frequency = 1
## xhat level trend
## 1992 375.2955 340.7778 34.51770
## 1993 386.4922 353.4656 33.02657
## 1994 400.4468 368.6397 31.80712
## 1995 407.7795 377.5373 30.24223
## 1996 428.5008 398.8674 29.63347
## 1997 415.2189 388.3294 26.88948
## 1998 428.4367 402.4214 26.01532
## 1999 448.0988 422.4897 25.60909
## 2000 472.1127 446.6056 25.50710
## 2001 487.5872 462.7216 24.86562
## 2002 521.5780 496.1289 25.44908
## 2003 547.7150 522.2220 25.49307
## 2004 592.6048 565.8715 26.73328
## 2005 605.4117 579.5689 25.84284
## 2006 643.8252 617.1786 26.64660
## 2007 646.5187 621.4036 25.11505
## 2008 681.3306 655.5955 25.73506
## 2009 709.8943 683.9784 25.91592
## 2010 726.3315 701.0217 25.30986
## 2011 630.1252 612.5850 17.54017
## 2012 700.4513 679.5360 20.91528
## 2013 745.1637 722.7269 22.43686
## 2014 760.2079 738.2437 21.96417
## 2015 732.3932 713.6119 18.78134
## 2016 738.9661 720.9654 18.00073
## 2017 761.6370 743.3376 18.29934
## 2018 793.1419 773.9982 19.14370
## 2019 840.3746 819.4349 20.93970
## 2020 909.5291 885.5066 24.02253# 预测
library(zoo)
library(forecast)
x_fore <- forecast(x_fit, h=4) # 4期预测
x_fore # 输出预测值，包含80%和95%的CI
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2021 940.2953 900.8494 979.7412 879.9680 1000.623
## 2022 964.7490 907.8597 1021.6384 877.7443 1051.754
## 2023 989.2028 917.5189 1060.8866 879.5718 1098.834
## 2024 1013.6565 928.3556 1098.9574 883.2000 1144.113
plot(
x_fore,
pch = "*", type = "o", lwd = 2, col = "red",
main = "", xlab = "时间(年)", ylab = "老年人口数（万人）"
)
lines(x_fit$fitted[ , 1], lwd = 2, col = "blue") # 添加拟合线
legend(
"topleft",
c("老年人口实际值", "指数平滑拟合值"),
col = c("red", "blue"),
lwd = c(2, 2),
pch = c("*", ""),
bty = "n"
) # 设置图例
# 保存数据
write.csv(x_fit$fitted, "data/Holt-fit.csv") # 保存拟合结果
write.csv(x_fore, "data/Holt-forecast.csv") # 保存预测结果
print("Work done!")
## [1] "Work done!"白噪声的三个性质：
一般不是很重要，因为即使均值为非零常数，也可以对时间序列数据中心化实现。
序列各项之间没有任何相关关系，序列在进行完全无序的随机波动。
序列的每个变量的方差都相等。若方差非齐性，则参数的显著性检验失去意义，模型拟合的精度会降低。
res = x - x_fit$fitted[ , 1] # 求出残差
# 时序图检验
plot(
res,
type = "o", pch = 16, col = "blue", lwd = 1.5,
main = "", xlab = "时间(年)", ylab = "残差"
) # 残差时序图
abline(h = 0, lty = 2) # 辅助线
par(mar = c(5, 5, 4, 3))
acf(res, lag = 30) # 自相关图
pacf(res, lag = 30) # 偏自相关图
ACF 与 PACF 均在两倍标准差范围内，说明残差序列平稳，且应该具有纯随机性。但从图标来判断序列是否是纯随机的具有一定的主观性，我们需要进行客观的统计检验。
Barlett 证明：如果一个时间序列是纯随机的，得到一个 \(n\) 期的观测序列 \(\{x_t, t=1,2,\cdots,n\}\)，那么该序列的延迟非零期的样本自相关系数将近似服从均值为零，方差为序列观察期数倒数的正态分布，即： \[ \hat{\rho_k}\sim N(0,\frac{1}{n}),\quad \forall k \neq 0 \]
于是，可以构造以下统计检验量来检验序列的纯随机性。
\[ \begin{alignat*}{2} & H_0: \rho_1 = \rho_2 = \cdots = \rho_m, \quad \forall m \geqslant 1 \\ & H_1: \exists \rho_k \neq 0, \quad \forall m\geqslant 1,k\leqslant m \end{alignat*} \]
Box 和 Pierce 推导出了 \(Q\) 统计量：
\[ Q=n \sum_{k=1}^{m} \hat{\rho_{k}}^{2} \sim \chi^2(m) \]
式中 \(n\) 为观测期数，\(m\) 为指定滞后期数，\(Q\) 统计量近似服从自由度为 \(m\) 的卡方分布。
\(Q\) 统计量在大样本场合下检验效果较好，但在小样本场合不太精确。Box 和 Ljung 推导出了适用于小样本场合下的 \(LB\) 统计量作为 \(Q\) 统计量的修正：
\[ LB=n(n+2) \sum_{k=1}^{m} (\frac{\hat{\rho_{k}}^{2}}{n-k}) \sim \chi^2(m) \]
# 纯随机性检验
Box.test(res, type = "Box-Pierce", lag = 20) # Q统计量
##
## Box-Pierce test
##
## data: res
## X-squared = 8.0796, df = 20, p-value = 0.9913
Box.test(res, type = "Ljung-Box", lag = 20) # LB统计量，更适用于小样本
##
## Box-Ljung test
##
## data: res
## X-squared = 12.671, df = 20, p-value = 0.891
# 多期滞后的随机性检验
for (i in 1:4) print(Box.test(res, type = "Ljung-Box", lag = 5*i))
##
## Box-Ljung test
##
## data: res
## X-squared = 4.3954, df = 5, p-value = 0.494
##
##
## Box-Ljung test
##
## data: res
## X-squared = 8.5199, df = 10, p-value = 0.5782
##
##
## Box-Ljung test
##
## data: res
## X-squared = 10.669, df = 15, p-value = 0.7757
##
##
## Box-Ljung test
##
## data: res
## X-squared = 12.671, df = 20, p-value = 0.891\(P\) 值均大于给定的显著性水平 \(\alpha\)，所以该序列不能拒绝纯随机性的原假设。但由残差时序图可以看出方差可能非齐，对于无自相关又异方差的序列，暂时没有较好的方法，故不做讨论。