平稳时间序列的模型的计算

目 录

摘要„„„„„„„„„„„„„„„„„„„„„„„1

第一章 绪论„„„„„„„„„„„„„„„„„„„2

1.1 时间序列模型的发展及其作用„„„„„„„„„„„2

1.2 什么是时间序列模型„„„„„„„„„„„„„„„2

1.3 本文研究的主要方法和手段„„„„„„„„„„„„2

1.4 本文主要研究思路及内容安排„„„„„„„„„„„2

第二章 ARMA模型„„„„„„„„„„„„„„„„„4

2.1 ARMA模型的基本原理„„„„„„„„„„„„„„„4

2.2 样本自协方差函数、自相关函数和偏相关函数„„„„4

2.3 ARMA模型识别方法„„„„„„„„„„„„„„„„5

2.4 模型参数估计„„„„„„„„„„„„„„„„„„6

第三章 实例分析„„„„„„„„„„„„„„„„„7

3.1 题目„„„„„„„„„„„„„„„„„„„„„„7

3.2 问题分析„„„„„„„„„„„„„„„„„„„„7

3.3 问题求解„„„„„„„„„„„„„„„„„„„„8

3.3.1数据的观测„„„„„„„„„„„„„„„„8

3.3.2数据处理„„„„„„„„„„„„„„„„„8

3.3.3求解自相关和偏相关函数„„„„„„„„„„8

3.4 模型的识别及求解„„„„„„„„„„„„„„„„9

3.5 结论„„„„„„„„„„„„„„„„„„„„„„11

参考文献„„„„„„„„„„„„„„„„„„„„„12 附录„„„„„„„„„„„„„„„„„„„„„„„12 评阅书„„„„„„„„„„„„„„„„„„„„„„15

《随机过程 》 课程设计任务书

摘 要

ARMA模型是研究时间序列的重要方法，由自回归模型（简称AR模型）与

滑动平均模型（简称MA模型）为基础“混合”构成。ARMA模型广泛应用在经济、工程等各个领域得益于其在具体预测方面的优势。在许多方面用该模型所作出的预测比其他传统经济计量方法更加精确。平稳时间序列模型主要有自回归模型(AR)、滑动平均模型(MA)和自回归滑动平均模型(ARMA)等，这些线性模型考虑因素较简单。自回归滑动平均模型(ARMA)计算简单，易于实时更新数据。

本文描述了ARMA模型的原理、自相关函数和偏相关函数的计算过程、模型

的识别方法以及ARMA模型的计算过程。并给出一组平稳时间序列的数据，对数据进行分析和处理，求出自相关系数和偏相关，并利用MATLAB软件画出自相关系数和偏相关图形，有图可知它们都是拖尾的，因此可以确定是ARMA(p,q)模型。接下来就是确定ARMA(p,q)的阶数，本文采用了AIC准则确定模型的阶数，在实际问题中，为使线性模型简单起见，通常p与q的数值被取得较小，却需都不为零。确定阶数后，就用我们学过的求解方法解出未知的参数，这样我们就得到了混合模型的表达式。

关键字：ARMA(p,q)模型，自相关函数，偏相关函数

第一章 绪 论

1.1 时间序列模型的发展及其作用

人们的一切活动，其根本目的无不在于认识和改造客观世界。时间序列分析

不仅可以从数量上揭示某一现象的发展变化规律或动态的角度刻画某一现象与其他现象之间的内在关系及其变化规律性，达到认识客观世界之目的，而且运用时间序列模型还可以预测和控制现象的未来行为，修正或重新设计系统以达到利用和改造客观之目的。时间序列分析在工程技术中有重要的作用，常用于做预报、控制等。为建立其随机线性模型，首先，我们应明白什么是时间序列：时间序列是随机序列，即参数离散的随机过程。由于工程中遇到的随机序列的参数经常为时间，故称随机序列为随机时间序列，简称时间序列。可以说，时间序列是随时间改变而随机变化的序列。平稳时间序列是平稳序列，它满足期望为零，且任意两个时刻的相关函数与时间t无关，仅与两个时刻的时间差相关。

1.2 什么是时间序列模型

所谓的时间序列模型就是用观测得到的若干数据，即样本数据的基础上，建

立符合事物发展规律的数学理论模型，就是根据时间序列Xt的一段样本观测数据Xt（1tN）所包含的信息，利用自协方差函数Rr和偏相关函数kk具有的“截尾”与“拖尾”性质，判定模型结构及阶数，时间序列的模型我们主要学习了AR(p)，MA(q)及混合模型ARMA(p,q)。 

1.3 本文研究的主要方法和手段

求解模型就是绘制出这些数据的图形观察图形的走向。通常该图形被分为两

种：平稳和非平稳图形。对于非平稳的图形，我们要用所学过的方法把这些非平稳的数据变成与之对应的平稳的数据，再用平稳的方法去研究它。对于平稳的图形，我们首先根据理论方法及数学软件，求解出自相关系数和偏相关系数，根据得出的自相关系数和偏相关系数分别绘制出图形，观察图形截尾还是拖尾，因此可以确定是哪类模型。接下来就是确定模型的阶数。确定阶数后即可求解出相关系数，这样就得到我们所要的模型表达式。

1.4 本文主要研究思路及内容安排

第一章介绍了时间序列模型的发展及其作用，什么是时间序列模型，课题研

究工作和论文的内容安排。

第二章介绍ARMA(p,q)的原理，自相关函数和偏相关函数，时间序列模型

模型的识别方法以及ARMA(p,q)的求解。

第三章利用具体问题，对问题分析，判断出符合ARMA(p,q)模型和并求解

ARMA(p,q)模型，并对文章进行总结。

第二章 ARMA模型

2.1 ARMA模型的基本原理

ARMA(p,q)模型：如果时间序列Xt满足：

Xt1Xt1pXtpat1at1qatq

则称时间序列为Xt服从(p,q)阶自回归滑动平均混合模型。或者记为：

(B)Xt(B)at，

特殊情况：q=0，模型即为AR(p)，p=0，模型即为MA(q)。

将预测指标随时间推移而形成的数据序列看作是一个随机序列，这一组随机过程所具有的依存关系体现着原始数据在时间上的延续性。一方面，影响因素的影响，另一方面，又有自身变动规律，假定影响因素为x1，x2，，xt由回归分析：

其中Y是预测对象的观测值， e为误差。作为预测对象Yt受到自身变化的

影响，其规律可由下式体现：

误差项在不同时期具有依存关系，由下式表示：

由此，获得ARMA模型表达式：

Yt01xt12xt2pxtp01t12t2qtqt

2.2 样本自协方差函数、自相关函数和偏相关函数

平稳序列,2,1,0,1,2， E(t)0，对于样本，定义自协方差函数：

1nkˆkjkj， nnj1ˆkˆk/ˆ0。同时为了保证ˆkk，ˆkk一般取则自相关函数为：

n50,kn/。常取4kn/10。 1k12k2nkn

为求出样本的偏相关函数，需要求解如下的Yule-Walker方程：

12k11k11111

22k22



1kkkk1211

ˆ1ˆ1ˆ2k1ˆˆˆ11k21ˆˆ12

ˆˆ2ˆkkk1ˆk1ˆ1ˆ111解得ˆ1ˆ2 ˆk

2.3 ARMA模型识别方法

通常平稳时间序列Zt，t0,1仅进行有限n次测量(n50)，得到一个样

1ntZt，本函数，且利用平稳序列各态历经性：Zj做变换，t1,n，nj1

将Z1,,Zn样本换算成为样本1,,n，然后再确定平稳时间序列{t,t0,1}的随机线性模型。

按样本自相关函数和样本偏相关函数的值分别作出点图，按“截尾”，“拖尾”

情况，查表确定模型的类别与阶数p，q。

理论上，线性模型的类别的确定可以根据kk和k的拖尾性和截尾性来判

别。但是在实际应用中，我们常用有一个样本算出k=k，kk=kk判别k，kk是拖尾还是截尾的。随机线性模型的三种形式：AR(p) ，MA(q)，ARMA(p,q)的判别分别如下：

（1）若k拖尾，kk截尾在kp处，则线性模型为AR(p)模型。k拖尾可

以用的点图判断，只要样本自相关函数的绝对值愈变愈小（k增大时）；但是，用样本偏相关函数判断kk 截尾在k=p处应该如何作呢？因为k>p时，样本偏相

关函数不为零，而偏相关函数kk0，这就给判断带来一定的困难。我们可以采

用一下方法解决：当kp时，平均20个样本偏相关函数中至多有一个使|kk| ≥

2/则认为kk截尾在kp处。

（2）若kk拖尾，k在kp处截尾，那么线性模型为MA(q)滑动平均模型。

^^^

kk拖尾可以根据样本偏相关函数的点图判断，只要|kk|愈变愈小（k增大时）。

但是，用样本自相关函数判断自相关函数k在kp处截尾可采用如下方法：当

kp时，若平均20个样本自相关函数中至多有一个使|k|≥

^^

（3)若样本自相关函数和样本偏相关函数都是拖尾的，则线性模型可以看成混和模型（或者，样本自相关函数和样本偏相关函数都不为截尾的，又被负指数型的数列所控制）。其具体的判别方法和上述一样。如下表所示：

2.4 模型参数估计

1、AR(p)模型参数估计：

2AR(p)模型有p2个参数：p,1,2,,p,。利用Yule-Walker方程，利

用Toeplitz矩阵求逆和作矩阵乘法的方法算样本偏相关函数kk。AR(p)模型的参数值不必作专门的计算，只要在样本偏相关函数计算的记录中取出样本参数值

0jj。即可。此时1,2,,p,都已经确定了，经过推理我们可以得到：

j12、MA(q)滑动平均模型参数估计： 2p

2ˆ2ˆ2ˆ2),k0ˆ(112q ˆk2ˆˆˆˆˆˆ(k1k+1qkq),1kq

ˆ,ˆˆ,ˆ2，即解这个非线性方程组。 可得q1个方程，求12q

3、ARMA(p,q)混和模型参数估计

对于满足一个条件：t1t1...ptpat1at1...patq采用先计算

2ˆ,ˆ,,ˆ，在计算ˆ,ˆˆ,ˆ的方法，具体如下：1）可利用Toeplitz矩阵和作矩阵12p12q

ˆ,ˆ,,ˆ。2）令...混和模型化为：乘法的方法求出12ptt1t1ptp'

t'at1at1...patq这是关于t'的MA(q)模型，用t'的样本协方差函数估计

2ˆˆ1,ˆ2ˆq,的值。

第三章 实例分析

3.1 题目

某简谐振动位移离平衡点的距离zt的观测数据如下所示，N=64，建立表示测

数据序列的时间序列模型。

表-1 原始数据

-0.17007 -2.05056 -4.18035 0.11153 2.17675 4.08507 0.24983 -4.92336 -8.15675 -4.42107 -2.64058 -0.52570 -0.22015 0.36515

2.61949 1.53203 -0.77529 -3.25439 -0.41904 3.92127 2.97641 1.05303 -3.07748 -3.41256 -5.37570 -3.53729 -0.79940 -1.07239

0.78588 -2.48208 -4.03418 0.16420 -0.39307 3.88025 6.03561

5.70750 1.45257 -2.97501 -3.16561 -3.15184 -5.66101 -4.55502 0.10608 0.04231 0.79915 -0.12247 -5.97528 -9.16058 -9.73835 -5.22347 0.62229 4.34356 1.32363 -1.25405 -3.90467 -3.02789

2.78487 7.61065 6.14580 3.51779 1.79635 0.74996 -0.19133 -1.40428

3.2 问题分析

本问题是让我们建立时间序列模型，题中给出了64个数据首先用所给数据

绘制图形。根据数据绘制出的图形如下（程序见附录1）：

[**************]

图（1）

由图形可以看出所给数据本身就是稳定的时间序列，接下来用MATLAB软件求解出自相关系数和偏相关系数，根据得出的自相关系数和偏相关系数分别绘制出图形，观察图形截尾还是拖尾，由图可知自相关系数和偏相关图形都拖尾，因此可以确定是ARMA(p,q)模型。接下来就是确定ARMA(p,q)的阶数。在实际问题中，为使线性模型简单起见，通常p与q的数值被取得较小，却需都不为零。确定阶数后，就用我们学过的求解方法解出未知的参数，这样我们就得到了ARMA(p,q)模型的表达式。

3.3 问题求解

3.3.1 数据的观测：

对一个时间序列作n次测量得到一个样本Z,Z,,Z，一般取n50。某简12n

谐振动位移离平衡点的距离zt的观测数据，得到n=64个样本值Z1,Z2,,Z60，（列于表1中）。表1的Z栏，共64个数据。将原始样本数据经过处理后变成时间t

序列W。 t

3.3.2 数据处理：

（1）求取误差时间序列Zt：

164

用MATLAB软件对原始数据进行处理，做变换WtZtZ，其中ZZt，代64t1入表1数据得：，1Z[(0.17007)(2.05056)(1.40428)]0.757464

（程序见附录2） WtZt(0.7574)，

3.3.3 数据处理：

ˆkˆk/ˆ0，（2）计算样本自协方差函数k，样本自方差函数k。 其中k0,1,2,3,4,5，

ˆk1k12k2nknn1nkjkj。由图-3数据可得：随着k的增大，k越nj1

来越小，具有拖尾性。

（3）接下来计算偏相关函数kk(k1)。利用Yule-Walker方程，

12k11k1111122k22



1kkkk1211

求得偏相关函数。

由机械振动位移图可看出这些数据画出的图形是平稳序列，由MATLAB软件

运行以上数据（程序见附录），得到自相关系数和偏相关系数如下表:

3.4 模型的识别及求解

ˆkk值，绘制ˆkk的曲线趋势来判断，ˆK值和偏相关函数ˆK和 根据自相关函数

分别画自相关系数和偏相关系数点图如下（程序分别见附录[3]，[4]）： 自相关函数图：

1

0.8

0.6

0.4

0.2

-0.2

-0.4

010

[1**********]0

图（3）

偏相关函数图如图：

0.80.6

0.40.20-0.2-0.4-0.6-0.8

[**************]

图（4）

Sample Autocorrelation Function (ACF)Sample Partial Autocorrelation Function

05

Lag

1015

Sample Partial Autocorrelations

Sample Autocorrelation

05

Lag

1015

图（5）

图形可以看出自相关函数和偏相关函数都拖尾，所以该数据符合ARMA(p,q)混合模型。用ARMA(p,q)模型对该数据进行拟合，即:

Xt1Xt1pXtpat1at1qatq

在计算机上，利用图（6）所示的AIC准则对候选模型参数进行计算和分析：（程序见附录5）

图（6）

根据计算知当p=2，q =2时的ARMA（2，2）模型，1 =1.637754，

2

2=0.7523818，1=-0.01194535，2=0.2609289，a=1.206746，AIC=30.7927

（最小）。根据模型定阶的AIC准则，比较上述计算结果可知，最能描述上述机械振动观测数据序列的模型结构应该是ARMA（2，2）模型，即：

Xt1.6378Xt10.7524Xt2at0.0119at10.26at2

3.5 结论

由以上问题的求解过程可知，所有的时间序列模型参数都是建立在样本统计参数的基础上，通过对数据处理，利用平稳时间序列的数据，得出的ARMA(p,q)模型，随机过程论目前已得到广泛的应用，在诸如天气预报、统计物理、天体物理、运筹决策、经济数学、安全科学、人口理论、可靠性及计算机科学等很多领域都要经常用到随机过程的理论来建立数学模型。时间序列分析也广泛应用到了解决许多问题上，如信号处理、金融预测、状态估计、控制和模式识别等.相信随着时间的发展，随机过程会越来越多的应用于我们的生活中，为我们带来更多的惊喜和发现。

参考文献

[2] 张善文 雷英杰 冯有前. MATLAB在时间序列分析中的应用 西安：西安电子科技大学出版社 2007 。

[2] 汪荣鑫编. 随机过程 西安：西安交通大学出版社 2006 。 [3] 刘卫国 Matlab程序设计与应用 北京 高等教育出版社 2006 。 [4] 田铮 秦超英 随机过程与应用 北京 科学出版社 2007 。

附录

附录1

z=[-0.17007 -2.05056 -4.18035 0.11153 2.17675 4.08507 0.24983 -4.92336 -8.15675 -4.42107 -2.64058 -0.52570 -0.22015 0.36515 2.61949 1.53203 -0.77529 -3.25439 -0.41904 3.92127 2.97641 1.05303 -3.07748 -3.41256 -5.37570 -3.53729 -0.79940 -1.07239 0.78588 -2.48208 -4.03418 0.16420 -0.39307 3.88025 6.03561 5.70750 1.45257 -2.97501 -3.16561 -3.15184 -5.66101 -4.55502 0.10608 0.04231 0.79915 -0.12247 -5.97528 -9.16058 -9.73835 -5.22347 0.62229 4.34356 1.32363 -1.25405 -3.90467 -3.02789 2.78487 7.61065 6.14580 3.51779 1.79635 0.74996 -0.19133 -1.40428];

>> x=[0:1:63]; plot(x,z) 附录2：

b=length(z);

m=sum(z)/b „„„„„„计算均值 w= z-m;

sumw=zeros(1,6); sumw1=0; for j=1:64

sumw1=sumw1+w(j)^2; „„„„计算0

sumw1 for k=0:5 for i=1:(b-k)

sumw(k+1)=sumw(k+1)+w(i)*w(i+k); .........计算k end

sumw(k+1) r=sumw/b; r0=sumw1/b;

p=r/r0 „„„„„„计算自相关函数k kk11=p(2) „„„„„„计算11 a2=[1,p(2);p(2),1]; >> a22=inv(a2);

>> kk2=a22*p(1,2:3)' „„„„„„„计算22

kk22=kk2(2,1); kk22=kk2(2,1)

a5=[1,p(2),p(3),p(4),p(5);p(2),1,p(2),p(3),p(4);p(3),p(2),1,p(2),p(3);p(4),p(3),p(2),1,p(2);p(5),p(4),p(3),p(2),1];

a55=inv(a5);

kk5=a55*p(1,2:6)';

kk55=kk5(5,1) „„„„„计算55

2

D=r0-kk11*r(2)-kk22*r(3) „„„„计算

附录3： zc=0;

for i=1:64 zc=zc+z(i); end

zp=zc/64; for k=1:64 w(k)=z(k)-zp; end w

rc=0; for k=1:64

rc=rc+1/64*w(k)*w(k); end

for k=1:62 x=0;

for j=1:64-k

x=x+1/64*w(j)*w(j+k); end r(k)=x;

p(k)=r(k)/rc; end p

m(1,1)=p(1); for k=1:61 x=0; y=0;

for j=1:k

x=x+p(k+1-j)*m(k,j); y=y+p(j)*m(k,j); end

m(k+1,k+1)=(p(k+1)-x)/(1-y); for j=1:k

m(k+1,j)=m(k,j)-m(k+1,k+1)*m(k,k-j+1); end end

(diag(m))' figure(1) x=1:62;

plot(x,p,x,0); figure(2) x=1:62;

plot(x,(diag(m))',x,0)

附录4：

[ACF]=autocorr(x,10)

[PartialACF]=parcorr(x,10) subplot(121);autocorr(x,10) subplot(122);parcorr(x,10)

附录5：

pppp=AICB(1,1); nna=1; nnb=1; for p=1:5 for q=1:5

if pppp>AICB(p,q) pppp=AICB(p,q); nna=p; nnb=q; end end end nna

nnb

评 阅 书

目 录

摘要„„„„„„„„„„„„„„„„„„„„„„„1

第一章 绪论„„„„„„„„„„„„„„„„„„„2

1.1 时间序列模型的发展及其作用„„„„„„„„„„„2

1.2 什么是时间序列模型„„„„„„„„„„„„„„„2

1.3 本文研究的主要方法和手段„„„„„„„„„„„„2

1.4 本文主要研究思路及内容安排„„„„„„„„„„„2

第二章 ARMA模型„„„„„„„„„„„„„„„„„4

2.1 ARMA模型的基本原理„„„„„„„„„„„„„„„4

2.2 样本自协方差函数、自相关函数和偏相关函数„„„„4

2.3 ARMA模型识别方法„„„„„„„„„„„„„„„„5

2.4 模型参数估计„„„„„„„„„„„„„„„„„„6

第三章 实例分析„„„„„„„„„„„„„„„„„7

3.1 题目„„„„„„„„„„„„„„„„„„„„„„7

3.2 问题分析„„„„„„„„„„„„„„„„„„„„7

3.3 问题求解„„„„„„„„„„„„„„„„„„„„8

3.3.1数据的观测„„„„„„„„„„„„„„„„8

3.3.2数据处理„„„„„„„„„„„„„„„„„8

3.3.3求解自相关和偏相关函数„„„„„„„„„„8

3.4 模型的识别及求解„„„„„„„„„„„„„„„„9

3.5 结论„„„„„„„„„„„„„„„„„„„„„„11

参考文献„„„„„„„„„„„„„„„„„„„„„12 附录„„„„„„„„„„„„„„„„„„„„„„„12 评阅书„„„„„„„„„„„„„„„„„„„„„„15

《随机过程 》 课程设计任务书

摘 要

ARMA模型是研究时间序列的重要方法，由自回归模型（简称AR模型）与

滑动平均模型（简称MA模型）为基础“混合”构成。ARMA模型广泛应用在经济、工程等各个领域得益于其在具体预测方面的优势。在许多方面用该模型所作出的预测比其他传统经济计量方法更加精确。平稳时间序列模型主要有自回归模型(AR)、滑动平均模型(MA)和自回归滑动平均模型(ARMA)等，这些线性模型考虑因素较简单。自回归滑动平均模型(ARMA)计算简单，易于实时更新数据。

本文描述了ARMA模型的原理、自相关函数和偏相关函数的计算过程、模型

的识别方法以及ARMA模型的计算过程。并给出一组平稳时间序列的数据，对数据进行分析和处理，求出自相关系数和偏相关，并利用MATLAB软件画出自相关系数和偏相关图形，有图可知它们都是拖尾的，因此可以确定是ARMA(p,q)模型。接下来就是确定ARMA(p,q)的阶数，本文采用了AIC准则确定模型的阶数，在实际问题中，为使线性模型简单起见，通常p与q的数值被取得较小，却需都不为零。确定阶数后，就用我们学过的求解方法解出未知的参数，这样我们就得到了混合模型的表达式。

关键字：ARMA(p,q)模型，自相关函数，偏相关函数

第一章 绪 论

1.1 时间序列模型的发展及其作用

人们的一切活动，其根本目的无不在于认识和改造客观世界。时间序列分析

不仅可以从数量上揭示某一现象的发展变化规律或动态的角度刻画某一现象与其他现象之间的内在关系及其变化规律性，达到认识客观世界之目的，而且运用时间序列模型还可以预测和控制现象的未来行为，修正或重新设计系统以达到利用和改造客观之目的。时间序列分析在工程技术中有重要的作用，常用于做预报、控制等。为建立其随机线性模型，首先，我们应明白什么是时间序列：时间序列是随机序列，即参数离散的随机过程。由于工程中遇到的随机序列的参数经常为时间，故称随机序列为随机时间序列，简称时间序列。可以说，时间序列是随时间改变而随机变化的序列。平稳时间序列是平稳序列，它满足期望为零，且任意两个时刻的相关函数与时间t无关，仅与两个时刻的时间差相关。

1.2 什么是时间序列模型

所谓的时间序列模型就是用观测得到的若干数据，即样本数据的基础上，建

立符合事物发展规律的数学理论模型，就是根据时间序列Xt的一段样本观测数据Xt（1tN）所包含的信息，利用自协方差函数Rr和偏相关函数kk具有的“截尾”与“拖尾”性质，判定模型结构及阶数，时间序列的模型我们主要学习了AR(p)，MA(q)及混合模型ARMA(p,q)。 

1.3 本文研究的主要方法和手段

求解模型就是绘制出这些数据的图形观察图形的走向。通常该图形被分为两

种：平稳和非平稳图形。对于非平稳的图形，我们要用所学过的方法把这些非平稳的数据变成与之对应的平稳的数据，再用平稳的方法去研究它。对于平稳的图形，我们首先根据理论方法及数学软件，求解出自相关系数和偏相关系数，根据得出的自相关系数和偏相关系数分别绘制出图形，观察图形截尾还是拖尾，因此可以确定是哪类模型。接下来就是确定模型的阶数。确定阶数后即可求解出相关系数，这样就得到我们所要的模型表达式。

1.4 本文主要研究思路及内容安排

第一章介绍了时间序列模型的发展及其作用，什么是时间序列模型，课题研

究工作和论文的内容安排。

第二章介绍ARMA(p,q)的原理，自相关函数和偏相关函数，时间序列模型

模型的识别方法以及ARMA(p,q)的求解。

第三章利用具体问题，对问题分析，判断出符合ARMA(p,q)模型和并求解

ARMA(p,q)模型，并对文章进行总结。

第二章 ARMA模型

2.1 ARMA模型的基本原理

ARMA(p,q)模型：如果时间序列Xt满足：

Xt1Xt1pXtpat1at1qatq

则称时间序列为Xt服从(p,q)阶自回归滑动平均混合模型。或者记为：

(B)Xt(B)at，

特殊情况：q=0，模型即为AR(p)，p=0，模型即为MA(q)。

将预测指标随时间推移而形成的数据序列看作是一个随机序列，这一组随机过程所具有的依存关系体现着原始数据在时间上的延续性。一方面，影响因素的影响，另一方面，又有自身变动规律，假定影响因素为x1，x2，，xt由回归分析：

其中Y是预测对象的观测值， e为误差。作为预测对象Yt受到自身变化的

影响，其规律可由下式体现：

误差项在不同时期具有依存关系，由下式表示：

由此，获得ARMA模型表达式：

Yt01xt12xt2pxtp01t12t2qtqt

2.2 样本自协方差函数、自相关函数和偏相关函数

平稳序列,2,1,0,1,2， E(t)0，对于样本，定义自协方差函数：

1nkˆkjkj， nnj1ˆkˆk/ˆ0。同时为了保证ˆkk，ˆkk一般取则自相关函数为：

n50,kn/。常取4kn/10。 1k12k2nkn

为求出样本的偏相关函数，需要求解如下的Yule-Walker方程：

12k11k11111

22k22



1kkkk1211

ˆ1ˆ1ˆ2k1ˆˆˆ11k21ˆˆ12

ˆˆ2ˆkkk1ˆk1ˆ1ˆ111解得ˆ1ˆ2 ˆk

2.3 ARMA模型识别方法

通常平稳时间序列Zt，t0,1仅进行有限n次测量(n50)，得到一个样

1ntZt，本函数，且利用平稳序列各态历经性：Zj做变换，t1,n，nj1

将Z1,,Zn样本换算成为样本1,,n，然后再确定平稳时间序列{t,t0,1}的随机线性模型。

按样本自相关函数和样本偏相关函数的值分别作出点图，按“截尾”，“拖尾”

情况，查表确定模型的类别与阶数p，q。

理论上，线性模型的类别的确定可以根据kk和k的拖尾性和截尾性来判

别。但是在实际应用中，我们常用有一个样本算出k=k，kk=kk判别k，kk是拖尾还是截尾的。随机线性模型的三种形式：AR(p) ，MA(q)，ARMA(p,q)的判别分别如下：

（1）若k拖尾，kk截尾在kp处，则线性模型为AR(p)模型。k拖尾可

以用的点图判断，只要样本自相关函数的绝对值愈变愈小（k增大时）；但是，用样本偏相关函数判断kk 截尾在k=p处应该如何作呢？因为k>p时，样本偏相

关函数不为零，而偏相关函数kk0，这就给判断带来一定的困难。我们可以采

用一下方法解决：当kp时，平均20个样本偏相关函数中至多有一个使|kk| ≥

2/则认为kk截尾在kp处。

（2）若kk拖尾，k在kp处截尾，那么线性模型为MA(q)滑动平均模型。

^^^

kk拖尾可以根据样本偏相关函数的点图判断，只要|kk|愈变愈小（k增大时）。

但是，用样本自相关函数判断自相关函数k在kp处截尾可采用如下方法：当

kp时，若平均20个样本自相关函数中至多有一个使|k|≥

^^

（3)若样本自相关函数和样本偏相关函数都是拖尾的，则线性模型可以看成混和模型（或者，样本自相关函数和样本偏相关函数都不为截尾的，又被负指数型的数列所控制）。其具体的判别方法和上述一样。如下表所示：

2.4 模型参数估计

1、AR(p)模型参数估计：

2AR(p)模型有p2个参数：p,1,2,,p,。利用Yule-Walker方程，利

用Toeplitz矩阵求逆和作矩阵乘法的方法算样本偏相关函数kk。AR(p)模型的参数值不必作专门的计算，只要在样本偏相关函数计算的记录中取出样本参数值

0jj。即可。此时1,2,,p,都已经确定了，经过推理我们可以得到：

j12、MA(q)滑动平均模型参数估计： 2p

2ˆ2ˆ2ˆ2),k0ˆ(112q ˆk2ˆˆˆˆˆˆ(k1k+1qkq),1kq

ˆ,ˆˆ,ˆ2，即解这个非线性方程组。 可得q1个方程，求12q

3、ARMA(p,q)混和模型参数估计

对于满足一个条件：t1t1...ptpat1at1...patq采用先计算

2ˆ,ˆ,,ˆ，在计算ˆ,ˆˆ,ˆ的方法，具体如下：1）可利用Toeplitz矩阵和作矩阵12p12q

ˆ,ˆ,,ˆ。2）令...混和模型化为：乘法的方法求出12ptt1t1ptp'

t'at1at1...patq这是关于t'的MA(q)模型，用t'的样本协方差函数估计

2ˆˆ1,ˆ2ˆq,的值。

第三章 实例分析

3.1 题目

某简谐振动位移离平衡点的距离zt的观测数据如下所示，N=64，建立表示测

数据序列的时间序列模型。

表-1 原始数据

-0.17007 -2.05056 -4.18035 0.11153 2.17675 4.08507 0.24983 -4.92336 -8.15675 -4.42107 -2.64058 -0.52570 -0.22015 0.36515

2.61949 1.53203 -0.77529 -3.25439 -0.41904 3.92127 2.97641 1.05303 -3.07748 -3.41256 -5.37570 -3.53729 -0.79940 -1.07239

0.78588 -2.48208 -4.03418 0.16420 -0.39307 3.88025 6.03561

5.70750 1.45257 -2.97501 -3.16561 -3.15184 -5.66101 -4.55502 0.10608 0.04231 0.79915 -0.12247 -5.97528 -9.16058 -9.73835 -5.22347 0.62229 4.34356 1.32363 -1.25405 -3.90467 -3.02789

2.78487 7.61065 6.14580 3.51779 1.79635 0.74996 -0.19133 -1.40428

3.2 问题分析

本问题是让我们建立时间序列模型，题中给出了64个数据首先用所给数据

绘制图形。根据数据绘制出的图形如下（程序见附录1）：

[**************]

图（1）

由图形可以看出所给数据本身就是稳定的时间序列，接下来用MATLAB软件求解出自相关系数和偏相关系数，根据得出的自相关系数和偏相关系数分别绘制出图形，观察图形截尾还是拖尾，由图可知自相关系数和偏相关图形都拖尾，因此可以确定是ARMA(p,q)模型。接下来就是确定ARMA(p,q)的阶数。在实际问题中，为使线性模型简单起见，通常p与q的数值被取得较小，却需都不为零。确定阶数后，就用我们学过的求解方法解出未知的参数，这样我们就得到了ARMA(p,q)模型的表达式。

3.3 问题求解

3.3.1 数据的观测：

对一个时间序列作n次测量得到一个样本Z,Z,,Z，一般取n50。某简12n

谐振动位移离平衡点的距离zt的观测数据，得到n=64个样本值Z1,Z2,,Z60，（列于表1中）。表1的Z栏，共64个数据。将原始样本数据经过处理后变成时间t

序列W。 t

3.3.2 数据处理：

（1）求取误差时间序列Zt：

164

用MATLAB软件对原始数据进行处理，做变换WtZtZ，其中ZZt，代64t1入表1数据得：，1Z[(0.17007)(2.05056)(1.40428)]0.757464

（程序见附录2） WtZt(0.7574)，

3.3.3 数据处理：

ˆkˆk/ˆ0，（2）计算样本自协方差函数k，样本自方差函数k。 其中k0,1,2,3,4,5，

ˆk1k12k2nknn1nkjkj。由图-3数据可得：随着k的增大，k越nj1

来越小，具有拖尾性。

（3）接下来计算偏相关函数kk(k1)。利用Yule-Walker方程，

12k11k1111122k22



1kkkk1211

求得偏相关函数。

由机械振动位移图可看出这些数据画出的图形是平稳序列，由MATLAB软件

运行以上数据（程序见附录），得到自相关系数和偏相关系数如下表:

3.4 模型的识别及求解

ˆkk值，绘制ˆkk的曲线趋势来判断，ˆK值和偏相关函数ˆK和 根据自相关函数

分别画自相关系数和偏相关系数点图如下（程序分别见附录[3]，[4]）： 自相关函数图：

1

0.8

0.6

0.4

0.2

-0.2

-0.4

010

[1**********]0

图（3）

偏相关函数图如图：

0.80.6

0.40.20-0.2-0.4-0.6-0.8

[**************]

图（4）

Sample Autocorrelation Function (ACF)Sample Partial Autocorrelation Function

05

Lag

1015

Sample Partial Autocorrelations

Sample Autocorrelation

05

Lag

1015

图（5）

图形可以看出自相关函数和偏相关函数都拖尾，所以该数据符合ARMA(p,q)混合模型。用ARMA(p,q)模型对该数据进行拟合，即:

Xt1Xt1pXtpat1at1qatq

在计算机上，利用图（6）所示的AIC准则对候选模型参数进行计算和分析：（程序见附录5）

图（6）

根据计算知当p=2，q =2时的ARMA（2，2）模型，1 =1.637754，

2

2=0.7523818，1=-0.01194535，2=0.2609289，a=1.206746，AIC=30.7927

（最小）。根据模型定阶的AIC准则，比较上述计算结果可知，最能描述上述机械振动观测数据序列的模型结构应该是ARMA（2，2）模型，即：

Xt1.6378Xt10.7524Xt2at0.0119at10.26at2

3.5 结论

由以上问题的求解过程可知，所有的时间序列模型参数都是建立在样本统计参数的基础上，通过对数据处理，利用平稳时间序列的数据，得出的ARMA(p,q)模型，随机过程论目前已得到广泛的应用，在诸如天气预报、统计物理、天体物理、运筹决策、经济数学、安全科学、人口理论、可靠性及计算机科学等很多领域都要经常用到随机过程的理论来建立数学模型。时间序列分析也广泛应用到了解决许多问题上，如信号处理、金融预测、状态估计、控制和模式识别等.相信随着时间的发展，随机过程会越来越多的应用于我们的生活中，为我们带来更多的惊喜和发现。

参考文献

[2] 张善文 雷英杰 冯有前. MATLAB在时间序列分析中的应用 西安：西安电子科技大学出版社 2007 。

[2] 汪荣鑫编. 随机过程 西安：西安交通大学出版社 2006 。 [3] 刘卫国 Matlab程序设计与应用 北京 高等教育出版社 2006 。 [4] 田铮 秦超英 随机过程与应用 北京 科学出版社 2007 。

附录

附录1

z=[-0.17007 -2.05056 -4.18035 0.11153 2.17675 4.08507 0.24983 -4.92336 -8.15675 -4.42107 -2.64058 -0.52570 -0.22015 0.36515 2.61949 1.53203 -0.77529 -3.25439 -0.41904 3.92127 2.97641 1.05303 -3.07748 -3.41256 -5.37570 -3.53729 -0.79940 -1.07239 0.78588 -2.48208 -4.03418 0.16420 -0.39307 3.88025 6.03561 5.70750 1.45257 -2.97501 -3.16561 -3.15184 -5.66101 -4.55502 0.10608 0.04231 0.79915 -0.12247 -5.97528 -9.16058 -9.73835 -5.22347 0.62229 4.34356 1.32363 -1.25405 -3.90467 -3.02789 2.78487 7.61065 6.14580 3.51779 1.79635 0.74996 -0.19133 -1.40428];

>> x=[0:1:63]; plot(x,z) 附录2：

b=length(z);

m=sum(z)/b „„„„„„计算均值 w= z-m;

sumw=zeros(1,6); sumw1=0; for j=1:64

sumw1=sumw1+w(j)^2; „„„„计算0

sumw1 for k=0:5 for i=1:(b-k)

sumw(k+1)=sumw(k+1)+w(i)*w(i+k); .........计算k end

sumw(k+1) r=sumw/b; r0=sumw1/b;

p=r/r0 „„„„„„计算自相关函数k kk11=p(2) „„„„„„计算11 a2=[1,p(2);p(2),1]; >> a22=inv(a2);

>> kk2=a22*p(1,2:3)' „„„„„„„计算22

kk22=kk2(2,1); kk22=kk2(2,1)

a5=[1,p(2),p(3),p(4),p(5);p(2),1,p(2),p(3),p(4);p(3),p(2),1,p(2),p(3);p(4),p(3),p(2),1,p(2);p(5),p(4),p(3),p(2),1];

a55=inv(a5);

kk5=a55*p(1,2:6)';

kk55=kk5(5,1) „„„„„计算55

2

D=r0-kk11*r(2)-kk22*r(3) „„„„计算

附录3： zc=0;

for i=1:64 zc=zc+z(i); end

zp=zc/64; for k=1:64 w(k)=z(k)-zp; end w

rc=0; for k=1:64

rc=rc+1/64*w(k)*w(k); end

for k=1:62 x=0;

for j=1:64-k

x=x+1/64*w(j)*w(j+k); end r(k)=x;

p(k)=r(k)/rc; end p

m(1,1)=p(1); for k=1:61 x=0; y=0;

for j=1:k

x=x+p(k+1-j)*m(k,j); y=y+p(j)*m(k,j); end

m(k+1,k+1)=(p(k+1)-x)/(1-y); for j=1:k

m(k+1,j)=m(k,j)-m(k+1,k+1)*m(k,k-j+1); end end

(diag(m))' figure(1) x=1:62;

plot(x,p,x,0); figure(2) x=1:62;

plot(x,(diag(m))',x,0)

附录4：

[ACF]=autocorr(x,10)

[PartialACF]=parcorr(x,10) subplot(121);autocorr(x,10) subplot(122);parcorr(x,10)

附录5：

pppp=AICB(1,1); nna=1; nnb=1; for p=1:5 for q=1:5

if pppp>AICB(p,q) pppp=AICB(p,q); nna=p; nnb=q; end end end nna

nnb

评 阅 书