前言
作为一种常用而有效的时间序列预测模型，指数平滑法（全称ExponenTial Smoothing，亦称作Error, Trend, Seasonality，简称ETS）模型所涉及的参数和算法较为简单，可以直接在数据库中通过SQL实现预测值的实时、批量输出；相比其他需要特定软件环境的复杂模型来说，ETS模型能够大大降低部署的难度和维护的复杂性。指数加权移动均值（exponentially weighted moving average, ewma）的计算则是ETS模型的实现基础，并且还可以用于计算统计指标、构建非时间序列模型的输入特征等。本文将介绍指数加权移动均值的计算方法，并提供其SQL和R语言实现代码；我将在另一篇文章中介绍如何在SQL中实现ETS模型预测值的输出。
计算方法
以\(x\)表示一条时间序列、\(t\)表示序列的时间点序号、\(w\)表示一个给定的权重、再以\(x_0\)表示一个给定的序列初始值，那么\(x\)在\(t\)时间点对应的指数加权移动平均数\(\bar{x}_t\)可以表示为：
\[ \begin{cases} \begin{align} \bar{x}_1 &= (1 - w) x_0 + w x_1, \quad &t = 1 \\ \bar{x}_{t} &= (1 - w) \bar{x}_{t-1} + w x_t, \quad &t = 2, ..., T \end{align} \end{cases} \quad (1) \]
若是递归地对算式右侧所含的指数加权移动平均项\((\bar{x}_{t-1}, \bar{x}_{t-2}, ...)\)进行分解和代入，那么\(\bar{x}_t\)也可以表示为：
\[ \bar{x}_{t} = (1 - w)^t x_0 + \sum_{i = 0}^{t - 1} w (1 - w)^i x_{t - i}, \quad t = 1, ..., T \quad (2) \]
假设\(x_0 = 0\)，那么\(\bar{x}_{t}\)是原序列\(x\)中前\(t\)个时间点的加权平均值，且各项权重能够组成一个以\(w\)为首项、\(1 - w\)为比、\(t\)为项数的等比数列\(a_t = w(1 - w)^{t - 1}\)。根据等比数列的求和公式，可算得\(a_t\)的和为\(S_t = 1 - (1 - w)^t\)；由于\(0 < 1 - w < 1\)，所以\(S_t < 1\)恒成立，且\(t\)越小，\(S_t\)与1的差值越大（例如，当\(1 - w = 0.75\)时，\(S_1 = 0.25\), \(S_2 = 0.44\), \(S_3 = 0.58\)；又例如，当\(1 - w = 0.25\)时，\(S_1 = 0.75\), \(S_2 = 0.94\), \(S_3 = 0.98\)）。这意味着，在按照上式计算序列\(x\)中前\(t\)个时间点的移动均值时，各项权重之和会小于1，导致结果的绝对值偏小（见图1）。为了对此偏差进行校正以得到更加符合直觉的结果，我们可以将\(\bar{x}_{t}\)除以\(S_t\)，以此保证权重之和恒等于1；如果以\(\bar{x}^{'}_{t}\)表示校正后的指数加权移动均值，则有：
\[ \bar{x}^{'}_{t} = \frac{\bar{x}_t}{1 - (1 - w)^t} \quad (3) \]
算式\((1)\)和\((2)\)中的\(x_0\)和\(w\)都属于算法参数，而ETS模型在拟合过程中会通过最小化误差自动选择最优的参数值（注意：ETS模型算法并不会对偏差进行校正）。但在某些应用场景中（例如，将计算结果当作一般回归模型或分类模型的输入特征），我们可能希望跳过不必要的拟合过程，人为设定参数；那么，在这种情况下应该如何决定它们的大小？以下补充一些对指数加权移动平均的理解，这些理解可能对人为设定参数取值有一定的参考作用。
首先，\(x_0\)的引入除了能使算式的逻辑变得严密之外，还能在一定程度上对指数加权移动均值的偏差进行校正。但是，相对于上面提到的校正方法，通过调整\(x_0\)进行校正会比较麻烦。因此，在多数情况下，我们可以直接令\(x_0 = 0\)，并同时用上面介绍的方法进行校正。
在默认\(x_0 = 0\)之后，再回顾算式\((2)\)，可知\(x_{t - i}\)（\(t - i > 0\)）的权重与\(x_t\)的权重之比为\((1 - w)^i\)。数学上可以证明，对于任意\(0 < w < 1\)，均有\((1 - w)^{\frac{1}{w}} < \frac{1}{e} \approx 0.37\)。如果我们把\(0.37\)视作一个临界值，且当\(x_{t-i}\)与\(x_t\)的权重之比小于该临界值时、\(t - i\)至\(1\)的观测值对总体移动均值的影响可以忽略；那么，因为\(i \geq \frac{1}{w}\)时\((1 - w)^i < 0.37\)恒成立，所以我们可以认为指数加权移动平均主要代表的就是最后\(\frac{1}{w}\)个时间点的观测值的集中程度。需要注意的是，这种理解缺乏严谨的理论依据，只能用作经验性的参考(Ng 2018)。
根据上面的理解，在具体应用中我们可以这样设置\(x_0\)和\(w\)的大小：对于前者，通常来说取\(x_0 = 0\)即可；而对于后者，如果业务需求主要考虑过去最近\(k\)个时间点观察值的影响，那么可以取\(w = \frac{1}{k}\)。以统计A企业季度收入的指数加权移动均值为例，假设由该企业在过去8个季度的收入值\((x_1, x_2, ..., x_8)\)（单位：亿元人民币）组成的序列为
\[ x_t = (11.0, 15.0, 22.0, 23.0, 25.0, 30.0, 37.0, 40.0) \]
在令\(x_0 = 0\)的基础上，如果着重考虑最近4个季度的情况，可以取\(w = 0.25\)；那么，校正前的各季度收入的指数加权移动均值\((\bar{x}_1, \bar{x}_2, ..., \bar{x}_8)\)为
\[ \bar{x}_t = (2.8, 5.8, 9.9, 13.1, 16.1, 19.6, 23.9, 28.0) \]
而校正后的结果\((\bar{x}^{'}_1, \bar{x}^{'}_2, ..., \bar{x}^{'}_8)\)为
\[ \bar{x}^{'}_t = (11.0, 13.3, 17.1, 19.2, 21.1, 23.8, 27.6, 31.1) \]
为了更直观地感受\(w\)的大小以及偏差校正对指数加权移动平均结果的影响（不考虑\(x_0\)的影响，并假设\(x_0 = 0\)），我们可以分别计算不同情况下A企业季度收入指数加权移动均值，并将其绘制成图1。其中，蓝色短划线、蓝色点线分别代表校正前、校正后的季度收入指数加权移动均值，而蓝色深浅则代表权重\(w\)的大小（图中只显示0.25、0.50、0.75三个取值的结果）；另外，图中还添加了代表季度收入原始值的黑色实线，以及代表季度收入累计均值的灰色实线。
从图中可以看出这些规律：1. 相对于校正后的结果，校正前的移动加权均值（的绝对值）序列在起始阶段明显偏小，但该偏差会随着\(w\)的增加或时间\(t\)的推移而减小；2. 在\(0 < w < 1\)范围内，指数加权移动均值序列会随着\(w\)的增大而更加接近原始值序列、随着\(w\)的减小而更加接近累计均值序列，且该情况在校正后的结果中尤为明显（事实上，对于校正后的结果而言，当\(w \rightarrow 1\)时会与原始值序列重合、而当\(w \rightarrow 0\)时会与累计均值序列重合）。
算法的SQL实现
考虑到实际应用中，往往需要同时计算多组序列的指数加权移动均值、且不同序列所含的时间点数量可能不一样（但都是连续的），这里我们再引入B企业在过去7个季度的收入序列\((13.0, 19.0, 20.0, 22.0, 26.0, 32.0, 34.0)\)（单位：亿元人民币），并给出利用SQL同时计算企业A和企业B季度收入移动均值的方法。
在关系数据库中（例如PostgreSQL；本节代码在PostgreSQL 14.2中测试通过），我们可以通过以下代码构建一个名为xt的表用于存储前述数据，其中的字段包括group_id（企业/分组标识）、t（时间点序号）、x（季度收入）。
CREATE TEMPORARY TABLE xt (
group_id VARCHAR,
t INT,
x DOUBLE PRECISION
);
INSERT INTO xt
VALUES ('A', 1, 11.0), ('B', 1, 13.0),
('A', 2, 15.0), ('B', 2, 19.0),
('A', 3, 22.0), ('B', 3, 20.0),
('A', 4, 23.0), ('B', 4, 22.0),
('A', 5, 25.0), ('B', 5, 26.0),
('A', 6, 30.0), ('B', 6, 32.0),
('A', 7, 37.0), ('B', 7, 34.0),
('A', 8, 40.0);以下是指数加权移动均值的计算方法和结果；其中的参数w为指数加权的权重、参数x0为序列的初始值、参数correct_bias代表是否校正偏差（设w = 0.25, x0 = 0, correct_bias = TRUE）：
DO $$
DECLARE
x0 DOUBLE PRECISION := 0;
w DOUBLE PRECISION := 0.25;
correct_bias BOOLEAN := TRUE;
BEGIN
-- 基于算式(1)和(3)
-- 需要利用递归查询，可移植性较低（有的数据库可能不支持）
-- 时间点（t）个数较多时，计算速度较快
-- 综合考虑，多数情况下不推荐使用
-- CREATE TEMPORARY TABLE result AS
-- WITH RECURSIVE ta AS (
-- SELECT CAST(NULL AS VARCHAR) AS group_id,
-- CAST(0 AS INT) AS t,
-- CAST(x0 AS DOUBLE PRECISION) AS x
-- UNION ALL
-- SELECT tb.group_id, tb.t AS t, (1 - w) * ta.x + w * tb.x AS x
-- FROM ta
-- INNER JOIN xt AS tb ON (ta.group_id IS NULL OR ta.group_id = tb.group_id) AND ta.t + 1 = tb.t
-- )
-- SELECT group_id, t, x / CASE WHEN correct_bias THEN (1 - (1 - w) ^ t) ELSE 1 END AS x
-- FROM ta
-- WHERE t >= 1
-- ORDER BY group_id, t;
-- 基于算式(2)和(3)
-- 无需利用递归查询，可移植性较高（大部分数据库都支持）
-- 分组（group_id）个数较多时，计算速度较快
-- 综合考虑，多数情况下推荐使用
CREATE TEMPORARY TABLE result AS
SELECT group_id, t, (sum(w * x * (1 - w) ^ i) + (1 - w) ^ t * x0) / CASE WHEN correct_bias THEN (1 - (1 - w) ^ t) ELSE 1 END AS x
FROM (SELECT ta.group_id, ta.t AS t, ta.t - tb.t AS i, tb.x
FROM xt AS ta
INNER JOIN xt AS tb ON ta.group_id = tb.group_id
WHERE ta.t >= tb.t) ta
GROUP BY group_id, t
ORDER BY group_id, t;
END $$;
SELECT * FROM result; group_id | t | x
----------+---+------
A | 1 | 11.0
A | 2 | 13.3
A | 3 | 17.1
A | 4 | 19.2
A | 5 | 21.1
A | 6 | 23.8
A | 7 | 27.6
A | 8 | 31.1
B | 1 | 13.0
B | 2 | 16.4
B | 3 | 18.0
B | 4 | 19.4
B | 5 | 21.6
B | 6 | 24.8
B | 7 | 27.4
(15 rows)算法的R语言实现
若只考虑A企业的季度收入，我们可以将其赋值给一个向量xt：
xt <- c(11.0, 15.0, 22.0, 23.0, 25.0, 30.0, 37.0, 40.0)下面定义的ewma()函数可以用于计算校正后的指数加权移动均值，其中参数x为原始序列数据、参数w为指数加权的权重、参数x0为序列的初始值、参数correct_bias代表是否校正偏差：
# 基于算式(1)和(3)
# 绝大多数情况下计算速度较快，推荐使用
ewma <- function(x, w = 0.25, x0 = 0, correct_bias = TRUE) {
result <- Reduce(function(xt_1, xt) {
(1 - w) * xt_1 + w * xt
}, x = x, accumulate = T, init = x0)[-1]
if (correct_bias) {
t <- seq_along(x)
result <- result / (1 - (1 - w) ^ t)
}
result
}
# 基于算式(2)和(3)
# 绝大多数情况下计算速度较慢，不推荐使用
# ewma <- function(x, w = 0.25, x0 = 0, correct_bias = TRUE) {
# result <- vapply(seq_along(x), function(t) {
# i <- seq_len(t)
# (1 - w) ^ t * x0 + sum(x[i] * w * (1 - w) ^ (t - i))
# }, double(1))
# if (correct_bias) {
# t <- seq_along(x)
# result <- result / (1 - (1 - w) ^ t)
# }
# result
# }以下是利用该函数算出的校正前和校正后的指数加权移动均值：
# 校正前
ewma(xt, w = 0.25, x0 = 0, correct_bias = FALSE)[1] 2.8 5.8 9.9 13.1 16.1 19.6 23.9 28.0# 校正后
ewma(xt, w = 0.25, x0 = 0, correct_bias = TRUE)[1] 11.0 13.3 17.1 19.2 21.1 23.8 27.6 31.1下面再利用R语言同时计算企业A和企业B季度收入的移动均值（设w = 0.25, x0 = 0, correct_bias = TRUE）。在R语言中，可以用数据框（dataframe）存储多组时间序列数据、以其中一个列表示不同组别的标识；需要特别注意的是，在计算前必须先后按组别（下例中的group_id列）和时间点序号（下例中的t列）对数据进行排序：
# 在实际应用中，需确保数据集中不存在缺失值
xt <- data.frame(
group_id = c("A", "A", "A", "A", "A", "A", "A", "A",
"B", "B", "B", "B", "B", "B", "B"),
t = c(1, 2, 3, 4, 5, 6, 7, 8,
1, 2, 3, 4, 5, 6, 7),
x = c(11.0, 15.0, 22.0, 23.0, 25.0, 30.0, 37.0, 40.0,
13.0, 19.0, 20.0, 22.0, 26.0, 32.0, 34.0)
)
# 先按组别、再按时间点排序
xt <- xt[order(xt$group_id, xt$t),]
# 分组计算指数加权移动均值
xt$x <- unlist(
tapply(xt$x, xt$group_id, function(x) {ewma(x, w = 0.25, x0 = 0, correct_bias = TRUE)})
)
# 显示结果
print(xt) group_id t x
1 A 1 11.0
2 A 2 13.3
3 A 3 17.1
4 A 4 19.2
5 A 5 21.1
6 A 6 23.8
7 A 7 27.6
8 A 8 31.1
9 B 1 13.0
10 B 2 16.4
11 B 3 18.0
12 B 4 19.4
13 B 5 21.6
14 B 6 24.8
15 B 7 27.4参考文献
Hyndman, Rob J, and George Athanasopoulos. 2021. Forecasting: Principles and Practice. 3rd edition. Melbourne, Australia: OTexts. https://OTexts.com/fpp3.
Ng, Andrew. 2018. “Optimization Algorithms - Understanding Exponentially Weighted Averages.” In Deep Learning Specialization - Improving Deep Neural Networks: Hyperparameter Tuning, Regularization and Optimization. Cousera.