时间序列分解：stl、prophet原理与实现
参考⽂章：
STL时序分解：
将时序分解为趋势项、季节项（周、⽉等）、余项。利⽤Lowess局部加权回归技术进⾏平滑；通过外循环设计体现鲁棒性。
分别⽤Yv， Tv，Sv，Rv分别代表数据，趋势项、季节项和余项，v的范围为0到N，那么Yv=Tv+Sv+Rv ，其中v=1,…,N （加法模型中，各项具有相同量纲、STL只能处理加法模型，可以先将数据取对数，进⾏STL分解后的各分量结果取指数即可）
T(k)v、Sv(k)为内循环中第k-1次pass结束时的趋势分量、周期分量，初始时T(k)v=0
每个周期相同位置的样本点组成⼀个⼦序列（subries），共有 n(pi)个⼦序列。
涉及参数：
·
n(i) 内层循环数，⼀般为1-2
· n(o) 外层循环数， ⼀般不超过10；没有异常值时，可以设为0
· n(pi) 为⼀个周期的样本数（p在()中会显⽰乱码，⽤pi表⽰），
· n(s) 为Step 2中LOESS平滑参数，窗⼝⼤⼩，越⼤则⼦序列越平滑，⼀般不⼩于7的奇数
· n(l) 为Step 3中LOESS平滑参数，窗⼝⼤⼩，通常为⼤于等于n(pi)的最⼩奇数
· n(t)为Step 6中LOESS平滑参数，窗⼝⼤⼩，通常n(pi)—2n(pi)的奇数
内循环步骤：
Step 1: 去趋势（Detrending），减去上⼀轮结果的趋势分量T(k)v，即Yv−T(k)v
Step 2: 周期⼦序列平滑，⽤LOESS (q=n(s), d=1 )对每个⼦序列做回归，并向前向后各延展⼀个周期（单个⼦序列前后各延⼀个点）；平滑结果组成temporary asonal ries，记为C(k+1)v
Step 3: 周期⼦序列的低通量过滤，对上⼀个步骤的结果序列C(k+1)v依次做长度为n(pi)、n(pi)、3的滑动平均（moving average），然后做LOESS (q=n(l), d=1)回归，得到结果序列L(k+1)v,v=1,…,N ；相
当于周期⼦序列的低通量L(k+1)v；
Step 4: 去除平滑周期⼦序列的低通量，得到季节项：S(k+1)v=C(k+1)v−L(k+1)v
Step 5: 去周期（Deasonalizing），减去季节项分量，Yv−S(k+1)v
Step 6 : 趋势平滑（Trend Smoothing），对于step5中去除周期之后的序列做LOESS (q=n(t), d=1)回归，得到趋势分量T(k+1)v
外循环步骤：
计算更新各样本点v的鲁棒权重ρv值，在对应内循环的Step 2与Step 6中做LOESS回归时，邻域权重需要乘以ρv。
ρv值的计算：ρv=B(|Rv|/h)，其中，h=6∗median(|Rv|)
数据点的余项Rv越⼤，对应ρv权重越⼩；外循环主要⽤于调节robustness weight。如果数据序列中有异常值，则余项会较⼤，可以减⼩影响.
总结：
outer loop:
计算robustness weight;
inner loop:
Step 1 去趋势：Yv−T(k)v;
Step 2 周期⼦序列平滑：C(k+1)v;
Step 3 周期⼦序列的低通量过滤：L(k+1)v;
Step 4 去除平滑周期⼦序列趋势：S(k+1)v=C(k+1)v−L(k+1)v;
Step 5 去周期：Yv−S(k+1)v;
Step 6 趋势平滑：T(k+1)v;
X-11、STL、prophet⽅式时序分解的python实现
1.导⼊相关模块
statsmodels模块实现X-11时序分解；
rpy2模块实现对R语⾔相关模块STL的调⽤，另需要安装R软件；
fbprophet 模块实现prophet分解。
"""
安装R软件后，需要定义相关环境变量，或者直接添加如下三⾏定义，路径为R软件安装位置
如果出现cannot load library R.dll的错误，尝试卸载rpy2模块，移步‘www.lfd.uci.edu/~gohlke/pythonlibs/’⽹页下载⾮官⽅rpy2模块Prophet安装失败可参考"/2129684.html"
holidays功能需要pandas version <1.1.0 (e.g. 1.0.5)
"""
import os
# os.environ['R_HOME']='D:/Program Files/R/R-4.0.3'
# os.environ['PATH'] += os.pathp + 'D:/Program Files/R/R-4.0.3/bin/X64/'
# os.environ['PATH'] += os.pathp + 'D:/Program Files/R/R-4.0.3/'
bjects as robjects
bjects import pandas2ri
bjects import r
import statsmodels.api as sm
from fbprophet import Prophet
import aborn as sns
import matplotlib.pyplot as plt
import math
import numpy as np
import pandas as pd
from pandas import DataFrame
2.利⽤statsmodels模块的asonal_decompo函数实现X-11时序分解（⾮STL分解）2.1算法实现函数decompo_x11如下：
def decompo_x11(dta, freq):
"""
算法⼆：
利⽤statsmodels模块的asonal_decompo函数实现X-11时序分解（⾮STL分解）
dta：预测数据，含'timestamp','value'两列
"""
dta = dta.t_index('timestamp')
dta['value']= dta['value']._numeric, errors='ignore')
dta.value.interpolate(inplace=True)
"""
asonal_decompo函数参数说明：
model："additive"为加法模型， "multiplicative"为乘法模型
freq：指定序列周期
"""
res = sm.tsa.asonal_decompo(dta.value, freq=freq, model="additive")
trend = d
asonal = res.asonal
residual = sid
# 分解模型展⽰
# res.plot()
# plt.show()
return trend, asonal, residual
2.2具体实例的相关说明，另两种算法实现参考本实例
实例数据说明：'report_date_level_3to8.csv’含两列：时间列和待预测时序列；
先进⾏周期7day的分解，对结果趋势项、随机项合并进⾏周期30day的⼆次分解；
add2mult函数：实现加法模型的乘法分解（将df数据取对数，进⾏加法分解，然后将分项取指数即可）。
"""
先对数据进⾏周期为7 day的分解，得到T_1、S_7、R_1项；
将T_1、R_1合并(乘法模型中的合并为乘)进⾏周期为30 day的分解，得到T_2、S_30、R_2项，即原始数据由T_2、S_7、S_30、R_2组成；
如何考虑节假⽇的影响，根据上述R_2值（历史节假⽇）进⾏修正。
具体步骤如下：
step1：待预测数据data取对数，得到log(data）
step2：对log(data）进⾏加法分解（频率为7），并将得到的三分项取指数运算，得到T_1、S_7、R_1
step3：T_1、R_1合并，取对数，进⾏加法分解（频率为30），并将得到的三分项取指数运算，得到T_2、S_30、R_2
step4：原始数据data由T_2、S_7、S_30、R_2的乘法模型组成，根据R_2值（历史节假⽇）进⾏修正
"""
# 获取purcha列的数据，并对列重命名
purcha_data = pd.read_csv('report_date_level_3to8.csv', ucols=[1,2])
purcha_data = ame(columns={'report_date':'timestamp','total_purcha_amtsum':'value'})
def add2mult(df, a_freq, decompo_type=decompo_x11):
"""
# 定义⼀个函数，实现乘法模型与加法模型的转变，即实现上述step1和step2
data：待处理数据，含'timestamp'，'value'两列
a_freq：时序分解周期
decompo_type：时序分解⽅式，x11或stl
返回分解之后的各项T、S、R，各项之间是相乘的关系，但是⽤到的分解模型是加法模型
"""
data = df.copy(deep=True)
"""step1: 待预测数据data取对数"""
data['value']= data['value'].apply(lambda x: math.log(x))
"""step2"""
trend, asonal, residual = decompo_type(data, a_freq)
# 将分解得到的三项由ries转化为dataframe数据类型
T_1 = pd.DataFrame({'timestamp': trend.index,'value': trend.values})
S_7 = pd.DataFrame({'timestamp': asonal.index,'value': asonal.values})
R_1 = pd.DataFrame({'timestamp': residual.index,'value': residual.values})
# ⾸位NAN的填充(X-11时序分解得到的trend、residual的前后半个周期均为空值，对空值⽤最近的⾮空值替代)
if decompo_type == decompo_x11:
m = T_1.shape[0]
T_1.iloc[:int(a_freq/2.0),1].fillna(T_1.iloc[int(a_freq/2.0)][1], inplace=True)
T_1.iloc[m-int(a_freq/2.0):,1].fillna(T_1.iloc[m-int(a_freq/2.0)-1]['value'], inplace=True)
R_1.iloc[:int(a_freq/2.0),1].fillna(R_1.iloc[int(a_freq/2.0)][1], inplace=True)
R_1.iloc[m-int(a_freq/2.0):,1].fillna(R_1.iloc[m-int(a_freq/2.0)-1]['value'], inplace=True)
# 取指数运算
T_1['value']= T_1['value'].apply(lambda x: p(x))
S_7['value']= S_7['value'].apply(lambda x: p(x))
R_1['value']= R_1['value'].apply(lambda x: p(x))
return T_1, S_7, R_1
2.3 具体x11时序分解实例代码
x11分解的结果，因取移动平均算法，trend和residual 分项会有前后半周期缺失，进⾏最近⾮nan值填充，在上述add2mult函数有体现
"""利⽤X-11算法对数据进⾏多周期时序分解(T、R⾸位半周期缺失)"""
# 周期为7的分解，得到T_1, S_7, R_1
T_1, S_7, R_1 = add2mult(purcha_data,7, decompo_type=decompo_x11)
# 对T_1、R_1合并（相乘），得到TR_1
TR_1 = py(deep=True)
TR_1['value']= T_1['value']* R_1['value']
# 周期为30的分解，得到T_2, S_30, R_2
T_2, S_30, R_2 = add2mult(TR_1,30, decompo_type=decompo_x11)
# 原始数据data由T_2、S_7、S_30、R_2的乘法模型组成，根据R_2值（历史节假⽇）进⾏修正plt.figure(figsize=(20,8))
plt.subplot(511)
plt.plot(purcha_data['value'], label="data")
plt.legend()
plt.subplot(512)
plt.plot(T_2['value'], label="trend")
plt.legend()
plt.subplot(513)
plt.plot(S_7['value'], label="asonal_7")
plt.legend()
plt.subplot(514)
plt.plot(S_30['value'], label="asonal_30")
plt.legend()
plt.subplot(515)
plt.plot(R_2['value'], label="residual")
plt.legend()
plt.show()
结果分解如下图：
3 利⽤rpy2模块实现对R语⾔stl函数的调⽤，进⾏STL时序分解3.1算法实现函数如下：
