《經(jīng)濟(jì)預(yù)測與決策技術(shù)及MATLAB實現(xiàn)》課件 第7-9章 時間序列預(yù)測法、干預(yù)模型預(yù)測法、馬爾可夫預(yù)測方法

經(jīng)濟(jì)預(yù)測與決策技術(shù)及MATLAB實現(xiàn)第7章時間序列預(yù)測法

7.1移動平均值預(yù)測法

7.2指數(shù)平滑預(yù)測法7.3季節(jié)指數(shù)預(yù)測法7.4時間序列分解法練習(xí)與提高（七）

7.5ARMA模型預(yù)測法7.6案例分析

7.1移動平均值預(yù)測法7.1.1一次移動平均法（1）一次移動平均法模型一次移動平均法是收集一組觀察值，計算這組觀察值的均值，利用這一均值作為下一期的預(yù)測值。其模型：其中，為t期的實際值；N為所選數(shù)據(jù)個數(shù)，為下一期（t＋1）的預(yù)測值?！纠?-1】我國2005--2020年的房地產(chǎn)開發(fā)投資數(shù)據(jù)如表，試用一次移動平均法預(yù)測2021年的房地產(chǎn)投資額（取N=3）。MATLAB程序clearX=[1.59 1.94 2.53 3.12 3.62 4.83 6.18 7.18...8.60 9.50 9.60 10.26 10.98 12.02 13.22 14.14];年份20052006200720082009201020112012投資額1.591.942.533.123.624.836.187.18年份20132014201520162017201820192020投資額8.609.509.6010.2610.9812.0213.2214.14N=3;fort=3:length(X)M1(t)=(X(t)+X(t-1)+X(t-2))/N;%一次移動平均值

X1(t+1)=M1(t);%下一期預(yù)測值endM1,X1t1=1:length(X);t2=4:length(X)+1plot(t1,X,'-+',t2,X1(4:end),'-O')xlabel('時間/年')ylabel('投資額/萬億元')legend('原始數(shù)據(jù)','預(yù)測值')X1=0 00 2.0200 2.5300 3.0900 3.8567 4.8767 6.0633 7.3200 8.4267 9.2333 9.7867 10.2800 11.0867 12.0733 13.12677.1.2二次移動平均法

（1）二次移動平均法的線性模型其中，為t期的實際值，為t＋T期的預(yù)測值，t為當(dāng)前的時期數(shù)，

T為由t至預(yù)測期的時期數(shù)?！纠?-2（續(xù)【例7-1】）利用二次移動平均法預(yù)測2021年房地產(chǎn)投資額（取N=3）。（1）先計算一次、二次移動平均值clearX=[1.59 1.94 2.53 3.12 3.62 4.83 6.18 7.18...8.60 9.50 9.60 10.26 10.98 12.02 13.22 14.14];N=3;fort=3:length(X)%t從3開始，前兩項用0占位置

M1(t)=(X(t)+X(t-1)+X(t-2))/N;%一次移動平均值endM1%t從5開始是因為M1的前2項為0M2(t)=(M1(t)+M1(t-1)+M1(t-2))/N;%二次移動平均值endM2（2）給出2010--2021年的預(yù)測值，并繪出預(yù)測圖。a=2*M1(5:end)-M2(5:end);b=2*(M1(5:end)-M2(5:end))/(N-1);T=1;y=a+b*Tt1=1:length(X);t2=6:length(X)+1;plot(t1,X,'-+',t2,y,'-O')xlabel('時間/年')ylabel('投資額/萬億元')legend('原始數(shù)據(jù)','預(yù)測值')y=4.17675.25226.74788.32569.786710.740011.046711.062211.306712.491113.926715.18897.2指數(shù)平滑預(yù)測法7.2.1一次指數(shù)平滑法（1）一次指數(shù)平滑法的基本模型

其中，為時間序列觀測值，為觀測值的指數(shù)平滑值，為平滑系數(shù)，?！纠?-3】（續(xù)例7-1）利用一次指數(shù)平滑法預(yù)測2021年的房地產(chǎn)投資額（=0.7、0.8、0.9）。clearX=[1.59 1.94 2.53 3.12 3.62 4.83 6.18 7.18...8.60 9.50 9.60 10.26 10.98 12.02 13.22 14.14];X0=X(1);X1=X(2:end);alpha=0.9;S0=X0;%初始值S1(1)=alpha*X1(1)+(1-alpha)*S0;%指數(shù)平滑值第一項fort=1:length(X1)-1S1(t+1)=alpha*X1(t+1)+(1-alpha)*S1(t);endS1%指數(shù)平滑值全部項S=[S0S1];MSE=sum((X1-S(1:length(X1))).^2)./length(X1)%均方誤差t1=1:length(X);t2=2:length(X)+1;plot(t1,X,'-+',t2,S,'-O')xlabel('時間/年');ylabel('投資額/萬億元')legend('原始數(shù)據(jù)','預(yù)測值')alpha=0.9，MSE=1.0020alpha=0.8，MSE=1.2294alpha=0.7，MSE=1.5482首頁

7.2.2二次指數(shù)平滑法（1）二次指數(shù)平滑法的線性模型為

【例7-4】（續(xù)例7-1）試用二次指數(shù)平滑法預(yù)測2021年房地產(chǎn)投資額（alpha＝0.9）。1）先計算一次、二次指數(shù)平滑值clearX=[1.59 1.94 2.53 3.12 3.62 4.83 6.18 7.18...8.60 9.50 9.60 10.26 10.98 12.02 13.22 14.14];X0=X(1);X1=X(2:end);alpha=0.9;S10=X0;%S1的初始值S1(1)=alpha*X1(1)+(1-alpha)*S10;%一次指數(shù)平滑值第一項fort=1:length(X1)-1S1(t+1)=alpha*X1(t+1)+(1-alpha)*S1(t);%一次指數(shù)平滑值第二項以后項endS1S20=X0;%S2的初始值S2(1)=alpha*S1(1)+(1-alpha)*S20;%二次指數(shù)平滑值第一項fort=1:length(S1)-1S2(t+1)=alpha*S1(t+1)+(1-alpha)*S2(t);%二次指數(shù)平滑值第二項以后項endS2（2）預(yù)測2021年及以前的全部預(yù)測值a=2*S1(1:end)-S2(1:end);b=alpha/(1-alpha)*(S1(1:end)-S2(1:end));T=1;y=a+b*TY=[X0,y];MSE=sum((X(2:end)-Y(1:end-1)).^2)./(length(X)-1)t1=1:length(X);t2=2:length(X)+1;plot(t1,X,'-+',t2,Y,'-O')xlabel('時間/年')ylabel('投資額/萬億元')legend('原始數(shù)據(jù)','預(yù)測值')為兩個平滑參數(shù)，取值范圍都在（0,1）區(qū)間;為第t期平滑值;為第t期趨勢值。7.2.4雙參數(shù)線性指數(shù)平滑法霍爾特雙參數(shù)線性指數(shù)平滑法的基本原理與二次指數(shù)平滑法相似，只不過它并不是進(jìn)行二次指數(shù)平滑，而是對序列的趨勢直接進(jìn)行平滑。雙參數(shù)線性指數(shù)平滑法的模型：【例7-6】（續(xù)【例7-1】）試用雙參數(shù)線性指數(shù)平滑法預(yù)測2021年房地產(chǎn)投資額（＝0.9）clearX=[1.59 1.94 2.53 3.12 3.62 4.83 6.18 7.18...8.60 9.50 9.60 10.26 10.98 12.02 13.22 14.14];alpha=0.9;gamma=0.9;S0=X(1);b0=(X(end)-X(1))/(length(X)-1);S(1)=alpha*X(1)+(1-alpha)*(S0+b0);b(1)=gamma*(S(1)-S0)+(1-gamma)*b0;fort=2:length(X)S(t)=alpha*X(t)+(1-alpha)*(S(t-1)+b(t-1));b(t)=gamma*(S(t)-S(t-1))+(1-gamma)*b(t-1);endfort=1:length(X)Y(t)=S(t)+b(t);endY%霍爾特雙參數(shù)線性指數(shù)平滑法預(yù)測值MSE=sum((X(2:end)-Y(1:end-1)).^2)./(length(X)-1)t1=1:length(X);t2=2:length(X)+1;plot(t1,X,'-+',t2,Y,'-o')7.3季節(jié)指數(shù)預(yù)測法首頁7.3.1季節(jié)性水平模型

如果時間序列沒有明顯的趨勢變動，而主要受季節(jié)變化和不規(guī)則變動影響時，可用季節(jié)性水平模型進(jìn)行預(yù)測。預(yù)測模型的方法：

（1）計算歷年同季的平均數(shù)（2）計算全季總平均數(shù)（3）計算各季的季節(jié)指數(shù)歷年同季的平均數(shù)與全時期的季平均數(shù)之比，即：若各季的季節(jié)指數(shù)之和不為4，季節(jié)指數(shù)需要調(diào)整為（4）利用季節(jié)指數(shù)法進(jìn)行預(yù)測【例7-7我國2015-2021年城鎮(zhèn)居民人均消費性支出28個季度數(shù)據(jù)如表7-2所示，并利用2021年第4季度數(shù)據(jù)作為基期數(shù)據(jù)，預(yù)測2022年四個季度居民人均消費性支出。

季度年份123420155534486752355756201659705215561262822017638755445915659920186749599662697098201971606405681476842020647860076762776020217495707174158326

（1）根據(jù)所給數(shù)據(jù)，畫出走勢圖，觀察季節(jié)性clearX=[5534 4867 5235 57565970 5215 5612 62826387 5544 591565996749 5996 6269 70987160 6405 681476846478 6007 6762 77607495 7071 7415 8326];Y=X';%將年度數(shù)據(jù)轉(zhuǎn)化為按列、季度數(shù)據(jù)按行排列Z=Y(:);%矩陣轉(zhuǎn)換為一列向量t=1:length(Z);plot(t,Z,'-o')xlabel('時間/季度')ylabel('人均消費支出/元')（2）計算季節(jié)指數(shù)并預(yù)測r=mean(X)%同季平均數(shù)R=mean(X(:))%全部季度平均數(shù)b=r./R%各季季節(jié)指數(shù)F=4/sum(b)*b%調(diào)整各季季節(jié)指數(shù)%以2021年第4季度作為基期I=X(end)*(F./F(4))%2022年第1~4季度預(yù)測值r=1.0e+03*6.53905.87216.28897.0721R=6.4430e+03b=1.01490.91140.97611.0976F=1.01490.91140.97611.0976I=1.0e+03*7.69836.91327.40388.32607.3.2季節(jié)性趨勢模型

當(dāng)時間序列既有季節(jié)性變動又有趨勢性變動時，先建立趨勢預(yù)測模型，在此基礎(chǔ)上求得季節(jié)指數(shù)，再建立預(yù)測模型。其過程如下：（1）計算歷年同季平均數(shù)r；（2）建立趨勢預(yù)測模型，求趨勢值；（3）計算出趨勢值后，再計算出歷年同季的平均數(shù)R；（4）計算趨勢季節(jié)指數(shù)（k）；用同季平均數(shù)r與趨勢值同季平均數(shù)R之比來計算；（5）對趨勢季節(jié)指數(shù)進(jìn)行修正；（6）求預(yù)測值。將預(yù)測期的趨勢值乘以該期的趨勢季節(jié)指數(shù)，即預(yù)測模型?！纠?-8】（續(xù)【例7-7】），試用季節(jié)性趨勢模型預(yù)測2022年1至4季度居民人均消費性支出。（1）根據(jù)所給數(shù)據(jù)，畫出走勢圖，觀察季節(jié)性和趨勢性。X=[55344867523557565970521556126282638755445915659967495996...6269 7098716064056814768464786007676277607495 707174158326];t=1:length(X);plot(t,X,'-o')%顯示結(jié)果見如圖7-7xlabel('時間');ylabel('銷售量')（2）計算各年同季平均數(shù)。r1=mean(X(1:4:length(X)));r2=mean(X(2:4:length(X)));r3=mean(X(3:4:length(X)));r4=mean(X(4:4:length(X)));r=[r1r2r3r4]%各年同季平均（3）計算趨勢預(yù)測值。p=polyfit(t,X,1);%擬合得長期趨勢參數(shù)T=polyval(p,t);%計算長期趨勢預(yù)測值（4）計算趨勢值各年同季平均。R1=mean(T(1:4:length(T)));R2=mean(T(2:4:length(T)));R3=mean(T(3:4:length(T)));R4=mean(T(4:4:length(T)));R=[R1R2R3R4]%趨勢值各年同季平均（5）計算并調(diào)整趨勢季節(jié)指數(shù)。k=r./R%趨勢季節(jié)指數(shù)K=4/sum(k)*k%調(diào)整趨勢季節(jié)指數(shù)（6）預(yù)測2022年四個季度居民人均消費性支出。t1=length(X)+1:length(X)+4;%2022年1至4季度時間T1=polyval(p,t1)%計算2022年趨勢預(yù)測值Y=K.*T1%計算2022年最終預(yù)測值plot(t,X,'-o',t1,Y,'-*',[t,t1],[T,T1],'-')xlabel('時間/季度')ylabel('人均消費支出/元')7.3.3季節(jié)性環(huán)比法模型環(huán)比法是指積累歷年（至少三年）各月或各季的歷史資料，逐期計算環(huán)比，加以平均，求出季節(jié)指數(shù)季節(jié)預(yù)測的方法。（1）求逐期環(huán)比：將本期實際值和前期實際值相比，即：第一期的環(huán)比不能計算（2）計算同季環(huán)比平均數(shù)

（3）計算各季連鎖指數(shù)以第一季度為固定基準(zhǔn)期，其連鎖指數(shù)為，后面各季平均環(huán)比逐期連乘，得各季連鎖指數(shù)：（4）根據(jù)趨勢變動修正連鎖指數(shù)如果沒有趨勢變動，基準(zhǔn)期的連鎖指數(shù)應(yīng)為1，若求出來的基準(zhǔn)期（第一季度）的連鎖指數(shù)不為1，則存在趨勢變動的影響，應(yīng)加以修正，其修正值為此時是第四季度的連鎖指數(shù)乘以第一季度的平均環(huán)比，即各季扣除d后的修正連鎖指數(shù)應(yīng)為：第一季度第二季度第三季度第四季度：（5）計算季節(jié)指數(shù)將各季修正連鎖指數(shù)，除以全部四個季度修正連鎖指數(shù)的平均數(shù)，得各季季節(jié)指數(shù)：（6）配合趨勢直線模型，計算趨勢值結(jié)合季節(jié)指數(shù)進(jìn)行預(yù)測，預(yù)測模型為：【例7-9】（續(xù)【例7-7】）試用季節(jié)性環(huán)比法預(yù)測2022年1至4季度居民人均消費性支出。（1）先畫出走勢圖X=[5534 486752355756597052155612628263875544591565996749 5996...6269 7098716064056814768464786007676277607495 707174158326];t=1:length(X);plot(t,X,'-o')%顯示結(jié)果如圖（2）計算季節(jié)指數(shù)h=x(2:end)./x(1:end-1)%各期環(huán)比h1=mean(h(4:4:end))%第1季度同季環(huán)比平均數(shù)h2=mean(h(1:4:end))%第2季度同季環(huán)比平均數(shù)h3=mean(h(2:4:end))%第3季度同季環(huán)比平均數(shù)h4=mean(h(3:4:end))%第4季度同季環(huán)比平均數(shù)H1=[h1h2h3h4]%四個季度同季環(huán)比平均數(shù)H2=[1h2h3h4]%第1季度基準(zhǔn)期為1四個季度

%同季環(huán)比平均數(shù)c=cumprod(H2)%四個季度連鎖指數(shù)c1=c(4)*H1(1)%第1季度連鎖指數(shù)d=(c1-1)/4%修正值C1=1%第1季度修正連鎖指數(shù)C2=c(2)-d%第2季度修正連鎖指數(shù)C3=c(3)-2*d%第3季度修正連鎖指數(shù)C4=c(4)-3*d%第4季度修正連鎖指數(shù)C=[C1C2C3C4]%匯總修正連鎖指數(shù)F=C./mean(C)%季節(jié)指數(shù)環(huán)比法季節(jié)指數(shù)季度年份12342015環(huán)比

0.87951.07561.09952016環(huán)比1.03720.87351.07611.11942017環(huán)比1.01670.86801.06691.11562018環(huán)比1.02270.88841.04551.13222019環(huán)比1.00870.89461.06391.12772020環(huán)比0.84310.92731.12571.14762021環(huán)比0.96590.94341.04861.1229同季環(huán)比平均數(shù)H10.98240.89641.07181.1236連鎖指數(shù)c1.00001.06040.89840.96071.0794修正值d0.0151修正連鎖指數(shù)C10.88130.93051.0341季節(jié)指數(shù)F1.04010.91660.96781.0756（3）求趨勢值p=polyfit(t,x,1)T=polyval(p,t)（4）求2016年1-4季度預(yù)測值1=length(X)+1:length(X)+4;T1=polyval(p,t1)Y=F.*T1%2022年預(yù)測值plot(t,X,'-o',t1,Y,'-+',[t,t1],[T,T1],'-k')7.4時間序列分解法（1）時間序列數(shù)據(jù)的影響因素主要有長期趨勢、季節(jié)變動、周期變動、不規(guī)則變動。（2）乘法分解模型為時間序列的全變動，為長期趨勢，為季節(jié)變動，

為循環(huán)變動，為不規(guī)則變動。（3）確定上述各個因素的步驟1）用MA=TC分析長期趨勢與循環(huán)變動；2）用X/MA=SI分析季節(jié)性與隨機(jī)性；3）用SI平均值分析季節(jié)性；4）用趨勢外推法分析長期趨勢T；5）用MA/T分析循環(huán)變動；6）將時間序列的T、S、C分解出來后，剩余的即為不規(guī)則變動，即在實際運(yùn)算時可以不考慮隨機(jī)因素，而直接用前三種因素來處理：即【例7-10】（續(xù)【例7-7】）試用序列分解模型預(yù)測2022年1至4季度居民人均消費性支出。（1）計算一次移動平均值MA，即得長期趨勢與循環(huán)變動。X=[5534 486752355756597052155612628263875544591565996749 5996...6269 7098716064056814768464786007676277607495 707174158326];fort=1:length(X)-3MA1(t)=(X(t)+X(t+1)+X(t+2)+X(t+3))/4;end%四項移動平均值所得出的結(jié)果應(yīng)放在第2.5季度的位置，若將其放在第三季度的位置，%需將一次移動平均值MA1的前后兩項再取平均數(shù)：MA=(MA1(1:end-1)+MA1(2:end))/2（2）季節(jié)與隨機(jī)性SI，即X與MA的比率。SI=100*X(3:end-2)./MA（3）用各年同季平均，去掉SI中的隨機(jī)性，得季節(jié)指數(shù)r，并修正季節(jié)指數(shù)R。r1=mean(SI(3:4:end));r2=mean(SI(4:4:end));r3=mean(SI(1:4:end));r4=mean(SI(2:4:end));r=[r1r2r3r4];R=r./mean(r)（4）擬合得長期趨勢T。m=1:length(X);p=polyfit(m,X,1)T=polyval(p,m);（5）計算循環(huán)變動CC=MA./T(3:end-2);（6）預(yù)測2022年1-4季度。n=length(X)+1:length(X)+4%2022年1至4季度時間T1=polyval(p,n)%2022年1至4季度趨勢項C1=mean(C)%取循環(huán)變動C的均值X1=T1.*C1.*R%2022年1至4季度預(yù)測值（3）用各年同季平均，去掉SI中的隨機(jī)性，得季節(jié)指數(shù)r，并修正季節(jié)指數(shù)R。r1=mean(SI(3:4:end));r2=mean(SI(4:4:end));r3=mean(SI(1:4:end));r4=mean(SI(2:4:end));r=[r1r2r3r4];R=r./mean(r)（4）擬合得長期趨勢T。m=1:length(X);p=polyfit(m,X,1)T=polyval(p,m);（5）計算循環(huán)變動CC=MA./T(3:end-2);（6）預(yù)測2022年1-4季度。n=length(X)+1:length(X)+4%2022年1至4季度時間T1=polyval(p,n)%2022年1至4季度趨勢項C1=mean(C)%取循環(huán)變動C的均值X1=T1.*C1.*R%2022年1至4季度預(yù)測值序號觀察值X移動平均值MA=TC比

率X/MA=SI長期趨勢T循環(huán)變動C15534

5261.27

24867

5348.81

352355402.5096.905436.350.9938457565500.50104.655523.880.9958559705591.13106.785611.420.9964652155704.0091.435698.961.0009756125821.8896.405786.501.0061862825915.13106.205874.041.0070963875994.13106.555961.571.00551055446071.6391.316049.111.00371159156156.5096.086136.651.00321265996258.25105.446224.191.00551367496359.00106.136311.731.00751459966465.6392.746399.271.01041562696579.3895.286486.801.01431670986681.88106.236574.341.01641771606801.13105.286661.881.02091864056942.5092.266749.421.02861968146930.5098.326836.961.01372076846795.50113.076924.500.98142164786739.2596.127012.030.96112260076742.2589.097099.570.94972367626878.8898.307187.110.95712477607139.00108.707274.650.98142574957353.63101.927362.190.99882670717506.0094.207449.731.0076277415

7537.26

288326

7624.80

7.5ARMA模型預(yù)測法7.5.1ARMA模型的基本形式（1）自回歸模型AR(p)其中，是獨立同分布的隨機(jī)變量序列稱時間序列服從p階自回歸模型AR(p)（2）移動平均模型MA(q)服從q階移動平均模型MA(q)（3）自回歸移動平均模型ARMA(p,q)服從（p,q）階自回歸移動平均模型ARMA(p,q)2．模型建立的條件及判定法時間序列的平穩(wěn)性自相關(guān)分析法它可以測定時間序列的隨機(jī)性和平穩(wěn)性，以及時間序列的季節(jié)性。7.5.2ARMA模型相關(guān)性分析及識別根據(jù)繪制的自相關(guān)分析圖和偏自相關(guān)分析圖，我們可以初步地識別平穩(wěn)序列的模型類型和模型階數(shù)。1．AR(p)模型（1）AR(p)的自相關(guān)函數(shù)滿足表明隨k的增加按指數(shù)形式衰減，呈“拖尾”狀。AR(1)模型AR(2)模型（2）AR(p)的偏相關(guān)函數(shù)可知偏相關(guān)函數(shù)具有“截尾”狀。2．MA(q)模型（1）MA(q)自相關(guān)函數(shù)（2）MA(q)偏相關(guān)函數(shù)由于任何一個可逆的MA過程都可以轉(zhuǎn)化為一個無限階的系數(shù)按幾何遞減的AR過程，所以MA過程的偏自相關(guān)函數(shù)同AR模型一樣呈緩慢衰減特征。3．ARMA(p,q)模型根據(jù)AR、MA模型可知ARMA模型的自相關(guān)函數(shù)和偏自相關(guān)函數(shù)也是無限延長的，其過程也是呈緩慢衰減，是拖尾的。三個基本模型的相關(guān)性特征根據(jù)相關(guān)性特征，可利用自相關(guān)函數(shù)與偏自相關(guān)函數(shù)的截尾性來識別模型類型。并利用偏相關(guān)函數(shù)PartialACF，確定AR模型的滯后階數(shù)；利用自相關(guān)函數(shù)ACF，確定MA模型的滯后階數(shù)。模型自相關(guān)函數(shù)偏自相關(guān)函數(shù)AR(p)拖尾p階截尾MA(q)q階截尾拖尾ARMA(p,q)拖尾拖尾4．自相關(guān)函數(shù)與偏相關(guān)函數(shù)的命令（2）計算并描繪時間序列的自相關(guān)函數(shù)格式：

autocorr(series,nLags,M,nSTDs)%繪出自相關(guān)函數(shù)圖

[ACF,Lags,Bounds]=autocorr(series,nLags,M,nSTDs)說明：series：時間序列

nLags：延遲數(shù)，默認(rèn)為20個ACF。

M：延遲階數(shù)，缺省時假設(shè)為高斯白噪聲。

nSTDs：表示計算出的相關(guān)函數(shù)ACF估計誤差的標(biāo)準(zhǔn)差；ACF：相關(guān)函數(shù)；Lags：對應(yīng)于ACF的延遲；Bounds：置信區(qū)間的近似上下限，假設(shè)序列是MA（M）模型（1）計算時間序列的相關(guān)系數(shù)格式：r=corrcoef(x1,x2)說明：計算兩時間序列x1,x2的相關(guān)系數(shù)r，其值在[0，1]之間?！纠?-11】x=randn(1000,1);%生成1000點的Gaussian白噪聲y=filter([1-11],1,x);%生成MA（2）過程autocorr(y,[],2)%繪出相關(guān)函數(shù)圖[ACF,Lags,Bounds]=autocorr(y,[],2)%計算95％置信度下的相關(guān)系數(shù)（3）計算并描繪時間序列的偏相關(guān)函數(shù)格式：parcorrr(series)[PACF,Lags,Bounds]=parcorr(series,nLags,R,nSTDs)說明：series：時間序列；

nLags：延遲數(shù)，缺省時計算在延遲點0，1，…，T

（T=min([20,length(series)-1])的PACF；

R：表示Lags延遲階數(shù)，缺省時假設(shè)為AR(R)過程；

nSTDs：表示計算出的相關(guān)函數(shù)PACF估計誤差的標(biāo)準(zhǔn)差PACF：相關(guān)函數(shù)；Lags：對應(yīng)于ACF的延遲；

Bounds：置信區(qū)間的近似上下限，假設(shè)序列是AR(R)過程?！纠?-12】x=randn(1000,1);%Gaussian白噪聲

y=filter(1,[1-0.60.08],x);%生成AR(2）過程

parcorr(y,[],2)%繪出偏相關(guān)函數(shù)圖，

[PACF,Lags,Bounds]=parcorr(y,[],2)%偏相關(guān)系數(shù)5．評價時間序列模型的準(zhǔn)則FPE準(zhǔn)則：是指最終預(yù)報誤差（FinalPredictionError）的定階準(zhǔn)則。主要用于AR模型、ARMA模型的階，其方法是以選用模型的一步誤差達(dá)到最小的相應(yīng)的階作為模型的階，用其預(yù)報效果的優(yōu)劣來確定該模型的階數(shù)。

7.5.3ARMA模型參數(shù)估計1、AR(p)模型參數(shù)矩估計Yule-Walker方程利用實際時間序列數(shù)據(jù)，首先求得自相關(guān)函數(shù)的估計值，代入Yule-Walker方程組，求得模型參數(shù)的估計值。2、MA(q)模型參數(shù)估計利用實際時間序列數(shù)據(jù)，求得自協(xié)方差函數(shù)的估計值求得模型參數(shù)的估計值。

3、ARMA(p,q)模型的參數(shù)估計先求得自相關(guān)函數(shù)的估計值，代入Yule-Walker方程組，求得模型參數(shù)的估計值，再改寫ARMA模型求解估計值

4、模型參數(shù)的MATLAB命令（1）AR模型參數(shù)估計格式：m=ar(y,n)[m,refl]=ar(y,n,approach,window)說明：y是數(shù)據(jù)結(jié)構(gòu)，由iddata函數(shù)得到：y=iddata(y)，后面y是給定的時間序列；

n是AR階次；

approach：估計時采用的方法：Approach＝’fb’：前向后；’ls’：最小二乘法；’yw’：Yule-Walker方法；’Burg’：基于Burg譜估計方法；

Window：處理Y中缺失值的方法，Window＝’now’：表示觀察值中沒有缺失值；Window＝’yw’：表示Yule-Walker方法處理缺失值；

m：AR模型的文字形式；

refl：AR模型的系數(shù)。（2）ARMAX模型參數(shù)估計自回歸移動平均各態(tài)歷經(jīng)ARMAX（AutoRegressiveMovingAverageeXogenous）模型，是考慮外部解釋變量X的模型。na,nb,nc是滯后多項式的階數(shù)，nk為延遲格式：Z=iddata(y)m=armax(Z,[nanbncnk])m=armax(Z,'na',na,'nb',nb,'nc',nc,'nk',nk)說明：y原始序列，Z是y的數(shù)據(jù)結(jié)構(gòu)；

na,nb,nc是滯后多項式的階數(shù)，nk為延遲（3）MA模型參數(shù)估計。用ARMAX模型可對MA模型進(jìn)行估計，只需在模型

A(q)＝1，B(q)＝0

格式：z=iddata(y)m=armax(z,’nc’,5)（4）ARMA模型參數(shù)估計用ARMAX模型可對ARMA模型進(jìn)行估計，只需在模型：B(q)＝0

格式：z=iddata(y);m=armax(z,[nanc]);（5）ARX模型參數(shù)估計

A(q)y(t)=B(q)u(t-nk)+e(t)格式：m=arx(data,[nanbnk])m=arx(data,'na',na,'nb',nb,'nk',nk)7.5.4ARMA模型的預(yù)測1．AR（p）模型的預(yù)測公式預(yù)測方差

GREEN函數(shù)2．MA(q)模型預(yù)測公式預(yù)測方差

若已知和新獲得的數(shù)據(jù)，則得的遞推公式：[,,…，]T

，初始值可取某個時刻3．ARMA(p,q)模型預(yù)測公式預(yù)測方差

的遞推公式：[,,…，]T

，當(dāng)時，上式最后一項為04．模型預(yù)測及誤差的MATLAB命令格式：yp=predict(m,y,k)說明：m表預(yù)測模型，y為實際輸出，k為預(yù)測區(qū)間；yp為預(yù)測輸出。當(dāng)k<inf,yp(t)為模型m與y（1,2，……t-k）的預(yù)測值；當(dāng)k=inf,yp(t)為模型m的純仿真值，默認(rèn)k=1；在計算AR模型預(yù)測時，k應(yīng)取1。格式：[yh,fit,x0]=compare(m,y,k)說明：Compare的預(yù)測原理與predict相同，但對預(yù)測進(jìn)行比較，并可繪出比較圖

格式：e=pe(m,data)%pe誤差計算，說明：采用yh=predict(m,data,1)進(jìn)行預(yù)測，然后計算誤差e=data-yh在無輸出情況下，繪出誤差圖，誤差曲線應(yīng)足夠小，黃色區(qū)域為99%的置信區(qū)間，誤差曲線在該區(qū)域內(nèi)表明通過檢驗。格式：[e,r]=resid(m,data,mode,lags)resid(r)計算并檢驗誤差。7.6案例分析7.6.1利用移動平均法預(yù)測GDP【例7-13】我國2008--2021年的國內(nèi)生產(chǎn)總值（GDP）數(shù)據(jù)如表，試用指數(shù)平滑法預(yù)測2022年的國內(nèi)生產(chǎn)總值。年份2008200920102011201220132014GDP31.9234.8541.2148.7953.8659.364.36年份2015201620172018201920202021GDP68.8974.6483.291.9398.65101.36113.35第一步，先編程查找二次指數(shù)平滑系數(shù)alpha.clearX=[31.9234.8541.2148.79 53.8659.364.36...68.89 74.64 83.291.9398.65101.36113.35];X0=X(1);X1=X(2:end);U=[];foralpha=0.1:0.1:0.9%在0.1至0.9之間查找alphaS0=X0;S1(1)=alpha*X1(1)+(1-alpha)*S0;fort=1:length(X1)-1S1(t+1)=alpha*X1(t+1)+(1-alpha)*S1(t);endS20=X0;S2(1)=alpha*S1(1)+(1-alpha)*S20;fort=1:length(S1)-1S2(t+1)=alpha*S1(t+1)+(1-alpha)*S2(t);end

a=2*S1(1:end)-S2(1:end);b=alpha/(1-alpha)*(S1(1:end)-S2(1:end));T=1;y=a+b*T;Y=[X0,y];MSE=sum((X(2:end)-Y(1:end-1)).^2)./(length(X)-1);U=[U,MSE];endR=find(U==min(U));%獲得均方誤差最小的最佳位置alpha=0.1:0.1:0.9;ALPHA=alpha(R)%獲得二次指數(shù)平滑系數(shù)alphaALPHA=0.8000第二步，對獲取的二次指數(shù)平滑系數(shù)alpha=0.8進(jìn)行預(yù)測。只需上面程序添加alpha=0.8，刪去U=[];foralpha=0.1:0.1:0.9……MSE=sum((X(2:end)-Y(1:end-1)).^2)./(length(X)-1);U=[U,MSE];endR=find(U==min(U)alpha=0.1:0.1:0.9;ALPHA=alpha(R)再添加繪圖t1=1:length(X);t2=2:length(X)+1;plot(t1,X,'-+',t2,Y,'-O')xlabel('時間/年')ylabel('GDP/萬億元')legend('原始數(shù)據(jù)','預(yù)測值')7.6.2利用ARMA模型預(yù)測股票價格【例7-14】平安銀行2021年10月11日至12月31日的交易日的收盤價數(shù)據(jù)如表7-7，試用ARMA模型預(yù)測未來五個交易日的股票價格。日期收盤價日期收盤價日期收盤價日期收盤價10月11日1411月1日19.3911月22日18.1212月13日18.2710月12日19.3511月2日18.1811月23日17.8812月14日17.5810月13日19.5811月3日18.0311月24日17.8712月15日17.5510月14日19.2111月4日17.8711月25日17.6812月16日17.7210月15日19.6611月5日17.6411月26日17.5812月17日17.5710月18日19.2911月8日17.4211月29日17.5112月20日17.5210月19日19.5711月9日17.5311月30日17.4412月21日17.5910月20日19.2411月10日17.4012月1日17.6412月22日17.3910月21日20.0011月11日18.3512月2日17.5912月23日17.3210月22日20.0411月12日18.2712月3日17.6512月24日17.3110月25日20.1211月15日18.4312月6日18.1012月27日17.2210月26日20.0511月16日18.2212月7日18.4612月28日17.1710月27日19.9211月17日18.1112月8日18.4712月29日16.7510月28日19.7211月18日17.8012月9日18.7012月30日16.8210月29日19.5011月19日18.1512月10日18.3812月31日16.48（1）輸入數(shù)據(jù)，畫圖觀察平穩(wěn)性clearX=[…];%將表格給定的60個收盤價數(shù)據(jù)用列向量X表示t=1:length(X);plot(t,X)%畫出時序圖xlabel('時間/天')ylabel('收盤價/元')運(yùn)行結(jié)果顯示如圖（2）從圖7-14可知，非平穩(wěn)，所以需要做差分處理化成平穩(wěn)序列Y=diff(X)%差分figure(2)plot(1:length(Y),Y)xlabel('時間/天')ylabel('差分值')（3）利用自相關(guān)函數(shù)圖與偏相關(guān)函數(shù)圖判斷模型的類型與階次figure(3)subplot(2,1,1)autocorr(Y)subplot(2,1,2)parcorr(Y)（4）從圖7-16可知，可確定為ARMA(2,2)模型Z=iddata(Y)m=armax(Z,[2,2])（5）利用模型進(jìn)行預(yù)測L=5%預(yù)測長度，五個交易日y=[Y;zeros(L,1)];%加初始窗p=iddata(y);P=predict(m,p,L)%預(yù)測但沒有顯示G=get(P)%獲取預(yù)測值的全部結(jié)構(gòu)內(nèi)容PT=G.OutputData{1,1}(length(Y)+1:length(Y)+L,1)%顯示未來長度L=5預(yù)測值D=[Y;PT]%全部差分值X1=cumsum([X(1);D])%全部差分值的還原值X2=X1(length(X)+1:end)%未來長度L=5預(yù)測值的還原值figure(4)compare(m,Z,L)%顯示預(yù)測值與原始數(shù)據(jù)比較圖X2=16.484216.481216.506316.471416.5064故2021年12月31日之后的五個交易日，平安銀行股票收盤價預(yù)測值分別為16.48元、16.48元、16.50元、16.47元和16.51元。在上述求解過程中，第（1）步和第（2）步在討論是否為平穩(wěn)序列時，采用的是觀察法。也可以使用平穩(wěn)性單位根（ADF）檢驗的函數(shù)命令：h1=adftest(X)%X原始序列h2=adftest(Y)%一次差分后的序列h1=0,X為非平穩(wěn)序列h2=1,Y為平穩(wěn)序列程序中的第（3）步確定模型類型和階次時，使用的是觀察自相關(guān)函數(shù)和偏相關(guān)函數(shù)圖，也可以使用AIC準(zhǔn)則來編程：Z=iddata(Y);%將X轉(zhuǎn)化為matlab接受的格式test=[];forp=1:6forq=1:6m=armax(Z,[p,q]);AIC=aic(m);%armax(p,q),選擇對應(yīng)FPE最小，AIC值最小的模型

test=[test;pqAIC];endendfork=1:size(test,1)iftest(k,3)==min(test(:,3))%選擇AIC值最小的模型

p_test=test(k,1);q_test=test(k,2);endendp_test,q_testp_test=4q_test=6ARMA(4,6)模型因此需將第（4）步求解模型的程序改為：m=armax(Z,[p_testq_test])經(jīng)濟(jì)預(yù)測與決策技術(shù)及MATLAB實現(xiàn)第8章干預(yù)分析模型預(yù)測法8.1干預(yù)分析模型練習(xí)與提高（八）

8.2案例分析

8.1干預(yù)分析模型基本形式（1）持續(xù)性干預(yù)變量。是指當(dāng)某一突發(fā)事件在T時刻發(fā)生以后，一直對所研究的問題有影響。這時我們可用階躍函數(shù)表示，其形式為：首頁（2）短暫性干預(yù)變量。指當(dāng)某一突發(fā)事件在T時刻發(fā)生，僅對該時刻有影響。這時我們可用脈沖函數(shù)表示，其形式為：

8.1.1干預(yù)分析模型基本變量8.1.2干預(yù)事件的形式對經(jīng)濟(jì)活動可能產(chǎn)生的經(jīng)濟(jì)效果和經(jīng)濟(jì)發(fā)展趨勢作出科學(xué)預(yù)見的學(xué)科。（1）干預(yù)事件的運(yùn)行突然開始，長期持續(xù)下去。

設(shè)干預(yù)對因變量的影響是固定的，從某一時刻T開始，但影響程度是未知的，即因變量的大小是未知的。干預(yù)模型為：w表示干預(yù)影響程度的未知參數(shù)干預(yù)事件發(fā)生后需延遲b個時期才能產(chǎn)生影響其最簡單模型為：（2）干預(yù)事件的影響逐漸開始，長期持續(xù)下去。（3）干預(yù)事件突然開始，產(chǎn)生暫時的影響。（4）干預(yù)事件逐漸開始，產(chǎn)生暫時的影響。8.1.3干預(yù)分析模型的基本形式設(shè)y(t)為經(jīng)濟(jì)過程的某個指標(biāo)值,x(t)表示無干預(yù)影響時的經(jīng)濟(jì)過程值，z(t)為干預(yù)影響過程值，則在干預(yù)影響下的經(jīng)濟(jì)過程值可表示為：其中，x(t)可用回歸法、季節(jié)指數(shù)以及ARMA等方法確定，z(t)由上述干預(yù)事件確定。首頁以ARMA模型說明建立干預(yù)模型的步驟：（1）利用無干預(yù)影響的序列（即的前半部分）建立ARMA模型：然后利用此模型進(jìn)行外推預(yù)測，得到的預(yù)測值作為不受干預(yù)影響的數(shù)值（2）將干預(yù)事件后的實際值（即y(t)的后半部分）減去預(yù)測值，得到受干預(yù)影響的效用值：利用這些結(jié)果估計干預(yù)影響部分首頁中的參數(shù)和若干預(yù)影響部分的形式為：則就是估計自回歸方程的參數(shù)。（3）計算凈化序列。凈化序列是指消除了干預(yù)影響的序列，它由實際的觀察值減去干預(yù)影響值得到，即：當(dāng)時，當(dāng)時，（4）利用凈化序列建立一個改進(jìn)的ARMA模型：(5)將第（4）步與第（2）步的結(jié)合，就得總的干預(yù)分析模型：（6）利用模型進(jìn)行預(yù)測當(dāng)時，當(dāng)時，2案例分析

8.2.1我國進(jìn)出口商品總額模型首頁選取2007年1月至2009年2月的海關(guān)進(jìn)出口總額數(shù)據(jù)，并以2008年9月金融危機(jī)爆發(fā)作為干預(yù)事件，因此將時間序列數(shù)據(jù)分為兩個時期：第一時期為2007年1月至2008年9月，第二時期為2008年10月至2009年2月。由于金融危機(jī)爆發(fā)并不是立刻產(chǎn)生完全的影響，而是隨著時間的推移，逐漸地感到這種影響的存在。因此干預(yù)影響可選為模型：海關(guān)進(jìn)出口總額首頁序號日期進(jìn)出口額調(diào)整(30日)額序號日期進(jìn)出口額調(diào)整(30日)額12007年1月1573.61522.8142008年2月16611779.622007年2月1404.51504.8152008年3月2044.51978.532007年3月1599.31547.7162008年4月2207.632207.642007年4月1779.91779.9172008年5月2209.42138.152007年5月1656.11602.7182008年6月2220.882220.962007年6月1795.91795.9192008年7月2479.412399.472007年7月1910.71849.1202008年8月2412.652334.882007年8月1977.131913.4212008年9月2437.532437.592007年9月2009.172009.2222008年10月2215.392143.9102007年10月1881.91821.2232008年11月1897.031897112007年11月2088.72088.7242008年12月1832.61773.5122007年12月2061.441994.9252009年1月14181372.3132008年1月1998.31933.8262009年2月1249.71339首頁（1）由于是月數(shù)據(jù)，考慮每月的天數(shù)不一致，所以都化為標(biāo)準(zhǔn)的30天，調(diào)整后的數(shù)據(jù)見表，并用y表示這26個數(shù)據(jù)組成的列向量。cleary=[1522.81504.81547.71779.91602.71795.91849.11913.4...2009.21821.22088.71994.91933.81779.61978.52207.62138.1...2220.92399.42334.82437.52143.918971773.51372.31339]';（2）根據(jù)2007年1月到2008年9月，即前21個數(shù)據(jù)，建立時間序列模型x1=y(1:21)t1=1:21;T1=[ones(21,1),t1'];formatlong[b1,bint1,r1,rint1,stats1]=regress(x1,T1)%建立回歸模型首頁（3）分離出干預(yù)影響的具體數(shù)據(jù)，估計干預(yù)模型的參數(shù)利用所得模型，對2008年10月以后的數(shù)據(jù)進(jìn)行預(yù)測，然后用實際值減去預(yù)測值，得到的差值就是金融危機(jī)產(chǎn)生的效應(yīng)值z。t2=22:26;x2=b1(1)+b1(2).*t2';z=y(22:end)-x2(1:end)利用z的數(shù)據(jù)，可估計出干預(yù)模型：u=z(2:end)v=z(1:end-1)T2=[ones(length(v),1),v];[b2,bint2,r2,rint2,stats2]=regress(u,T2)首頁（4）計算凈化序列。T＝21，t＝1,2,…,26X(1:21)=y(1:21);fort=22:26X(t)=y(t)+b2(2)*(X(t-1)-y(t-1))-b2(1);endX首頁（5）對凈化序列建立擬合模型t3=1:26;T3=[ones(26,1),t3'];[b3,bint3,r3,rint3,stats3]=regress(X',T3)（6）組建干預(yù)分析模型將（3）、（5）步驟組合，即得干預(yù)分析模型：首頁（7）運(yùn)用模型進(jìn)行預(yù)測fort=1:21X0(t)=b3(1)+b3(2)*t;y0=X0;%2008年9月及以前的預(yù)測值endfort=22:26X0(t)=b3(1)+b3(2)*t;y0(t)=X0(t)+b2(2)*(y0(t-1)-X0(t-1))+b2(1);%2008年10月以后的預(yù)測值

endformaty0首頁（8）畫圖比較plot(t3,y,'+',t3,y0,'o')經(jīng)濟(jì)預(yù)測與決策技術(shù)及MATLAB實現(xiàn)

第9章馬爾可夫預(yù)測方法9.1馬爾可夫鏈基本理論9.2.1市場占有率預(yù)測9.2案例分析9.2.2股票價格走勢預(yù)測9.2.3加權(quán)馬氏鏈法預(yù)測證券指數(shù)走勢9.2.4期望利潤預(yù)測

練習(xí)與提高（九）9.1.1馬爾可夫鏈基本概念（1）馬爾可夫鏈?zhǔn)醉?.1馬爾可夫鏈基本理論注:而與以前的狀態(tài)(2)一步轉(zhuǎn)移概率馬氏鏈在時刻n處于狀態(tài)i

的條件下，到時刻n+1轉(zhuǎn)移到狀態(tài)j

的條件概率，即稱為在時刻n的一步轉(zhuǎn)移概率，首頁注:由于概率是非負(fù)的，且過程從一狀態(tài)出發(fā)，經(jīng)過一步轉(zhuǎn)移后，必到達(dá)狀態(tài)空間中的某個狀態(tài)一步轉(zhuǎn)移概率滿足(3)一步轉(zhuǎn)移矩陣稱為在時刻n的一步轉(zhuǎn)移矩陣首頁即有有限馬氏鏈狀態(tài)空間I={0，1，2，…，k}首頁(4)齊次馬氏鏈即則稱此馬氏鏈為齊次馬氏鏈（即關(guān)于時間為齊次）(5)初始分布首頁注馬氏鏈在初始時刻有可能處于I中任意狀態(tài)，初始分布就是馬氏鏈在初始時刻的概率分布。(6)絕對分布概率分布稱為馬氏鏈的絕對分布或稱絕對概率定態(tài)分布即首頁在馬氏鏈的研究中，須研究“從已知狀態(tài)i出發(fā)，經(jīng)過n次轉(zhuǎn)移后，系統(tǒng)將處于狀態(tài)j”的概率.(7)n步轉(zhuǎn)移矩陣1)n步轉(zhuǎn)移概率系統(tǒng)在時刻m從狀態(tài)i經(jīng)過n步轉(zhuǎn)移后處于狀態(tài)j的概率稱為n步轉(zhuǎn)移概率由于馬氏鏈?zhǔn)驱R次的，這個概率與m無關(guān)首頁顯然有2)n步轉(zhuǎn)移矩陣稱為n步轉(zhuǎn)移矩陣規(guī)定首頁注（1）用一步轉(zhuǎn)移概率表示多步轉(zhuǎn)移概率首頁注I={1，2，…，N}由矩陣的乘法規(guī)則，得表示：在時刻n，各狀態(tài)的概率等于其初始狀態(tài)的概率與n步轉(zhuǎn)移概率矩陣之積。若鏈?zhǔn)驱R次的，則有首頁(8)遍歷性定義1使得則稱此馬氏鏈具有遍歷性馬氏鏈的遍歷性表明不論從哪一個狀態(tài)i出發(fā)，當(dāng)轉(zhuǎn)移的步數(shù)n充分大時，轉(zhuǎn)移到狀態(tài)j的概率都接近于正常數(shù)首頁(9)平穩(wěn)分布有給轉(zhuǎn)移概率P，若存在一個概率分布中的j=0，1，2，…，s則稱為平穩(wěn)分布。9.1.2馬爾可夫鏈預(yù)測原理1．馬氏鏈近期預(yù)測原理定理1設(shè){Xn}為一個齊次馬氏鏈，其狀態(tài)空間為I，絕對概率為

n步轉(zhuǎn)移概率為，則有：即表明馬氏鏈的絕對概率由其初始分布和n步轉(zhuǎn)移概率完全確定定理2C-K方程（Chapman-Kolmogorov方程）設(shè){Xn}為一個齊次馬氏鏈，其狀態(tài)空間為I，

n步轉(zhuǎn)移概率為，則有：

注：（1）C-K方程的矩陣形式為：

（2）n步轉(zhuǎn)移矩陣與一步轉(zhuǎn)移矩陣之間的關(guān)系

（3）定理1與定理2的結(jié)合即為馬爾可夫鏈預(yù)測模型：

2．馬氏鏈穩(wěn)態(tài)概率分布預(yù)測原理則此馬氏鏈?zhǔn)潜闅v的，且中的是方程組j=0，1，2，…，s的滿足條件的唯一解注1定理表明不論從鏈中哪一狀態(tài)i出發(fā)，都能以正概率經(jīng)有限次轉(zhuǎn)移到達(dá)鏈中預(yù)先指定的其它任一狀態(tài)。定理給出了求平穩(wěn)分布的方法。注23.馬爾可夫鏈預(yù)測基本步驟：（1）劃分狀態(tài)區(qū)間，確定狀態(tài)空間I=[1,2，…，m]；（2）按步驟（1）所劃分狀態(tài)區(qū)間，確定資料序列中各時段指標(biāo)值所對應(yīng)的狀態(tài)；（3）對步驟（2）所得的結(jié)果進(jìn)行統(tǒng)計計算，得馬爾可夫鏈的一步轉(zhuǎn)移頻數(shù)矩陣N及概率矩陣P?！榜R氏性”檢驗方法:（4）進(jìn)