%-----AR(3)--------
clear;
clc;
close all;
N=1000; %步长
rand('seed',10); %改变种子改变数据
a=-[ -1.2000 0.2900 -0.0180]; %poly([0.1 0.2 0.9])
e=0.6*randn(1,N+1); %白噪声序列
y(1)=e(1); %初值y(1),y(2),y(3)
y(2)=a(1)*y(1)+e(2);
y(3)=a(1)*y(2)+a(2)*y(1)+e(3);
for t=4:N+1
y(t)=a(1)*y(t-1)+a(2)*y(t-2)+a(3)*y(t-3)+e(t);%AR(3)的RLS算法实现
end
n=3;
lamda=1;%遗忘因子
sita(:,3)=zeros(n,1);%初值,均为0
P(:,:,3)=10^5*eye(n);%P的初值,选取足够大的单位阵
for t=3:N
fai(:,t+1)=[y(t) y(t-1) y(t-2)]';%最小二乘格式的fai(t)的表达式
%参数估值
sita(:,t+1)=sita(:,t)+P(:,:,t)*fai(:,t+1)/(lamda+fai(:,t+1)'*P(:,:,t)*fai(:,t+1))*...
[y(t+1)-fai(:,t+1)'*sita(:,t)];
P(:,:,t+1)=1/lamda*(P(:,:,t)-P(:,:,t)*fai(:,t+1)/(lamda+fai(:,t+1)'*P(:,:,t)*fai(:,t+1))*...
fai(:,t+1)'*P(:,:,t));
end
ee(1)=y(1);
taoe(1)=ee(1)^2;%y(0)的初值为0
for t=2:N
%对噪声方差的估计
ee(t)=y(t)-fai(:,t)'*sita(:,t-1);
taoe(t)=taoe(t-1)+1/t*[ee(t)^2-taoe(t-1)];
end
t=1:N;
subplot(2,2,1);
plot(t,sita(1,t),'r');
line([0,N],[a(1),a(1)])
subplot(2,2,2);
plot(t,sita(2,t),'r');
line([0,N],[a(2),a(2)])
subplot(2,2,3);
plot(t,sita(3,t),'r');
line([0,N],[a(3),a(3)])
subplot(2,2,4);
plot(t,taoe,'r');
line([0,N],[0.36,0.36])
其中,e代表,为白噪声序列
lamda代表,为遗忘因子
sita代表,估值参数之一
fai代表,为估值参数之一
ee代表
taoe代表