自回归移动平均程
自回归移动平均程
————————————————————————————————作者：————————————————————————————————日期：
A . 自回归移动平均过程(),ARMA p q 理论部分
1．基本概念
(),ARMA p q 表达式为：
112211.......t t t p t p t t q t q Y c Y Y Y φφφεθεθε-----=++++++++ （1）
写成滞后算子的形式为：
()()2
1
2
1
1....1...p
q
p
t
q
t
L L L Y c L L φφφθθε----=++++ （2）
两侧同时除以()2121....p p L L L φφφ----，从而得到
()t t Y L μψε=+ （3）
其中
()()
()
1
2
1
2
1...1....q
q
p
p L L L L L
L
θθψφφφ+++=
----
()12/1....p c μφφφ=----
j
j ψ
∞
=<∞∑
从而可以发现，(),ARMA p q 过程的平稳性完全取决于回归参数()12,,...,p φφφ而与移动平均参数无关。即(),ARMA p q 过程的平稳性条件为特征方程：
2121....0p p z z z φφφ----=
的根在单位圆外。
（1）变形：
()()()112211.......t t t p t p t t q t q Y Y Y Y μφμφμφμεθεθε------=-+-++-++++ （4） 两边同时乘以()t j Y μ--，求期望得到自协方差。当j q >时，结果方程的形式p 阶自协方差形式：
1122....j j j p j p γφγφγφγ---=+++ 1,2,.....j q q =++ （5）
从而解为
1122....j j j j p p h h h γλλλ=+++ （6）
j q ≤时的自协方差函数比较复杂，并且不具有应用意义。不过(),ARMA p q 过程
的自相关函数都具有拖尾特征。
(),ARMA p q 过程容易出现的两个问题：
1）过度参数化问题。例如一个白噪声过程t t Y ε=也可以用()()11t t
L Y L ρρε-=-表示。此时无论ρ取何值，利用()()11t t L Y L ρρε-=-都能够很好的拟合数据，因此造成估计的困难。
2）(),ARMA p q 过程的表达式（54）的滞后多项式进行因式分解得到
()()()()()()()121211....111...1p t q t L L L Y L L L λλλμηηηε----=--- （7）
假设自回归算子()2121....p p L L L φφφ----和移动平均算子()11...q q L L θθ+++存在共同根（公因子），同时除以公因子，得到的过程()1,1ARMA p q --和原来的
(),ARMA p q 过程相同。
表1 时间序列模型性质表
模型 性质
AR(p)
MA(q)
ARMA(p,q)
模型方程 ()p t t L Y φε= ()t q t Y L θε= ()()p t q t L Y L φθε= 平稳条件
()0p z φ=的根在
单位圆外
无条件平稳
()0p z φ=的根在
单位圆外
可逆条件 无条件可逆 ()0q z θ=的根在
单位圆外 ()0q z θ=的根在
单位圆外 ACF 自相关 拖尾 在q 截尾 拖尾 PACF 偏自相关 在p 截尾
拖尾
拖尾
2.(),ARMA p q 的预测
2.1．预测原理（基于条件的预测）：
定义1：均方误差
对于任何预测都存在误差，我们需要给出一个损失函数来度量预测偏离一个特定的量的程度。假定一个二次损失函数，选择*1t t
Y +，使得
(
)2
*11t t t
E Y Y
++- （8）
最小。表达式（8）称为预测值*1t t
Y
+的均方误差，记做()()2
**1
11t t t
t t
MSE Y
E Y
Y
+++=-。
定理1：最小均方误差预测就是t X 条件下1t Y +的期望。 证明：
假定*1t t Y +为基于条件期望以外的其他函数()t g X 的预测()*1t t t
Y g X +=，其MSE 为：
()()
()()()()()()()()(){}
()()(){}
2
12
1112
2
1111112
2
1111 22t t t t t t t t t t t t t t t t t t t t t t t t t t t E Y g X E Y E Y X E Y X g X E Y E Y X E E Y X g X E Y E Y X Y X g X E Y E Y X E E Y X g X E η++++++++++++++-??=-+-??
????=-+-????
????+--????
????=-+-+???? （9）
因为在t X 的条件下，()1t t E Y X +与()t g X 都是常数，因此
(){}()(){}
11110
t t t t t t t t t E X E Y E Y X X E Y X g X η++++????
????=--????
= （10） 根据迭代期望法则，（10）的期望就是无条件期望，即
()(){}
110t t t t E E E X X ηη++??==?? （11）
从而，（9）变为
()()()()()2
2
2
1111t t t t t t t t E Y g X E Y E Y X E E Y X g X ++++????-=-+-???? （12）
右边第一项为常数，因此如果希望均方误差最小，只有：
()()1t t t E Y g +=X X （13）
定理得证。
定义2：线性投影
假设预测*1t t Y +为t X 的线性函数，即*1t t
Y +'=t a X 。如果存在一个α，使得预测误差()1t t Y X α+'-与t X ，即
()10t t t E Y X X α+'''-=???? （14）
则预测t X α'称为1t Y +关于t X 的线性投影。
定理2：在线性预测族中，线性投影具有最小均方误差。 证明和定理1相似。
线性投影是随机过程总体特征的归纳；而OLS 回归是对样本观察值的归纳。 定理3 多重投影定理
如果2t Y +的1t +期的预测是t 期信息的投影，则结果为2t Y +的t 期最小均方误差预测。
例1.()1,1ARMA 过程预测
解：()1,1ARMA 过程()()()11t t L Y L φμθε--=+，当1,1φθ<<时，满足平稳性和可逆性。因此预测为：
(
)()11?11t t s t s L L
Y Y L L L θφμμφθ++??+-=+-??-+?? 其中
()()()()()
()
22221221111...1...1...1s s s
s s s
s L L L L L L L L L L L L L
θφφφθφφφθφφφφθφφ+
+--??
+??-????
++++++??=+????=+++++=
- 从而预测为
()()
()
()111?111s
s s
s t t t s t
L
Y Y Y L
L
L
φθφφθφφμμμμφθθ--+++-=+
-=+
--++
对于2,3,...s =，预测服从递归算法：
()
1t s t t s t Y Y μφμ++-=+-))
即在一期以后，预测按几何方式以速度φ收敛于无条件均值μ。前一期的预测为：
()()11t t t t t Y Y Y L
φθμμφμθεθ++=+-=-++))
其中
()()111t t t t t t t Y Y Y Y εμφμθε---=----=-)
))
例2：(),ARMA p q 过程预测 解：对于(),ARMA p q 过程
()()()2
2
1
2
1
2
1...1...p
q
p
t
q
L L L Y L L L φφφμθθθ-----=++++
1期预测为
()()()
1211121... ...t t p t p t t t t q t q
Y Y Y Y μφμφμφμθεθεθε--+--=+-+-++-++++)
)))
其中1t t t t Y Y ε-=-))
。前s 期预测为：
()()()(
)
(
)
(
)
1212111212...1,2,..., ......1,2,...
p t s t t s t t s p t s t s t q t q
t s t
p t s t
t s t t s p t Y Y Y s q
Y Y Y Y s q q φμφμφμθεθεθεμφμφμφμ
+-+-+-+--++-+-+-?-+-++-?
=?
++++-=??-+-++-=++??)))
)))))
))
当s q >时，预测为由自回归系数决定的p 阶差分方程。 3. 沃尔德分解和ARMA 建模 3.1．沃尔德分解
定理4 Wold 分解定理：
任何零均值协方差平稳过程t Y 可表示成如下形式
0t j t j t j Y k ψε∞
-==+∑ （15）
其中01ψ=，且20j j ψ∞
=<∞∑。t ε是白噪声（新生量），表示以Y 的滞后项预测t Y 产生的误差：
()12
?,,...t t t t t Y E Y Y Y ε--≡- （16） ()12
?,,...t t t t k E k Y Y --= （17） 对于任意的j ，t k 的值与t j ε-无关。t k 由Y 的过去值确定，称为t Y 的线性确定性分量。0j t j j ψε∞
-=∑称为线性非确定性分量，若0t k =，该过程为纯线性不确定的。
Wold 分解定理仅依赖于Y 的稳定的二阶矩。因此描述了Y 的最优线性预测。
3.2．Box-Jenkins 建模 3.2.1.建模基本思想
将某个时间序列的SACF 和SPACF 的行为与各种理论ACF 和PACF 的行为匹配起来，挑选最佳匹配（或一组匹配的集合），估计模型的未知参数()2,,i i φθσ，并检查从模型拟和得到的残差，已发现可能的模型错误。
步骤：
1）变换数据，是数据满足协方差平稳性假设（单位根检验和季节调整）。 2）对序列的(),ARMA p q 过程的参数,p q 做一个初始的较小值猜测。 3）估计()L φ和()L θ的系数。
4）初步诊断分析。保证所得模型和数据特征相符。 3.2.2．样本自相关SACF
0/j j ργγ=)
)
)
（18）
其中
()()11T
j t t j t j y y y y T γ-=+=--∑)
0,1,2,...,1j T =- （19）
1
1T
t t y y T ==∑ （20）
由于实际上假定了协方差平稳性，因此当j →∞，总体自协方差趋向于零。
SACF 的检验统计量为：
()*
21k
i i Q k T r ==∑ （21）
其渐进分布服从自由度为k 的卡方X 分布，即()*
2~a
k Q k X 。
3.2.3．偏自相关函数SPACF
m 阶偏自相关系数的估计是y 关于常数项和最近m 个值的OLS 回归的最
末一个系数：
()()()11211?...m m m t t t m t m t y c y y y e ααα+--+=+++++)))) （22）
其中t e )
代表OLS 回归的残差。
3.2.4．选择模型的标准：（存在多个行为匹配的模型） 1） AIC 标准：（Akaike 信息标准）
()()21?,log 2AIC p q p q T σ
-=++ （23） 2） BIC 标准
()()21?,log log BIC p q p q T T σ
-=++ （24） 3）首先设定()B φ和()B θ的阶数上限，max p 和max q ，并规定{}max 0,1,....,p p =和
{}max 0,1,...,q q =，则选择的阶数1p 和1q 由法则确定：
()()11,min ,AIC p q AIC p q =或 ()()11,min ,BIC p q BIC p q = （25） 3.3(),ARMA p q 在Eviews 中的实现
(1) 通过自相关分析图判断平稳性：如果序列的自相关系数很快地趋于零，即落入随机区间，则时序是平稳的，否则是非平稳的。
(2)自相关图的实现：主菜单中选择quick/series Statisttics/correlogram ，在对话框中输入分析的序列名称。如index ，点击OK 弹出相关图定义。选择之后，点击OK ，从而得到时间序列的自相关和偏自相关分析图。 (3)根据相关图和偏自相关图判断自回归和移动平均的阶数。
(4)模型参数的估计方法： 在主窗口选择Quick/Estimate/Equation ，输入index ar(1) ar(2) …ar(p) ma(1) ma(2)…ma(q) 点击OK 进入。
(5)结果中要求AIC 和BIC 越小越好。而且最后两行的数值落在单位圆内。 (6) 模型的检验：
1）对模型的残差序列进行白噪声检验。检验残差序列的样本自相关系数是否为零。检验统计量为卡方检验。残差序列的自相关函数为
()121n
t t k
t k k n t
t e e r e e
-=+==
∑∑
1,2,....,k m = （26）
m 为最大滞后期。一般取/4n 。检验统计量为
()()
212m
k k r e Q n n n k
==+-∑ （27）
在零假设下，Q 服从卡方分布。给定置信度，如果()2
Q X m p q α≤--，则不能拒
绝残差序列相互独立的原假设，通过检验。否则拒绝原假设。直接对残差序列的检验，分析残差序列的自相关图。
2) 检查是否过度拟合。用高阶的模型进行拟合，并与原模型比较。
4.极大似然估计 4.1引言
考虑ARMA 模型：
112211.......t t t p t p t t q t q Y c Y Y Y φφφεθεθε-----=++++++++ （28）
其中()20,t WN εσ:。前面我们假定知道总体参数()211,,...,,,...,,p q c φφθθσ，此时利用过程（28）进行预测。
本章我们要研究在仅能观测到Y 的情况下，如何估计()2
1
1
,,...,,,...,,p
q
c φφθθσ。
估计方法为极大似然估计。令()2
1
1
,,...,,,...,,p
q
c φφθθσ=θ表示总体参数向量。假定我们观察到一个样本量为T 的样本()12,,...,T y y y 。计算所实现样本的联合概率密度函数：
()11,,...,11,,...,T T Y Y Y T T f y y y θ-- （29）
这可以看作是观察到样本发生的概率。使得“概率”最大的θ值就是最优估计。这种思想就是极大似然估计的思想。极大似然估计需要设定白噪声的分布。如果
t ε是高斯白噪声，则得到的函数为高斯似然函数。
极大似然估计的步骤： 1）计算似然函数（29）。
2）利用求极大值方法求使得函数值最大的θ值。 4.2 高斯(),ARMA p q 过程的似然函数
对于高斯(),ARMA p q 过程
112211.......t t t p t p t t q t q Y c Y Y Y φφφεθεθε-----=++++++++ （30）
其中()20,t iidN εσ:。总体参数向量为()21212,,,...,,,,...,,p q c φφφθθθσ=θ。
自回归过程的似然函数的近似以y 的初始值为条件，移动平均过程似然函数的近似以ε的初始值为条件。(),ARMA p q 过程以y 和ε的初始值为条件。
假设初始值()0011,,...,p y y y --+'=y 和()0011,,...,q εεε--+'=ε给定，则利用实现
{}12,,..,T y y y ，迭代得到：
11221122......t t t t p t p t t q t q y c y y y εφφφθεθεθε------=--------- （31）
可得1,2,....,t T =的序列{}12,,..,T εεε。则条件似然函数为：
()()
()()11001100,,...,,22
2
1ln ,,...,,;ln 2ln 222T T T T Y Y Y T
t t L f y y y T T θεπσσ
--===---∑Y εY εθ （32）
4.3极大似然估计的统计推断 4.3.1. 极大似然估计的标准差
如果样本量T 足够大，则极大似然估计θ)
近似表示为：
()11,N T ?--≈0θθ% （33）
其中0θ代表真实参数向量。矩阵?称为信息矩阵。信息矩阵的二阶导数估计为
()2
1
L T
θθ
θ?θθ-=?=-'
??)
)
（34）
其中()L θ为对数似然函数：
()()111ln ;T t T
t t Y t L f y θ--ψ==ψ∑θ （35）
1t -ψ表示t 时刻y 的所有历史观测值。利用数值方法可以计算出对数似然函数的
二阶导数。（34）代入（35），得到
(
)(
)
()1
2
0L E θθθθθθθθθ-???'--?-=??'????
)
)
) （36） 4.3.2．似然比检验
假设原假设：参数向量θ中存在m 个限制（例如某些系数等于零）。首先求
出无限制极大似然估计；在求出存在限制情况下的极大似然估计。令()
L θ)
表示
无限制对数似然函数。()L θ
%表示限制对数似然函数。明显()()
L L θθ≥)
%，检验统计量为：
()()
()22L L m θθχ??-??
)
%: （37） 利用显著性检验法和置信区间法可以对原假设进行检验。 4.3.3．拉格朗日乘子检验
标准差检验需要计算无限制极大似然估计θ)
。似然比检验既要计算有限制极大似然估计量，又要计算无限制极大似然估计量。而如果计算有限制极大似然估
计量比较简单，则可以利用拉格朗日乘子检验。拉格朗日乘子检验是从有限制极大似然估计这一原假设出发，即原假设：有限制极大似然估计量为真。
令θ为一个()1a ?向量，令θ)
为有m 个限制的极大似然估计量。令
()1,...;t t f y y θ-为第t 个观察值的条件密度，令(),t h Y θ表示对数条件密度对限制估计θ的导数组成的()1a ?向量：
(
)
()
1ln ,...,,t t t f y y h Y θθ
θθθ
-=?=
?)
)
（38）
此时拉格朗日乘子统计量为：
()
()
()1
12
1
1,,T
T t
t t t T h Y h Y m θ?θχ--=='????≈????????
∑∑)) （39） 其中?为信息矩阵，其估计量()2
1
L T
θθ
θ?θθ-=?=-'
??)
)
。
B. 实验部分（MATLAB ）
105页例2.2: 利用平稳序列的前5个自协方差建立ARMA(2,2)模型
1. 根据延伸的Yule-Walker 方程
31213
224a a γγγγγγ??????= ? ? ??????? 计算出1a 和2a ,得到)6265.0,0894.0(),(21-=a a ;
⒉ 利用公式
121121(),||2
k k k T y k k k k k k k a a k γγγγγγγγγγ++-+--??
?=≤ ? ?
??
计算出((0),(1),(2))y y y γγγ,得到
3.2729) 2.4287- 7.1278())2(),1(),0((=y y y γγγ;
⒊ 利用公式
222
1
()b A C γσ
=
-∏ 20T C C σγ=-∏ 其中
1
lim T k k
k
k -→∞∏=ΩΓΩ，0100A ??
= ???
， 10C ??= ???
，1
212k k k k
γγγγ+???
Ω= ???L L
计算出0119.42=σ和)8158.0,3334.0(),(21-=b b 。
⒋ 所要求的模型为
2
12
1*8158.0*3334.0*6265.0*0894.0----+-+-=t t t t t t X X X εεε,
t Z ∈
其中{}t ε是)0119.4,0(WN 。
k 6 12 20 30 40 >=51
1b
-0.3171 -0.3281 -0.3333 -0.3333 -0.3334 -0.3334 2b
0.7375 0.7986 0.8122 0.8154 0.8157 0.8158 2σ
4.4379
4.0981
4.0299
4.0139
4.0122
4.0119
程序1：
gamma1=[0.29 -1.06;0.69 0.29]; r=[0.69 -0.12]; a1=inv(gamma1)*r' a1 =
0.0894 -0.6265 程序2： a=[-1; a1];
gammaY=zeros(1,3);
for k=0:2
gamma2=[4.61 -1.06 0.29 0.69 -0.12]; gamma3=zeros(3,3);
gamma3(1,1)=gamma2(1,abs(k)+1); gamma3(1,2)=gamma2(1,abs(k+1)+1); gamma3(1,3)=gamma2(1,abs(k+2)+1); gamma3(2,1)=gamma2(1,abs(k-1)+1); gamma3(2,2)=gamma2(1,abs(k)+1); gamma3(2,3)=gamma2(1,abs(k+1)+1); gamma3(3,1)=gamma2(1,abs(k-2)+1); gamma3(3,2)=gamma2(1,abs(k-1)+1);
gamma3(3,3)=gamma2(1,abs(k)+1);
gammaY(k+1)=a'*gamma3*a;
end
gammaY
gammaY =
7.1278 -2.4287 3.2729
程序3。4：
A=[0 1;0 0;];C=[1;0];
gamma=[gammaY(1,2);gammaY(1,3)];
k=6;
Omega=zeros(2,k);
Omega(1,1)=gammaY(1,2);Omega(2,1)=gammaY(1,3);Omega(1,2)=gammaY(1,3); Gamma=zeros(k,k);
for i=1:k
Gamma(i,i)=gammaY(1,1);
end
for i=2:k
Gamma(i,i-1)=gammaY(1,2);Gamma(i-1,i)=gammaY(1,2);
end
for i=3:k
Gamma(i,i-2)=gammaY(1,3);Gamma(i-2,i)=gammaY(1,3);
end
Pai=Omega*inv(Gamma)*(Omega)';
sigma=gammaY(1,1)-C'*Pai*C;
b2=(gamma-A*Pai*C)/sigma;
sigma
b2
(4) %零点%
A=[0.6265,-0.0894,1]; roots(A)
ans =
0.0713 + 1.2614i
0.0713 - 1.2614i
B=[0.8158, -0.3334, 1]; roots(B)
ans =
0.2043 + 1.0881i
0.2043 - 1.0881i
0.0713 + 1.2614i A(z) 0.0713 - 1.2614i
1.2634
0.2043 + 1.0881i B(z) 0.2043 - 1.0881i
1.1072
C . 实验 ARMA 模型建模与预测（Eviews ）
一、实验目的
学会通过各种手段检验序列的平稳性；学会根据自相关系数和偏自相关系数来初步判断ARMA 模型的阶数p 和q ，学会利用最小二乘法等方法对ARMA 模型进行估计，学会利用信息准则对估计的ARMA 模型进行诊断，以及掌握利用ARMA 模型进行预测。掌握在实证研究中如何运用Eviews 软件进行ARMA 模型的识别、诊断、估计和预测和相关具体操作。
二、基本概念
宽平稳：序列的统计性质不随时间发生改变，只与时间间隔有关。
AR 模型：AR 模型也称为自回归模型。它的预测方式是通过过去的观测值和现在的干扰值的线性组合预测，自回归模型的数学公式为:
1122t t t p t p t y y y y φφφε---=++++L
式中: p 为自回归模型的阶数i φ（i=1，2， K ，p ）为模型的待定系数，t ε为误差， t y 为一个平稳时间序列。
MA 模型：MA 模型也称为滑动平均模型。它的预测方式是通过过去的干扰值和现在的干扰值的线性组合预测。滑动平均模型的数学公式为:
1122t t t t q t q y εθεθεθε---=----L
式中: q 为模型的阶数； j θ（j=1，2，K ，q ）为模型的待定系数；t ε为误差；
t y 为平稳时间序列。
ARMA 模型：自回归模型和滑动平均模型的组合，便构成了用于描述平稳随机过程的自回归滑动平均模型ARMA ， 数学公式为:
11221122t t t p t p t t t q t q y y y y φφφεθεθεθε------=++++----L L
三、实验内容及要求 1、实验内容：
（1）根据时序图判断序列的平稳性；
（2）观察相关图，初步确定移动平均阶数q 和自回归阶数p ；
（3）运用经典B-J 方法对某企业201个连续生产数据建立合适的ARMA （,p q ）模型，并能够利用此模型进行短期预测。 2、实验要求：
（1）深刻理解平稳性的要求以及ARMA 模型的建模思想；
（2）如何通过观察自相关，偏自相关系数及其图形，利用最小二乘法，以及信息准则建立合适的ARMA 模型；如何利用ARMA 模型进行预测； （3）熟练掌握相关Eviews 操作，读懂模型参数估计结果。
四、实验指导
1、模型识别
（1）数据录入
打开Eviews软件，选择“File”菜单中的“New--Workfile”选项，在“Workfile structure type”栏选择“Unstructured /Undated”，在“Date range”栏中输入数据个数201，点击ok，见图1，这样就建立了一个工作文件。
图1 建立工作文件窗口
点击File/Import，找到相应的Excel数据集，打开数据集，出现图2的窗口，在“Data order”选项中选择“By observation”即按照观察值顺序录入，第一个数据是从a2开始的，所以在“Upper-left data cell”中输入a2，本例只有一列数据，在“Names for series or number if named in file”中输入序列的名字production或1，点击ok，则录入了数据。
图2
（2）绘制序列时序图
双击序列production，点击view/Graph/line，则出现图3的序列时序图，
时序图看出201个连续生产的数据是平稳的，这个判断比较粗糙，需要用统计方
法进一步验证。
76
80
84
88
92
255075100125150175200
PRODUCTI ON
图3
（3）绘制序列相关图
双击序列production，点击view/Correlogram，出现图4，我们对原始数据序列做相关图，因此在“Correlogram of”对话框中选择“Level”即表示对原始序列做相关，在滞后阶数中选择14（201
??
??），点击ok，即出现相关图5。
图4
从相关图看出，自相关系数迅速衰减为0，说明序列平稳，但最后一列白噪声检验的Q统计量和相应的伴随概率表明序列存在相关性，因此序列为平稳非白噪声序列。我们可以对序列采用B-J方法建模研究。
图5
（4）ADF检验序列的平稳性
通过时序图和相关图判断序列是平稳的，我们通过统计检验来进一步证实这个结论，双击序列production，点击view/unit root test，出现图6的对话框，我们对序列本身进行检验，序列不存在明显的趋势，所以选择对常数项，不带趋势的模型进行检验，其他采用默认设置，点击ok，出现图7的检验结果，表明拒绝存在一个单位根的原假设，序列平稳。
图6
图7
（5）模型定阶
由图5看出，偏自相关系数在k=3后很快趋于0即3阶截尾，尝试拟合AR （3）；自相关系数在k=1处显著不为0，当k=2时在2倍标准差的置信带边缘，可以考虑拟合MA（1）或MA（2）；同时可以考虑ARMA（3，1）模型等。
在序列工作文件窗口点击View/Descriptive Statistics/Histogram and States对原序列做描述统计分析见图8，可见序列均值非0，我们通常对0均值平稳序列做建模分析，所以需要在原序列基础上生成一个新的0均值序列。点击主菜单Quick/Generate Series，在对话框中输入赋值语句Series x=production-84.11940，点击ok则生成新序列x，这个序列是0均值的平稳非