(整理)ARMA算法整理.

信息通信网络时序指标动态阈值选取方法研究
整篇文章分为三部分，第一部分是点预测，第二部分是阈值d 的选取，第三部分是结合两者进行区间预测。

其中第一部分是重点，分为两个小部分，分别为前期预处理检验和模型建立，模型建立部分又分别由6个小部分组成。

一、建立模型进行中心点预测
思路：根据给出的数据序列，利用自相关系数，偏相关系数的性质，选择合适的模型进行模拟，如AR 模型（AR(p)），MA 模型（MA(q)），ARMA 模型（ARMA(p ，q)），并确定它们的阶数。

然后估计模型中未知参数的值，并利用AIC 准则来进行模型优化，从而可以对未来数据进行预测。

注：）（p AR 定义：⎪⎪⎩⎪⎪⎨⎧∀=≠===≠+++++=---t s x E t s E Var E x x x x t s s t t t p
t
p t p t t t ,0(,0)(,)(,0)(02
22110）
εεεσεεφεφφφφε )(q MA 定义：⎪⎩
⎪
⎨⎧≠===≠----+=---t
s E Var E x s t t t q q t q t t t t ,0)(,)(,0)(022211εεσεεθεθεθεθεμε
),(q p A R M A 定义：⎪⎪⎩⎪⎪⎨⎧∀=≠===≠≠---+++++=-----t
s x E t s E Var E x x x x t s s t t t q
p
q t q t t p t p t t t ,0)(,0)(,)(,0)(0,02
1122110εεεσεεθφεθεθεφφφφε 建模： 1. 前提准备
得到数据之后（比如移动公司一个月的通话时长记录），我们要对数据进行一个预处理，判定给出的数据满足为平稳非白噪声序列，
才可以利用上述几种模型对该数据序列进行建模。

（这是一个前提条件，所以，这就要求我们在选取数据时要有意识的控制）
(1) 平稳性检验
这里我们利用时序图检验的方法进行平稳性检验。

所谓时序图，就是一个平面二维坐标图，通常横轴表示时间，纵轴表示序列取值。

平稳序列的时序图应该显示出该序列始终在一个常数值附近随机波动，而且波动的范围有界的特点。

如果观察序列的时序图显示出明显的趋势性或周期性，那它通常不是平稳序列。

这样，我们根据时序图就可以判断是否是平稳的。

(2) 纯随机性检验
纯随机性序列：如果序列值彼此之间没有任何相关性，那就意味着该序列是一个没有记忆的序列，过去的行为对将来的发展没有丝毫影
响，这种序列称之为纯随机序列。

白噪声序列是典型的纯随机序列。

这里我们要验证我们所要研究的数据序列不是纯随机序列，即过去的值对现在有影响，才能建立ARMA 模型，从而进行预测。

方法：利用LB 统计量
∑=-+=m
k k k
n n n LB 12)ˆ(
)2(ρ，式中n 为序列观测期数，m 为指定延迟期数，k ρ
ˆ为自相关系数（当前x 与k 期前x 的相关系数）。

且
）（m k
n n n LB m
k k 2
12~)ˆ(
)2(χρ⋅-+=∑=。

一般只要计算出来延迟6期和延迟12期
的LB 及所对应的P 值就可以判断序列的随机性。

如果计算结果P 值很小，基本上以0.05为标准，只要小于0.05即可断定该序列不是纯随机序列，属于非白噪声序列。

其中P 值计算方法：k 自由度的）（
k 2
χ函数的密度函数为⎪⎪
⎩⎪⎪⎨⎧Γ=--其他
,00,)2(21
)(2122 x x k x f x
k k
，对)(x f 进行积分，代入之前计算的LB 计
算值，得到⎰-=LB
dx x f p 0
)(1。

2. 建模步骤
假如我们的观察值序列通过序列预处理，可以判定为平稳非白噪声序列，那么我们就可以利用模型对该序列建模。

建模的基本步骤如图所示：
（1）求出该观察值序列的样本自相关系数（ACF ）和样本偏自相关系数（PACF ）的值。

（2）根据样本自相关系数和偏自相关系数的性质，选择阶数适当的
),(q p ARMA 模型进行拟合，即确定p,q 的值。

（3）估计模型中未知参数的值。

（4）检验模型的有效性。

如果拟合模型通不过检验，专享步骤（2），重新选择模型再拟合。

（5）模型优化。

充分考虑各种可能，建立多个拟合模型，从所有通过检验的拟合模型中选择最优模型。

（6）利用拟合模型，预测序列的将来走势。

N
模型检验
Y
平稳，非白噪声序列
计算样本自相关系数（ACF ）和偏自相关系数（PACF ） ARMA 模型识别
估计模型中未知参数的值
模型优化
2.1 计算样本自相关系数和偏自相关系数
因为我们是通过考察平稳序列样本自相关系数和偏自相关系数的性质来选择合适的模型拟合观察值序列，所以模型拟合的第一步是
要根据观察值序列的取值求出该序列的样本自相关系数n}k ,0ˆ{k ρ
和样本偏自相关系数n}k ,0ˆ{kk φ的值。

样本自相关系数可以根据以下公式求得：
n k x x
x x x x
n
t t
k
n t k t t
k 0,)()
)((ˆ1
2
1
∀---=∑∑=-=+ρ
样本偏自相关系数可以利用样本自相关系数的值，根据以下公式求
得：n k D
D k kk 0,ˆˆˆ∀=φ
式中，1
ˆˆˆ1ˆˆˆ1ˆ212111
----=k k k k D
ρρρ
ρρ
ρ，k k k k
D ρ
ρρ
ρρρ
ρˆˆˆˆ1ˆˆˆ1ˆ212111
--=
其中，k D ˆ是将D ˆ的第k 列变为⎪⎪⎪⎪⎪⎭
⎫ ⎝⎛k ρρρ
ˆˆˆ21 。

2.2 模型识别（计算p, q ）
计算出样本自相关系数和偏自相关系数的值之后，就要根据他们表现出来的性质，选择适当的ARMA 模型拟合观察值序列。

这个过程实际上就是要根据样本自相关系数和偏自相关系数的性质估计自
预测序列将来的走势
相关阶数p
ˆ和移动平均阶数q ˆ，因此，模型识别过程也称为模型定阶过程。

理论依据：（ARMA 模型定阶的基本原则）
k ρ
ˆ kk φˆ 模型定阶
拖尾
q 阶截尾
拖尾 p 阶截尾
拖尾 拖尾
)(p AR 模型 )(q MA 模型
),(q p ARMA 模型
方法：下面我们考虑，在实际操作中，怎样判定是截尾或拖尾。

即当样本自相关系数或偏自相关系数在延迟若干阶之后衰减为小值波动时，什么情况下该看做相关系数截尾，什么情况下该看做相关系数在延迟若干阶之后正常衰减到零值附近做拖尾波动。

由于当样本容量n 充分大时，样本自相关系数近似服从正态分布：
)1,0(~ˆn
N k ⋅ρ
样本偏自相关系数也同样近似服从这个正态分布：)1
,
0(~ˆn
N kk ⋅φ 根据正态分布的性质，有
95.0)2ˆ2Pr(95.0)2
ˆ2Pr(≥≤≤-≥≤≤-
n n n
n kk k φρ
注：
)
21ˆ2Pr()2ˆ2
Pr(≤≤
-=≤≤-
n
k n
k n
ρρ
95
.04773.022
222
21=⨯=⎰--
=
dx x π
所以可以利用2倍标准差范围
）（n
n 2,2-辅助判断。

如果样本自相关系数或偏自相关系数在最初的d 阶明显大于2倍标准差范围，而后几乎95%的自相关系数都落在2倍标准差的范围以内，而且由非零相关系数衰减为小值波动的过程非常突然，这时通常视为相关系数截尾，且截尾阶数为d 。

如果有超过5%的样本自相关系数落入2倍标准差范围之外，或者是由显著非零的相关系数衰减为小值波动的过程比较缓慢或者非常连续，这时通常视为相关系数不截尾，即拖尾。

这样，我们就可以根据最初的理论依据和计算的比较结果来选择合适的模型进行模拟。

（1）若观察序列的自相关系数拖尾，偏自相关系数截尾，且截尾阶数为p ，则选）（p AR 模型。

（2）若观察序列的自相关系数截尾，且截尾阶数为q ，偏自相关系数拖尾，则选）（q MA 模型。

（3）若观察序列的自相关系数和偏自相关系数均拖尾，则选
),(q p ARMA 模型。

2.3 参数估计
选择好拟合模型之后，下一步就是要利用序列的观察值确定该模型的口径，即估计模型中未知参数的值。

比如：)2(AR 模型t t t t x x x εφφ++=--2211中未知参数21φφ，的估计为
212
1221121ˆ1ˆˆˆ,ˆˆ1ˆ1ˆρ
ρρφρρρφ--=--=。

)1(MA 模型11--=t t t x εθε中未知参数1θ的估计为1211
ˆ2ˆ411ˆρ
ρθ-+-=。

)11(，ARMA 模型1111---+=t t t t x x εθεφ中未知参数，1φ1θ的估计为121ˆˆˆρ
ρφ=，11221221ˆˆˆ2ˆ1,2,2
4
2,24
ˆρφρφθ---=⎪⎪⎩
⎪⎪⎨⎧≥---≤-+
=c c c c c c c 当然，这种低阶的ARMA 模型可以根据公式直接带入求参数估计即可，但当阶数比较大时，参数估计值会非常复杂，这时候我们采用最大似然估计或最小二乘估计方法。

极大似然估计：
q t q t t p t p t t x x x -------+++=εθεθεφφ 1111
记 ),(~
1n x x x = '
11,,~
）
，，（q p θθφφβ = 最终，求解似然方程组
⎪⎪⎩⎪⎪⎨
⎧=∂∂+∂Ω∂=∂∂=-=∂∂0~2)~(1~ln 21)~;~(~02)~(2)~
;~(24
22ββσβββ
σβσβσεεεε
S x l S n x l 式中，]~'~[21
ln 2
1)ln(2)2ln(2)~
;~
(12
2x x n n x l -Ω-Ω---=ε
εσσπβ x x S ~'~
)~
(1-Ω=β ⎪⎪⎪⎪
⎪⎭
⎫
⎝⎛=Ω∑∑∑∑∞
=∞=-+∞
--+∞-02
1010
2
i i i n i i i n i i i i
G G G G G G
，i G 为Green 系数. 但是，由于)~
(βS 和Ωln 都不是β~
的显式表达式，因而上述似然方程组
实际上是由1++q p 个超越方程构成，通常需要经过复杂的迭代算法才能求出未知参数的极大似然估计值。

最小二乘估计：
q
t q t t p t p t t x x x -------+++=εθεθεφφ 1111
记 '11,,~
）
，，（q p θθφφβ =
q t q t p t p t t x x F -------++=εθεθφφβ 1111)~
( 残差项为：)~
(βεt t t F x -=
残差平方和为：∑==n
t t Q 12)~
(εβ
∑=----+++---=n
t q t q t q p t p t t x x x 1
2111)(εθεθφφ
使残差平方和达到最小的那组参数即为β~
的最小二乘估计值。

同极大似然估计一样，由于)~(βQ 不是β~
的显性函数，未知参数的最小二乘估计值通常也得借助迭代法求出。

网上有极大似然法和最小二乘法的程序代码，有的是matlab 的。

2.4 模型检验
模型的检验主要是检验模型的有效性，一个好的拟合模型应该能够提取观察值序列中几乎所有的样本相关信息。

换言之，拟合残差项
（真实值与拟合值的差）中将不再蕴含任何相关信息，即残差序列应该为白噪声序列。

这样，我们就通过检验残差序列是否为白噪声序列来说明所建模型是否通过检验。

方法与上面的序列随机性检验方法一样，利用LB 统计量，
）（m k
n n n LB m
k k 2
12~)ˆ(
)2(χρ⋅-+=∑=，不过这里要验证是属于白噪声序列，所
以最后计算结果P 值要大于0.05才可断定残差序列是纯随机序列，进而说明该拟合模型通过了检验。

2.5 模型优化
对于给定的观察值序列，我们能够建立多个拟合模型，模型优化就是从所有通过检验的拟合模型中选择最优模型。

在这里，我们利用AIC 准则（最小信息量准则）来选择最优模型。

AIC 准则的思想是认为一个拟合模型的好坏可以从两方面去考察：一方面是用来衡量拟合程度的似然函数值；另一方面是模型中未知参数的个数。

一个好的拟合模型应该是一个拟合精度和未知参数个数的综合最优配置，使得AIC 最小的ARMA 模型为最优模型。

),(q p ARMA 模型的AIC 函数为：
)2(2)ˆln(2+++=q p n AIC εσ
，其中)(ˆ2t Var εσε=n
n
i i
∑=-=1
2
)(εε。

2.6序列预测
通过以上所有的步骤，我们得到了最优的拟合模型，然后可以对未来的某一点进行预测。

二、阈值的选取
第一部分的模型得到的预测结果是一个点值，然而在现实生活中，由于环境的复杂性和因素的不可控制性，预测结果应该是一个以预测值为中心点上下浮动的区间更为合适。

这里，我们利用聚类分析的思想来确定区间的大小，即阈值的选取。

方法：先利用K-均值聚类方法将同一时间点的不同y值分为合适的k类，然后再根据最小的类均值
最大的类均值
阈值-
=来确定阈值的大小。

K-均值聚类步骤：
（1）在n个样品中随机选择k个样品作为初始凝聚点，（或者将所有样品分成k个初始类，然后将这k个类的均值作为初始凝聚点）。

（2）对除凝聚点之外的所有样品样品逐个归类，将每个样品归入离它最近的凝聚点所在的类（采用欧氏距离计算），然后将该类的凝聚点更新为这一类目前的均值。

（3）重复步骤（2），直到每个聚类不再发生变化或满足某个终止条件为止。

一般的迭代终止条件为误差平方和最小或限定聚类次数。

误差平方和为∑
=-
=
n
i
k
i
i x
J
1
2
)(
μ，其中有n个数据，分为k类，)(i
k
μ为i x
所在类的均值。

关于K-均值聚类，网上有程序代码。

三、模型结合动态阈值进行更准确预测
若要实现将来某一时刻的区间预测，由第一步的中间点预测加上第二步计算的所在时间点的上下浮动区间，即可得到该时刻的区间预测。

总结：第一步的整体思路和作者一样，斜率概念没用到；第二步没有用遗忘曲线知识，直接对数据进行了聚类。

另外，实际中移动电话时长记录数据应该具有很强的周期性，利用非平稳序列（ARIMA模型）先消除季节性变化影响更好一些，但却不好建模，不知道作者怎么考虑的。

这些可以与客户再进一步沟通。