时间序列 - 论文笔记本

前言：

关于时间序列：我做了很多摸索探究，这里做一个简单的总结记录

---

预测质量评价指标

R 平方: 可决系数，取值范围为 $[0,+\infty )$ ，其值越大，表示拟合效果越好 调用接口为 sklearn.metrics.r2_score，计算公式如下：

$$R^2=1-\frac{S{S}_{res}}{S{S}_{tot}}$$

---

平均绝对误差: 即所有单个观测值与算术平均值的偏差的绝对值的平均。这是一个可解释的指标，因为它与初始系列具有相同的计量单位。取值范围为 $[0,+\infty )$ ，调用接口为 sklearn.metrics.mean_absolute_error ，计算公式如下：

$$MAE=\frac{\sum _{i=1}^n|y_i-{\hat{y}}_i|}{n}$$

---

中值绝对误差: 与平均绝对误差类似，即所有单个观测值与算术平均值的偏差的绝对值的中值。而且它对异常值是不敏感。取值范围为 $[0,+\infty )$ ，调用接口为 sklearn.metrics.median_absolute_error ,计算公式如下：$$MedAE=median(|y_1-{\hat{y}}_1|,...,|y_n-{\hat{y}}_n|)$$

---

均方差：最常用的度量标准，对大偏差给予较高的惩罚，反之亦然，取值范围为 $[0,+\infty )$ ，调用接口为 sklearn.metrics.mean_squared_error ,计算公式如下：$$MSE=\frac{1}{n}\sum _{i=1}^n(y_i-{\hat{y}}_i{)}^2$$

---

均方对数误差: 这个与均方差类似，通过对均方差取对数而得到。因此，该评价指标也更重视小偏差。这指标通常用在呈指数趋势的数据。取值范围为 $[0,+\infty )$ ，调用接口为 sklearn.metrics.mean_squared_log_error ,计算公式如下：
$$MSLE=\frac{1}{n}\sum _{i=1}^n(log(1+y_i)-log(1+{\hat{y}}_i){)}^2$$

---

平均绝对百分比误差：这与 MAE 相同，但是是以百分比计算的。取值范围为 $[0,+\infty )$ ，计算公式如下：$$MAPE=\frac{100}{n}\sum _{i=1}^n\frac{|y_i-{\hat{y}}_i|}{y_i}$$

---

平均绝对百分比误差的实现如下：

def mean_absolute_percentage_error(y_true, y_pred):
    return np.mean(np.abs((y_true - y_pred) / y_true)) * 100

导入上述的评价指标:

from sklearn.metrics import (mean_absolute_error, mean_squared_error,
                             mean_squared_log_error, median_absolute_error,
                             r2_score)

差分

---

为什么要差分？

自回归模型在概念上类似于线性回归，后者所做的假设在这里也成立。

时间序列数据必须是静止的，以消除与过去数据的任何明显相关性和共线性。

• 也就是，在固定时间序列数据中，样本观察的属性或值不取决于观察它的时间戳。

例如，给定一个地区的年度人口的假设数据集，如果观察到人口每年增加两倍或增加固定数量，则该数据是非平稳的。
任何给定的观察都高度依赖于年份，因为人口价值将取决于它与任意过去一年的差距。在使用时间序列数据训练模型时，这种依赖性会导致不正确的偏差。

• 为了消除这种相关性，ARIMA 使用差分使数据平稳。最简单的差分涉及取两个相邻数据点的差值。

例如，上图左图显示了谷歌 200 天的股价。右边的图是第一张图的不同版本——这意味着它显示了谷歌股票 200 天的变化。在第一张图中可以观察到一种模式（每过100天就会上涨），这些趋势是非平稳时间序列数据的标志。

然而，在第二张图中没有观察到趋势或季节性，也没有观察到增加的方差。因此，我们可以说差分版本是平稳的。

平滑

---

一般情况下，处理时间序列的核心任务就是根据历史数据来对未来进行预测。这可以通过许多模型来完成。先来介绍一个最老也是最简单的模型：移动平均。

在移动平均中，假设 ${\hat{y}}_t$ 仅仅依赖 $k$ 个最相近的值，对这 $k$ 个值求平均值得到 ${\hat{y}}_t$ 。公式如下式所示：$${\hat{y}}_t=\frac{1}{k}\sum _{n=1}^k{y}_{t-n}$$

很明显，这种方法不能预测未来很久的数据。因为，为了预测下一个的值，就需要实际观察到之前的值。但这种方法可以对原始数据进行平滑。在进行平滑时，窗口越宽，也就是 k 的值越大，趋势越平滑。对于波动非常大的数据，这种处理可以使其更易于分析。

---

we try something！

通过举个例子来进行说明，这里有一份真实的手机游戏数据，记录的是用户每小时观看的广告和每天游戏货币的支出

• 画图函数

def plotMovingAverage(
    series, window, plot_intervals=False, scale=1.96, plot_anomalies=False
):
    """
    series - 时间序列
    window - 滑动窗口尺寸 
    plot_intervals -置信区间
    plot_anomalies - 显示异常值 
    """
    rolling_mean = series.rolling(window=window).mean()

    plt.figure(figsize=(15, 4))
    plt.title("Moving average\n window size = {}".format(window))
    plt.plot(rolling_mean, "g", label="Rolling mean trend")

    # 画出置信区间
    if plot_intervals:
        mae = mean_absolute_error(series[window:], rolling_mean[window:])
        deviation = np.std(series[window:] - rolling_mean[window:])
        lower_bond = rolling_mean - (mae + scale * deviation)
        upper_bond = rolling_mean + (mae + scale * deviation)
        plt.plot(upper_bond, "r--", label="Upper Bond / Lower Bond")
        plt.plot(lower_bond, "r--")

        # 画出奇异值，upper_bond:上界 ，lowwer_bond下界
        if plot_anomalies:
            anomalies = pd.DataFrame(index=series.index, columns=series.columns)
            anomalies[series < lower_bond] = series[series < lower_bond]
            anomalies[series > upper_bond] = series[series > upper_bond]
            plt.plot(anomalies, "ro", markersize=10)

    plt.plot(series[window:], label="Actual values")
    plt.legend(loc="upper left")
    plt.grid(True)

• 窗口平滑（hour=4）
  
• 窗口平滑（hour=12）
  
• 窗口平滑（hour=24）
  

平滑的目的不是预测的有多准确，而是得到变化趋势，上图所示，是24小时（每天的趋势）

当对时间数据进行平滑时，可以清楚的看到整个用户查看广告的动态过程。在整个周末期间（2017-09-16），整个值变得很高，这是因为周末许多人都会有更多的时间。

• 给出置信区间

置信区间：【rolling_mean - (mae + scale * deviation)，rolling_mean + (mae + scale * deviation)】

• 找出异常值

在这里的检测结果出乎意料，从图中可以看到该方法的缺点：它没有捕获数据中的每月季节性，并将几乎所有的峰值标记为异常。

指数平滑

与加权平滑不同，加权平滑只是加权时间序列最后的 k 个值，而 指数平滑 则是一开始加权所有可用的观测值，而当每一步向后移动窗口时，进行指数地减小权重，这个过程可以使用下面的公式进行表达。$${\hat{y}}_t=\alpha \cdot y_t+(1-\alpha )\cdot {\hat{y}}_{t-1}$$

这里，预测值是当前真实值和先前预测值之间的加权平均值。 𝛼 权重称为平滑因子。它定义了预测将「忘记」最后一次可用的真实观测值的速度。 𝛼 越小，先前观测值的影响越大，系列越平滑。

定义指数平滑的函数：

def exponential_smoothing(series, alpha):
    """
    series - 时间序列
    alpha - 平滑因子
    """
    result = [series[0]]  # 第一个值是一样的，存储平滑的结果
    for n in range(1, len(series)):
        result.append(alpha * series[n] + (1 - alpha) * result[n - 1])
    return result

• 画图函数

def plotExponentialSmoothing(series, alphas):
    """
    画出不同 alpha 值得平滑结果
    series - 时间序列
    alphas - 平滑因子
    """
    with plt.style.context("seaborn-white"):
        plt.figure(figsize=(15, 4))
        for alpha in alphas:
            plt.plot(
                exponential_smoothing(series, alpha), label="Alpha {}".format(alpha)
            )
        plt.plot(series.values, "c", label="Actual")
        plt.legend(loc="best")
        plt.axis("tight")
        plt.title("Exponential Smoothing")
        plt.grid(True)

$\alpha =0.3/0.05$

双指数平滑

为了提取更多细节信息，将时间序列数据分解为两个分量：截距（即水平） ℓ 和斜率（即趋势） 𝑏 。对于截距，可以使用之前所述的平滑方法进行平滑，对于趋势，则假设时间序列的趋势(未来方向变化)取决于先前值的加权变化，对趋势应用相同的指数平滑。因此，可以使用以下函数来表示
在上面的式子中，

• 第一个函数是截距计算公式，其中的第一项取决于序列的当前值；第二项现在分为以前的水平和趋势值。
• 第二个函数描述趋势，它取决于当前时间点的水平变化和趋势和先前值。在这种情况下， 𝛽 系数是指数平滑的权重。
• 第三个函数描述的最终预测是截距和趋势的总和。

def double_exponential_smoothing(series, alpha, beta):
    """
    series - 序列数据
    alpha - 取值范围为 [0.0, 1.0], 水平平滑参数
    beta - 取值范围为 [0.0, 1.0], 趋势平滑参数
    """
    # 原始序列和平滑序列的第一个值相等
    result = [series[0]]
    for n in range(1, len(series) + 1):
        if n == 1:
            level, trend = series[0], series[1] - series[0]
        if n >= len(series):  # 预测
            value = result[-1]
        else:
            value = series[n]
        last_level, level = level, alpha * value + (1 - alpha) * (level + trend)
        trend = beta * (level - last_level) + (1 - beta) * trend
        result.append(level + trend)
    return result

结果：

从上面的图可以看出，可以通过调整两个参数 𝛼 和 𝛽 来得到不同的结果。 𝛼 负责序列的平滑， 𝛽 负责趋势的平滑。值越大，最近观测点的权重越大，模型系列的平滑程度越小。某些参数组合可能会产生奇怪的结果，尤其是手动设置时。

三指数平滑：

你可能已经猜到了，这个想法是增加第三个组成部分：季节性成分。这意味着如果时间序列预计没有季节性，就不应该使用这种方法。模型中的季节性成分将解释围绕截距和趋势的重复变化，并且将通过季节的长度来指定，换言之，通过变化重复的时段来指定。对于季节中的每个观测值，都有一个单独的组成部分，例如，如果季节的长度为 7 天（每周季节性），则有 7 个季节性成分，也即是一周中的每一天。

截距现在取决于序列的当前值减去相应的季节性成分。而趋势保持不变，季节性分量取决于系列的当前值减去截距和其先前值。考虑到所有可用季节的成分都是平滑的，例如，如果我们有星期一季节性成分，那么它只会与其他星期一平均。

你可以阅读更多关于 平滑参考资料如何工作以及如何完成趋势和季节性成分的初始近似。有了季节性因素，模型不仅可以预测未来的一两步，而且还可以预测未来任意的 𝑚 步，这听起来挺nb的。

三次指数平滑模型也称为 Holt-Winters 模型，得名于发明人的姓氏 Charles Holt 和他的学生 Peter Winters。此外，模型中还引入了 Brutlag 方法，用来创建置信区间，其公式描述如下：

其中 $T$ 是季节的长度，$d$ 是预测的偏差。其他参数取自三指数平滑。想要了解更多细节的内容，可以参考 这篇论文。

• 实现起来非常复杂，我实现不出来了，CV了一下大佬的代码

class HoltWinters:


    def __init__(self, series, slen, alpha, beta, gamma, n_preds, scaling_factor=1.96):
        self.series = series
        self.slen = slen
        self.alpha = alpha
        self.beta = beta
        self.gamma = gamma
        self.n_preds = n_preds
        self.scaling_factor = scaling_factor

    def initial_trend(self):
        sum = 0.0
        for i in range(self.slen):
            sum += float(self.series[i + self.slen] - self.series[i]) / self.slen
        return sum / self.slen

    def initial_seasonal_components(self):
        seasonals = {}
        season_averages = []
        n_seasons = int(len(self.series) / self.slen)
        # 计算季节平均值
        for j in range(n_seasons):
            season_averages.append(
                sum(self.series[self.slen * j : self.slen * j + self.slen])
                / float(self.slen)
            )
        # 计算初始值
        for i in range(self.slen):
            sum_of_vals_over_avg = 0.0
            for j in range(n_seasons):
                sum_of_vals_over_avg += (
                    self.series[self.slen * j + i] - season_averages[j]
                )
            seasonals[i] = sum_of_vals_over_avg / n_seasons
        return seasonals

    def triple_exponential_smoothing(self):
        self.result = []
        self.Smooth = []
        self.Season = []
        self.Trend = []
        self.PredictedDeviation = []
        self.UpperBond = []
        self.LowerBond = []

        seasonals = self.initial_seasonal_components()

        for i in range(len(self.series) + self.n_preds):
            if i == 0:  # 季节性成分初始化
                smooth = self.series[0]
                trend = self.initial_trend()
                self.result.append(self.series[0])
                self.Smooth.append(smooth)
                self.Trend.append(trend)
                self.Season.append(seasonals[i % self.slen])

                self.PredictedDeviation.append(0)

                self.UpperBond.append(
                    self.result[0] + self.scaling_factor * self.PredictedDeviation[0]
                )

                self.LowerBond.append(
                    self.result[0] - self.scaling_factor * self.PredictedDeviation[0]
                )
                continue

            if i >= len(self.series):  # 预测
                m = i - len(self.series) + 1
                self.result.append((smooth + m * trend) + seasonals[i % self.slen])

                # 在预测时，增加一些不确定性
                self.PredictedDeviation.append(self.PredictedDeviation[-1] * 1.01)

            else:
                val = self.series[i]
                last_smooth, smooth = (
                    smooth,
                    self.alpha * (val - seasonals[i % self.slen])
                    + (1 - self.alpha) * (smooth + trend),
                )
                trend = self.beta * (smooth - last_smooth) + (1 - self.beta) * trend
                seasonals[i % self.slen] = (
                    self.gamma * (val - smooth)
                    + (1 - self.gamma) * seasonals[i % self.slen]
                )
                self.result.append(smooth + trend + seasonals[i % self.slen])

                # 偏差按Brutlag算法计算
                self.PredictedDeviation.append(
                    self.gamma * np.abs(self.series[i] - self.result[i])
                    + (1 - self.gamma) * self.PredictedDeviation[-1]
                )

            self.UpperBond.append(
                self.result[-1] + self.scaling_factor * self.PredictedDeviation[-1]
            )

            self.LowerBond.append(
                self.result[-1] - self.scaling_factor * self.PredictedDeviation[-1]
            )

            self.Smooth.append(smooth)
            self.Trend.append(trend)
            self.Season.append(seasonals[i % self.slen])

参数：

• series - 时间序列
• slen - 季节长度
• alpha, beta, gamma - 三指数模型系数
• n_preds -预测
• scaling_factor -置信值设置

时间序列交叉验证很简单，我们在时间序列的一小部分上训练模型，从开始到时间节点 𝑡 ，对下一个 𝑡+𝑛 时刻进行预测，并计算误差。然后，将训练样本扩展到 𝑡+𝑛 值，从 𝑡+𝑛 做出预测，直到 𝑡+2∗𝑛 ，并继续移动时间序列的测试段，直到达到最后一次可用的观察。因此，在初始训练样本和最后一次观察之间，可以得到与 𝑛 一样多的数据份数。

• 也就是说，一段一段的移动，一段一段的数据作为测试集

data = ads.Ads[:-20]  # 留下一些数据用于测试

# 初始化模型的参数
x = [0, 0, 0]

# 最小化损失函数
opt = minimize(
    timeseriesCVscore,
    x0=x,
    args=(data, mean_squared_log_error),
    method="TNC",
    bounds=((0, 1), (0, 1), (0, 1)),
)

# 取出最优的参数
alpha_final, beta_final, gamma_final = opt.x
print(alpha_final, beta_final, gamma_final)

# 使用最优参数训练模型
model = HoltWinters(
    data,
    slen=24,
    alpha=alpha_final,
    beta=beta_final,
    gamma=gamma_final,
    n_preds=50,
    scaling_factor=3,
)
model.triple_exponential_smoothing()

从图中可以看出，模型能够成功地拟合原始的时间序列，并捕捉到了每日季节性和整体向下趋势，甚至一些异常也检测到了。如果你看模型的偏差，你可以清楚地看到模型对序列结构的变化作出了相当大的反应，但随后迅速将偏差返回到正常值，基本上是「遗忘」掉过去。该模型的这一特性使我们能够快速构建异常检测系统，即使对于嘈杂的系列数据，也不需要花费太多时间和精力来准备数据和训练模型

时序的特征提取

• 时序的时差
• 窗口统计：
  • 序列在一个窗口中的最大最小值
  • 一个窗口中的中值和平均值
  • 窗口的方差
  • 等等
• 时间特征：
  • 年月日，时分秒
  • 假期、特殊的节日
• 目标编码
• 其他模型的预测结果

画图函数：

def plotModelResults(
    model, X_train=X_train, X_test=X_test, plot_intervals=False, plot_anomalies=False
):
    """
    画出模型、事实值、预测区间和异常值的图表
    """
    prediction = model.predict(X_test)

    plt.figure(figsize=(15, 4))
    plt.plot(prediction, "g", label="prediction", linewidth=2.0)
    plt.plot(y_test.values, label="actual", linewidth=2.0)

    if plot_intervals:
        cv = cross_val_score(
            model, X_train, y_train, cv=tscv, scoring="neg_mean_absolute_error"
        )
        mae = cv.mean() * (-1)
        deviation = cv.std()

        scale = 1.96
        lower = prediction - (mae + scale * deviation)
        upper = prediction + (mae + scale * deviation)

        plt.plot(lower, "r--", label="upper bond / lower bond", alpha=0.5)
        plt.plot(upper, "r--", alpha=0.5)

        if plot_anomalies:
            anomalies = np.array([np.NaN] * len(y_test))
            anomalies[y_test < lower] = y_test[y_test < lower]
            anomalies[y_test > upper] = y_test[y_test > upper]
            plt.plot(anomalies, "o", markersize=10, label="Anomalies")

    error = mean_absolute_percentage_error(prediction, y_test)
    plt.title("Mean absolute percentage error {0:.2f}%".format(error))
    plt.legend(loc="best")
    plt.tight_layout()
    plt.grid(True)

特征重要度可视化：

def plotCoefficients(model):
    """绘制模型的排序系数值
    """

    coefs = pd.DataFrame(model.coef_, X_train.columns)
    coefs.columns = ["coef"]
    coefs["abs"] = coefs.coef.apply(np.abs)
    coefs = coefs.sort_values(by="abs", ascending=False).drop(["abs"], axis=1)

    plt.figure(figsize=(15, 4))
    coefs.coef.plot(kind="bar")
    plt.grid(True, axis="y")
    plt.hlines(y=0, xmin=0, xmax=len(coefs), linestyles="dashed")

时差

• 直接对时间序列进行差分

时间特征

• 现在提取更多的特征，添加小时、日期、是否是周末等特征值。

---

目标编码

• 一种编码类别特征的版本: 按平均值编码。按平均值编码指的是用目标变量的平均值进行编码
• 日和、周和、月和…

特征选择

用模型来回归选择

• 岭回归
• LASSO
• XGBOOST…

深度学习

1. RNN
2. LSTM
3. etc…

原理你自己get吧，我用tensorflow简单的玩一会lstm

#导入深度学习包
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import LSTM , Dense ,Dropout
from sklearn.preprocessing import MinMaxScaler
from keras.wrappers.scikit_learn import KerasRegressor
test_split = round(len(df1)*0.20)
test_split

#切分数据，切记不可随机打乱，而是有顺序的切分
test_split = round(len(df1)*0.20)
df_for_training = df1[:-312]
df_for_testing = df1[-312:]

#缩放数据，归一化
scaler = MinMaxScaler(feature_range=(0,1))
df_for_training_scaler = scaler.fit_transform(df_for_training)
df_for_testing_scaler = scaler.transform(df_for_testing)
df_for_training_scaler


#构造时间步，返回3维数据集（样本数，时间步，属性维度）
def createXY(dataset,n_past):
    dataX=[]
    dataY=[]
    for i in range(n_past,len(dataset)):
        dataX.append(dataset[i-n_past:i,0:dataset.shape[1]])
        dataY.append(dataset[i,0])
    return np.array(dataX),np.array(dataY)

trainX,trainY = createXY(df_for_training_scaler,3)
testX,testY = createXY(df_for_testing_scaler,3)

from sklearn.model_selection import GridSearchCV
#定义模型
def build_model(optimizer):
    grid_model = Sequential()
    grid_model.add(LSTM(50,return_sequences=True,input_shape=(3,13)))  #第一层50个神经元
    #中间的参数自己设置，我这仅供参考，神经元不是越多越好，估量一下参数的数量哦！
    grid_model.add(LSTM(50))#第二层50个神经元
    grid_model.add(Dropout(0.2))#损失降到0.2一下退出
    grid_model.add(Dense(1)) #输出一个时间步

    grid_model.compile(loss='mse',optimizer=optimizer)#评估指标与优化器
    return grid_model


grid_model = KerasRegressor(build_fn=build_model,verbose=1,validation_data=(testX,testY))



#定义评价指标
def mape(true,predictd):
    inside_sum = np.abs(predictd-true)/true
    return round(100*np.sum(inside_sum)/inside_sum.size,2)
#r2_score
from sklearn.metrics import r2_score,mean_squared_error

• 78轮训练，每轮迭代10次，损失收敛
  

• 预测结果
  
• LSTM model MAE is 4.09%
• LSTM model R2_score is 96.87%

那么我根据具体需求更改一下时间步呢？

• 这是t+n的逐步输出，细节化

注意点

• 结果不好不要怪模型，模型很棒，是你的问题好吧！！！
• 数据和特征工程决定模型的高下，别怪模型+1
• 数据预处理应该占用你50%的时间，有因就有果，中间漏了啥，自己积累一些经验
• 机器学习，可不一定输给深度学习，没有免费的午餐定理可以解释

XGBOOST在时序问题的表现，一点不弱呀!! 关键是还具有可解释性嘞

• 准确率score= 0.961624695892833
• 均方根误差rmse= 3.086622697986724
• 均方误差mse= 9.527239679726847

---

论文精读

2020 ICLR）N-BEATS: Neural Basis Expansion Analysis For Interpretable Time Series Forecasting

我们专注于使用深度学习解决单变量时间序列点预测问题。我们提出了一种基于前向和后向残差链接以及非常深的全连接层堆栈的深度神经架构。该架构具有许多理想的特性，包括可解释性、适用于各种目标领域且训练速度快。我们在几个知名的数据集上测试了所提出的架构，包括M3、M4和TOURISM竞赛数据集，其中包含来自不同领域的时间序列。我们展示了对于两种N-BEATS配置，我们的模型在所有数据集上均表现出最先进的性能，相比统计基准提高了11%的预测准确性，并且比去年M4竞赛的获胜者——神经网络和统计时间序列模型之间的领域调整手工混合模型——提高了3%的预测准确性。我们的模型的第一种配置不使用任何时间序列特定组件，其在异构数据集上的表现强烈表明，与常规智慧相反，深度学习基元（如残差块）本身足以解决广泛的预测问题。最后，我们展示了如何增强所提出的架构以提供可解释的输出，而准确性的损失不大。

本文集中于解决

• 单变量（一维序列，输入和输出都是一维的）
• 多步（一次输出多个未来时刻的预测值，而不是自回归）
• 点预测（输出是一个预测值而不是一个预测范围）的问题。预测长度为H，过去输入的长度为nH（n是可以自己设置的常数）。

Doubly Residual Stacking 双残差的堆叠

每个block有两个输出，一个是backcast、一个是forecast，可以把它们理解为，分别是过去的信息和未来的信息。未来的信息可以理解，直接加在一起最终得到预测结果，那过去的信息怎么理解呢？实际上，这就是这篇文章用残差连接的一个技巧，我觉得有点像boosting的思想。可以看模型结构Stack那一块（模型图里中间那一块），经过每个block后，下一个block的输入就会是该block的输入再减去该block的backcast输出。相当于是，stack input中可能有很多不同杂糅的信息，第一个block用stack input中的一部分来预测，然后就可以把该部分（backcast）减掉，将剩下的再输入到第二个block，继续上述操作。这样就相当于，每个block只负责stack input中一部分的预测。同样利用boosting的思想是回归提升树，上个树预测一个值，下一个树预测它离GT还差多少，预测这个残差，然后下下个树再预测剩下的残差。NBeats也利用了这种残差，即将所有block的输出加一起是最终的输出，此外，NBeats还利用了另一种残差，就是是减输入，一开始很复杂的stack input，经过一个block减去backcast，下一个block的输入就更简单，相当于逐渐把输入简化。所以，论文中提到是双残差结构（DOUBLY RESIDUAL STACKING）。

可解释性

每个block中，输出backcast和forecast并不是直接输出Backcast和Forecast，而是输出系数
${\theta }^b$和 ${\theta }^f$，然后再将系数输入到函数 $g^b(.)$和 $g^f(.)$中（函数可以是设定好的，也可以是可学习的） 。

为什么要这样做呢，为什么不直接FC输出Backcast和Forecast呢，而非要先输出系数再将系数输入到函数中得到Backcast和Forecast呢？

• 因为这样可以有一些可解释性，我们知道了系数 。

基（Basis）这个概念在机器学习中很经典，所以在很多领域的DL的方法，都会用到它的思想，有些用深度字典（Dictionary）的方法也可以将其理解为基。

• 你可以把它理解为，将基根据系数来加权组合，就可以从系数中了解到，哪个基对当前输出更重要，也可以自己设计基。本文就是这样，预测系数 ${\theta }^b$和 ${\theta }^f$ ，然后函数 $g^b(.)$和 $g^f(.)$设计为根据系数来对基进行加权。

在这里，基实际上就是一些序列，比如可以将一个正弦曲线设置为基，那这个基可以表示一些季节性的分量，再比如可以将一个线性直线或者二次曲线（多项式）设置为基，这个基可以表示一些趋势性的分量。

• 把基加权线性相加，即可得到预测Forecast或者是Backcast。

集成

1. 训练多个网络，不同的网络采用不同的正则化方法，例如dropout，l2-norm正则化
2. 不同的损失函数：sMAPE、MASE、MAPE
3. 对于每一个时间步H，每一个模型的时间不都不一样，2H、3H、4H…7H，多尺度的思想
4. 采用bagging策略，每个模型随机初始化训练

Neural basis expansion analysis with exogenous variables:Forecasting electricity prices with NBEATSx

• NBEATS的输入和输出是一维序列，是序列的值本身，不掺杂其他的变量维度。

但是在时间序列预测中，协变量（covariate）（也叫外生变量 exogenous variables）是很关键的，比如很多Transformer-based models输入都会用到协变量。在这里简要介绍一下协变量是什么意思，感兴趣的也可以看一下TFT这篇文章，它对不同种类的协变量分别做了处理。

• 协变量一般分为几种，静态协变量，时变协变量，时变协变量又分为未来不可知和未来可知。

静态协变量的一个例子是，我有几个风力发电站，我想预测每个风力发电站的发电量，可能风电站的位置对预测也很重要，所以它的位置坐标可以当成协变量，之所以称为静态的，是因为对每个发电站来说，位置永远是不变的。

未来可知时变协变量的一个例子是，我有周一到周五的数据，想预测周末的，那可以将日期视为协变量。在输入的时候，除了过去的数据本身，还可以将每个时间点属于周几、属于哪个月、甚至是否是国庆节等法定节假日这些日期相关的变量当成协变量给输入进去。为什么叫时变呢，因为这些协变量肯定随时间变化而变化，为什么叫未来可知呢，你想预测周末，你肯定知道你要预测的未来点它究竟属于哪个日期，是不是节假日，也就是说这种协变量实际上是不用预测就已经知道的。

未来不可知时变协变量的一个例子是，我要预测某地未来的温度，我不仅知道过去的温度，我还知道过去的降水量，由于温度可能和降水量有关，所以降水量可以当成协变量输入来辅助预测温度。降水量肯定在未来不知道，这和前面提到的未来可知时变协变量形成了鲜明的对比，也就是说，这种协变量在预测时，你只知道它过去的信息，无法得到它未来的信息，要是真知道了不就跟作弊一样。

• 所以这篇文章就想把NBEATS改进一下，使其也能够输入并处理协变量，同时将协变量对预测的影响做成可解释的，模型结构如下

• 在NBEATS模型中，Stack根据基的不同可以分为建模趋势和建模周期的，即下面两种基：season和trend

• NBEATSx模型则引入了下面这组基，加了一种Stack，用来建模协变量：
  

（Google）Long-term Time Series Forecasting with TiDE: Time-series Dense Encoder

摘要：

最近的研究表明，在长期时间序列预测中，简单的线性模型可以胜过几种基于Transformer的方法。出于这个原因，我们提出了一种基于多层感知器（MLP）的编码器-解码器模型，称为时间序列密集编码器（TiDE），用于长期时间序列预测。该模型具有线性模型的简单性和速度，同时也能处理协变量和非线性依赖关系。从理论上讲，我们证明了我们模型的最简线性模拟可以在某些假设下实现线性动力系统（LDS）的近似最优误差率。从实证上看，我们展示了我们的方法可以在流行的长期时间序列预测基准测试中匹配或胜过之前的方法，同时比最佳的基于Transformer的模型快5-10倍。

TSF的模型发展：

本文介绍了一种简单而有效的深度学习架构，用于预测长期时间序列预测基准测试中比现有的基于神经网络的模型获得更优异的性能。我们的基于多层感知器（MLP）的模型非常简单，没有任何自我注意力、循环或卷积机制。因此，与许多基于Transformer的解决方案相比，它在上下文和预测长度方面具有线性计算比例。

主要贡献如下：

• 我们提出了Time-series Dense Encoder (TiDE)模型架构，用于长期时间序列预测。TiDE使用密集型MLP编码时间序列的过去和协变量，然后使用密集型MLP解码时间序列的未来和协变量。
• 我们分析了我们模型的最简线性模型，并证明当LDS（线性动态系统）[Kal63]的设计矩阵的最大奇异值远离1时，该线性模型可以在LDS中实现近乎最优的误差率。我们在模拟数据集上进行了实证验证，发现线性模型的表现优于LSTMs和Transformers。
• 在流行的实际长期预测基准测试中，我们的模型与以前的基于神经网络的基线相比表现更好或类似（最大数据集上的均方误差降低了>10%）。同时，与最佳的Transformer基于模型相比，TiDE的推断速度快5倍，训练速度快10倍以上。

研究背景

• 长期预测模型可以广义地分为多元模型和单元模型。

多元模型接受所有相关时间序列变量的过去，并预测所有时间序列的未来作为这些过去的联合函数。这包括经典的VAR模型[ZW06]。

单元模型将主要关注基于神经网络的长期预测的先前工作。

• LongTrans 使用具有LogSparse设计的注意层来捕获近似线性空间和计算复杂度的本地信息。

• Informer使用ProbSparse自注意机制来实现对上下文长度的次二次依赖。

• Autoformer 使用趋势和季节分解以及次二次自注意机制。

• FEDFormer 使用频率增强结构

• Pyraformer 使用具有线性复杂度并能够关注不同粒度的金字塔自注意力。

另一方面，单变量模型仅将时间序列变量的未来预测作为该时间序列的过去和协变量特征的函数，即其他时间序列的过去在推理期间不作为输入的一部分。

• 单变量模型分为两种：局部和全局。

局部单变量模型通常针对每个时间序列变量进行训练，并且推理也是针对每个时间序列进行的。经典模型，如AR、ARIMA、指数平滑模型和Box-Jenkins方法属于这一类。

全局单变量模型包括一个共享模型，它接收一个时间序列的过去数据（以及协变量），以预测其未来。然而，该模型的权重是在整个数据集上联合训练的。这个类别主要包括基于深度学习的架构。

在长期预测的背景下，最近观察到一个简单的线性全局单变量模型可以优于基于Transformer的多变量方法。

• Dlinear学习一个从上下文到时间范围的线性映射，指出了对自注意机制的次二次逼近的不足之处。
• 实际上，最近的一个模型PatchTST已经表明，将连续的时间序列片段作为标记输入到基础的自注意机制中，可以在长期预测基准测试中击败DLinear的性能。
  • 请注意，如果它们在相同的测试集上进行相同的任务评估，所有模型类别都可以在多变量长期预测任务上进行公平比较。

例如，它们可以是全局协变量（适用于所有时间序列），例如星期几、假期等，或者针对时间序列的特定因素，例如在需求预测用例中某一天某一产品的折扣。

我们还可以有时间序列的静态属性，用a(i)表示，例如在零售需求预测中不随时间变化的产品特征。在许多应用中，这些协变量对于准确预测至关重要，良好的模型架构应该具备处理它们的能力。

---

最近，观察到简单的线性模型在多个长期预测基准测试中能够胜过基于Transformer的模型。

另一方面，当未来对过去的依赖关系存在固有的非线性时，线性模型将不足。

此外，线性模型无法建模预测与协变量之间的依赖关系，因为没有使用时间协变量，因为它们会损害性能。

---

我们介绍了一种基于简单且高效的MLP的体系结构，用于长期时间序列预测。

在我们的模型中，我们以一种能够处理过去数据和协变量的方式添加MLP形式的非线性。

该模型被称为TiDE（时间序列密集编码器），因为它使用密集的MLP对时间序列的过去和协变量进行编码，然后解码编码后的时间序列以及未来的协变量。

我们的模型以通道独立的方式应用，即模型的输入是一个时间序列的过去和协变量y(i)1:L，x(i)1:L，a(i)，并将其映射到该时间序列的预测值ˆy(i)L+1:L+H。
映射到

• 请注意，模型的权重是使用整个数据集进行全局训练的。我们模型中的一个关键组件是图中右侧的MLP残差块。

模型

残差块：我们在我们的架构中使用残差块作为基本层。它是一个具有ReLU激活的隐藏层的MLP。它还具有完全线性的跳过连接。我们在将隐藏层映射到输出的线性层上使用了dropout，并在输出处使用了层归一化。

• 这两个dropout和层归一化可以通过模型级超参数开启或关闭。

• 其实就是加了Dropout、ReLU非线性、残差连接、LayerNorm的MLP

模型代码

import torch
import torch.nn as nn
import torch.nn.functional as F
from einops import rearrange


# B: Batchsize
# L: Lookback
# H: Horizon
# N: the number of series
# r: the number of covariates for each series
# r_hat: temporalWidth in the paper, i.e., \hat{r} << r
# p: decoderOutputDim in the paper
# hidden_dim: hiddenSize in the paper


class ResidualBlock(nn.Module):
    def __init__(self, in_dim, hidden_dim, out_dim, dropout_rate=0.2):
        super(ResidualBlock, self).__init__()
        self.linear_1 = nn.Linear(in_dim, hidden_dim)
        self.linear_2 = nn.Linear(hidden_dim, out_dim)
        self.linear_res = nn.Linear(in_dim, out_dim)
        self.dropout = nn.Dropout(dropout_rate)
        self.layernorm = nn.LayerNorm(out_dim)

    def forward(self, x):
        # x: [B,L,in_dim] or [B,in_dim]
        h = F.relu(self.linear_1(x))  # [B,L,in_dim] -> [B,L,hidden_dim] or [B,in_dim] -> [B,hidden_dim]
        h = self.dropout(self.linear_2(h))  # [B,L,hidden_dim] -> [B,L,out_dim] or [B,hidden_dim] -> [B,out_dim]
        res = self.linear_res(x)  # [B,L,in_dim] -> [B,L,out_dim] or [B,in_dim] -> [B,out_dim]
        out = self.layernorm(h+res)  # [B,L,out_dim] or [B,out_dim] 

        # out: [B,L,out_dim] or [B,out_dim]
        return out


class Encoder(nn.Module):
    def __init__(self, layer_num, hidden_dim, r, r_hat, L, H, featureProjectionHidden):
        super(Encoder, self).__init__()
        self.encoder_layer_num = layer_num
        self.horizon = H
        self.feature_projection = ResidualBlock(r, featureProjectionHidden, r_hat)
        self.first_encoder_layer = ResidualBlock(L + 1 + (L + H) * r_hat, hidden_dim, hidden_dim)
        self.other_encoder_layers = nn.ModuleList([
            ResidualBlock(hidden_dim, hidden_dim, hidden_dim) for _ in range(layer_num-1)
            ])

    def forward(self, x, covariates, attributes):
        # x: [B*N,L], covariates: [B*N,1], attributes: [B*N,L+H,r]

        # Feature Projection
        covariates = self.feature_projection(covariates)  # [B*N,L+H,r] -> [B*N,L+H,r_hat]
        covariates_future = covariates[:, -self.horizon:, :]  # [B*N,H,r_hat]

        # Flatten
        covariates_flat = rearrange(covariates, 'b l r -> b (l r)')  # [B*N,L+H,r_hat] -> [B*N,(L+H)*r_hat]

        # Concat
        e = torch.cat([x, attributes, covariates_flat], dim=1)  # [B*N,L+1+(L+H)*r_hat]

        # Dense Encoder
        e = self.first_encoder_layer(e)  # [B*N,L+1+(L+H)*r_hat] -> [B*N,hidden_dim]
        for i in range(self.encoder_layer_num-1):
            e = self.other_encoder_layers[i](e)  # [B*N,hidden_dim] -> [B*N,hidden_dim]

        # e: [B*N,hidden_dim], covariates_future: [B*N,H,r_hat]
        return e, covariates_future


class Decoder(nn.Module):
    def __init__(self, layer_num, hidden_dim, r_hat, H, p, temporalDecoderHidden):
        super(Decoder, self).__init__()
        self.decoder_layer_num = layer_num
        self.horizon = H
        self.last_decoder_layer = ResidualBlock(hidden_dim, hidden_dim, p * H)
        self.other_decoder_layers = nn.ModuleList([
                ResidualBlock(hidden_dim, hidden_dim, hidden_dim) for _ in range(layer_num-1)
            ])
        self.temporaldecoder = ResidualBlock(p + r_hat, temporalDecoderHidden, 1)

    def forward(self, e, covariates_future):
        # e: [B*N,hidden_dim], covariates_future: [B*N,H,r_hat]

        # Dense Decoder
        for i in range(self.decoder_layer_num-1):
            e = self.other_decoder_layers[i](e)  # [B*N,hidden_dim] -> [B*N,hidden_dim]
        g = self.last_decoder_layer(e)  # [B*N,hidden_dim] -> [B*N,p*H]

        # Unflatten
        matrixD = rearrange(g, 'b (h p) -> b h p', h=self.horizon)  # [B*N,p*H] -> [B*N,H,p]

        # Stack
        out = torch.cat([matrixD, covariates_future], dim=-1)  # [B*N,H,p+r_hat]

        # Temporal Decoder
        out = self.temporaldecoder(out)  # [B*N,H,p+r_hat] -> [B*N,H,1]
        
        # out: [B*N,H,1]
        return out


class TiDE(nn.Module):
    def __init__(
            self,
            L,
            H,
            r,
            r_hat,
            p,
            hidden_dim,
            encoder_layer_num,
            decoder_layer_num,
            featureProjectionHidden,
            temporalDecoderHidden,
        ):
        super(TiDE, self).__init__()
        self.encoder = Encoder(encoder_layer_num, hidden_dim, r, r_hat, L, H, featureProjectionHidden)
        self.decoder = Decoder(decoder_layer_num, hidden_dim, r_hat, H, p, temporalDecoderHidden)
        self.residual = nn.Linear(L, H)

    def forward(self, x, covariates, attributes):
        # x: [B,L,N], covariates: [B,L+H,N,r], attributes: [B,N,1]
        batch_size = x.size(0)
        
        # Channel Independence: Convert Multivariate series to Univariate series
        x = rearrange(x, 'b l n -> (b n) l')  # [B,L,N] -> [B*N,L]
        covariates = rearrange(covariates, 'b l n r -> (b n) l r')  # [B,L+H,N,r] -> [B*N,L+H,r]
        attributes = rearrange(attributes, 'b n 1 -> (b n) 1')  # [B,N,1] -> [B*N,1]
        
        # Encoder
        e, covariates_future = self.encoder(x, covariates, attributes)

        # Decoder
        out = self.decoder(e, covariates_future)  # out: [B*N,H,1]

        # Global Residual
        prediction = out.squeeze(-1) + self.residual(x)  # prediction: [B*N,H]

        # Reshape
        prediction = rearrange(prediction, '(b n) h -> b h n', b=batch_size)  # [B*N,H] -> [B,H,N]

        # prediction: [B,H,N]
        return prediction

TimesNet

同于Autoformer只集中于时间序列预测，本文提出的TimesNet是一个通用的时间序列神经网络骨干，可处理各种不同的时间序列任务，如最常见的任务：预测、分类、异常检测等等。

• 其实几乎所有的时间序列预测模型也可以当做是通用骨干，比如Autoformer，Informer，FEDformer，Preformer这些Transformer-based模型中只采用Encoder就相当于是一个时间序列的特征提取器，区别在于它们捕获时序依赖性的方式不同。

比如Autoformer是用Auto-Correlation，Informer中的概率稀疏Attention，FEDformer的频域Attention，Preformer中的Multi-Scale Segment-Correlation。

• 还有那些MLP-based模型比如DLinear也可以当做是通用骨干，它是直接采用线性层权重来表示时序依赖性。

1D变2D
这是本文的核心。大部分现有方法都是作用于时间序列的时间维度，捕获时序依赖性。

实际上，现实时间序列一般都有多种模式，比如不同的周期，各种趋势，这些模式混杂在一起。

如果直接对原始序列的时间维度来建模，真正的时序关系很可能隐藏在这些混杂的模式中，无法被捕获。

• 考虑到：现实世界的时间序列通常具有多周期性，比如每天周期、每周周期、每月周期；
• 而且，每个周期内部的时间点是有依赖关系的（比如今天1点和2点），不同的相邻周期内的时间点也是有依赖关系的（比如今天1点和明天1点），作者提出将1D的时间维度reshape成2D的，示意图如下。
• 下图左侧的时间序列具有三个比较显著的周期性（Period 1、Period 2、Period 3），将其reshape成三种不同的2D-variations，2D-variations的每一列包含一个时间段（周期）内的时间点，每一行包含不同时间段（周期）内同一阶段的时间点。变成2D-variations之后，就可以采用2D卷积等方式来同时捕获时间段内部依赖和相邻时间段依赖。

• 每个周期都包括周期内变化和周期间变化。我们基于多个周期将原始的1D时间序列转换为一组2D张量，可以统一周期内和周期间的变化。

那么怎么确定时间序列中的周期性呢？

采用傅里叶变换。给时间序列做傅里叶变换后，主要的周期会呈现对应的高幅值的频率分量。设定超参数k，然后只取top k个最大的幅值对应的频率分量，即可得到top k个主要的周期，这和Autoformer中的处理类似。具体操作如下图，左侧是确定top k个周期，在此只画了三个，然后将1D的时间序列reshape成3种不同的2D-variations（不能整除的可以用padding），对这三种2D-variations用2D卷积进行处理之后再聚合结果即可。

• 单变量示例，说明时间序列中的2D结构。通过发现周期性，我们可以将原始的1D时间序列转换为结构化的2D张量，这可以方便地通过2D卷积核进行处理。通过对时间序列的所有变量进行相同的重塑操作，我们可以将上述过程扩展到多变量时间序列。

TimeBlock

得到k个2D-variations之后该怎么处理呢？

本文提出了TimesBlock，每层TimesBlock又分为两步。首先是要先对这些2D-variations分别用2D卷积（可以是ResNet、ConvNeXt等）或者其他的视觉骨干网络（比如Swin，Vit）处理；其次将k个处理后的结果再聚合起来。

对于第一步，本文采用了一种参数高效的Inception block。Inception block是GoogleNet中的模块，包含多个尺度的2D卷积核。如下图左侧蓝色区域，处理k个2D-variations的Inception block是参数共享的。

• 因此，模型整体的参数量不会随着超参数k的增大而增大，因此本文将其称为参数高效的Inception block（Parameter-efficient Inception block）。

AutoFormer

LightCTS

时序异常检测综述

• 根据输入数据的类型可以分为单变量和多变量。在下文我也称之为单通道和多通道。

• 根据异常的类型可以分为点级别异常、子序列级别异常、整个序列级别异常。

  • 对于整个序列异常，这种情况只能发生在多变量序列中。比如多个商品的销售量时序数据，我们希望找到其中哪几个商品的销售量不太正常，这就是希望检测出某些单变量序列是否异常。

单变量、点级别异常

最直觉的做法就是，如果某个点的观察值和期望值差距过大，大于某个阈值，则说明该点是异常的

• 这种技术也被称为基于模型的方法，是最常见的技术。

所以，基于模型的方法需要估计序列在某个点的期望值。

如果该方法估计某个点的期望值时用了过去、现在、未来的数据，则该方法是基于估计模型的方法；如果该方法估计某个点的期望值时只用了过去的数据，则该方法是基于预测模型的方法。前者不可用于实时在线异常检测场景，后者则可以。

在线异常检测场景中，可以将新来的数据加到旧数据中重新训练模型，也可以使用增量学习的方式更新模型。

基于密度的方法

还有一种可以想象到的做法，找到处于低密度区域的离群点就是异常点。比如，有一个离群点，其他点与该点的距离小于某个距离阈值的点的数量小于某个数量阈值：

• 这种方式在非时序的数据异常检测中非常常见，但是时序数据是有顺序的，因此距离的概念更加复杂。有些方法限制在一个局部滑动窗口内找距离近的点的数量。

基于直方图的方法

使用多个桶来拟合时间序列，形成一个最优近似直方图如下。如果去掉某些点，使得桶的数量更少的情况下，拟合误差还更小，则说明这些被去掉的点是异常值。下图很形象：

多变量、点级别异常

如果将上述单变量异常检测算法应用于多变量时序数据，则没有考虑变量（通道）之间的依赖性，所以效果不佳。

为了克服这一问题，有些方法采用降维的方式，将多变量序列转化为更少维度的变量序列，而且迫使降维后的序列变量之间尽可能相互独立、不相关，这样的话，就依然可以在降维后的多变量序列中采用单变量异常检测算法。如下图所示：

• 还可以采用基于模型的方法，比如重建，用自编码器来重建。异常点重建误差更大一些。或者预测模型，异常点预测误差更大一些。

• 还可以采用基于差异的方法，计算多元点的期望值和实际值之间或它们的表示之间的两两不相似性。

单变量、子序列级别异常

多变量、子序列级别异常

整个序列级别异常

这个有点像为序列分类，分为异常或者正常。只不过，样本分布可能是高度不平均的，因为异常序列要远远低于正常序列。

• 一种方法是对每个单变量时间序列抽取低维特征，一个变量的序列可以用一个低维特征表示。然后用降维后的特征点来聚类找出离群点，即可找到对应的异常序列。

• 另一种方法是基于差异，分析时间序列之间的成对差异。也可以用一些聚类的技术。