z(k)-z(k-1)+-z(k-3)=u(k-1)+u(k-2)-u(k-3)+v(k)-0.2v(k-1) 其中v(k)是服从N(0,1)分布的不相关随机噪声，对该系统用最小二乘递推法对参数进行辨识，并用理论值进行比较，绘制输入输出和辨识参数曲线、辨识误差曲线，给出详细具体的matlab代码

以下是对该系统进行最小二乘递推辨识的Matlab代码，包括数据预处理、递推最小二乘辨识和性能分析：

%% 数据预处理
load('data.mat'); % 加载数据文件
u = data.u; % 输入信号
y = data.y; % 输出信号
v = 1*randn(size(u)); % 生成服从N(0,1)分布的不相关随机噪声

% 数据分割
u1 = u(1:3); % 前3个数据
y1 = y(1:3);
u2 = u(4:end); % 后面的数据
y2 = y(4:end);
v2 = v(4:end);

%% 递推最小二乘辨识
N = length(u2); % 数据长度
P = 1e6*eye(5); % 参数误差协方差矩阵
theta = zeros(5,1); % 参数向量
epsilon = 0; % 误差估计值

for k = 1:N
% 计算当前时刻的输出信号
yk = theta(1)*y2(k-1) + theta(2)*y2(k-3) + theta(3)*u2(k-1) + theta(4)*u2(k-2) - theta(5)*u2(k-3) + v2(k) - 0.2*v2(k-1);

% 更新参数误差协方差矩阵
P = P - (P*[y2(k-1);y2(k-3);u2(k-1);u2(k-2);-u2(k-3)]*[y2(k-1),y2(k-3),u2(k-1),u2(k-2),-u2(k-3)]*P)/(1+[y2(k-1),y2(k-3),u2(k-1),u2(k-2),-u2(k-3)]*P*[y2(k-1);y2(k-3);u2(k-1);u2(k-2);-u2(k-3)]);

% 更新参数向量
theta = theta + P*[y2(k-1);y2(k-3);u2(k-1);u2(k-2);-u2(k-3)]*(yk - y2(k-1)*theta(1) - y2(k-3)*theta(2) - u2(k-1)*theta(3) - u2(k-2)*theta(4) + u2(k-3)*theta(5));

% 更新误差估计值
epsilon = epsilon + (yk - y2(k-1)*theta(1) - y2(k-3)*theta(2) - u2(k-1)*theta(3) - u2(k-2)*theta(4) + u2(k-3)*theta(5))^2;
end

%% 性能分析
% 绘制输入输出曲线和辨识输出曲线
figure;
plot(u2,'b'); hold on;
plot(y2,'r'); hold on;
y_ident = zeros(size(y2));
for k = 4:N
y_ident(k) = theta(1)*y2(k-1) + theta(2)*y2(k-3) + theta(3)*u2(k-1) + theta(4)*u2(k-2) - theta(5)*u2(k-3) + v2(k) - 0.2*v2(k-1);
end
plot(y_ident,'g');
legend('输入信号','输出信号','辨识输出信号');
xlabel('时间');

% 绘制辨识参数曲线和理论参数曲线
figure;
plot([1:N],ones(1,N)*1,'k--'); hold on;
plot([1:N],ones(1,N)*(-1),'k--'); hold on;
plot([1:N],theta(1)*ones(1,N),'r'); hold on;
plot([1:N],theta(2)*ones(1,N),'b'); hold on;
plot([1:N],theta(3)*ones(1,N),'g'); hold on;
plot([1:N],theta(4)*ones(1,N),'m'); hold on;
plot([1:N],theta(5)*ones(1,N),'c');
legend('理论参数范围','a1','a2','b1','b2','b3');
xlabel('时间');

% 计算辨识误差曲线
err = zeros(size(y2));
for k = 4:N
err(k) = y2(k) - theta(1)*y2(k-1) - theta(2)*y2(k-3) - theta(3)*u2(k-1) - theta(4)*u2(k-2) + theta(5)*u2(k-3) - v2(k) + 0.2*v2(k-1);
end

% 绘制辨识误差曲线
figure;
plot(err);
xlabel('时间');
ylabel('误差');