CAR(2)模型的最小二参数估计的MATALB实现

%-----CAR(2)--------
clear;
clc;
close all;
N=4000;                             %步长
rand('seed',10);                    %改变种子改变数据
e=0.6*randn(1,N+1);                 %白噪声序列
a=-[-0.5000    0.0400];             %poly([0.1 0.4])
b=[3 1];
rand('seed',4);                     %改变种子改变数据
u=2*rand(1,N+1)-1;                  %u(t)在[-1,1]之间服从均匀分布

y(1)=e(1);                          %初值y(1),y(2),y(3)
y(2)=a(1)*y(1)+b(1)*u(1)+e(2);

for t=3:N+1
    y(t)=a(1)*y(t-1)+a(2)*y(t-2)+b(1)*u(t-1)+b(2)*u(t-2)+e(t);%AR(3)的RLS算法实现
end

n=2;
lamda=0.998;%遗忘因子
sita(:,2)=zeros(2*n,1);%初值,均为0
P(:,:,2)=10^5*eye(2*n);%P的初值,选取足够大的单位阵

for t=2:N
    fai(:,t+1)=[y(t) y(t-1) u(t)  u(t-1)]';%最小二乘格式的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(3,2,1);
plot(t,sita(1,t),'r');
line([0,N],[a(1),a(1)])
axis([0 N 0 1]);
title('a(1)');
subplot(3,2,2);
plot(t,sita(2,t),'r');
line([0,N],[a(2),a(2)])
axis([0 N -1 1]);
title('a(2)');
subplot(3,2,3);
plot(t,sita(3,t),'r');
line([0,N],[b(1),b(1)])
axis([0 N 0 5]);
title('b(1)');
subplot(3,2,4);
plot(t,sita(4,t),'r');
line([0,N],[b(2),b(2)])
axis([0 N 0 3]);
title('a(2)');
subplot(3,2,5);
plot(t,taoe,'r');
axis([0 N 0 1]);
line([0,N],[0.36,0.36])
title('taoe');
subplot(3,2,6);
plot(t,u(t),'r');
title('u(t)');

其中,e代表,为白噪声序列

lamda代表,为遗忘因子

sita代表,估值参数之一

fai代表,为估值参数之一

ee代表

taoe代表

其中一次的测试结果如下图所示: