求教时间序列分析

leadersheng

1. #灰色预测模型G(1,1)
   
2. GM11<-function(x0,t){ #x0为输入学列，t为预测个数
   
3.   x1<-cumsum(x0) #一次累加生成序列1-AG0序列(累加序列)
   
4.   b<-numeric(length(x0)-1)# 初始化长度为length(x0)-1的整数部分个numeric类型且值为0的数据集b
   
5.   n<-length(x0)-1# n 为length(x0)-1长度 因为需要生成MEAN（紧邻均值）生成序列 其长度少1
   
6.   for(i in 1:n){ #生成x1的紧邻均值生成序列
   
7.     b[i]<--(x1[i]+x1[i+1])/2
   
8.     b} #得序列b，即为x1的紧邻均值生成序列
   
9.   D<-numeric(length(x0)-1)
   
10.   D[]<-1
    
11.   B<-cbind(b,D)#作B矩阵
    
12.   BT<-t(B)#B矩阵转置
    
13.   M<-solve(BT%*%B)#求BT*B得逆
    
14.   YN<-numeric(length(x0)-1)
    
15.   YN<-x0[2:length(x0)]
    
16.   alpha<-M%*%BT%*%YN  #模型的最小二乘估计参数列满足alpha尖
    
17.   alpha2<-matrix(alpha,ncol=1)# 将结果变成一列
    
18.   # 得到方程的两个系数
    
19.   a<-alpha2[1]
    
20.   u<-alpha2[2]
    
21.   
    
22.   # 下面为结果输出
    
23.   # 输出参数估计值及模拟值
    
24.   cat("GM(1,1)参数估计值：",'\n',"发展系数-a=",-a,"  ","灰色作用量u=",u,'\n','\n') #利用最小二乘法求得参数估计值a,u
    
25.   y<-numeric(length(c(1:t)))# t为给定的预测个数
    
26.   y[1]<-x1[1] # 第一个数不变
    
27.   for(w in 1:(t-1)){  #将a,u的估计值代入时间响应序列函数计算x1拟合序列y
    
28.     y[w+1]<-(x1[1]-u/a)*exp(-a*w)+u/a
    
29.   }
    
30.   cat("x(1)的模拟值：",'\n',y,'\n')
    
31.   xy<-numeric(length(y))
    
32.   xy[1]<-y[1]
    
33.   for(o in 2:t){ #运用后减运算还原得模型输入序列x0预测序列
    
34.     xy[o]<-y[o]-y[o-1]
    
35.   }
    
36.   cat("x(0)的模拟值：",'\n',xy,'\n','\n')                       
    
37.   
    
38.   # 计算残差e
    
39.   e<-numeric(length(x0))
    
40.   for(l in 1:length(x0)){
    
41.     e[l]<-x0[l]-xy[l] #得残差序列（未取绝对值）
    
42.   }
    
43.   cat("绝对残差：",'\n',e,'\n')
    
44.   #计算相对误差
    
45.   e2<-numeric(length(x0))
    
46.   for(s in 1:length(x0)){
    
47.     e2[s]<-(abs(e[s])/x0[s]) #得相对误差
    
48.   }
    
49.   cat("相对残差：",'\n',e2,'\n','\n')
    
50.   cat("残差平方和=",sum(e^2),'\n')
    
51.   cat("平均相对误差=",sum(e2)/(length(e2)-1)*100,"%",'\n')
    
52.   cat("相对精度=",(1-(sum(e2)/(length(e2)-1)))*100,"%",'\n','\n')
    
53.   
    
54.   #后验差比值检验
    
55.   avge<-mean(abs(e));esum<-sum((abs(e)-avge)^2);evar=esum/(length(e)-1);se=sqrt(evar)  #计算残差的均方差se
    
56.   avgx0<-mean(x0);x0sum<-sum((x0-avgx0)^2);x0var=x0sum/(length(x0));sx=sqrt(x0var)  #计算原序列x0的方差sx
    
57.   cv<-se/sx  #得验差比值（方差比）
    
58.   cat("后验差比值检验:",'\n',"C值=",cv,'\n')#对后验差比值进行检验，与一般标准进行比较判断预测结果好坏。
    
59.   #计算小残差概率
    
60.   P<-sum((abs(e)-avge)<0.6745*sx)/length(e)
    
61.   cat("小残差概率:",'\n',"P值=",P,'\n')
    
62.   
    
63.   if(cv < 0.35 && P>0.95){     
    
64.     cat("C<0.35, P>0.95,GM(1,1)预测精度等级为：好",'\n','\n')
    
65.   }else{
    
66.     if(cv<0.5 && P>0.80){
    
67.       cat("C值属于[0.35,0.5), P>0.80,GM(1,1)模型预测精度等级为：合格",'\n','\n')
    
68.     }else{
    
69.       if(cv<0.65 && P>0.70){
    
70.         cat("C值属于[0.5,0.65), P>0.70,GM(1,1)模型预测精度等级为：勉强合格",'\n','\n')
    
71.       }else{
    
72.         cat("C值>=0.65, GM(1,1)模型预测精度等级为：不合格",'\n','\n')
    
73.       }
    
74.     }
    
75.   }
    
76.   #画出输入序列x0的预测序列及x0的比较图像
    
77.   plot(xy,col='blue',type='b',pch=16,xlab='时间序列',ylab='值')
    
78.   points(x0,col='red',type='b',pch=4)
    
79.   legend("topleft",c('预测','原始'),title = "预测时序与原始时序对比",pch=c(16,4),lty=l,col=c('blue','red'))
    
80. }
    
81. x<-c(-6.05，1.95，-4.04，-17.02，3.08，-0.88，-7.10，-9.96，-9.77，-2.34，1.37，0.45，-6.22，6.57，0.14，0.93，-2.50，-5.57，-1.27，3.60，2.52，0.79，-0.63，2.51，-2.83，2.69，-1.99，1.34，0.09，-3.22，4.24，2.54，9.27，-1.77，5.40，0.83，4.84，0.36，10.60，8.88，-0.82，-3.64，-2.18，7.38，11.73，7.75，-0.89，0.22，4.28，-2.05 )
    
82. GM11(x,length(x)+2)
    
83. 
    
84. 
    
85. 
    
86. -6.05，1.95，-4.04，-17.02，3.08，-0.88，-7.10，-9.96，-9.77，-2.34，1.37，0.45，-6.22，6.57，0.14，0.93，-2.50，-5.57，-1.27，3.60，2.52，0.79，-0.63，2.51，-2.83，2.69，-1.99，1.34，0.09，-3.22，4.24，2.54，9.27，-1.77，5.40，0.83，4.84，0.36，10.60，8.88，-0.82，-3.64，-2.18，7.38，11.73，7.75，-0.89，0.22，4.28，-2.05
    
87. 
    
88. 
    
89. 
    
90. 
    
91. 
    
92. 
    
93. 
    
94. 
    
95. 
    
96. 
    
97. 
    
98. 
    
99. 
    
100. 
     
101. 
     
102. 
     
103. 
     
104. 
     
105. 
     
106. 
     
107. 
     
108. 
     
109. 
     
110. 
     
111. 
     
112. 这是灰色预测的MATLAB程序

复制代码

leadersheng

时间序列模型
   

leadersheng

1. function [yucezhi]=shijianxulie(data,cishu);
   
2. [hang lie]=size(data);
   
3. yucezhi=zeros(hang,lie);
   
4. z=zeros(hang,cishu);
   
5. for nt=1:hang
   
6.     x=data(nt,:);
   
7. xx=x;
   
8. x=[0 diff(x)];
   
9. [NnN Nn]=size(x);
   
10. %%%%%%用AR模型
    
11. for n=1:sqrt(Nn)
    
12.     x2=x((n+1):Nn);
    
13.     for j=1:(Nn-n)
    
14.         for ji=1:n
    
15.             x1(j,ji)=x(n-ji+j);
    
16.         end
    
17.     end
    
18. dd=pinv((x1'*x1))*x1'*x2';
    
19.    
    
20.     p=dd;%%%%%%%%p为行向量
    
21.     for j=n+1:Nn
    
22.         for ii=1:n
    
23.             ax(ii)=p(ii)*x(j-ii);%%%%拟合出的xt
    
24.         end
    
25.         ct(j)=x(j)-sum(ax);%%%%%%%%白噪声序列
    
26.     end
    
27.     t=[];
    
28.     t=ct(n+1:Nn).^2;
    
29.     b(n)=sum(t)/(Nn-n);
    
30.     AIC(n)=log(b(n))+2*n/Nn;%%%%%%%%AIC准则
    
31.     x1=[];x2=[];
    
32. end
    
33. %%%%%%%%%%%%%%%%%求AIC的最小值
    
34. pp=AIC(1);
    
35. nn=1;
    
36. for i=1:sqrt(Nn)
    
37.     if pp>AIC(i)
    
38.         pp=AIC(i);
    
39.         nn=i;
    
40.     end
    
41. end
    
42. %%%%%%%%%%%%%%%AIC中最小值为pp，回归节数为nn
    
43. %%%%%%%%%求白噪声序列%%%%%%
    
44. for j=nn+1:Nn
    
45.     for ii=1:nn
    
46.         ax(ii)=p(ii)*x(j-ii);
    
47.     end
    
48.     ct(j)=x(j)-sum(ax);%%%%%%%%%白噪声序列%%%%%%
    
49. end
    
50. at0(nn+1:Nn)=ct(nn+1:Nn);
    
51. %%%%%%%%%%%%%估计参数aa bb
    
52. jies=12;
    
53. if jies>=Nn
    
54.     jies=Nn-1;
    
55. end
    
56. for p=1:jies
    
57.     for q=1:jies
    
58.         hh=[p,q,nn];
    
59.         L=max(hh);
    
60.         
    
61.         X=[];e1=[];YY=[];Y=[];m=[];n=[];aabb=[];pa=[];pb=[];
    
62.         for j=1:Nn-L
    
63.             for km=1:p
    
64.                 X(j,km)=x(L+j-km);
    
65.             end
    
66.         end
    
67.         
    
68.         
    
69.         for j=1:Nn-L
    
70.             for jm=1:q
    
71.                 e1(j,jm)=at0(L-jm+j);
    
72.             end
    
73.         end
    
74.         
    
75.        for i=1:Nn-L
    
76.            YY(i)=x(L+i);
    
77.        end
    
78.        Y=YY';
    
79.      XX=[X e1];
    
80.       
    
81.       
    
82.      aabb=pinv(XX)*Y;
    
83.         
    
84.         
    
85.       
    
86.         pa=aabb(1:p,1);%%%%%%%%%%%pa为自回归系数
    
87.         pb=aabb(p+1:p+q,1);%%%%%%%%%%%pb为移动平均系数
    
88.     %%%%%%%%%%%%计算AICB（p,q)
    
89.         for j=L+1:Nn
    
90.             for ii=1:p
    
91.                 ax(ii)=pa(ii)*x(j-ii);
    
92.             end
    
93.             for kk=1:q
    
94.                 bx(kk)=pb(kk)*x(j-kk);
    
95.             end
    
96.             abtt(j)=x(j)-sum(ax)-sum(bx);%%%%%%白噪声序列%%%%%%
    
97.         end
    
98.         t=[];
    
99.         t=abtt(L+1:Nn).^2;
    
100.         ab(p,q)=sum(t)/(Nn-L);
     
101.         AICB(p,q)=log(ab(p,q))+2*(p+q)/Nn;
     
102.     end
     
103. end
     
104. %%%%%%%%%%%%%%%%求AICB(p,q)的最小值
     
105. pppp=AICB(1,1);
     
106. nna=1;
     
107. nnb=1;
     
108. for ik=1:jies
     
109.     for jk=1:jies
     
110.         if pppp>AICB(ik,jk)
     
111.             pppp=AICB(ik,jk);
     
112.             nna=ik;
     
113.             nnb=jk;
     
114.         end
     
115.     end
     
116. end
     
117. %%%%%%%%%%%%%%%%确定模型为ARMA（p,q)%%%%%%%%%%求自回归系数和移动平均系数
     
118. 
     
119. X=[];e1=[];
     
120. p=nna;
     
121. q=nnb;
     
122. hh=[p,q];
     
123. L=max(hh);
     
124.   
     
125. 
     
126. for j=1:Nn-L
     
127.             for km=1:p
     
128.                 X(j,km)=x(L+j-km);
     
129.             end
     
130. end   
     
131. 
     
132. for j=1:Nn-L
     
133.             for jm=1:q
     
134.                 e1(j,jm)=at0(L-jm+j);
     
135.             end
     
136. end
     
137. 
     
138.    for i=1:Nn-L
     
139.            YY(i)=x(L+i);
     
140.        end
     
141.        Y=YY';
     
142.       
     
143.        XX=[X e1];
     
144.         aabb=pinv(XX)*Y;
     
145.         
     
146.       
     
147.       abc=aabb;
     
148.         pma=[];pmb=[];
     
149.         pma=abc(1:p,1)';%%%%%%%%%自回归系数（行向量
     
150.         pmb=abc(p+1:p+q,1)';%%%%%%%%%%移动平均系数（行向量
     
151.         %%%%%%%%%%%模型ARMA（p,q)%%%%%%自回归系数pma%%%%%%%移动平均系数pmb
     
152.         %%%
     
153.      
     
154.    %%signa=sum( Y-X*pma'-e1*pmb').^2/(Nn-L)
     
155.   %%%%%%%%%%%%%%%预测%%%%%%%%%%%%%%%
     
156.   
     
157.   yyuce=[x(1:L)';X*pma'+e1*pmb'];
     
158.   
     
159.   
     
160.    ee=(Y-X*pma'-e1*pmb');
     
161.     ee=[zeros(L,1);ee; zeros(cishu,1)]';%%%%%%%白噪声（L+1:Nn)（行向量
     
162.    
     
163.     c=zeros(1,p);
     
164.     cc=zeros(1,q);
     
165.     for i=1:p
     
166.         c(i)=x(Nn-i+1);
     
167.     end
     
168.     for i=1:q
     
169.         cc(i)=ee(Nn-i+1);
     
170.     end
     
171.    z1=c*pma'+cc*pmb';
     
172.    z(nt,1)=z1;
     
173.   for k=2:cishu
     
174.       if p==1
     
175.           z2=z1;
     
176.           z1=z2*pma'+cc*pmb';
     
177.       end
     
178.    if p>1      
     
179.    for i=1:p-1
     
180.         c(p-i+1)=c(p-i);
     
181.    end
     
182.     c(1)=z1;
     
183.     for i=1:q
     
184.         cc(i)=ee(Nn-i+k);
     
185.     end
     
186.         z1=c*pma'+cc*pmb';
     
187.    end
     
188.    z(nt,k)=z1;   
     
189.   end
     
190. 
     
191.    z(nt,1)=z(nt,1)+xx(Nn);
     
192.    for k=2:cishu
     
193.        z(nt,k)=z(nt,k)+z(nt,k-1);
     
194.    end
     
195.    yyuce=(yyuce)';
     
196.    yuce(1)=xx(1);
     
197.    for i=1:Nn-1
     
198.        yuce(i+1)=yyuce(i)+xx(i+1);
     
199.    end
     
200.    
     
201.    yucezhi(nt,:)=yuce;
     
202. 
     
203. end
     
204.    yucezhi=[yucezhi z];
     
205. 
     
206. end

复制代码