一种基于压缩感知的时间序列预测方法

张亚宁，马军海

（天津大学 管理与经济学部，天津 300072）

0 引言

对于时间序列数据来说已知数据和未知数据间通常都存在某种线性或非线性的关系，时间序列预测方法就是通过将它们间的这种关系估计出来而进行预测的。在现有的预测方法中以参数模型估计法最为流行，其中又以AR模型、ARMA和ARIMA模型的研究最为广泛。这三种模型都假定时间序列是由某个白噪声E={ei}，i=1，2…，N激励而来的。参数模型估计法的任务就是在噪声未知的前提下，估计出模型参数。以AR模型为例，AR模型的参数估计往往需要假定观测数据为平稳时间序列，然后根据自相关函数建立起Yule-Walker方程(或 Wiener-Hopf方程)，再通过 Levinson递推算法(或最小二乘)估计出参数[1]。但对于很多经济数据来说，采用AR模型预测存在一些问题，如经济数据往往呈现出递增、递减或周期性的变化趋势，通常都不能满足平稳时间序列均值和方差恒定的前提，再比如所能获得的经济数据量经常有限，而自相关函数往往需要大量的数据才能较为准确的估计出来。除此之外，p阶AR模型只能考虑到距离预测值最近的p个已知数据的贡献，而忽略其他数据。ARMA和ARIMA模型都以AR模型为基础，因此也存在类似的问题。为解决这些问题，本文提出了基于压缩感知[2]的参数估计方法，该方法并不局限于前p个观测值对预测值的贡献，而是根据数据的特点在一定范围内找到最优的选择；并且该方法对时间序列的平稳性没有要求，只要求数据间确实存在某种可以表达的关系。

1 数学模型

若预测值为它的前期值和随机项的线性函数，则时间序列数据预测模型可由式(1)描述：

其中，xn+1为预测数据，xi，i=1，2，…，n为前期值，ai为模型的待估计参数，en+1为随机项，服从相互独立的均值为0、方差为σ2的正态分布。如果每一组观测数据，均可由同一组模型参数表达，则式(1)表达为矩阵形式：

为便于表达，本文将式（2）表达为矩阵相乘的形式：p=Ru+e。矩阵向量p、R、u和e的定义见式(2)，其中R是一个M+1行、N+1列矩阵。根据极大似然估计法，似然函数L为：

由于随机项向量e未知，模型参数向量u无法通过直接解方程的方法求出。由于似然函数L取极大值等价于，因此模型参数向量u的可通过式(4)求解：

下面利用正交匹配追踪法实现式(4)所示的优化问题的求解。

2 算法原理与实现

AR模型和ARMA模型在经济数据预测中的成功应用表明，未知的经济数据可通过有限个已知的经济数据比较好的预测出来。据此，本文做如下假设：在式(2)中，当M和N都比较大时，向量为稀疏向量。即，向量u中非零元素的个数远小于向量u中元素的总数，但这些非零元素的个数、大小和位置仍未知。至此，本节将时间序列预测问题归结为：在u为稀疏向量的前提下求解式(4)的问题。压缩感知理论是求解稀疏信号有力工具，并且已经在信号处理领域当中取得了巨大的成功，其基本思想是在向量u为稀疏向量的前提下，式(2)存在唯一解，且该解可求。

2.1 算法原理

本节选用压缩感知理论中经典的正交匹配追踪(Orthogonal Matching Pursuit,OMP)为基础来求解上述问题。OMP算法属于贪婪算法的范畴，它的基本思想是，在对问题求解时，总是做出在当前看来是最好的选择，其数学原理见文献[2]。本文在OMP算法的基础上，针对式(4),提出的参数模型求解流程如下：

（1）初始化残差r=p、迭代次序k=0、支持集Λ为空集。

（2）计算残差r与矩阵R各列的相关系数，并找出相关系数最大的列对应的列标ik。

（3）首先将列标ik存入支持集中Λ=Λ∪{ik}，再将矩阵R中列标不属于支持集Λ的元素置为零得到RΛ，最后在新的系统矩阵RΛ上根据（5）式估计出最小二乘解ut：

（4）更新残差 r=p-RΛut。

（5）判断终止条件是否成立。成立，则转步骤6。不成立，迭代次序自增1：k=k+1，并转步骤2。

（7）选择置信度η，判断去均值后的残差数据r在置信度为η的条件下是否服从0均值方差为的正态分布。是，则输出模型参数向量ut；否，选择其他参数重新计算。

显然，OMP算法的基本思想是从矩阵R中选择与预测向量p最接近的列来逼近预测向量p并求出残差r，然后再从投影矩阵R中选择其他列向量进一步消除残差。根据式(5)可得：

可见，更新后观测矩阵RΛ的列向量始终与残差r=p-RΛut保持正交。残差r的能量必定不包含观测矩阵RΛ列向量的贡献。随着矩阵RΛ维数不断的扩大，残差r必然不断衰减。可见，该方法最大的特点是在残差最小的前提下，自适应的选择参数向量u中的非零元素，而不是像AR模型那样，事先指定非零元素的个数和位置。步骤5中的终止条件有多种选择方式，本文仅给出两种选择方式作为参考：

方式2：多次实验确定最佳终止条件，k≥P时停止，P通过多次实验确定。

具体选择何种方式可根据实际情况自行选择。除此之外，影响预测结果的还有矩阵R的行数M+1和列数N+1。其中，N主要影响非零元素的搜索范围，因此N应当在条件允许的情况下，选择一个足够大的数值。参数M主要影响方程组的个数，如果M选择过小，则模型参数较易受误差的影响；如果M选择过大，则同样可能导致模型参数的不精确。这主要是因为式(2)暗含一个假设，即所有方程均满足同一组模型参数u，如果M选择过大，这个假设可能不成立。在实际应用中，参数M应当根据经验在一定范围内逐次尝试选出最优值。

2.2 预测

为了验证本文方法的有效性，笔者在软件MATLAB(2010b)上实现了上述方法，并应用于对上证指数和我国第三产业值的预测。第一组实验选用1994年3月1日到2013年2月28日的日收盘价对上证指数进行预测。该实验主要用于说明本文算法的执行过程。为了给予参数向量u较大的选取空间，本文式（2）中的参数M和N分别选为49和199，对上证指数进行预测。预测结果与真实值的对比见图1（b）。关于步骤5的终止条件，本文选择第二种方式来确定。具体来说就是，针对观测到的数据进行预测，并逐次尝试迭代次数k的所有可能取值，选定预测值平均绝对百分比（Mean Absolute Percentage Error,MAPE）最低的迭代次数来进行预测。MAPE在每次迭代后的值见图1.（a）。可见，，随着k的增加，MAPE则表现为先下降后上升，在k=10时最小，因此对上证指数数据的预测中，本文推荐将步骤5中的迭代终止条件设定为迭代次序大于等于10。据此，对上证指数中的序列进行预测，预测结果与真实值见图2（b），此时的平均绝对百分比为2.00%。经检验，残差通过了显著性水平为0.05的高斯分布检验。

图1 对上证指数的预测

第二组实验选取我国第三产业产值1952～2011年的数据来进行预测。本文式（2）中的参数M取为6，N取为7。关于步骤5的终止条件，本文选择第二种方式来确定。具体来说就是逐次尝试迭代次数k的所有可能取值，选定预测值平均绝对百分比（Mean Absolute Percentage Error,MAPE）最低的迭代次数来进行预测。由图2（a）可知，随着迭代次数k的增加，MAPE与之呈现出正相关的关系，因此步骤5的终止条件为P=1，此时MAPE=0.0573，残差通过了显著性水平为0.05的高斯分布检验，预测结果与真实值见图2（b）。

图2 对我国第三产业产值的预测

在2007～2011年数据未知的情况下，表1给出了对我国第三产业产值样本数据进行预测得到的预测值与真实值的对比，经检验，残差同样也通过了显著性水平为0.05的高斯分布检验。另外还给出了本方法与ARMA（1,1）模型利用同样样本数据得到的2007～2011年的预测值和真实值及误差百分比。通过以上对比可知本文方法在预测精度上是优于ARMA模型的。最后，用本文方法对2012～2016年的我国第三产业产值进行预测，预测结果见表2。

表1 本文方法与ARMA模型预测结果对比 （单位：亿元）

表2 对2012年到2016年我国第三产业产值的预测 （单位：亿元）

3 结论

本文提出了一种基于压缩感知的时间序列预测方法，首先根据极大似然估计理论建立了数学模型，再在压缩感知的理论框架下将参数估计归结为一个压缩感知中的经典问题，进而采用正交匹配追踪法求解出参数，最后得到时间序列的预测值。利用该方法对我国第三产业产值的预测结果表明，本文方法与ARMA模型相比具有更高的预测精度，预测结果有很高的实用参考价值。将本文方法和AR模型做对比不难发现，如果将向量u的非零元素设定在前p个，则该方法将退化为p阶AR模型。也就是说，AR模型是本文方法的一个特例。本文方法的预测精度高于ARMA模型，另外由于对数据的平稳性，季节性等没有要求，因此也有着更广的适用范围。

[1]李子奈,叶阿忠.高级计量经济学[M].北京：清华大学出版社,2003.

[2]Tropp J,Gilbert A.Signal Recovery from Random Measurements Via Orthogonal Matching Pursuit[J].IEEE Transactions on Information Theory,2007,53(12).