时序数据分解

时序数据

  时序数据作为与时间强相关数据，有着独特的特点，但是也有很多通用的数据的性质。

1.通过数学期望与协方差进行特征相关性的计算；
2.平稳性检验

  定义上的平稳性指的是固定时间和位置的概率分布与所有时间和位置的概率分布相同的随机过程。其数学期望和方差这些参数也不随着时间和位置变化。
  即固定时间长度，所得到的时序滑窗中分布和概率都是相同的，这在真实数据中完全不可能达到。
  因此提出弱平稳的定义，它要求时间序列满足均值平稳性和二阶平稳性。
需要满足两个条件：1). 对于所有时刻的t，数学期望都是一个常数；2).对于所有的时刻t和任意的间隔k，值之间的协方差与时间t无关，它仅仅以来于间隔k，满足方差平稳性。

  时序数据可分为两种：连续性预测（数值预测，范围估计），离散型预测（事件预测）等。
  与普通回归模型不同的是，时间序列模型依赖于数值在时间上的先后顺序，并依赖从中学习到的抽象特征进行输出。同样大小的值改变顺序后输入模型产生的结果是不同的。
  时序数据可以分为平稳序列，即存在某种周期，季节性及趋势的方差和均值不随时间而变化的序列，和非平稳序列。

几种常见的时序分析方法：

1.传统参数法：

  首先需要对观测值序列进行平稳性检测，如果不平稳，则对其进行差分运算直到差分后的数据平稳；
  在数据平稳后则对其进行白噪声检验，白噪声是指零均值常方差的随机平稳序列；
  如果是平稳非白噪声序列就计算ACF（自相关系数）、PACF（偏自相关系数），进行ARMA等模型识别，对已识别好的模型，确定模型参数，最后应用预测并进行误差分析。
  该方法需要一定的统计学知识，且对大数据的数据挖掘场景中比较难适用，因为需要大量的参数化建模。比如，有一个连锁店的销售数据，需要预测每个门店的未来销量，就需要对每个门店都进行单独建模，使用起来较为繁琐，门槛较高。

2.时间序列分解

  时间序列分解法是数年来一直非常有效的方法，一个时间序列往往是一下几种变化形式的叠加或耦合：

• 长期趋势(Secular trend, T)：长期趋势指现象在较长时期内持续发展变化的一种趋向或状态。
• 季节变动(Seasonal Variation, S)：季节波动是由于季节的变化引起的现象发展水平的规则变动
• 循环波动(Cyclical Variation, C)：循环波动指以若干年为期限，不具严格规则的周期性连续变动
• 不规则波动(Irregular Variation, I): 不规则波动指由于众多偶然因素对时间序列造成的影响

3.常用的时间序列分解模型，可以分为以下几种形式：

• 加法形式：
  X t = T t + S t + C t + I t X_t= T_t + S_t + C_t + I_t Xt​=Tt​+St​+Ct​+It​
• 乘法形式：
  X t = T t ∗ S t ∗ C t ∗ I t X_t= T_t * S_t * C_t * I_t Xt​=Tt​∗St​∗Ct​∗It​
• 加乘混合模式：以上两种方式的混合

4.时间序列的长期趋势分析

• 平均移动法
  在原时间序列内一次求连续若干其的平均数作为其中某一期的趋势值，如此逐项递移得到一系列的移动平均数，形成一个平均数时间序列；

• 时间回归法
  使用回归分析中的最小二乘法，以时间t或t的函数为自变量拟合趋势方程。常见的趋势方程如下：
  一阶线性方程
  二次（多次）方程曲线
  指数曲线

5.时间序列季节变动性

• 乘法模型
  乘法模型中的季节成分通过季节指数来反映。常用的方法称为移动平均趋势剔除法。步骤如下：
  1）计算一动平均值
  2）从序列中剔除移动平均值
  3）季节指数（S）= 同季平均数/总季平均数 *100%

6.时序长期看存在着某一个循环往复，通常通过剩余法来计算循环变动成分C：

  如果有季节成分，计算季节指数，得到季节调整后的数据TCI
  根据趋势方程从季节调整后的数据中消除长期趋势，得到序列CI
  对消去季节成分和趋势值的序列CI进行移动平均以消除不规则波动，得到循环变动成分C

二、 时序分解方法：

fast robust stl：

1.去噪：

  使用双边滤波对时间序列数据进行去噪（来自图像处理中的边缘保留滤波器的思想）。双边滤波的主要思想是使用具有相似值的邻居来平滑时间序列。应用于时间序列数据时，可以完全保留趋势的突变。
滤波公式为：
y ′ = ∑ j ∈ J w j t y j , J = t , t ± 1 , t ± 2 , . . . , t ± H y' = \sum_{j\in J} w_j^t y_j,J =t,t±1,t±2,...,t±H y′=j∈J∑​wjt​yj​,J=t,t±1,t±2,...,t±H
  其中J表示长度为2H+1的滤波窗口，滤波权重通过两个高斯函数给出，
w j t = 1 z e − ∣ j − t ∣ 2 2 δ d 2 − ∣ y j − y t ∣ 2 2 δ d 2 w_j^t = \frac{1}{z} e^{-\frac{|j-t|^2}{2\delta^2_d} - \frac{|y_j-y_t|^2}{2\delta^2_d}} wjt​=z1​e−2δd2​∣j−t∣2​−2δd2​∣yj​−yt​∣2​

其中1/z是归一化因子， δ d δ_d δd​和 δ i δ_i δi​是两个控制平滑度的参数。可以看出双边滤波（bilateral filtering）类似于平滑数据中的核平滑技术，它使用权重函数或者核函数来平均观测值。

def denoise_step(sample, H=3, dn1=1., dn2=1.):
    def get_denoise_value(idx):
        start_idx, end_idx = get_neighbor_idx(len(sample), idx, H)
        idxs = np.arange(start_idx, end_idx)
        weight_sample = sample[idxs]

        weights = np.array(list(map(lambda j: bilateral_filter(
            j, idx, sample[j], sample[idx], dn1, dn2), idxs)))
        return np.sum(weight_sample * weights)/np.sum(weights)

    idx_list = np.arange(len(sample))
    denoise_sample = np.array(list(map(get_denoise_value, idx_list)))
    return denoise_sample

2.趋势性提取：

去除噪声之后，在提取趋势分量时，仍然存在季节性分量和异常值的干扰，。文中假设季节性分量变化缓慢，首先对去噪信号进行季节性差分操作以减轻季节性影响

之后通过gt恢复趋势信号的一阶差分

def trend_extraction(sample, season_len, reg1=10., reg2=0.5):
    sample_len = len(sample)
    season_diff = sample[season_len:] - sample[:-season_len]
    g = season_diff[:, None]

    assert len(season_diff) == (sample_len - season_len)
    q = np.concatenate([season_diff, np.zeros([sample_len*2-3])])
    q = np.reshape(q, [len(q), 1])

    M = get_toeplitz([sample_len-season_len, sample_len-1],
                     np.ones([season_len]))
    D = get_toeplitz([sample_len-2, sample_len-1], np.array([1, -1]))
    P = np.concatenate([M, reg1*np.eye(sample_len-1), reg2*D], axis=0)
    '''
    l1 approximation solver
    '''
    q = matrix(q)
    P = matrix(P)
    delta_trends = l1(P, q)

    relative_trends = get_relative_trends(delta_trends)

    return sample-relative_trends, relative_trends

3.季节性提取

去除相对趋势成分之后，y‘’可以视为被污染的季节性。

def seasonality_extraction(sample, season_lens, alphas, z, K=2, H=5, ds1=50., ds2=1.):
    sample_len = len(sample)
    idx_list = np.arange(sample_len)

    def get_season_value(idx, season_len):
        idxs = get_season_idx(sample_len, idx, season_len, K, H)
        if idxs.size == 0:
            return sample[idx]

        weight_sample = sample[idxs]
        weights = np.array(list(map(lambda j: bilateral_filter(
            j, idx, sample[j], sample[idx], ds1, ds2), idxs)))
        season_value = np.sum(weight_sample * weights)/np.sum(weights)
        return season_value

    seasons_tilda = 0.
    for season_len, alpha in zip(season_lens, alphas):
        seasons_tilda += alpha * np.array(
            [get_season_value(idx, season_len) for idx in idx_list])
    seasons_tilda /= z

    return seasons_tilda

4.多季节分解

def adjustment(sample, relative_trends, seasons_tilda, season_len):
    num_season = int(len(sample)/season_len)

    trend_init = np.mean(seasons_tilda[:season_len*num_season])

    trends_hat = relative_trends + trend_init
    seasons_hat = seasons_tilda - trend_init
    remainders_hat = sample - trends_hat - seasons_hat
    return [trends_hat, seasons_hat, remainders_hat]