利用ARMA、AR、MA模型,以及周期图等进行系统参数估计

利用ARMA、AR、MA模型，以及周期图等进行系统参数估计
[Copy to clipboard][ - ]
CODE:
N=456;
B1=[1 0.3544 0.3508 0.1736 0.2401];
A1=[1 -1.3817 1.5632 -0.8843 0.4096];
w=linspace(0,pi,512);
H1=freqz(B1,A1,w);%产生信号的频域响应
Ps1=abs(H1).^2;
SPy11=0;%20次AR(4)
SPy12=0;%20次AR(8)
SPy13=0;%20次AR周期图
SPy14=0;%20次ARMA(4,4)
SPy15=0;%20次ARMA(8,8)
VSPy11=0;%20次AR(4)
VSPy12=0;%20次AR(8)
VSPy13=0;%20次AR周期图
VSPy14=0;%20次ARMA(4,4)
VSPy15=0;%20次ARMA(8,8)
for k=1:20
%采用自协方差法对AR模型参数进行估计%
%gA1：AR模型的参数；gE1：激励白噪声的方差%
y1=filter(B1,A1,randn(1,N)).*[zeros(1,200),ones(1,256)];
[Py11,F]=pcov(y1,4,512,1);%AR(4)的估计%
[Py12,F]=pcov(y1,8,512,1);%AR(8)的估计%
[Py13,F]=periodogram(y1,[],512,1);
SPy11=SPy11+Py11;
SPy12=SPy12+Py12;
SPy13=SPy13+Py13;
VSPy11=VSPy11+abs(Py11).^2;
VSPy12=VSPy12+abs(Py12).^2;
VSPy13=VSPy13+abs(Py13).^2;
figure(1)
plot(w./(2*pi),Ps1,F,Py11);
legend('真实功率谱','20次AR(4)估计图');
hold on;
figure(2)
plot(w./(2*pi),Ps1,F,Py12);
legend('真实功率谱','20次AR(8)估计图');
hold on;
figure(3)
plot(w./(2*pi),Ps1,F,Py13);
legend('真实功率谱','20次周期图法估计图');
hold on;
%------------ARMA模型---------------%
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);%ARMA(4,4)的R
r4=zeros(8,1);%ARMA(4,4)的r
for i=1:8
r4(i,1)=-Ry(260+i);
for j=1:4
R4(i,j)=Ry(260+i-j);
end
end
R4%R矩阵
r4%r矩阵
a4=inv(R4'*R4)*R4'*r4%利用最小二乘法得到的a的估计参数
%----------------ARMA(4,4)对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模型产生数据%---因为(q=4)<<L<<N=256,所以选取L=32---%
[Ama4,Ema4]=arburg(y24,32),%利用数据y2估计AR(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;%估计谱
%if Py14>200
% PPy14=200;
%elseif Py14<200
% PPy14=Py14;
%end
SPy14=SPy14+Py14;
figure(4)
plot(w./(2*pi),Ps1,w./(2*pi),Py14);
legend('真实功率谱','20次ARMA(4,4)的估计图');
hold on;
R8=zeros(16,8);%ARMA(8,8)的R
r8=zeros(16,1);%ARMA(8,8)的r
for i=1:16
r8(i,1)=-Ry(264+i);
for j=1:8
R8(i,j)=Ry(264+i-j);
end
end
R8%R矩阵
r8%r矩阵
a8=inv(R8'*R8)*R8'*r8%利用最小二乘法得到的a的估计参数
%----------------ARMA(8,8)对MA的参数b(1)-b(8)进行估计----------------------%
A1
A18=[1,a8']%AR的参数a(1)-a(8)的估计值
B18=fliplr(conv(fliplr(B1),fliplr(A18)))%MA模型的分子
y28=filter(B18,A1,randn(1,N));%.*[zeros(1,200),ones(1,256)];%由估计出的MA模型产生数据%---因为(q=8)<<L<<N=256,所以选取L=48---%
[Ama8,Ema8]=arburg(y28,48),%利用数据y2估计AR(32)的参数
B1
b8=arburg(Ama8,8)%求出MA模型的参数
%---求功率谱---%
%w=linspace(0,pi,512);
%H1=freqz(B1,A1,w)
H18=freqz(b8,A14,w);%产生信号的频域响应
%Ps1=abs(H1).^2;%真实谱
Py15=abs(H18).^2;%估计谱
%if Py15>200
% PPy15=200;
%elseif Py15<200
% PPy15=Py15;
%end
SPy15=SPy15+Py15;
VSPy15=VSPy15+abs(Py15).^2;
figure(5)
plot(w./(2*pi),Ps1,w./(2*pi),Py15);
legend('真实功率谱','20次ARMA(8,8)的估计图');
hold on;
%-----------------------------------%
end
V2=VSPy12/20-abs(SPy12/20).^2;
V3=VSPy13/20-abs(SPy13/20).^2;
V4=VSPy14/20-abs(SPy14/20).^2;
V5=VSPy15/20-abs(SPy15/20).^2;
figure(6)
plot(w./(2*pi),Ps1,F,SPy11/20);
legend('真实功率谱','20次AR(4)估计的均值');
figure(7)
plot(w./(2*pi),Ps1,F,SPy12/20);
legend('真实功率谱','20次AR(8)估计的均值');
figure(8)
plot(w./(2*pi),Ps1,F,SPy13/20);
legend('真实功率谱','20次周期图估计的均值');
figure(9)
plot(w./(2*pi),Ps1,w./(2*pi),SPy14/20);
legend('真实功率谱','20次ARMA(4,4)估计的平均值');
figure(10)
plot(w./(2*pi),Ps1,w./(2*pi),SPy15/20);
legend('真实功率谱','20次ARMA(8,8)估计的平均值');
figure(12)
plot(F,V1,F,V2,F,V3,w./(2*pi),V4,w./(2*pi),V5);
legend('AR(4)方差','AR(8)方差','周期图方差','ARMA(4,4)方差','ARMA(8,8)方差');
axis([0 0.5 0 2000]);
请问happy教授:
1.以上这个网址怎么打不开呀?本人急需参考与此有关的内容( 有关如何用matlab设计ARMA模型的问题),请帮忙提供一下,谢谢!
2.我对一时间序列建模,但不知怎么来判断到底建模是否合适呢？对ARMA模型来说其残差、标准差达到多少数量级算这个模型最好呢？象我采取的ARMA模型残差是不是挺大的？是否要改为其他的模型如AR或MA模型等，请问您做过非平稳时间序列ARMA的建模问题吗？
Discrete-time IDPOL Y model: A(q)y(t) = C(q)e(t)
A(q) = 1 - 0.135 q^-1 - 0.8649 q^-2
C(q) = 1 + 0.2315 q^-1
Estimated using ARMAX from data set mydata_y
Loss function 115.824 and FPE 115.835
Discrete-time IDPOL Y model: A(q)y(t) = C(q)e(t)
A(q) = 1 - 1.556 q^-1 + 0.07357 q^-2 + 1.031 q^-3 - 0.5484 q^-4
C(q) = 1 - 1.171 q^-1 + 0.5629 q^-2 - 0.0342 q^-3
Estimated using ARMAX from data set mydata_y
Loss function 90.7377 and FPE 90.7571
对于模型具体阶次的判断是否需要求出各模型系数的95%置信区间再通过F检验准则来判断呢，请问模型系数的置信区间又是怎么求呢？是否需用最小二乘？那又怎么用呢？请帮忙,不胜感激!谢谢！
3.有相关的例子可以给我参考一下吗?。