模型的識(shí)別與預(yù)測VIP

模型的識(shí)別與預(yù)測實(shí)驗(yàn)內(nèi)容依照某AR模型生成一段數(shù)據(jù)（1000），同時(shí)用另一MA模型生成一段數(shù)據(jù)（200），合成一段1200長度的數(shù)據(jù)1）依賴于這1200個(gè)數(shù)據(jù)的前800個(gè)數(shù)據(jù)，識(shí)別這段數(shù)據(jù)背后的AR模型。2）在1）的基礎(chǔ)上對(duì)新數(shù)據(jù)進(jìn)行預(yù)測，并通過后續(xù)的400個(gè)數(shù)據(jù)進(jìn)行判別（數(shù)據(jù)模型是否匹配）或者模型的修正（修正只需要提供思路和方法）。理論基礎(chǔ)1.時(shí)間序列模型介紹時(shí)間序列是隨時(shí)間改變而隨機(jī)地變化的序列。時(shí)間序列分析的目的是找出它的變化規(guī)律，即線性模型，主要有三種：AR模型（自回歸模型）、MA模型（滑動(dòng)平均模型）和ARMA模型（自回歸滑動(dòng)平均模型或混合模型）。設(shè){Xt}為零均值的實(shí)平穩(wěn)時(shí)間序列，階數(shù)為p的AR模型定義為其其中{}成為自回歸系數(shù)，白噪聲序列{}成為新信息序列；階數(shù)為q的MA模型定義為其中{}稱為滑動(dòng)平均系數(shù)；P階自回歸q階ARMA模型定義為記為ARMA（p，q）。模型的識(shí)別根據(jù)教材對(duì)平穩(wěn)時(shí)間序列的特性分析，對(duì)初步識(shí)別平穩(wěn)時(shí)間序列的類型提供了依據(jù)，如表1所示：表1各時(shí)間序列模型的特性類別模型AR（p）MA（q）ARMA（p，q）自相關(guān)函數(shù)拖尾截尾k=q處拖尾偏相關(guān)函數(shù)截尾k=p處拖尾拖尾模型階數(shù)的確定樣本自相關(guān)函數(shù)和樣本偏相關(guān)函數(shù)設(shè)有零均值平穩(wěn)時(shí)間序列{}的一段樣本觀測值，樣本協(xié)方差函數(shù)估計(jì)式為同理樣本自相關(guān)函數(shù)定義為MA模型階數(shù)的確定設(shè){}是正態(tài)的零均值平穩(wěn)MA（q）序列，而對(duì)于充分大的N，可以認(rèn)為的分布近似于正態(tài)分布，從而，的截尾性判斷如下：首先計(jì)算（?。?，因?yàn)閝值未知，故令q值從小到大，分別檢驗(yàn)滿足或的比例是否占總個(gè)數(shù)M的68.3%或95.5%。第一個(gè)滿足上述條件的q就是的截尾處，及MA（q）模型的階數(shù)。模型參數(shù)的估計(jì)當(dāng)選定模型及確定階數(shù)后，進(jìn)一步的問題就是估計(jì)出模型的位置參數(shù)。參數(shù)估計(jì)的方法有矩法，最小二乘法及極大似然估計(jì)等。在工程計(jì)算中，只要用樣本自協(xié)方差函數(shù)或樣本自相關(guān)函數(shù)中的一部分?jǐn)?shù)值。這里不詳細(xì)介紹。模型的檢驗(yàn)由樣本觀測序列{}，經(jīng)過模型的識(shí)別，階數(shù)的確定和參數(shù)估計(jì)，可以初步建立{}的模型。這樣建立的模型一般還需要進(jìn)行統(tǒng)計(jì)檢驗(yàn)，只有經(jīng)檢驗(yàn)確認(rèn)模型基本上能反映{}的統(tǒng)計(jì)特性時(shí)，用它進(jìn)行預(yù)測才能獲得良好的效果。模型的所謂自相關(guān)函數(shù)檢驗(yàn)法，其基本思想是，如果模型是正確的，則模型的估計(jì)值與實(shí)際觀測值所產(chǎn)生的殘差序列應(yīng)是隨機(jī)干擾產(chǎn)生的誤差，即{}應(yīng)是白噪聲序列。否則，模型不正確。模型的預(yù)報(bào)根據(jù)時(shí)間序列觀測數(shù)據(jù)，建立了一個(gè)與實(shí)際問題相適應(yīng)的模型后，就可以利用過去與現(xiàn)在的觀測值，對(duì)該序列未來時(shí)刻的取值進(jìn)行估計(jì)，即預(yù)報(bào)。關(guān)于預(yù)報(bào)效果優(yōu)劣的標(biāo)準(zhǔn)，下面采用的是在平穩(wěn)性最小方差意義下的預(yù)報(bào)。具體步驟從AR模型中提取1000個(gè)樣本點(diǎn)以及MA模型中提取200個(gè)樣本點(diǎn)，分別存儲(chǔ)于x2和x1中；取x2的前800個(gè)數(shù)據(jù)作為模型識(shí)別的樣本觀測值，存儲(chǔ)于X2中，計(jì)算樣本自相關(guān)函數(shù)和樣本偏相關(guān)函數(shù)，并通過它們的特性識(shí)別模型；根據(jù)自相關(guān)或者偏相關(guān)確定模型階數(shù)，用最小二乘法求模型系數(shù)；提取1000個(gè)數(shù)據(jù)中后200個(gè)數(shù)樣本據(jù)檢驗(yàn)所得模型是否正確，用卡方假設(shè)檢驗(yàn)，其自由度為M，顯著性水平為0.05；預(yù)測模型的新數(shù)據(jù)，并與原始數(shù)據(jù)相比較；通過后續(xù)的400個(gè)數(shù)據(jù)進(jìn)行判別（數(shù)據(jù)模型是否匹配）結(jié)論程序基本滿足實(shí)驗(yàn)要求，能通過一組樣本數(shù)據(jù)通過計(jì)算自相關(guān)函數(shù)和偏相關(guān)函數(shù)來識(shí)別其背后的模型，并確定模型的階數(shù)以及模型的系數(shù)，得到一個(gè)完整的時(shí)間序列。在此基礎(chǔ)上，可以通過該數(shù)據(jù)模型預(yù)報(bào)一組新的數(shù)據(jù)，與原始數(shù)據(jù)相比，能有較好的預(yù)報(bào)性，還能判別一組數(shù)據(jù)與這個(gè)模型是否匹配。我關(guān)于模型修正的思路：為將這后續(xù)400個(gè)數(shù)據(jù)與再從原始的800個(gè)數(shù)據(jù)提取出的400個(gè)數(shù)據(jù)合并，重新進(jìn)行一次模型的建立，理論上能得到較好的修正模型。分析總結(jié)盡管程序已基本完成實(shí)驗(yàn)內(nèi)容，但還是有一些不盡人意的地方，還存在如下問題：已給AR模型為2階，程序中未考慮到更高階的情況，即使識(shí)別到的是更高階模型，也作二階來處理，會(huì)存在誤差，這點(diǎn)有待完善；在確定模型階數(shù)時(shí)，偏相關(guān)函數(shù)的計(jì)算不夠準(zhǔn)確也可能導(dǎo)致階數(shù)的判定出錯(cuò)，但在后續(xù)加了檢驗(yàn)程序可在一定程度上解決該問題；對(duì)于一階情況也未作處理，樣本數(shù)據(jù)被判別為別的模型時(shí)也未處理，這點(diǎn)可以后續(xù)完善。附錄clearall;clc;%從MA模型X(t)=a(t)-1.3a(t-1)+0.4a(t-2)中提取200個(gè)樣本點(diǎn)K=200;a1=randn(1,K);x1=zeros(1,K);x1(1)=1;x1(2)=1;fori=3:Kx1(i)=a1(i)-1.3*a1(i-1)+0.4*a1(i-2);endX1=x1;%從AR模型X(t)+0.8X(t-1)=a(t)中提取1000個(gè)樣本點(diǎn)N=1000;a2=randn(1,N);x2=zeros(1,N);x2(1)=0.1;x2(2)=0.1;x2(3)=1;fori=3:Nx2(i)=1.5*x2(i-1)-0.5*x2(i-2)+a2(i);end%取x2的前800個(gè)數(shù)據(jù)作為模型識(shí)別的樣本觀測值L=800;X2=x2(1,1:L);%計(jì)算樣本自相關(guān)函數(shù)和樣本偏相關(guān)函數(shù)r=zeros(1,L);%樣本協(xié)方差矩陣fork=0:L-1%求樣本協(xié)方差rk=0;fori=1:L-krk=rk+X2(i)*X2(i+k);endr(k+1)=(1/L)*rk;endp=zeros(1,L);fori=1:L%求樣本自相關(guān)函數(shù)p(i)=r(i)/r(1);endM=L/10;pp=zeros(M,M);%樣本偏相關(guān)系數(shù)矩陣pp(1,1)=p(2);p(1)=[];fork=1:M-1%求樣本偏相關(guān)函數(shù)s1=0;s2=0;forj=1:ks1=s1+p(k+1-j)*pp(k,j);s2=s2+p(j)*pp(k,j);endpp(k+1,k+1)=(p(k+1)-s1)/(1-s2);forj=1:kpp(k+1,j)=pp(k,j)-pp(k+1,k+1)*pp(k,k+1-j);endendpp1=diag(pp);%提取偏相關(guān)系數(shù)subplot(2,1,1)plot(p);title('自相關(guān)圖');subplot(2,1,2)plot(pp1);title('偏相關(guān)圖');%識(shí)別模型，返回模型的階數(shù)%取M=L/10個(gè)自相關(guān)和偏相關(guān)系數(shù)count=0;forf1=1:Mifabs(p(f1))<=1/sqrt(L)count=count+1;endifcount/M>=0.683break;endendcount_1=0;forf2=1:Mifabs(pp1(f2))<=1/sqrt(L)count_1=count_1+1;endif(count_1)/M>=0.683breakendendiff1==M&&f2<Mdisp('建立AR模型');%確定模型階數(shù)fori=1:Mifabs(pp1(i))<=1/sqrt(L)jieshu=i-1;breakendendfprintf('其階數(shù)為:%d\n',jieshu);%最小二乘法求模型系數(shù)forn=1:jieshuY=X2';Y(1:n)=[];m=L-n;X=[];fori=1:mforj=1:nX(i,j)=X2(n+i-j);endendxishu=inv(X'*X)*(X'*Y);end%用后200個(gè)數(shù)據(jù)檢驗(yàn)所得模型是否正確X3=x2(1,L+1:end);%提取后200個(gè)樣本點(diǎn)T=N-L;a3=zeros(1,T);%由觀測樣本點(diǎn)計(jì)算得到的白噪聲點(diǎn)ra=ones(1,T+1);%白噪聲的協(xié)方差a3(1)=X3(1);a3(2)=X3(2)-xishu(1)*X3(1);ifjieshu>1fori=1:T-2%用模型估計(jì)值與樣本觀測值產(chǎn)生白噪聲序列a3(i+2)=X3(i+2)-xishu(1)*X3(i+1)-xishu(2)*X3(i);end%求白噪聲協(xié)方差M=T/10;fori=0:Ms3=0;forj=1:T-is3=s3+a3(j)*a3(j+i);endra(i+1)=s3/T;end%求白噪聲自相關(guān)rra=ones(1,M);fori=1:Mrra(i)=ra(i+1)/ra(1);end%卡方假設(shè)檢驗(yàn)，自由度為M，顯著性水平為0.05Qm=0;fori=1:MQm=Qm+(rra(i))^2;endQm=T*Qm;kafang=chi2inv(0.95,M);%求得卡方值ifQm<kafangdisp(['統(tǒng)計(jì)量','Qm=',num2str(Qm),'小于給定水平的卡方值，故所估計(jì)模型合適']);elsedisp(['統(tǒng)計(jì)量','Qm=',num2str(Qm),'大于給定水平的卡方值，故所估計(jì)模型不合適,請(qǐng)繼續(xù)執(zhí)行...']);endendendiff1==M&&f2==Mdisp('建立ARMA模型');endiff1<M&&f2<Mdisp('初始數(shù)據(jù)處理出錯(cuò)，重新執(zhí)行');endiff1<M&&f2==Mdisp('建立MA模型');fori=1:Mif(abs(p(i)))<=1/sqrt(L)jieshu=i-1;breakendendfprintf('其階數(shù)為:%d\n',jieshu);end%預(yù)測模型的新數(shù)據(jù)new_X=ones(1,200);new_X(1)=X2(end);new_X(2)=xishu(2)*X2(end)+xishu(1)*new_X(1);fori=1:198new_X(i+2)=xishu(1)*new_X(i+1)+xishu(2)*new_X(i);endfiguresubplot(2,1,1)plot(1:N,x2,'r.')title('原始數(shù)據(jù)')subplot(2,1,2)plot(1:200,new_X,'b*')holdplot(1:200,X3(1:200),'r.')legend('預(yù)測點(diǎn)','后200點(diǎn)')title('預(yù)測數(shù)據(jù)與模型生成后200樣本點(diǎn)比較')%通過后續(xù)的400個(gè)數(shù)據(jù)進(jìn)行判別（數(shù)據(jù)模型是否匹配）X4(1:200)=X1;X4(201:400)=X3;%提取后續(xù)400個(gè)數(shù)據(jù)U=400;a4=zeros(1,U);%由觀測樣本點(diǎn)計(jì)算得到的白噪聲點(diǎn)ra1=ones(1,U+1);%白噪聲的協(xié)方差a4(1)=X4(1);a4(2)=X4(2)-xishu(1)*X4(1);ifjieshu>1fori=1:U-2%用模型估計(jì)值與樣本觀測值產(chǎn)生白噪聲序列a4(i+2)=X4(i+2)-xishu(1)*X4(i+1)-xishu(2)*X4(i);end%求白噪聲協(xié)方差M=U/10;fori=0:Ms4=0;forj=1:U-is4=s4+a4(j)*a4(j+i);endra1(i+1)=s4/U;end%求白噪聲自相關(guān)rra1=ones(1,M);fori=1