多周期时间序列分解算法——MSTL原理

1.时间序列分解和STL算法

时间序列算法分解一般是指把一个时间序列分解成：趋势序列，周期序列，残差序列。

时间序列分解算法最广为人知的可能是STL算法。它只能分解出一个周期序列。有很多博客和文章叙述了STL分解的原理，例如博客 时间序列分解算法：STL。其中也有原论文的链接，看原论文也较容易理解。

2.MSTL分解算法

MSTL可以理解成multiple-STL。它和STL的区别是：STL分解的输入包含一个周期数，而这个算法可以输入多个周期数。MSTL可以分离出多个周期序列。
　　　STL：Y(t)=Trend(t)+Season(t)+Remainder(t)
　　MSTL：Y(t)=Trend(t)+Season1(t)+Season2(t)+…+SeasonN(t)+Remainder(t)

它的原理简述如下：
————————————————————————————————————
Step0：初始化。
　　　　原始序列设为x，空缺值用插值补全；
　　　　将给定的所有周期数按从小到大排序，记为数组msts；
　　　　各个周期成分seas初始值均设为0，seas可以理解成一个二维数组，第一维的长度为
　　　　周期的个数；
　　　　给定算法内部循环次数iterate；
　　　　去周期成分设为deseas，初始值为原始序列deseas=x；
Step1：对下述{}内的过程循环iterate次后停止，特别的第j次循环如下。
　　　　{从小到大遍历周期数组msts中的周期，运行下述(1.a)~(1.d)，特别的，对第i个周期
　　　　msts[i]，执行下述操作：
　　　　　　（1.a）更新去周期成分：deseas=deseas+seas[i]（seas[i]为第i个周期对应的序
　　　　　　　列,初始设为0向量)；
　　　　　　（1.b）对deseas用stl分解算法以第i个周期（msts[i]）进行分解，分解模型记为
　　　　　　　fit；
　　　　　　（1.c）更新第i个周期对应的周期序列为fit的周期序列，即
　　　　　　　　　　　　　　　　　seas[i]=seasonal(fit)；
　　　　　　（1.d）更新去周期成分，deseas=deseas–seas[i]。
　　　　　遍历msts中的每个周期数之后，计算原序列的其他成分：
　　　　（2.a）计算原始序列的趋势成分，trend=trendcycle(fit)。
　　　　（3.a）计算原始序列分解的残差成分，remainder=deseas-trend。
　　　　}
Step2：输出原始序列的分解结果：趋势序列trend，若干周期序列seas，残差序列　　
　　　　remainder。
————————————————————————————————————
可以看到：
（1）第一次循环，从小到大遍历周期时，用STL算法把各个周期序列“剥离”出来，存到seas
　　　里。然后计算一个什么周期都没有trend和remainder。注意这时候deseas等于　　　
　　　trend+remainder。
（2）第二次循环，从小到大遍历周期时，这时STL分解的是trend+reminder+seas[i]，没有其
　　　他周期项的干扰，这时分解出第i个周期对应的周期序列应该更好。用它来更新seas[i]。
　　　然后再次更新trend和remainder。注意这时候deseas一直等于trend+remainder。
（3）依次遍历，不断地优化所有的子序列。

3.参考资料

上述算法是从R的mstl()源码解析出来的，可能有误差，详情请参考R函数mstl()源码地址 。