Kaggle -Linear Regression with Time Series

时间序列的回归拟合

使用线性回归方法预测时间序列问题。过程使用时间序列特有的两个特性：滞后性（lags）和时间步长（time steps）。

时间序列问题：在预测应用中，通常使用固定频率记录，比如每天，每月

时间步长

时间步长特征是我们可以直接从时间索引推导出来的特征。最基本的时间步长特征是time dummy，他计算序列从开始到结束的时间步长。如果一个序列的值，可以从他出现的时间开始预测，那么他就是时间相关的。

滞后性

为了实现时间的延迟特性，可以调整观测值，使其看起来出现在较晚的时间。例如1步滞后：

Date	values	lag_1
2016-01-01	25	Nan
2016-01-01	65	25
2016-01-01	21	65
2016-01-01	3	21

上图就是滞后指标（Lag_1， 就是前一天的值）和当天指标有正相关关系。第一天的值很大，就可以预测第二天的值也会比较大。
更一般地说，延迟特性可以建立具有串行依赖关系的模型。当一个观测值可以从以前的观测值中预测时，时间序列具有序列依赖性。

一般来说，一个好的和时间相关的模型通常包含一些时间不长和滞后特性的组合

1. 趋势 – Trend

趋势反映一个序列的平均值的长期且持续的变化， 是一个序列中变化最慢的变量，也是最重要的时间尺度。下面将分析平均值的趋势。注意，序列中任何一个缓慢变换的量都可以作为一个趋势。

趋势除了作为复杂模型的基线外，还可以用来作为那些不能自动学习到趋势特征的算法的特征，例如随机森林或XGBoost

1.1 平均移动 – Moving Average Plot

为了计算时间序列的趋势，使用平均移动图。使用平均移动图可以平滑锯齿，使图像只反映长期变化。

1.2 工程的趋势 – Engineering Trend

当确定了趋势的形状后就可以尝试使用趋势对问题建模。在线性回归问题中，采用了

$target=a\ast time+b$

那么在了解了趋势的形状以后，就可以根据趋势的形状进行建模。如果趋势是抛物线，就可以使用

$target=a\ast time\ast \ast 2+b+time+c$

1.3 季节性

当一个时间序列的均值具有周期性变化时，便成这个时间序列有季节性。
下面接触两种模式下的季节性指标。

1. indicators（指标）：适合较少的观察，例如周期为一周
2. Fourier features（傅里叶特征）：适合观察对象很多时，例如周期为一年

1.3.1 季节性曲线（Seasonal Plots）

类似于Average Moving Plots，可以使用Seasonal Plots观察季节。季节图显示的是时间序列的片段（segments），每个片段对应于某个共同的时期（时期），这个时期就是想要观察的季节。

1.3.2 季节性指标（Seasonal Indicators）

如果将一个季节周期视为一个分类的特征，那么应该将其设置为one - hot向量，而不是0代表一个季节，1代表一个季节，2代表一个季节…因为特征的数值的大小反应其重要程度。假设春0夏1秋2冬3，不能说明冬天要比春天重要。
另外，在令季节作为特征时，如果舍弃了一个季节属性，那么是有利于回归问题的，如下舍弃了周一。

添加季节性指标可以帮助模型区分季节性期间的平均值。季节性指标的作用类似于一个开关。

1.4 傅里叶特征和周期图

当观察long seasons的时候，更适合使用傅里叶特征，因为单一的indicator是不切实际的，我们尝试使用多个特征来表征季节曲线的整体形状（而不是单一的indicator）。在捕捉季节变化的频率时，我们希望建立的模型的形状和季节的形状一样，因此会用到sin 或 cos一类曲线。
傅里叶特征是一对正弦和余弦曲线，从最长的周期开始，每个潜在频率对应一对。傅里叶对建模年度周期性的频率：每年一次、每年两次、每年三次，等。

如果将一组正弦曲线作为特征加入训练集当中，线性回归算法将计算出适合目标序列中的季节性成分（傅里叶特征）的权重。通过傅里叶特征只对季节性的“主要影响”进行建模，会减少向训练集中添加的特征，这会减少计算时间和过拟合风险。

1.4.1 用周期图选择傅里叶特征

通过周期图可以决定我们应当在特征集中包含多少个傅里叶对。
在对季度进行傅里叶建模时就可以忽略每周的周期性变化，因为它又可以用indicator来建模。例如使用indicator对每周的season进行建模，使用傅里叶对每季度的season进行建模。

1.4.2 计算傅里叶特征

下方是使用自定义的代码段计算傅里叶特征

def fourier_features(index, freq, order):
    '''
    根据输入的index和周期freq，返回order阶数的傅里叶特征
    :index  : 输入的数组，（我们只关心它的长度）
    :freq   : 周期
    :order  : 阶数
    :return : 
    '''
    time = np.arange(len (index), dtype = np.float32)
    k = 2 * np.pi * (1 / freq) * time  # 最后求的是 sin(k) = sin(2Πω * x)
    features = {}
    for Order in range(1, order + 1):
        # 建立1，2，3...,order阶傅里叶
        features.update({
            f'sin_{freq}_{Order}':np.sin(Order * k),
            f'cos_{freq}_{Order}':np.cos(Order * k)
        })
    return pd.DataFrame(features, index = index)
fourier_features(np.arange(1000), 365//4, 4)

在程序中使用stasmodels.tsa.deterministic库创建indicators和Fourier特征

2. 自相关

首先引入相关性的概念。

2.1 相关性

相关性是两个随机变量$X_i和Y_j$之间线性关系的强度和方向。

相关系数越接近于 1，表明二者越具有正相关性；越接近 -1，越具有负相关性。

2.1.1 皮尔逊相关系数

计算随即变量之间相关性的函数有很多，其中最常用的便是皮尔逊相关系数。
$${\rho }_{X,Y}=\frac{cov(X,Y)}{{\sigma }_X{\sigma }_Y}=\frac{E((X-{\mu }_X)(Y-{\mu }_Y))}{{\sigma }_X{\sigma }_Y}$$

2.1.2 平均移动曲线

平均移动曲线一类又称为滚动统计特征。通常使用GBDT算法，如XGBoost算法时，是一个经常被选择的优秀特征。但要注意在作为特征使用时要避免超前，例如将center设置为False，lag设置为1

y_ = y.shift(1)
y_ = y_.rolling(
    window = 8,
    center = False
).mean()
y_ = y_.dropna()

2.2 自相关 autocorrelation

又称作序列相关（serial correlation），指一个时间序列这些值前后自己相关，也就是$X_{i-1}和X_i$之间的相关系数。

一个时间序列的自相关系数又叫做自相关函数（ACF）。

在选择**滞后（lag）**作为特征加入特征集时，无需将所有延迟都加进来。因为2位滞后的数据所含的信息往往是1阶滞后信息的衰减信息，所以如果滞后2中没有产生新的信息，我们就没有理由再去添加它了。

2.3 偏相关 partial autocorrelation

所谓滞后k偏自相关系数就是指，在给定中间k-1个随机变量$X_{t-1},X_{t-2},...X_{t-k+1}$的条件下，即剔除了中间k-1个随机变量的干扰后，$X_{t-k}$对$X_t$影响的相关度量。
$${\rho }_{X_t,X_{t-k}|X_{t-1}...X_{t-k+1}}=\frac{E[(X_t-\hat{E}{X}_t)(X_{t-k}-\hat{E}{X}_{t-k})]}{E(X_{t-k}-\hat{E}{X}_{t-k}{)}^2}$$

• 概括地讲：偏相关系数反应本次滞后带来的新的信息的数量。

从上图可以得到的信息有，滞后1-6是有用的，而从滞后7开始，新的增益在置信区间外，可以视为无关特征。至于滞后11所显示出的有用增益，需要认为是假阳性。但要注意的是，并不是所有的间隔着无关特征的滞后（如滞后11）都需要认为是假阳性，有可能当天数据就与12天前的数据相关。

上图称为相关图，本质上就代表着傅里叶特征。因为真实世界中存在着大量的非线性相关数据，所以在采用之后特征之前需要先进行相关图的绘制和观察。

使用df.corr()计算两个数组的相关系数

2.3 平稳性 stationarity

• 指一个时间序列的概率分布不随时间变化的性质。

一般通过差分使数据变得平稳。但是做了差分操作后会丢失部分数据的性质。非平稳的数据不可以直接做回归，例如国家GDP的随机变量X在增长，身高随机变量Y也在增长，但不能说身高的增长导致了GDP的增长，原因就是这两个时间序列都是不平稳的。

2.4 时间依赖性和序列依赖性

时间依赖性指周期性明显的依赖于时间，随时间的变化而呈现出周期性。一般直观理解的具有周期性的数据都是具有时间依赖性的。本节着重介绍序列依赖性

序列依赖性

不同于时间依赖性，尽管序列依赖性也呈现出与时间相关的周期性，但其每个点的实际取值和所在时间段，或者说受时间步长的影响不大，更加直接的反应其周期性的是，某一时刻的值如何依赖前一时刻的值。

如图所示，流感（上图数据）趋势数据显示的是不规则的周期（区别于没有周期），而非正常的周期性：

1. 流感的高峰期往往发生在新年前后， 但有时会提前或滞后
2. 有时会更大或更小。

使用滞后特征来对问题建模使得可以根据不断变化的条件做出反应，而不像时间依赖性数据那样被固定在特定时间上。

import pandas as pd
def make_lags(ts, lags):
    return pd.concat({
        f'y_{lags}_{i}':ts.shift(i)
        for i in range(1, lags + 1)
    },
    axis = 1)

X = make_lags(flu_tends.Fluvists, lags=4)
X = X.fillna(0,0)

! 待解决的问题

使用滞后特征的模型在预测时，因为需要预测前一天的结果，就不得不使用测试集中的数据。至于如何解决此问题，将在以后作答

2.5 去季节化 Deseasonalize

对一个并没有很好的季节性和趋势的时间序列进行建模时， 趋势和季节性都会在相关图和滞后图中产生序列依赖性。为了去除纯粹的季节性行为，需要首先对数据进行去季节化处理（Deseasonalize）。

去季节化采用的操作通常是使用季节性特征进行预测，然后用原始数据减去预测结果得到去季节化后的数据。

fourier = CalendarFourier(freq = 'M', order = 4)

y = pd...
dp = DeterministicProcess(
    constant = True,
    order = 1,
    seasonal = True,
    drop = True,
    index = y.index,
    additional_items = [fourier]
)
X = dp.in_sample()

model = LinearRegression()
model.fit(X, y)
pred = model.predict(X)

y_deseason = y - pred

3. 混合模型

线性回归可以很好的推断趋势；XGBoost模型可以很好的反应交互，结合二者的优点可以做到更好。

3.1 组件

“趋势，时间依赖和序列依赖” 加上一些不可避免的误差就构成了时间序列的组件。记公式为
$$series=trend+seasons+cycles+errors$$
其中每个成员都叫做一个组件（component）

模型的残差是模型训练的目标和模型预测差值。根据特征绘制残差，就可以得到模型从未从特征中了解的信息。

将学习时间序列的组件作为一个迭代的过程。首先学习趋势，减去它的级数，然后学习周期（时间依赖和季节依赖）并减去周期，最后就只剩下那些不可预知的错误。

然后将所有学习到的组件加在一起，就得到了完整的模型。

以上就是线性回归要做的内容

3.2 残差 Residuals

3.1节中使用单一的算法——线性回归一次性学习了所有的组件，但是实际应用上可能对不同组件使用不同模型，为每个组件选择其最好的算法。为此，我们使用一种算法去拟合原始序列，另一种算法去拟合残差序列。

# 使用第一种模型训练原始序列
model1.fit(X_train_1, y_train_1)
y_pred = model1.pred(X_train)  # 为了得到残差，这里在训练集上预测

# 使用第二种模型训练残差
model_2.fit(X_train_2, y_train - y_pred_1)
y_pred_2 = model2.pred(X_train_2)

# 最后得到总的预测结果
y_pred = y_pred_1 + y_pred_2

这里X_train_1和X_train_2的维度不同但个数相同。原因是，因为我们希望两次训练模型训练的内容不同，例如采用第一个模型了解趋势，那么第二个模型训练时X_train_2通常不需要趋势特征。

3.3 混合模型

构造混合模型常见的策略是，一个简单的学习算法（通常是线性模型）后跟着一个复杂的、非线性的模型如GBDT或深度神经网络。

3.3.1 回归问题

回归模型通常分为两种预测策略，一种是转换特征，一种是转换目标

• 回归——转换特征 :
  特征转换学习算法是学习一些特征函数，将特征作为输入，然后组合和转换他们，以产生符合目标值的输出。常见的有线性回归和神经网络。

• 回归——转换目标 :
  目标转换算法是利用特征将targets进行分组，然后通过平均目标值进行预测。一组特征仅仅表明要平均哪个组，代表性的算法有决策树和k邻近。

3.3.2 决策树不能预测趋势

特征转换算法通常可以推断出训练集之外的结果，但是会在一定程度上受限于训练集。

如果给出了一段time dummy，线性回归会在训练集外按照相同的time dummy继续绘制趋势曲线。

但是决策树一类算法只会依据训练集的最后一步推断未来。也就是说决策树类算法不能用于趋势回归，同样的，基于决策树的XGBoost和GBDT算法也不能用于趋势回归问题。

红色曲线是决策树预测的效果

3.3.3 堆叠的混合模型 stacked hybrids

混合模型设计的初衷就是为了弥补决策树的这种缺点。使用线性回归推断趋势，转换目标以消除趋势，然后将XGBoost用于非趋势残差。为了混合一个神经网络，可以将另一个模型的预测作为一个特征。

残差拟合的方法和梯度增强（gradient boosting）算法的方法是一样的，所以称为混合增强算法，使用预测作为特征的方法被称作“叠加”，统称***堆叠的混合模型（stacked hybrids）***

3.4 大致过程

1. 首先使用线性回归模型了解每个系列的趋势。

1. 构建趋势特征

dp = DeterministicProcess(
    order  = 2,
    constant = True,
    index = df.index,
    drop = True
)

X = dp.in_sample()

如果想要预测多组数据：线性回归可以实现多输出回归，但XGBoost不能。所以使用XGboost时需要转换为长格式（类别按行索引），然后设置一个标签使XGBoost区分种类

2. 从日期中提取月份（或周）创建季节性特征

3. 把得到的趋势预测结果转换为长格式，然后从原始序列中减去他们。这将得到XGBoost需要的残差序列。

这样一来，XGBoost就弥补了其无法预测长期趋势的缺点。