matlab实习实验报告

实验一时间序列ARMA学习一、基本形式ARMA模型分为以下三种：1.自回归模型（AR）如果时间序列满足其中是独立同分布的随机变量序列，且满足：以及E()=0，则称时间序列为服从p阶的自回归模型。AR模型是用过去的观测值和现在的干扰值的线性组合来预测未来值，式子中的p为阶数，β为待定参数，y为一个平稳序列。自回归模型的平稳条件：滞后算子多项式的根均在单位圆外，即φ(B)=0的根大于1。2.移动平均模型（MA）如果时间序列满足则称时间序列为服从q阶移动平均模型；MA模型是用过去的干扰值和现在的干扰值的线性组合预测未来的值，式中q为阶数，α为待定参数，y为一个平稳序列。移动平均模型平稳条件：任何条件下都平稳。3.自回归滑动平均模型（ARMA）如果时间序列满足：则称时间序列为服从(p,q)阶自回归滑动平均混合模型。模型求解流程图：二、实例分析太阳的光球表面有时会出现一些暗的区域，它是磁场聚集的地方，这就是太阳黑子。黑子是太阳表面可以看到的最突出的现象。一个中等大小的黑子大概和地球的大小差不多。长期的观测发现，黑子多的时候，其他太阳活动现象也会比较频繁。从外文网站查询到太阳黑子数量的变化，从而预测接下来数量的走势。96.7126.5150116.739.299.573.328.341.781.771.7104.3125.8166.772.538.76653.315.785.5120.571.7116.7264.3142.375.547.5130.776.223.566.277.380.592.8142171.79473.348.863.335.354.27573.3141.7122.2152101.258.345.26043.7107.873.378139.2126.5109.584.583.377.752.85055.864.578.3158148.7105.5110.5118.362.736.763.562.7104.281.7110.5147.2125.799.798.866.76521.386.762.883.3这里只给出部分数据，原表格附在文档最后，此案例中选取最近50多年的数据局进行分析。1.AR模型求解首先整理数据，进行平稳性判断，使用matlab语言中的adftest()函数即可。数据为平稳数据，因而无需做平稳化处理，使用AR模型对原数据作出预测。Matlab代码为：clccleara=csvread('heizi.csv');b=a(2700:3234,[1,4]);c=b(:,2);diff(c);adftest(c)%figure(2);subplot(2,1,1)autocorr(c,500)%²set(gca,'Xlim',[0500]);subplot(2,1,2)parcorr(c,50);set(gca,'Xlim',[050]);z=iddata(c);%×m=armax(z(1:end),'na',20,'nc',0);%yp=predict(m,c,1);figure(3)plot(c,'-.');holdonplot(yp,'r')gridlegend('³õʼֵ','Ô¤²âÖµAR')bzc=sum((yp-c).^2)结果如下：2.MA模型求解Matlab代码：clccleara=csvread('heizi.csv');b=a(2700:3234,[1,4]);c=b(:,2);diff(c);adftest(c)figure(2);subplot(2,1,1)autocorr(c,500)set(gca,'Xlim',[0500]);subplot(2,1,2)parcorr(c,50);set(gca,'Xlim',[050]);z=iddata(c);m=armax(z(1:end),'na',0,'nc',4);yp=predict(m,c,1);figure(3)plot(c,'-.');holdonplot(yp,'r')gridlegend('³õʼֵ','Ô¤²âÖµMA')bzc=sum((yp-c).^2)结果如下：3.ARMA模型求解对于ARMA模型，需要做出阶数判断，用ACF和PACF确定p、q阶数,在使用ARMA模型求解。Matlab代码：clccleara=csvread('heizi.csv');b=a(2700:3234,[1,4]);c=b(:,2);diff(c);adftest(c)figure(2);subplot(2,1,1)autocorr(c,500)%²自相关set(gca,'Xlim',[0500]);subplot(2,1,2)parcorr(c,50);%偏自相关set(gca,'Xlim',[050]);结果如下：由图可以看出ACF和PACF都是拖尾，故可用ARMA模型求解，模拟数据matlab代码如下：a=csvread('heizi.csv');b=a(2700:3234,[1,4]);c=b(:,2);diff(c);adftest(c)%figure(2);subplot(2,1,1)autocorr(c,500)%set(gca,'Xlim',[0500]);subplot(2,1,2)parcorr(c,50);%set(gca,'Xlim',[050]);z=iddata(c);m=armax(z(1:end),'na',19,'nc',19);yp=predict(m,c,1);figure(3)plot(c,'-.');holdonplot(yp,'r')gridlegend('³õʼֵ','Ô¤²âÖµMA')bzc=sum((yp-c).^2)结果如下：比较图像可看出，ARMA数据的整体误差小于AR和MA数据的误差，因而ARMA方法的提出是有用的。再用误差平方和的比较定量的显示模型的差距。bzc=sum((yp-c).^2)用matlab计算得出数据：ARMAARMA3.1636e+057.5764e+052.4427e+05比较数据可得ARMA模型更加精确。三、改进方法通过讨论，我们决定对ARMA模型做改进，主要是在定阶和模型的误差方面作出改进。传统的ARMA模型定阶方法是比较繁琐的，对于有些不太好处理的数据，很难判断出符合的阶数，误差也比较大。对于上述问题，我们做了合理的改进，主要改进方法如下：1.定阶的方法（确定模型p,q）在传统的ARMA模型中，定阶方法有观察自相关偏相关图像确定阶数，aic最小原则确定阶数，等等。但是在实际操作中，我们发现这些步骤不好做，对于上述例子来说，自相关图像就不好观察，无法确定确切的拖尾截尾的点，因而无法判断。使用aic最小原则通常比较繁琐，数据量大的时候计算又耗费时间，在例题中为了得到aic最小的pq,我们做了实际计算，pq均在20以内时，计算完需要一分钟以上，效率不高，如果有更多的数据，将无法处理。对于这个问题，我们想到用神经网络训练让其判断数据的p，q。需要大量的原始数据来训练网络，用R语言生成符合时间序列阶数的大量数据，操作如下：R语言生成数据：arima.sim()分别取p,q为0,1,2,3,4,5生成符合的平稳数据，各生成1100个数据，1000个用来训练网络，100用来检测。Python代码如下：from__future__importprint_functionimportnumpyasnpimportrandomasrfromkeras.layers.coreimportDense,Dropout,Activationfromkeras.modelsimportSequentialfromkeras.modelsimportload_modeldefgetdata(p,q):'''生成数据的函数输入阶数p，q，输出对应的平稳序列'''ifp==0andq==0:a=r.random()returnnp.array([aforiinrange(20)])elifp==0andq!=0:elps,e=np.random.random(size=20),np.random.normal(size=20)e=e/np.max(e)cita=np.random.normal(0,1,q)cita=cita/np.max(cita)y=np.zeros(20)foriinrange(q,20):y[i]=elps[i]-np.dot(cita,e[i-q:i])y=y/np.max(y)returnyelifp!=0andq==0:elps,e=np.random.random(size=20),np.random.normal(size=20)e=e/np.max(e)fai=np.random.random(size=p)ifp=q:y=np.append(fai[0:p],np.zeros(20-p))foriinrange(p,20):y[i]=-np.dot(fai,y[i-p:i])+elps[i]y=y/np.max(y)returnyelse:elps,e=np.random.random(size=20),np.random.normal(size=20)e=e/np.max(e)fai,cita=np.random.random(size=p),np.random.normal(0,1,q)cita=cita/np.max(cita)ifp=q:y=np.append(fai[0:p],np.zeros(20-p))foriinrange(p,20):y[i]=-np.dot(fai,y[i-p:i])+elps[i]-\np.dot(cita,e[i-q:i])else:y=np.append(fai[0:p],np.zeros(20-p))foriinrange(q,20):y[i]=-np.dot(fai,y[i-p:i])+elps[i]-\np.dot(cita,e[i-q:i])y=y/np.max(y)returnydefgetlabel(p,q):'''对标签进行独热编码'''y=np.zeros(36)y[p*6+q]=1returny#生成数据并归一化x_train,y_train,x_test,y_test=[],[],[],[]x,y,x1=[],[],[]withopen('data1.txt','r')asf:lines=f.readlines()forlineinlines:x1.append(float(line))foriinrange(500):x_train.append(np.array(x1[i*20:i*20+20]))y_train.append(getlabel(0,2))x,y,x1=[],[],[]withopen('data2.txt','r')asf:lines=f.readlines()forlineinlines:x1.append(float(line))foriinrange(500):x_train.append(np.array(x1[i*20:i*20+20]))y_train.append(getlabel(2,0))ps,qs=[0,1,2,3,4,5],[0,1,2,3,4,5]forpinps:forqinqs:foriinrange(1000):x_train.append(getdata(p,q))y_train.append(getlabel(p,q))x_train=np.array(x_train)y_train=np.array(y_train)forpinps:forqinqs:foriinrange(100):x_test.append(getdata(p,q))y_test.append(getlabel(p,q))x_test=np.array(x_test)y_test=np.array(y_test)#设置初始随机数种