ARMA模型自回归滑动平均模型

% 本程序的目的是模拟一个ARMA模型，然后进行时频归并。考察归并前后模型的变化。

% 这个ARMA模型的一般形式用黑盒子模型表示为A（q）y（t）=C（q）e（t）。q是滞后算子。
% 或者是：（1+a1*q^（-1）+a2*q^（-2）+a3^（-3）+a4*q^（-4））y（t）=（1+c1*q^（-1）+c2*q^（-2）+c3^（-3）+c4*q^（-4））e（t） 
% 这里多项式A和C都只写出4阶，因为一般的经济时间序列阶数都不高。
clear
tic
% ====================第一步，模拟一个ARMA模型并绘制ACF，PACF图========================
%s首先设定ARMA模型的多项式系数。ARMA模型中只有多项式A（q）和C（q），
%把A（q）的系数都设为0就得到MA模型，把C（q）的系数都设为0就得到AR模型。
a1 = -（0.5）^（1/3）;
a2 = （0.5）^（2/3）;
a3 = 0;
a4 = 0;
c1 = 0;
c2 = 0;
c3 = 0;
c4 = 0; 
obv = 3000;                       %obv是模拟的观测数目。 
A = [1 a1 a2 a3 a4];
B = [];                           %因为ARMA模型没有输入，因此多项式B是空的。
C = [1 c1 c2 c3 c4];
D = [];                           %把D也设为空的。
F = [];                           %ARMA模型里的F多项式也是空的。
m = idpoly（ABCDF11）         %这样就生成了ARMA模型，把它存储在m中。NoiseVariance被设定为1，1也是默认值。抽样间隔Ts设为1。 
error = randn（obv1）;             %生成一个obv*1的正态随机序列。准备用作模型的误差项。
e = iddata（[]error1）;           %用randn函数生成一个噪声序列。存储在e中。抽样间隔是1秒。
%u = [];                          %因为是ARMA模型，没有输出。所以把u设为空的。这句可省略。
y = sim（me）; 
get（y）                            %使用get函数来查看动态系统的所有性质。
r=y.OutputData;                   %把y.OutputData的全部值赋给变量r，r就是一个obv*1的向量。 
figure（1）
plot（r）                           %绘出y随时间变化的曲线。或者写成plot（y）结果是一样的。 

%==========================第二步，绘制ARMA序列r的ACF和PACF图=======================
%如果直接用函数autocorr和parcorr画图，总是会把滞后一阶的ACF和PACF也画出来，不好看，所以用下述方法画。
% figure（2）
% subplot（211）
% n=100;
% [ACFLagsBounds]=autocorr（rn2）;
% x=Lags（2:n）;
% y=ACF（2:n）;                            %注意这里的y和前面y的完全不同。
% h=stem（xy‘fill‘‘-‘）;
% set（h（1）‘Marker‘‘.‘）
% hold on
% ylim（[-1 1]）; 
% a=Bounds（11）*ones（1n-1）;
% line（‘XData‘x‘YData‘a‘Color‘‘red‘‘linestyle‘‘--‘）
% line（‘XData‘x‘YData‘-a‘Color‘‘red‘‘linestyle‘‘--‘）
% xlabel（‘lags‘）
% ylabel（‘ACF‘）  
% title（‘ACF with  2stds Bounds（i.e. approximate 95% confidence interval）‘）
% 
% subplot（212）
% [PACFLags