智能控制系统实验报告ARMA 模型ARMA(p,q)是一个线性时间序列预测模型，适用于平稳的时间序列，即对于任何时刻t ，都有()t E Z μ=,E(a t )=0.协方差矩阵(')t t E a a ∑=,对于任意0l ≠有(')0t t l E a a +=。

AR 模型11t t p t p t Z Z Z νφφε--=++++(0.1)当误差项t ε自相关时，可以被有限阶滑动平均表示11t t t q t q a a a ε--=+Θ++Θ(0.2)这里t a 是零均白噪声，协方差矩阵a ∑非奇异。

结合AR 过程和MA 误差项，得到ARMA 过程：111111t t p t p tt p t p t t q t qZ Z Z Z Z a a a νφφενφφ------=++++=+++++Θ++Θ (0.3)对于一个很大的阶数n ，AR(n)接近ARMA(p,q)1()()nt t i t i i Z n Z a n -==+∑∏(0.4)从(0.4)得到残差的估计值1ˆˆ()()nt t i t ii a n Z n Z -==-∏∑ (0.5)其中ˆ()i n ∏利用多变量最小二乘法求解，然后使用估计值ˆ()t a n 建立多变量回归模型1111ˆˆt t p t p t t q t q Z Z Z a aa φφ----=++++Θ++Θ (0.6)1111[,,:,.]()ˆ()ˆ()t t p t p q t t t q Z Z Z a n a n a n φφ----⎡⎤⎢⎥⎢⎥⎢⎥=ΘΘ+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦(0.7)(1:)0[,]T Z Y A =ΦΘ+ (0.8) 01[,,]T A a a = (0.9)000'ˆˆˆ()a A A n T∑= (0.10)最小二乘法求解公式，以AR(p)为例。

111121[,,,,]t t p t p tt p tt p Z Z Z a Z a Z νφφνφφφ----=++++⎡⎤⎢⎥⎢⎥=+⎢⎥⎢⎥⎢⎥⎣⎦112111112111[,,][,,,]T T k Tp p p T p Z Z Z Z Z Z Z Z νφφ---⨯---⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦ (0.11)记11(1)101(1),:[,,]:[,,,],[1,',,']':[,,]1T k T T k Tp k kp t t t p T kp TZ ,,Z a a v Y Z Z Y Y φφ⨯⨯⨯+-+-+⨯Φ==== Z =Y +AZ:=[]A Y φ (0.12)()()k z I uvec ββ=⊗+=Y φ (0.13)最小化目标函数 11()'[(')]'()[(')]u k T a k S u u z I I z I βββ--=∑=-⊗⊗∑-⊗Y Y(0.14)解得 111ˆ'(')(+)'(')'(')---===+ZY YY Y A Y YY AY YY φφφ (0.15)以下是仿真例题：例题1：使用AR 模型预测德国西部经济数据。

该数据是三元时间序列，包括投资(Investment)，可支配输入(income)和消费支出(consumption)三个变量。

现使用AR 模型预测INCOME ，将1960年到1978年的75组用于训练，其余用于测试，利用最小二乘回归估计系数矩阵12[,,,,]p v φφφ 代码如下：clear clcdata1=load('e1.txt'); % data1=data1'; data2=log(data1);data3=data2(:,2:end)-data2(:,1:en d-1); %差分数据[K,num]=size(data3); data=data3; T=73;recordaic=100; for p=2:-1:2 clear Z; clear Zc;for t=p+1:T+pY(:,t-p)=data(:,t); %构造矩阵Ymiddle1=data(:,(t-p):(t-1)); middle=fliplr(middle1);Z(:,t-p)=[1;middle(:)]; %构造矩阵Zend% ---------------直接求伪逆%B1=Y*pinv(Z)%------------------y=Y(:);middleb=inv(Z*Z')*Z;b=kron(middleb,eye(K))*y;%------------------最小二乘函数b1=regress(Y(:),kron(Z',eye(K)));B1=reshape(b,K,K*p+1) %权值矩阵%-------------B=reshape(b,K,K*p+1); %权值矩阵%B=[-0.017,-0.32,0.146,0.961,-0.1 61,0.115,0.934;0.016,0.044,-0.153,0.2 89,0.050,...%0.019,-0.010;0.013,-0.002,0.225,-0.26 4,0.034,0.355,-0.022];for t=p+1:num % for t=T+p+1:numO(:,t-p)=data(:,t);middle2=data(:,(t-p):(t-1));middlec=fliplr(middle2);Zc(:,t-p)=[1;middlec(:)];endYc=B*Zc;%----------直接伪逆Yc1=B1*Zc;data41=Yc1+data2(:,p+1:end-1); %第三个到第八十九个数dataf1=exp(data41); %反对数rmse1=(sum((data1(3,3:91)-dataf1( 3,:)).^2)/(89-1))^0.5 %计算均方根误差%-----------data4=Yc+data2(:,p+1:end-1); %%反差分第三个到第八十九个数dataf=exp(data4); %反对数U=dataf(:,1:T)-data1(:,p+1:T+p);sigmaU=(1/T)*U*U';aic=log(det(sigmaU))+2*p*K*K/T;if aic<recordaicrecordaic=aic;m=p;Yc_b=Yc;dataf_b=dataf;endendsubplot(2,1,1)plot(1:92-1-m ,Yc_b(3,:),'b',1:92 -1-m ,data3(3,m+1:91),'r'); %差分数据画图subplot(2,1,2)plot(1:92-1-m ,data1(3,m+1:91),'b ',1:92-1-m ,dataf_b(3,:),'r'); %实际数据画图rmse=(sum((data1(3,m+1:91)-dataf( 3,:)).^2)/(91-m-1))^0.5 %计算均方根误差运行结果如下：实验分析：效果图如下图所示上面的子图是对数差分数据预测效果，下面的子图是反差分反对数的预测效果图，预测均方根误差为25.1096.可以看出，预测效果很好。

例题2：课件ARMA_AR_MA_RepresentGranger Analysis,P34的例题AR(p)模型转换成MA(∞)模型。

代码如下：function [psi]=ARchangeMA(AR_fi,num_psi) %输入AR_fi和num_psi，求MA模型中的psi psi=[];p=size(AR_fi,2); %AR_fi的列数if (p<=0)errordlg('model is not corrected');endif(p>0)k=size(AR_fi,1); %AR_fi的行数p=p/k; %求pD=eye(k);I=eye(k);O= zeros(k,k*(p-1));J=[I,O]; %求Jif(p==1)xi=AR_fi;J=I;elsexi=[AR_fi;[eye(k*(p-1)),zeros(k,k)]]; %求xi,xi为kp×kp阶矩阵end%判断是否平稳lamta=eig(xi);for i=1:p*kif(abs(lamta(i))>1) %判断此模型成立条件,lamta为xi的特征值，lamta=1/zreturnendend%求psipsi=eye(k);for i=1:num_psi-1psi=[psi,J*xi^i*J'];endend运行结果：例题3：卡尔曼滤波仿真中所选的模型的状态方程是AR(2)模型，基于老师所给数据。

代码如下：clc;clear all;H=diag([1 1 1 0 0 0]);G=zeros(6,6);v=[-0.017,0.016,0.013,-0.017,0.016,0.013]';phai_1=[-0.320,0.146,0.961;0.044,-0.153,0.289;-0.002,0.225,-0.264];phai_2=[-0.161,0.115,0.934;0.050,0.019,-0.010;0.034,0.355,-0.022];B=zeros(6,6);F=diag(ones(1,6));B(1:3,1:3)=phai_1;B(1:3,4:6)=phai_2;B(4:6,1:3)=diag(ones(1,3));Original_data=zeros(92,3);Original_data(:,1)=[180;179;185;192;211;202;207;214;231;229;234;237;206;250;259 ;263;264;280;282;292;286;302;304;307;317;314;306;304;292;275;273;301;280;289;30 3;322;315;339;364;371;375;432;453;460;475;496;494;498;526;519;516;531;573;551;5 38;532;558;524;525;519;526;510;519;538;549;570;559;584;611;597;603;619;635;658; 675;700;692;759;782;816;844;830;853;852;833;860;870;830;801;824;831;830;]; Original_data(:,2)=[451;465;485;493;509;520;521;540;548;558;574;583;591;599;610 ;627;642;653;660;694;709;734;751;763;766;779;808;785;794;799;799;812;837;853;87 6;897;922;949;979;988;1025;1063;1104;1131;1137;1178;1211;1256;1290;1314;1346;1385;1416;1436;1462;1493;1516;1557;1613;1642;1690;1759;1756;1780;1807;1831;1873;1 897;1910;1943;1976;2018;2040;2070;2121;2132;2199;2253;2276;2318;2369;2423;2457; 2470;2521;2545;2580;2620;2639;2618;2628;2651;];Original_data(:,3)=[415;421;434;448;459;458;479;487;497;510;516;525;529;538;546 ;555;574;574;586;602;617;639;653;668;679;686;697;688;704;699;709;715;724;746;75 8;779;798;816;837;858;881;905;934;968;983;1013;1034;1064;1101;1102;1145;1173;12 16;1229;1242;1267;1295;1317;1355;1371;1402;1452;1485;1516;1549;1567;1588;1631;1 650;1685;1722;1752;1774;1807;1831;1842;1890;1958;1948;1994;2061;2056;2102;2121; 2145;2164;2206;2225;2235;2237;2250;2271;];Sample_Original_data=Original_data(4:92,:)';data=zeros(92,3);data(:,1)=log10(Original_data(:,1));data(:,2)=log10(Original_data(:,2));data(:,3)=log10(Original_data(:,3));Difference_data=zeros(91,3);for i=2:92Difference_data(i-1,:)=data(i,:)-data(i-1,:);endStudy_data=zeros(6,89);%跟踪的信号Output_data=zeros(6,89);Estimate_data_t_t_1=zeros(6,89);Estimate_data_t_t=zeros(6,89);Estimate_Smoothing=zeros(6,89);a=1.0e-003 *[0.3922 0.0154 0.0236; 0.0154 0.0274 0.0122; 0.0236 0.0122 0.0168]; Esti_sigma=diag(ones(1,6));Esti_sigma(1:3,1:3)=a;Esti_sigma_y=zeros(6,6);Esti_sigma_y(1:3,1:3)=(0.1)*a;Esti_sigma_y(4:6,4:6)=a;Esti_sigma_out=zeros(6,6);Esti_sigma_t_1=zeros(6,6);Esti_sigma_t=0.004*diag(ones(1,6));Esti_eror=zeros(3,89);kalman_p=zeros(6,6);%相关参数for k=1:89Study_data(1:3,k)=Difference_data(k+2,:)';Study_data(4:6,k)=[0 0 0]';%Study_data(4:6,k)=Difference_data(k+1,:)';%Estimate_data_t_t_1(:,k)=B*Study_data(:,k)+F*v;endfor t=1:89if t==1Estimate_data_t_t_1(:,t)=v;%Esti_sigma_t_1=Esti_sigma;elseEstimate_data_t_t_1(:,t)=B*Estimate_data_t_t(:,t-1)+v;%Esti_sigmat_1=B*Esti_sigma_t+Esti_sigma;endEsti_sigma_t_1=B*Esti_sigma_t+Esti_sigma;Output_data(:,t)=H*Estimate_data_t_t_1(:,t);Esti_sigma_out=H*Esti_sigma_t_1*H'+Esti_sigma_y+0.05*Esti_sigma;%Esti_sigma_out=H*Esti_sigma_t_1*H'+Esti_sigma_y+Esti_sigma;kalman_p=Esti_sigma_t_1*(H')*(inv(Esti_sigma_out));Estimate_data_t_t(:,t)=Estimate_data_t_t_1(:,t)+kalman_p*(Study_data(:,t)-Outpu t_data(:,t));Esti_sigma_t=Esti_sigma_t_1-kalman_p*Esti_sigma*kalman_p';Esti_error(:,t)=Study_data(1:3,t)-Output_data(1:3,t);endkalman_pfiguresubplot(221)stairs(Output_data(1,:),'r');hold onstairs(Study_data(1,:),'b');grid on;legend('预报值','真实值');xlabel('投资');subplot(222)stairs(Output_data(2,:),'r');hold on;stairs(Study_data(2,:),'b');grid on;legend('预报值','真实值');xlabel('收入');subplot(223)stairs(Output_data(3,:),'r');hold on;stairs(Study_data(3,:),'b');grid on;legend('预报值','真实值');xlabel('消费');subplot(224)stairs(Esti_error(1,:),'r');hold on;stairs(Esti_error(2,:),'b');hold on;stairs(Esti_error(3,:),'b:');grid on;legend('投资误差','输入误差','消费误差');运行结果：训练出的卡尔曼滤波的系数矩阵：实验分析：由于此仿真中的输出方程为随意选定，输出方程的随机干扰噪声与实际偏差较大，故此仿真的跟踪效果不是很理想。