《随机过程》
课程设计（论文）
题 目: 平 稳 时 间 序 列 的 MA(q) 模 型 的 预 报
学 院： 理 学 院 专 业： 数 学 与 应 用 数 学 班 级： 数 学 09-1 班 学 生 姓 名： 学 生 学 号： 指 导 教 师：
年 月 日
目 录 目 录
任务书 ...................................................................................... 错误！未定义书签。 摘 要 ....................................................................................... 错误！未定义书签。
第1章 绪 论 ............................................................................................................ 1 第2章 MA（q）基本原理 ........................................................................................ 2
2.1 MA（q）模型 ................................................................................................ 3 2.2模型的识别 .................................................................................................. 3 2.3 MA模型功率谱估计设计思路 ..................................................................... 4 第3章 Matlab的软件实现.6 3.1 matlab的运行程序
3.2程序流程图 ................................................................................................. 第4章具体应用 ..........................................................................................................
4.1平稳时间序列的标准 ....................................................................................
4.2求w的自相关函数....................................................................................................... 参考文献………………………………………………………………………………………
附 录 ........................................................................................................................... 评阅书……………………………………………
随机过程 课程设计任务书
摘 要
在实际生活中，我们经常要遇到时间序列的分析问题。如某地逐月的降水量、市场上的小麦的逐年的平均价格、病人的心电图的采样信号等等都是时间问题的例子近年来由于时间序列的分析理论和方法的迅速发展，在工业、国防等的领域里人们常常对一系列观察数据分析研究。由于受各种因素的影响，表现随机性，但又有某种统计上的关系，按时间先后次序排列的研究的集合叫时间序列
时间序列分析，就是根据有序随机变量或者观测得到有血数据之间的相互依赖所包含的信息，用概率统计的方法定量建立一个合适的数学模型，并且根据这个模型对相应的所反应的序列过程或进行控制，在自然科学中，工程技术及社会、经济学的建模分析中起着重要作用的平稳时间序列---滑动平均模型，简称MA(q)模型。
MA(q)模型的基本思路是：将预测对象随时间的推移而形成的数据序列视为一个随机序列，用一定的数学模型来近似描述这个序列。这个模型一旦被识别就可以从时间序列的过去值及现在的值来预测未来的值。现在统计方法、计算经济模型，在某种程度上已经能够帮助企业对未来进行预测。
通过一组平稳时间序列的数据，求出自相关函数和偏向关函数，画出图判别平稳的时间符合MA(q)模型，在参数估计求出MA(q)模型。
根据自信关函数的截尾性和篇相关函数的来初步判断MA(q)模型的q也就是根据正态分布的性质来取介数。得到了拟合模型然后求出参数用迭代法的方法来解出方程的解在检验拟合的效果。
关键字：时间序列分析 自相关函数截尾 迭代法
利用平稳时间序列的MA(q)模型的预报
第1章 绪 论
现代信号分析中，对于常见的具有各态历经的平稳随机信号，不可能用清楚的数学关系式来描述，但可以利用给定的N个样本数据估计一个平稳随机信号的功率谱密度叫做功率谱估计(PSD)。它是数字信号处理的重要研究内容之一。功率谱估计可以分为经典功率谱估计(非参数估 和现代功率谱估计(参数估计)。功率谱估计在实际工程中有重要应用价值，如在语音信号识别、雷达杂波分析、波达方向估计、地震勘探信号处理、水声信号处理、系统辨识中非线性系统识别、物理光学中透镜干涉、流体力学的内波分析、太阳黑子活动周期研究等许多领域．发挥了重要作用。
谱估计分为两大类：非参数化方法和参数化方法。非参数化谱估计又叫做经典谱估计，其主要缺陷是频率分辨率低；而参数化谱估计又叫做现代谱估计，它具有频率分辨率高的优点。因此，有时又把它叫做高分辨率谱估计。经典的谱估计方法也是直接按定义用有限长度数据来估计。有两条基本途径：(1)先估计相关函数，再经傅氏变换得功率谱估计。(2)把功率谱和信号幅频特性的平方联系起来。不论采用哪一条途径，存在的共同问题是估计的方差特性不好，而且估计值沿频率轴的起伏甚烈。数据越长，这一现象越严重。因此需要采用一些措施来改进估计的性能。
经典功率谱估计是将数据工作区外的未知数据假设为零，相当于数据加窗。经典功率谱估计方法分为：相关函数法(BT法)、周期图法以及两种改进的周期图估计法即平均周期图法和平滑平均周期图法，其中周期图法应用较多，具有代表性。 现代功率谱估计
现代功率谱估计即参数谱估计方法是通过观测数据估计参数模型再按照求参数模型输出功率的方法估计信号功率谱。主要是针对经典谱估计的分辨率低和方差性能不好等问题提出的。主要方法有最大熵谱分析法(AR模型法)、Pisarenko谐波分解法、Prony提取极点法、Prony谱线分解法以及Capon最大似然法等。
第
2章 MA（q）基本原理
2.1 MA（q）模型
设{Xt}为零均值的是平稳时间序列，阶数为q的MA(q)模型定义为
其中
xt=at-θt- θqat-q (2.1)
{θ
k
,k
=1,2, q}称为滑动平均系数，并简记模型为MA(Q.)满足
MA(Q)模型的随机序列称为MA(Q)序列延迟算子表示，以写成
X
t
=θ(B)at，其中，
θ(B)=1-θ1B- -θqBq
(2.2)
对于MA(q)模型如果满足条件：θ(B)=0的根全在单位园外，即所有的模于1，就称条件为MA(Q)模型的可逆性条件。当模型满足可逆性条件时，θ-1(B)存在，此时(2.2)的式可以写成
at=θ-1(B)Xt
它称为逆转形式。模型xt可以看做是白噪声序列{at}输入线性系统中得输出MA模型
就是滑动平均模型。MA模型的系统H(z)表示为
H(z)=G(1+∑biz-i)
1
可以看出，MA模型的系统函数H(z)只有零点，没有极点，故MA模型又称为全零点模型。模型的阶由多项式的阶q决定，q阶MA模型表示为MA(q)。如果在白噪声w(n)激励下模型的输出为x(n),则模型输入、输出关系的时域表达式为
x(n)=G[w(n)+∑biw(n-i)]=G∑biw(n-i),其中b0=1
i=1
i=0
2
此式为MA模型的差分方程。由该方程可以看出，MA模型现在的输出x(n)是现在的输入w(n)和过去的q个输入w(n-i)(i=1,2,…q)的加权和，即输出时输入的滑动平均。所以MA模型叫滑动平均模型或动均模型。相应地，将白噪声w(n)激励MA模型产生的输出x(n)叫做MA过程。由于MA模型是全零点模型，所以MA过程的功率谱具有深谷而无尖锐的峰，具有这一特点的随机信号适合选择MA模型。
由于MA(q)过程的功率谱：
2 qqq
22j2πf2-j2πfk2*
Sx(f)MA=σB(e)=σ∑bke=σ∑∑bkble-j2π(l-k)f
k=0k=0l=0
3
令m=l-k，l=k＋m，可得 q0
⎛q-m*⎫-j2πfm⎛q*⎫-j2πfm⎤2⎡S(f)=σbbe+bb ⎪ ∑⎢∑ ∑kk+m⎪⎥xMA ∑kk+m⎪⎪e
m=-q⎝k=-m⎭⎭⎣m=0⎝k=0⎦
4
又MA过程的自相关函数是：
q-m
⎫2*
σ∑bkbk+m=Rx(m)m=0,1,⋅⋅⋅,q⎪
⎪k=0
⎬q-mq
2*2**
及σ∑bkbk+m=σbl-mbl=Rx(-m)=Rx(m)m=-q, ,-1⎪ ⎪k=-ml=0⎭
q
Sx(f)MA=Rx(m)e-j2πfm因此功率谱
m=-q
∑
想运用平稳时间序列模型对实际用生活中的问题进行预报和控制，首先我们得知道是哪类时间序列模型，然后才能运用此模型进行相关分析。因此，如何判别时间序列模型以及怎样确定模型的参数成为解决本问题的关键。
2.2 模型的识别
MA(q)序列的自相关函数
用Xt-k乘以上式两边，再取均值(由于序列的值为零，故自相关函数与协方 函数相同)为了不混淆，记所有得协方差函数为λk，有
λk=E[XtXt-k]
=E[atat-k]-∑θjE[atat-k-j]-∑θiE[at-iat-k]
j=1
i=1
=E[(at-θ1at-1- -θqat-q)(at-k-θiat-k-1- θqat-k-q)] +∑∑θiθjE[at-iat-j]
i=1j=1q
q
利用
2
,t=s⎧σa
E[asat]=⎨
,t≠s0⎩
显然上式 第二项对一切k都为零，其余的各项k.. (1)当k=0时，有
2
λ0=E[a]+∑θE[a]=σ+∑θi2σa
2
t
2i
2t-i
2a
i=1
i=1
q
q
（2）当1≤k≤q时。有
λk=-θkE[a]+2t-k
i=k+1
∑θθ
q
ii-k
E[a]=-θkσ=-θkσ+2t-i2a2a
i=k+1
∑θθ
q
ii-k
2 σa
（3）当k>q时。右边的四项都为0，此时λk=0 用λ0除以λk得标准化自相关函数ρk=λk序列的协方差函数λk和自相关函数
*
σ2∑bkbk+m=Rx(m)
q-mk=0
q
综上所述可得MA(Q)0简称它为自相关函数。
ρk为
及σ2
k=-m
∑bb
q-m
*
kk+m
=σ2
b
l=0
*
l-ml*
b=Rx(-m)=Rx(m)
⎫
m=0,1,⋅⋅⋅,q⎪
⎪⎬
m=-q, ,-1⎪
⎪⎭
2.3MA模型功率谱估计设计思路
由Kolmogorov定理，即任何ARMA(p,q)过程或MA(q)过程都能用无限阶的AR(∞)过程表示；同样，任何ARMA(p,q)过程或AR(p)过程也可用一个MA(∞)过程表示。即：MA(q)模型可等效成AR(∞)模型。
11
H(z)=B(z)==
MA模型D(z)B(z)7
∞
x(n)=-dkx(n-k)+u(n)求逆Z变换，得
k=1
8
故可等效成AR(∞)模型。
相对于MA模型的参数估计来说，AR模型的参数估计比较简单，而且现在对于AR模型的参数估计的理论研究的比较深入，基于这点，本次算法原理是由N点数据x(n)建立一个p阶的AR模型，从而求出p阶AR系数a，再利用求解AR系数建立线性预测，等效q阶MA模型(p>>q),再利用求解AR模型求出MA模型的参数b，从而实现MA功率谱估计。
∑
第3章 Matlab的软件实现
3.1 matlab的运行程序 N=456;
B1=[1 0.3544 0.3508 0.1736 0.2401]; A1=[1];
w=linspace(0,pi,512);
H1=freqz(B1,A1,w);%产生信号的频域响应 Ps1=abs(H1).^2;
SPy11=0;%20次AR(4) SPy14=0;%20次MA(4) VSPy11=0;%20次AR(4) VSPy14=0;%20次MA(4) for k=1:20
%采用自协方差法对AR模型参数进行估计%
y1=filter(B1,A1,randn(1,N)).*[zeros(1,200),ones(1,256)]; [Py11,F]=pcov(y1,4,512,1);%AR(4)的估计% [Py13,F]=periodogram(y1,[],512,1); SPy11=SPy11+Py11;
VSPy11=VSPy11+abs(Py11).^2; %------------MA模型---------------% y=zeros(1,256); for i=1:256 y(i)=y1(200+i); end
ny=[0:255];
z=fliplr(y);nz=-fliplr(ny);
nb=ny(1)+nz(1);ne=ny(length(y))+nz(length(z)); n=[nb:ne]; Ry=conv(y,z); R4=zeros(8,4); r4=zeros(8,1); for i=1:8
r4(i,1)=-Ry(260+i); for j=1:4
R4(i,j)=Ry(260+i-j); end end R4 r4
a4=inv(R4'*R4)*R4'*r4%利用最小二乘法得到的估计参数 %对MA的参数b(1)-b(4)进行估计% A1
A14=[1,a4']%AR的参数a(1)-a(4)的估计值
B14=fliplr(conv(fliplr(B1),fliplr(A14)));%MA模型的分子
y24=filter(B14,A1,randn(1,N));%.*[zeros(1,200),ones(1,256)];%由估计出的MA模型产生数据 [Ama4,Ema4]=arburg(y24,32), B1
b4=arburg(Ama4,4)%求出MA模型的参数 %---求功率谱---% w=linspace(0,pi,512); %H1=freqz(B1,A1,w)
H14=freqz(b4,A14,w);%产生信号的频域响应 %Ps1=abs(H1).^2;%真实谱 Py14=abs(H14).^2;%估计谱 SPy14=SPy14+Py14;
VSPy14=VSPy14+abs(Py14).^2; end
figure(1)
plot(w./(2*pi),Ps1,w./(2*pi),SPy14/20);
legend('真实功率谱','20次MA(4)估计的平均值');
3.2程序流程图
图1 程序流程图
第4章 具体应用
4.1、平稳时间序列的标准
z=[1196.8 1181.3 1222.6 1229.3 1221.5 1148.4 1250.2 1174.4 1234.5 1209.7 1206.5 1204.0 1234.1 1146.0 1304.9 1221.9 1244.1 1194.4 1281.5 1277.3 1238.9 1267.5 1200.9 1245.5 1249.9 1220.1 1267.4 1182.3 1221.7 1187.1 1261.6 1274.5 1196.4 1222.6 1174.7 1212.6 1215.0 1191.0 1179.0 1224.0 1183.0 1288.0 1274.0 1218.0 1263.0 1205.0 1210.0 1243.0 1266.0
1200.0 1306.0 1209.0 1248.0 1208.0 1231.0 1244.0 1296.0 1221.0 1287.0 1191.0];
命令：（1） %计算z的平均值
m=mean(z)
m =
1.2268e+003
（2） %z的标准化
for k=1:60
w(ｋ)=z(ｋ)-m; %验证股票的价格数据是平
end
w
Columns 1 through 5
-29.9850 -45.4850 -4.1850 2.5150 -5.2850
Columns 6 through 10
-78.3850 23.4150 -52.3850 7.7150 -17.0850
Columns 11 through 15
-20.2850 -22.7850 7.3150 -80.7850 78.1150
Columns 16 through 20
-4.8850 17.3150 -32.3850 54.7150 50.5150
Columns 21 through 25
12.1150 40.7150 -25.8850 18.7150 23.1150
Columns 26 through 30
-6.6850 40.6150 -44.4850 -5.0850 -39.6850
Columns 31 through 35
34.8150 47.7150 -30.3850 -4.1850 -52.0850
Columns 36 through 40
-14.1850 -11.7850 -35.7850 -47.7850 -2.7850
Columns 41 through 45
-43.7850 61.2150 47.2150 -8.7850 36.2150
Columns 46 through 50
-21.7850 -16.7850 16.2150 39.2150 -26.7850
Columns 51 through 55
79.2150 -17.7850 21.2150 -18.7850 4.2150
Columns 56 through 60
17.2150 69.2150 -5.7850 60.2150 -35.785
稳时间序列
plot(X)
4.2求w的自相关函数
[ACF]=autocorr(X,20)
ACF = Columns 1 through 5
1.0000 -0.1800 0.2775 -0.0793 0.0518
Columns 6 through 10
0.0594 0.2017 -0.1074 0.1029 -0.1273
Columns 11 through 15
-0.0001 -0.0201 0.0774 -0.1903 -0.1265
Columns 16 through 20
-0.0743 -0.0659 0.0272 -0.2209 -0.1177
Column 21
-0.0807
%画出自相关函数图形
（2） autocorr(X,20)
4.3根据自相关函数的截尾性求出模型
根据自相关函数的截尾性取M=≈7,当K=1时,
^^12+2ρ)=(1+2⨯0.685)=0.81在ρk
(k=2,3,4,5,6,7,8)中满足1
ˆ1+ρˆ2)]=0.192根据拟合MA(5)的残0.的仅占+2(ρ81当k=2时差方差最小拟合后发现θ5和θ4MA(3)，的系数比较好因而模型
为MA(3)。模型结果如图表
7参考文献：
[1] 张贤达 现代信号处理，清华大学出版社，北京，2002
[2] 丁玉美 数字信号处理一时域离散随机信号处理 ,西安电子科技大学出版社2002
[3] 陈怀琛 吴大正．MATLAB及在电子信息课程中的应用(第二版),电子工业出版社，2004
[5] 程卫东 冯新华 基于参数化功率谱估计的正弦频率估计
[6] 赵红怡 数字信号处理及其Matlab实现
评 阅 书
《随机过程》
课程设计（论文）
题 目: 平 稳 时 间 序 列 的 MA(q) 模 型 的 预 报
学 院： 理 学 院 专 业： 数 学 与 应 用 数 学 班 级： 数 学 09-1 班 学 生 姓 名： 学 生 学 号： 指 导 教 师：
年 月 日
目 录 目 录
任务书 ...................................................................................... 错误！未定义书签。 摘 要 ....................................................................................... 错误！未定义书签。
第1章 绪 论 ............................................................................................................ 1 第2章 MA（q）基本原理 ........................................................................................ 2
2.1 MA（q）模型 ................................................................................................ 3 2.2模型的识别 .................................................................................................. 3 2.3 MA模型功率谱估计设计思路 ..................................................................... 4 第3章 Matlab的软件实现.6 3.1 matlab的运行程序
3.2程序流程图 ................................................................................................. 第4章具体应用 ..........................................................................................................
4.1平稳时间序列的标准 ....................................................................................
4.2求w的自相关函数....................................................................................................... 参考文献………………………………………………………………………………………
附 录 ........................................................................................................................... 评阅书……………………………………………
随机过程 课程设计任务书
摘 要
在实际生活中，我们经常要遇到时间序列的分析问题。如某地逐月的降水量、市场上的小麦的逐年的平均价格、病人的心电图的采样信号等等都是时间问题的例子近年来由于时间序列的分析理论和方法的迅速发展，在工业、国防等的领域里人们常常对一系列观察数据分析研究。由于受各种因素的影响，表现随机性，但又有某种统计上的关系，按时间先后次序排列的研究的集合叫时间序列
时间序列分析，就是根据有序随机变量或者观测得到有血数据之间的相互依赖所包含的信息，用概率统计的方法定量建立一个合适的数学模型，并且根据这个模型对相应的所反应的序列过程或进行控制，在自然科学中，工程技术及社会、经济学的建模分析中起着重要作用的平稳时间序列---滑动平均模型，简称MA(q)模型。
MA(q)模型的基本思路是：将预测对象随时间的推移而形成的数据序列视为一个随机序列，用一定的数学模型来近似描述这个序列。这个模型一旦被识别就可以从时间序列的过去值及现在的值来预测未来的值。现在统计方法、计算经济模型，在某种程度上已经能够帮助企业对未来进行预测。
通过一组平稳时间序列的数据，求出自相关函数和偏向关函数，画出图判别平稳的时间符合MA(q)模型，在参数估计求出MA(q)模型。
根据自信关函数的截尾性和篇相关函数的来初步判断MA(q)模型的q也就是根据正态分布的性质来取介数。得到了拟合模型然后求出参数用迭代法的方法来解出方程的解在检验拟合的效果。
关键字：时间序列分析 自相关函数截尾 迭代法
利用平稳时间序列的MA(q)模型的预报
第1章 绪 论
现代信号分析中，对于常见的具有各态历经的平稳随机信号，不可能用清楚的数学关系式来描述，但可以利用给定的N个样本数据估计一个平稳随机信号的功率谱密度叫做功率谱估计(PSD)。它是数字信号处理的重要研究内容之一。功率谱估计可以分为经典功率谱估计(非参数估 和现代功率谱估计(参数估计)。功率谱估计在实际工程中有重要应用价值，如在语音信号识别、雷达杂波分析、波达方向估计、地震勘探信号处理、水声信号处理、系统辨识中非线性系统识别、物理光学中透镜干涉、流体力学的内波分析、太阳黑子活动周期研究等许多领域．发挥了重要作用。
谱估计分为两大类：非参数化方法和参数化方法。非参数化谱估计又叫做经典谱估计，其主要缺陷是频率分辨率低；而参数化谱估计又叫做现代谱估计，它具有频率分辨率高的优点。因此，有时又把它叫做高分辨率谱估计。经典的谱估计方法也是直接按定义用有限长度数据来估计。有两条基本途径：(1)先估计相关函数，再经傅氏变换得功率谱估计。(2)把功率谱和信号幅频特性的平方联系起来。不论采用哪一条途径，存在的共同问题是估计的方差特性不好，而且估计值沿频率轴的起伏甚烈。数据越长，这一现象越严重。因此需要采用一些措施来改进估计的性能。
经典功率谱估计是将数据工作区外的未知数据假设为零，相当于数据加窗。经典功率谱估计方法分为：相关函数法(BT法)、周期图法以及两种改进的周期图估计法即平均周期图法和平滑平均周期图法，其中周期图法应用较多，具有代表性。 现代功率谱估计
现代功率谱估计即参数谱估计方法是通过观测数据估计参数模型再按照求参数模型输出功率的方法估计信号功率谱。主要是针对经典谱估计的分辨率低和方差性能不好等问题提出的。主要方法有最大熵谱分析法(AR模型法)、Pisarenko谐波分解法、Prony提取极点法、Prony谱线分解法以及Capon最大似然法等。
第
2章 MA（q）基本原理
2.1 MA（q）模型
设{Xt}为零均值的是平稳时间序列，阶数为q的MA(q)模型定义为
其中
xt=at-θt- θqat-q (2.1)
{θ
k
,k
=1,2, q}称为滑动平均系数，并简记模型为MA(Q.)满足
MA(Q)模型的随机序列称为MA(Q)序列延迟算子表示，以写成
X
t
=θ(B)at，其中，
θ(B)=1-θ1B- -θqBq
(2.2)
对于MA(q)模型如果满足条件：θ(B)=0的根全在单位园外，即所有的模于1，就称条件为MA(Q)模型的可逆性条件。当模型满足可逆性条件时，θ-1(B)存在，此时(2.2)的式可以写成
at=θ-1(B)Xt
它称为逆转形式。模型xt可以看做是白噪声序列{at}输入线性系统中得输出MA模型
就是滑动平均模型。MA模型的系统H(z)表示为
H(z)=G(1+∑biz-i)
1
可以看出，MA模型的系统函数H(z)只有零点，没有极点，故MA模型又称为全零点模型。模型的阶由多项式的阶q决定，q阶MA模型表示为MA(q)。如果在白噪声w(n)激励下模型的输出为x(n),则模型输入、输出关系的时域表达式为
x(n)=G[w(n)+∑biw(n-i)]=G∑biw(n-i),其中b0=1
i=1
i=0
2
此式为MA模型的差分方程。由该方程可以看出，MA模型现在的输出x(n)是现在的输入w(n)和过去的q个输入w(n-i)(i=1,2,…q)的加权和，即输出时输入的滑动平均。所以MA模型叫滑动平均模型或动均模型。相应地，将白噪声w(n)激励MA模型产生的输出x(n)叫做MA过程。由于MA模型是全零点模型，所以MA过程的功率谱具有深谷而无尖锐的峰，具有这一特点的随机信号适合选择MA模型。
由于MA(q)过程的功率谱：
2 qqq
22j2πf2-j2πfk2*
Sx(f)MA=σB(e)=σ∑bke=σ∑∑bkble-j2π(l-k)f
k=0k=0l=0
3
令m=l-k，l=k＋m，可得 q0
⎛q-m*⎫-j2πfm⎛q*⎫-j2πfm⎤2⎡S(f)=σbbe+bb ⎪ ∑⎢∑ ∑kk+m⎪⎥xMA ∑kk+m⎪⎪e
m=-q⎝k=-m⎭⎭⎣m=0⎝k=0⎦
4
又MA过程的自相关函数是：
q-m
⎫2*
σ∑bkbk+m=Rx(m)m=0,1,⋅⋅⋅,q⎪
⎪k=0
⎬q-mq
2*2**
及σ∑bkbk+m=σbl-mbl=Rx(-m)=Rx(m)m=-q, ,-1⎪ ⎪k=-ml=0⎭
q
Sx(f)MA=Rx(m)e-j2πfm因此功率谱
m=-q
∑
想运用平稳时间序列模型对实际用生活中的问题进行预报和控制，首先我们得知道是哪类时间序列模型，然后才能运用此模型进行相关分析。因此，如何判别时间序列模型以及怎样确定模型的参数成为解决本问题的关键。
2.2 模型的识别
MA(q)序列的自相关函数
用Xt-k乘以上式两边，再取均值(由于序列的值为零，故自相关函数与协方 函数相同)为了不混淆，记所有得协方差函数为λk，有
λk=E[XtXt-k]
=E[atat-k]-∑θjE[atat-k-j]-∑θiE[at-iat-k]
j=1
i=1
=E[(at-θ1at-1- -θqat-q)(at-k-θiat-k-1- θqat-k-q)] +∑∑θiθjE[at-iat-j]
i=1j=1q
q
利用
2
,t=s⎧σa
E[asat]=⎨
,t≠s0⎩
显然上式 第二项对一切k都为零，其余的各项k.. (1)当k=0时，有
2
λ0=E[a]+∑θE[a]=σ+∑θi2σa
2
t
2i
2t-i
2a
i=1
i=1
q
q
（2）当1≤k≤q时。有
λk=-θkE[a]+2t-k
i=k+1
∑θθ
q
ii-k
E[a]=-θkσ=-θkσ+2t-i2a2a
i=k+1
∑θθ
q
ii-k
2 σa
（3）当k>q时。右边的四项都为0，此时λk=0 用λ0除以λk得标准化自相关函数ρk=λk序列的协方差函数λk和自相关函数
*
σ2∑bkbk+m=Rx(m)
q-mk=0
q
综上所述可得MA(Q)0简称它为自相关函数。
ρk为
及σ2
k=-m
∑bb
q-m
*
kk+m
=σ2
b
l=0
*
l-ml*
b=Rx(-m)=Rx(m)
⎫
m=0,1,⋅⋅⋅,q⎪
⎪⎬
m=-q, ,-1⎪
⎪⎭
2.3MA模型功率谱估计设计思路
由Kolmogorov定理，即任何ARMA(p,q)过程或MA(q)过程都能用无限阶的AR(∞)过程表示；同样，任何ARMA(p,q)过程或AR(p)过程也可用一个MA(∞)过程表示。即：MA(q)模型可等效成AR(∞)模型。
11
H(z)=B(z)==
MA模型D(z)B(z)7
∞
x(n)=-dkx(n-k)+u(n)求逆Z变换，得
k=1
8
故可等效成AR(∞)模型。
相对于MA模型的参数估计来说，AR模型的参数估计比较简单，而且现在对于AR模型的参数估计的理论研究的比较深入，基于这点，本次算法原理是由N点数据x(n)建立一个p阶的AR模型，从而求出p阶AR系数a，再利用求解AR系数建立线性预测，等效q阶MA模型(p>>q),再利用求解AR模型求出MA模型的参数b，从而实现MA功率谱估计。
∑
第3章 Matlab的软件实现
3.1 matlab的运行程序 N=456;
B1=[1 0.3544 0.3508 0.1736 0.2401]; A1=[1];
w=linspace(0,pi,512);
H1=freqz(B1,A1,w);%产生信号的频域响应 Ps1=abs(H1).^2;
SPy11=0;%20次AR(4) SPy14=0;%20次MA(4) VSPy11=0;%20次AR(4) VSPy14=0;%20次MA(4) for k=1:20
%采用自协方差法对AR模型参数进行估计%
y1=filter(B1,A1,randn(1,N)).*[zeros(1,200),ones(1,256)]; [Py11,F]=pcov(y1,4,512,1);%AR(4)的估计% [Py13,F]=periodogram(y1,[],512,1); SPy11=SPy11+Py11;
VSPy11=VSPy11+abs(Py11).^2; %------------MA模型---------------% y=zeros(1,256); for i=1:256 y(i)=y1(200+i); end
ny=[0:255];
z=fliplr(y);nz=-fliplr(ny);
nb=ny(1)+nz(1);ne=ny(length(y))+nz(length(z)); n=[nb:ne]; Ry=conv(y,z); R4=zeros(8,4); r4=zeros(8,1); for i=1:8
r4(i,1)=-Ry(260+i); for j=1:4
R4(i,j)=Ry(260+i-j); end end R4 r4
a4=inv(R4'*R4)*R4'*r4%利用最小二乘法得到的估计参数 %对MA的参数b(1)-b(4)进行估计% A1
A14=[1,a4']%AR的参数a(1)-a(4)的估计值
B14=fliplr(conv(fliplr(B1),fliplr(A14)));%MA模型的分子
y24=filter(B14,A1,randn(1,N));%.*[zeros(1,200),ones(1,256)];%由估计出的MA模型产生数据 [Ama4,Ema4]=arburg(y24,32), B1
b4=arburg(Ama4,4)%求出MA模型的参数 %---求功率谱---% w=linspace(0,pi,512); %H1=freqz(B1,A1,w)
H14=freqz(b4,A14,w);%产生信号的频域响应 %Ps1=abs(H1).^2;%真实谱 Py14=abs(H14).^2;%估计谱 SPy14=SPy14+Py14;
VSPy14=VSPy14+abs(Py14).^2; end
figure(1)
plot(w./(2*pi),Ps1,w./(2*pi),SPy14/20);
legend('真实功率谱','20次MA(4)估计的平均值');
3.2程序流程图
图1 程序流程图
第4章 具体应用
4.1、平稳时间序列的标准
z=[1196.8 1181.3 1222.6 1229.3 1221.5 1148.4 1250.2 1174.4 1234.5 1209.7 1206.5 1204.0 1234.1 1146.0 1304.9 1221.9 1244.1 1194.4 1281.5 1277.3 1238.9 1267.5 1200.9 1245.5 1249.9 1220.1 1267.4 1182.3 1221.7 1187.1 1261.6 1274.5 1196.4 1222.6 1174.7 1212.6 1215.0 1191.0 1179.0 1224.0 1183.0 1288.0 1274.0 1218.0 1263.0 1205.0 1210.0 1243.0 1266.0
1200.0 1306.0 1209.0 1248.0 1208.0 1231.0 1244.0 1296.0 1221.0 1287.0 1191.0];
命令：（1） %计算z的平均值
m=mean(z)
m =
1.2268e+003
（2） %z的标准化
for k=1:60
w(ｋ)=z(ｋ)-m; %验证股票的价格数据是平
end
w
Columns 1 through 5
-29.9850 -45.4850 -4.1850 2.5150 -5.2850
Columns 6 through 10
-78.3850 23.4150 -52.3850 7.7150 -17.0850
Columns 11 through 15
-20.2850 -22.7850 7.3150 -80.7850 78.1150
Columns 16 through 20
-4.8850 17.3150 -32.3850 54.7150 50.5150
Columns 21 through 25
12.1150 40.7150 -25.8850 18.7150 23.1150
Columns 26 through 30
-6.6850 40.6150 -44.4850 -5.0850 -39.6850
Columns 31 through 35
34.8150 47.7150 -30.3850 -4.1850 -52.0850
Columns 36 through 40
-14.1850 -11.7850 -35.7850 -47.7850 -2.7850
Columns 41 through 45
-43.7850 61.2150 47.2150 -8.7850 36.2150
Columns 46 through 50
-21.7850 -16.7850 16.2150 39.2150 -26.7850
Columns 51 through 55
79.2150 -17.7850 21.2150 -18.7850 4.2150
Columns 56 through 60
17.2150 69.2150 -5.7850 60.2150 -35.785
稳时间序列
plot(X)
4.2求w的自相关函数
[ACF]=autocorr(X,20)
ACF = Columns 1 through 5
1.0000 -0.1800 0.2775 -0.0793 0.0518
Columns 6 through 10
0.0594 0.2017 -0.1074 0.1029 -0.1273
Columns 11 through 15
-0.0001 -0.0201 0.0774 -0.1903 -0.1265
Columns 16 through 20
-0.0743 -0.0659 0.0272 -0.2209 -0.1177
Column 21
-0.0807
%画出自相关函数图形
（2） autocorr(X,20)
4.3根据自相关函数的截尾性求出模型
根据自相关函数的截尾性取M=≈7,当K=1时,
^^12+2ρ)=(1+2⨯0.685)=0.81在ρk
(k=2,3,4,5,6,7,8)中满足1
ˆ1+ρˆ2)]=0.192根据拟合MA(5)的残0.的仅占+2(ρ81当k=2时差方差最小拟合后发现θ5和θ4MA(3)，的系数比较好因而模型
为MA(3)。模型结果如图表
7参考文献：
[1] 张贤达 现代信号处理，清华大学出版社，北京，2002
[2] 丁玉美 数字信号处理一时域离散随机信号处理 ,西安电子科技大学出版社2002
[3] 陈怀琛 吴大正．MATLAB及在电子信息课程中的应用(第二版),电子工业出版社，2004
[5] 程卫东 冯新华 基于参数化功率谱估计的正弦频率估计
[6] 赵红怡 数字信号处理及其Matlab实现
评 阅 书