ARMA.txt
上传用户:zjjxqc
上传日期:2014-11-15
资源大小:2k
文件大小:4k
- % 本程序的目的是模拟一个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(A,B,C,D,F,1,1) %这样就生成了ARMA模型,把它存储在m中。NoiseVariance被设定为1,1也是默认值。抽样间隔Ts设为1。
- error = randn(obv,1); %生成一个obv*1的正态随机序列。准备用作模型的误差项。
- e = iddata([],error,1); %用randn函数生成一个噪声序列。存储在e中。抽样间隔是1秒。
- %u = []; %因为是ARMA模型,没有输出。所以把u设为空的。这句可省略。
- y = sim(m,e);
- 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(2,1,1)
- % n=100;
- % [ACF,Lags,Bounds]=autocorr(r,n,2);
- % x=Lags(2:n);
- % y=ACF(2:n); %注意这里的y和前面y的完全不同。
- % h=stem(x,y,'fill','-');
- % set(h(1),'Marker','.')
- % hold on
- % ylim([-1 1]);
- % a=Bounds(1,1)*ones(1,n-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(2,1,2)
- % [PACF,Lags,Bounds]=parcorr(r,n,2);
- % x=Lags(2:n);
- % y=PACF(2:n);
- % h=stem(x,y,'fill','-');
- % set(h(1),'Marker','.')
- % hold on
- % ylim([-1 1]);
- % b=Bounds(1,1)*ones(1,n-1);
- % line('XData',x,'YData',b,'Color','red','linestyle','--')
- % line('XData',x,'YData',-b,'Color','red','linestyle','--')
- % xlabel('lags')
- % ylabel('PACF')
- % title('PACF with 2stds Bounds(i.e., approximate 95% confidence interval)')
- %%============================第三步,对r进行m阶时频归并========================
- %%这一步对前面得到的ARMA序列r进行m阶时频归并。归并的方法是把每m个数据项加起来。
- %%m的取值根据实际情况来决定。
- m = 3; % m=3对应于把月度数据归并为季度数据。
- R = reshape(r,m,obv/m); % 把向量 r 变形成m*(obv/m)的矩阵R.如果obv=3000,m=3,则R为3*1000的矩阵。
- aggregatedr = sum(R); % sum(R)计算矩阵R每一列的和。得到的1*(obv/m)行向量aggregatedr就是时频归并后得到的序列。
- dlmwrite('E:Temporal AggregationEssays and ProgramsMatlab codeoutput.txt',aggregatedr','delimiter','t','precision',6,'newline','pc');
- % 至此完成了对r的时频归并。
- %%=====================第四步,画出aggregatedr的ACF和PACF图以判断其模型============
- figure(3)
- subplot(2,1,1)
- n=100;
- bound = 1;
- [ACF,Lags,Bounds]=autocorr(aggregatedr,n,2);
- x=Lags(2:n);
- y=ACF(2:n); %注意这里的y和前面y的完全不同。
- h=stem(x,y,'fill','-');
- set(h(1),'Marker','.')
- hold on
- ylim([-bound bound]);
- a=Bounds(1,1)*ones(1,n-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(2,1,2)
- [PACF,Lags,Bounds]=parcorr(aggregatedr,n,2);
- x=Lags(2:n);
- y=PACF(2:n);
- h=stem(x,y,'fill','-');
- set(h(1),'Marker','.')
- hold on
- ylim([-bound bound]);
- b=Bounds(1,1)*ones(1,n-1);
- line('XData',x,'YData',b,'Color','red','linestyle','--')
- line('XData',x,'YData',-b,'Color','red','linestyle','--')
- xlabel('lags')
- ylabel('PACF')
- title('PACF with 2stds Bounds(i.e., approximate 95% confidence interval)')
- t=toc