1、國內生產總值序列分析一問題的提出選取1978-2006歷年國內生產總值數據如下，試對該時間序列進行建模并預測。二問題分析與模型建立首先畫出數據的走勢圖，這一時間序列是具有明顯趨勢且不含有周期性變化經濟波動序列，即為非平穩的時間序列，對此序列進行建模預測需要用上面介紹的非平穩時間序列分析方法。采用模型：Xt=t+Yt其中t表示Xt中隨時間變化的趨勢值，Yt是Xt中剔除t后剩余部分。歷年國內生產總值時間序列圖三模型求解1確定性趨勢確定趨勢是按指數趨勢發展的t=abt-àlnt=lna+tlnb線性回歸分析程序：t=1978:2006;x=3624.10 4038.20 4517.80 4

2、862.40 5294.70 5934.50 7171.00 8964.40 10202.20 11962.50 14928.30 16909.20 18547.90 21617.80 26638.10 34634.40 46759.40 58478.10 67884.60 74462.60 78345.20 82067.46 89468.10 97314.80 105172.34 116898.40 136515.00 182321.00 209407.00;X=ones(29,1) t' %回歸的資料矩陣y=log(x)' %線性化B,BINT,R,RINT,STATS =

3、 regress(y,X) %回歸y2= exp(B(1)+B(2).*t) %預測值plot(t,x,t,y2,'+'); %回歸效果圖>>B =-290.48640.1510STATS = 1.0e+003 * 0.0010 2.1838 0 0.0000原始數據與指數回歸數據對比圖得到B=-290.4864,0.1510STATS=1.0e+003*0.0010,2.1838,0即由上圖可知僅用指數回歸的效果較差。2.隨機性趨勢（1）殘差序列Yt=Xt-tr=x-y2; %殘差數列plot(t,r,'O'); %殘差散點圖參差序列散點圖觀察殘差

4、序列的散點圖可知，該序列有很大的波動性，可認為是非平穩的，應該經過多次差分使其平穩。（2）二次差分后序列Yt=Yt-2Yt-1+Yt-2r1=diff(r); %殘差的一階差分r11=0 r1; %補數列差分后的項為0plot(t,r11,'o'); %一階差分散點圖r2=diff(r1); %二階差分r21=0 0 r2; %補數列差分后的項為0plot(t,r21,'o'); %二階差分散點圖一階差分散點圖二階差分散點圖（3）wt的時間序列分析A. 將序列r2t零均值化,序列wt的樣本自相關函數pk程序如下：w=r2-mean(r2); %零均值化gamao

5、=var(w); %求方差for j=1:27gama(j)=w(j+1:end)*w(1:end-j)'/27;endrho=gama/gamao %樣本自相關系數bar(rho) %條狀圖自相關系數條形圖B. 樣本偏相關函數kk程序如下：f(1,1)=rho(1);for k=2:27s1=rho(k);s2=1; %計算的初始值for j=1:k-1s1=s1-rho(k-j)*f(k-1,j);s2=s2-rho(j)*f(k-1,j);endf(k,k)=s1/s2; %對角上的樣本偏相關系數for j=1:k-1f(k,j)=f(k-1,j)-f(k,k)*f(k-1,k-

6、j); %不在對角上的樣本偏相關系數endendpcorr=diag(f)' %提取偏相關函數bar(pcorr) %條形圖偏自相關函數C. 模型定階的程序：for i=0:3for j=0:3spec= garchset('R',i,'M',j,'Display','off'); %指定模型的結構coeffX,errorsX,LLFX = garchfit(spec,w); %擬合參數num=garchcount(coeffX); %計算擬合參數的個數aic,bic=aicbic(LLFX,num,27);fprintf

7、('R=%d,M=%d,AIC=%f,BIC=%fn',i,j,aic,bic); %顯示計算結果endend結果如下：R=0,M=0,AIC=554.744695,BIC=557.336369R=0,M=1,AIC=548.981658,BIC=552.869169R=0,M=2,AIC=548.671841,BIC=553.855188R=0,M=3,AIC=550.112192,BIC=556.591376R=1,M=0,AIC=550.968125,BIC=554.855636R=1,M=1,AIC=550.239945,BIC=555.423293R=1,M=2,AI

8、C=551.360349,BIC=557.839534R=1,M=3,AIC=546.975261,BIC=554.750283R=2,M=0,AIC=552.918590,BIC=558.101938R=2,M=1,AIC=559.057147,BIC=565.536332R=2,M=2,AIC=551.171163,BIC=558.946184R=2,M=3,AIC=552.530372,BIC=561.601230R=3,M=0,AIC=553.140182,BIC=559.619367R=3,M=1,AIC=553.153087,BIC=560.928109R=3,M=2,AIC=56

9、1.426269,BIC=570.497127R=3,M=3,AIC=553.156375,BIC=563.523070得到結果顯示，可以認為是ARMA（1,3）模型D. 對模型wt=c+1wt-1+t+1t-1+2t-2+3t-3進行參數估計程序： spec = garchset('R',1,'M',3,'Display','off'); %指定模型的結構 coeffX,errorsX,LLFX = garchfit(spec,w) %擬合參數 運行結果如下:coeffX = Comment: 'Mean: ARMAX(

10、1,3,0); Variance: GARCH(0,0)' Distribution: 'Gaussian' R: 1 M: 3 C: 35.8925 AR: -0.6250 MA: -0.0387 0.0387 -1.0000 VarianceModel: 'GARCH' K: 3.5044e+007 Display: 'off'于是ARMA（1,3）模型為wt=35.8925-0.6250wt-1+t-0.0387t-1+0.0387t-2-1.0000t-3E. 模型的檢驗和預測程序： spec= garchset('R',1,'M',3); %指定模型的結構 coeff,errors,LLF,innovations,sigmas,summary = garchfit(spec,w) %擬合參數 h