ARMA.txt
上传用户:zjjxqc
上传日期:2014-11-15
资源大小:2k
文件大小:4k
源码类别:

书籍源码

开发平台:

Matlab

  1. % 本程序的目的是模拟一个ARMA模型,然后进行时频归并。考察归并前后模型的变化。
  2. % 这个ARMA模型的一般形式用黑盒子模型表示为A(q)y(t)=C(q)e(t)。q是滞后算子。
  3. % 或者是:(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) 
  4. % 这里多项式A和C都只写出4阶,因为一般的经济时间序列阶数都不高。
  5. clear
  6. tic
  7. % ====================第一步,模拟一个ARMA模型并绘制ACF,PACF图========================
  8. %s首先设定ARMA模型的多项式系数。ARMA模型中只有多项式A(q)和C(q),
  9. %把A(q)的系数都设为0就得到MA模型,把C(q)的系数都设为0就得到AR模型。
  10. a1 = -(0.5)^(1/3);
  11. a2 = (0.5)^(2/3);
  12. a3 = 0;
  13. a4 = 0;
  14. c1 = 0;
  15. c2 = 0;
  16. c3 = 0;
  17. c4 = 0; 
  18. obv = 3000;                       %obv是模拟的观测数目。 
  19. A = [1 a1 a2 a3 a4];
  20. B = [];                           %因为ARMA模型没有输入,因此多项式B是空的。
  21. C = [1 c1 c2 c3 c4];
  22. D = [];                           %把D也设为空的。
  23. F = [];                           %ARMA模型里的F多项式也是空的。
  24. m = idpoly(A,B,C,D,F,1,1)         %这样就生成了ARMA模型,把它存储在m中。NoiseVariance被设定为1,1也是默认值。抽样间隔Ts设为1。 
  25. error = randn(obv,1);             %生成一个obv*1的正态随机序列。准备用作模型的误差项。
  26. e = iddata([],error,1);           %用randn函数生成一个噪声序列。存储在e中。抽样间隔是1秒。
  27. %u = [];                          %因为是ARMA模型,没有输出。所以把u设为空的。这句可省略。
  28. y = sim(m,e); 
  29. get(y)                            %使用get函数来查看动态系统的所有性质。
  30. r=y.OutputData;                   %把y.OutputData的全部值赋给变量r,r就是一个obv*1的向量。 
  31. figure(1)
  32. plot(r)                           %绘出y随时间变化的曲线。或者写成plot(y)结果是一样的。 
  33. %==========================第二步,绘制ARMA序列r的ACF和PACF图=======================
  34. %如果直接用函数autocorr和parcorr画图,总是会把滞后一阶的ACF和PACF也画出来,不好看,所以用下述方法画。
  35. % figure(2)
  36. % subplot(2,1,1)
  37. % n=100;
  38. % [ACF,Lags,Bounds]=autocorr(r,n,2);
  39. % x=Lags(2:n);
  40. % y=ACF(2:n);                            %注意这里的y和前面y的完全不同。
  41. % h=stem(x,y,'fill','-');
  42. % set(h(1),'Marker','.')
  43. % hold on
  44. % ylim([-1 1]); 
  45. % a=Bounds(1,1)*ones(1,n-1);
  46. % line('XData',x,'YData',a,'Color','red','linestyle','--')
  47. % line('XData',x,'YData',-a,'Color','red','linestyle','--')
  48. % xlabel('lags')
  49. % ylabel('ACF')  
  50. % title('ACF with  2stds Bounds(i.e., approximate 95% confidence interval)')
  51. % subplot(2,1,2)
  52. % [PACF,Lags,Bounds]=parcorr(r,n,2);
  53. % x=Lags(2:n);
  54. % y=PACF(2:n);
  55. % h=stem(x,y,'fill','-');
  56. % set(h(1),'Marker','.')
  57. % hold on
  58. % ylim([-1 1]);
  59. % b=Bounds(1,1)*ones(1,n-1);
  60. % line('XData',x,'YData',b,'Color','red','linestyle','--')
  61. % line('XData',x,'YData',-b,'Color','red','linestyle','--')
  62. % xlabel('lags')
  63. % ylabel('PACF')  
  64. % title('PACF with  2stds Bounds(i.e., approximate 95% confidence interval)') 
  65. %%============================第三步,对r进行m阶时频归并========================
  66. %%这一步对前面得到的ARMA序列r进行m阶时频归并。归并的方法是把每m个数据项加起来。
  67. %%m的取值根据实际情况来决定。 
  68. m = 3;                            % m=3对应于把月度数据归并为季度数据。
  69. R = reshape(r,m,obv/m);           % 把向量 r 变形成m*(obv/m)的矩阵R.如果obv=3000,m=3,则R为3*1000的矩阵。
  70. aggregatedr = sum(R);             % sum(R)计算矩阵R每一列的和。得到的1*(obv/m)行向量aggregatedr就是时频归并后得到的序列。
  71. dlmwrite('E:Temporal AggregationEssays and ProgramsMatlab codeoutput.txt',aggregatedr','delimiter','t','precision',6,'newline','pc');
  72. % 至此完成了对r的时频归并。 
  73. %%=====================第四步,画出aggregatedr的ACF和PACF图以判断其模型============
  74. figure(3)
  75. subplot(2,1,1) 
  76. n=100;
  77. bound = 1;
  78. [ACF,Lags,Bounds]=autocorr(aggregatedr,n,2);
  79. x=Lags(2:n);
  80. y=ACF(2:n);                            %注意这里的y和前面y的完全不同。
  81. h=stem(x,y,'fill','-');
  82. set(h(1),'Marker','.')
  83. hold on
  84. ylim([-bound bound]); 
  85. a=Bounds(1,1)*ones(1,n-1);
  86. line('XData',x,'YData',a,'Color','red','linestyle','--')
  87. line('XData',x,'YData',-a,'Color','red','linestyle','--')
  88. xlabel('lags')
  89. ylabel('ACF')  
  90. title('ACF with  2stds Bounds(i.e., approximate 95% confidence interval)') 
  91. subplot(2,1,2)
  92. [PACF,Lags,Bounds]=parcorr(aggregatedr,n,2);
  93. x=Lags(2:n);
  94. y=PACF(2:n);
  95. h=stem(x,y,'fill','-');
  96. set(h(1),'Marker','.')
  97. hold on
  98. ylim([-bound bound]);
  99. b=Bounds(1,1)*ones(1,n-1);
  100. line('XData',x,'YData',b,'Color','red','linestyle','--')
  101. line('XData',x,'YData',-b,'Color','red','linestyle','--')
  102. xlabel('lags')
  103. ylabel('PACF')  
  104. title('PACF with  2stds Bounds(i.e., approximate 95% confidence interval)') 
  105. t=toc