-280- 第二十四章 時間序列模型 時間序列是按時間順序排列的、隨時間變化且相互關(guān)聯(lián)的數(shù)據(jù)序列。分析時間序列的方法構(gòu)成數(shù)據(jù)分析的一個重要領(lǐng)域，即時間序列分析。 時間序列根據(jù)所研究的依據(jù)不同，可有不同的分類。 1按所研究的對象的多少分，有一元時間序列和多元時間序列。 2按時間的連續(xù)性可將時間序列分為離散時間序列和連續(xù)時間序列兩種。 3按序列的統(tǒng)計特性分，有平穩(wěn)時間序列和非平穩(wěn)時間序列。如果一個時間序列的概率分布與時間 t 無關(guān)，則稱該序列為嚴格的（狹義的）平穩(wěn)時間序列。如果序列的一、二階矩存在，而且對任意時刻 t 滿足： （ 1）均值為常數(shù) （ 2）協(xié)方差為時間間隔 的函數(shù)。 則稱該序列為寬平穩(wěn)時間序列，也叫廣義平穩(wěn)時間序列。我們以后所研究的時間序列主要是寬平穩(wěn)時間序列。 4按時間序列的分布規(guī)律來分，有高斯型時間序列和非高斯型時間序列。 1 確定性時間序列分析方法概述 時間序列預測技術(shù)就是通過對預測目標自身時間序列的處理，來研究其變化趨勢的。一個時間序列往往是以下幾類變化形式的疊加或耦合。 （ 1）長期趨勢變動。它是指時間序列朝著一定的方向持續(xù)上升或下降，或停留在某一水平上的傾向，它反映了客觀事物的主要變化趨勢。 （ 2）季節(jié)變動。 （ 3）循環(huán)變動。通常是指周期為一年以上，由非季節(jié)因素引起的漲落起伏波形相似的波動。 （ 4）不規(guī)則變動。通常它分為突然變動和隨機變動。 通常用tT 表示長期趨勢項，tS 表示季節(jié)變動趨勢項，tC 表示循環(huán)變動趨勢項，tR表示隨機干擾項。常見的確定性時間序列模型有以下幾種類型： （ 1）加法模型 tttttRCSTy += （ 2）乘法模型 tttttRCSTy = （ 3）混合模型 ttttRSTy += tttttRCTSy += 其中ty 是觀測目標的觀測記錄， 0)( =tRE ，22)( =tRE 。 如果在預測時間范圍以內(nèi)，無突然變動且隨機變動的方差2 較小，并且有理由認為過去和現(xiàn)在的演變趨勢將繼續(xù)發(fā)展到未來時，可用一些經(jīng)驗方法進行預測。 2 移動平均法 移動平均法是根據(jù)時間序列資料逐漸推移，依次計算包含一定項數(shù)的時序平均數(shù)，以反映長期趨勢的方法。當時間序列的數(shù)值由于受周期變動和不規(guī)則變動的影響，起伏較大，不易顯示出發(fā)展趨勢時，可用移動平均法，消除這些因素的影響，分析、預測序 -281-列的長期趨勢。 移動平均法有簡單移動平均法，加權(quán)移動平均法，趨勢移動平均法等。 2.1 簡單移動平均法 設(shè)觀測序列為Tyy ,1L，取移動平均的項數(shù) TN 。一次簡單移動平均值計算公式為： )(111)1(+=NttttyyyNML )(1)(1)(1)1(11 NtttNttNttyyNMyyNyyN+=+=L （ 1） 當預測目標的基本趨勢是在某一水平上下波動時，可用一次簡單移動平均方法建立預測模型： )(11)1(1 +=NttttyyNMyL，L,1, += NNt ， （ 2） 其預測標準誤差為： NTyySTNttt=+= 12)(， （ 3） 最近 N 期序列值的平均值作為未來各期的預測結(jié)果。一般 N 取值范圍：2005 N 。當歷史序列的基本趨勢變化不大且序列中隨機變動成分較多時， N 的取值應(yīng)較大一些。否則 N 的取值應(yīng)小一些。在有確定的季節(jié)變動周期的資料中，移動平均的項數(shù)應(yīng)取周期長度。選擇最佳 N 值的一個有效方法是，比較若干模型的預測誤差。預測標準誤差最小者為好。 例 1 某企業(yè) 1 月 11 月份的銷售收入時間序列如表 1 示。試用一次簡單滑動平均法預測第 12 月份的銷售收入。 表 1 企業(yè)銷售收入 月份 t 1 2 3 4 5 6 銷售收入ty 533.8 574.6 606.9 649.8 705.1 772.0 月份 t 7 8 9 10 11 銷售收入ty 816.4 892.7 963.9 1015.1 1102.7 解： 分別取 5,4 = NN 的預測公式 4321)1(1+=tttttyyyyy ， 11,5,4L=t 54321)2(1+=ttttttyyyyyy ， 11,5L=t 當 4=N 時，預測值 993.6)1(12=y ，預測的標準誤差為 150.5411)(1152)1(1=tttyyS 當 5=N 時，預測值 182.4)2(12=y ，預測的標準誤差為 -282- 958.2511)(1162)2(2=tttyyS 計算結(jié)果表明， 4=N 時，預測的標準誤差較小，所以選取 4=N 。預測第 12 月份的銷售收入為 993.6。 計算的 Matlab 程序如下： clc,clear y=533.8 574.6 606.9 649.8 705.1 772.0 816.4 892.7 963.9 1015.1 1102.7； m=length(y)； n=4,5； %n 為移動平均的項數(shù) for i=1:length(n) %由于 n 的取值不同， yhat 的長度不一致，下面使用了細胞數(shù)組 for j=1:m-n(i)+1 yhati(j)=sum(y(j:j+n(i)-1)/n(i)； end y12(i)=yhati(end)； s(i)=sqrt(mean(y(n(i)+1:m)-yhati(1:end-1).2)； end y12,s 簡單移動平均法只適合做近期預測，而且是預測目標的發(fā)展趨勢變化不大的情況。如果目標的發(fā)展趨勢存在其它的變化， 采用簡單移動平均法就會產(chǎn)生較大的預測偏差和滯后。 2.2 加權(quán)移動平均法 在簡單移動平均公式中，每期數(shù)據(jù)在求平均時的作用是等同的。但是，每期數(shù)據(jù)所包含的信息量不一樣，近期數(shù)據(jù)包含著更多關(guān)于未來情況的信心。因此，把各期數(shù)據(jù)等同看待是不盡合理的，應(yīng)考慮各期數(shù)據(jù)的重要性，對近期數(shù)據(jù)給予較大的權(quán)重，這就是加權(quán)移動平均法的基本思想。 設(shè)時間序列為LL,21 tyyy ；加權(quán)移動平均公式為 NNtNttwwwwywywywM+=+LL211221， Nt （ 4） 式中twM 為 t 期加權(quán)移動平均數(shù)；iw 為1+ity 的權(quán)數(shù)，它體現(xiàn)了相應(yīng)的ty 在加權(quán)平均數(shù)中的重要性。 利用加權(quán)移動平均數(shù)來做預測，其預測公式為 twtMy =+1 （ 5） 即以第 t 期加權(quán)移動平均數(shù)作為第 1+t 期的預測值。 例 2 我國 1979 1988 年原煤產(chǎn)量如表 2 所示，試用加權(quán)移動平均法預測 1989 年的產(chǎn)量。 表 2 我國原煤產(chǎn)量統(tǒng)計數(shù)據(jù)及加權(quán)移動平均預測值表 年份 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 原煤產(chǎn)量ty 6.35 6.20 6.22 6.66 7.15 7.89 8.72 8.94 9.28 9.8 三年加權(quán)移動平均預測值 6.235 6.4367 6.8317 7.4383 8.1817 8.6917 9.0733 -283-相對誤差（） 6.38 9.98 13.41 14.7 8.48 6.34 7.41 解 取 1,2,3321= www ，按預測公式 12323211+=+ttttyyyy 計算三年加權(quán)移動平均預測值，其結(jié)果列于表 2 中。 1989 年我國原煤產(chǎn)量的預測值為（億噸） 48.9694.828.928.931989=+=y 這個預測值偏低，可以修正。其方法是：先計算各年預測值與實際值的相對誤差，例如1982 年為 %38.666.6235.666.6=將相對誤差列于表 2 中，再計算總的平均相對誤差。 %5.9%100)44.5889.521(%1001 =ttyy由于總預測值的平均值比實際值低 %5.9 ，所以可將 1989 年的預測值修正為 4788.10%5.9148.9=計算的 MATLAB 程序如下： y=6.35 6.20 6.22 6.66 7.15 7.89 8.72 8.94 9.28 9.8； w=1/6；2/6；3/6； m=length(y)；n=3； for i=1:m-n+1 yhat(i)=y(i:i+n-1)*w； end yhat err=abs(y(n+1:m)-yhat(1:end-1)./y(n+1:m) T_err=1-sum(yhat(1:end-1)/sum(y(n+1:m) y1989=yhat(end)/(1-T_err) 在加權(quán)移動平均法中，tw 的選擇，同樣具有一定的經(jīng)驗性。一般的原則是：近期數(shù)據(jù)的權(quán)數(shù)大，遠期數(shù)據(jù)的權(quán)數(shù)小。至于大到什么程度和小到什么程度，則需要按照預測者對序列的了解和分析來確定。 2.3 趨勢移動平均法 簡單移動平均法和加權(quán)移動平均法，在時間序列沒有明顯的趨勢變動時，能夠準確反映實際情況。但當時間序列出現(xiàn)直線增加或減少的變動趨勢時，用簡單移動平均法和加權(quán)移動平均法來預測就會出現(xiàn)滯后偏差。因此，需要進行修正，修正的方法是作二次移動平均，利用移動平均滯后偏差的規(guī)律來建立直線趨勢的預測模型。這就是趨勢移動平均法。 一次移動的平均數(shù)為 -284- )(111)1(+=NttttyyyNML 在一次移動平均的基礎(chǔ)上再進行一次移動平均就是二次移動平均，其計算公式為 )(1)(1)1()1()2(1)1(1)1()2(NtttNtttMMNMMMNM+=+=L （ 6） 下面討論如何利用移動平均的滯后偏差建立直線趨勢預測模型。 設(shè)時間序列 ty 從某時期開始具有直線趨勢，且認為未來時期也按此直線趨勢變化，則可設(shè)此直線趨勢預測模型為 TbayttTt+=+ ，L,2,1=T （ 7） 其中 t 為當前時期數(shù)； T 為由 t 至預測期的時期數(shù)；ta 為截距；tb 為斜率。兩者又稱為平滑系數(shù)。 現(xiàn)在，我們根據(jù)移動平均值來確定平滑系數(shù)。由模型（ 7）可知 ttya = tttbyy =1tttbyy 22= ttNtbNyy )1(1=+所以 tttttttttNttttbNyNbNNyNbNybyyNyyyM21)1(21)1()(11)1(=+=+=+=+LLL因此 tttbNMy21)1(= （ 8） 由式（ 7） ，類似式（ 8）的推導，可得 tttbNMy21)1(11=（ 9） 所以 tttttbMMyy =)1(1)1(1（ 10） 類似式（ 8）的推導，可得 tttbNMM21)2()1(= （ 11） 于是，由式（ 8）和式（ 11）可得平滑系數(shù)的計算公式為 =)(122)2()1()2()1(ttttttMMNbMMa（ 12） 例 3 我國 1965 1985 年的發(fā)電總量如表 3 所示，試預測 1986 年和 1987 年的發(fā)電總量。 -285-表 3 我國發(fā)電量及一、二次移動平均值計算表 年份 t 發(fā)電總量 yt 一次移動平均， N 6 二次移動平均， N 61965 1 676 1966 2 825 1967 3 774 1968 4 716 1969 5 940 1970 6 1159 848.3 1971 7 1384 966.3 1972 8 1524 1082.8 1973 9 1668 1231.8 1974 10 1688 1393.8 1975 11 1958 1563.5 1181.1 1976 12 2031 1708.8 1324.5 1977 13 2234 1850.5 1471.9 1978 14 2566 2024.2 1628.8 1979 15 2820 2216.2 1792.8 1980 16 3006 2435.8 1966.5 1981 17 3093 2625 2143.4 1982 18 3277 2832.7 2330.7 1983 19 3514 3046 2530 1984 20 3770 3246.7 2733.7 1985 21 4107 3461.2 2941.2 解 由散點圖 1 可以看出，發(fā)電總量基本呈直線上升趨勢，可用趨勢移動平均法來預測。 0 5 10 15 20 2550010001500200025003000350040004500圖 1 原始數(shù)據(jù)散點圖 取 6=N ，分別計算一次和二次移動平均值并列于表 3 中。 2.3461)1(21=M ， 2.2941)2(21=M 再由公式（ 12） ，得 3981.12)2(21)1(2121= MMa 208)(162)2(21)1(2121= MMb 于是，得 21=t 時直線趨勢預測模型為 TyT2081.398121+=+預測 1986 年和 1987 年的發(fā)電總量為 1.4192121221986=+yyy 1.4397221231987=+yyy 計算的 MATLAB 程序如下： -286- clc,clear load y.txt %把原始數(shù)據(jù)保存在純文本文件 y.txt 中 m1=length(y)； n=6； %n 為移動平均的項數(shù) for i=1:m1-n+1 yhat1(i)=sum(y(i:i+n-1)/n； end yhat1 m2=length(yhat1)； for i=1:m2-n+1 yhat2(i)=sum(yhat1(i:i+n-1)/n； end yhat2 plot(1:21,y,*) a21=2*yhat1(end)-yhat2(end) b21=2*(yhat1(end)-yhat2(end)/(n-1) y1986=a21+b21 y1987=a21+2*b21 趨勢移動平均法對于同時存在直線趨勢與周期波動的序列， 是一種既能反映趨勢變化，又可以有效地分離出來周期變動的方法。 3 指數(shù)平滑法 一次移動平均實際上認為最近 N 期數(shù)據(jù)對未來值影響相同，都加權(quán)N1；而 N 期以前的數(shù)據(jù)對未來值沒有影響，加權(quán)為 0。但是，二次及更高次移動平均數(shù)的權(quán)數(shù)卻不是N1，且次數(shù)越高，權(quán)數(shù)的結(jié)構(gòu)越復雜，但永遠保持對稱的權(quán)數(shù)，即兩端項權(quán)數(shù)小，中間項權(quán)數(shù)大，不符合一般系統(tǒng)的動態(tài)性。一般說來歷史數(shù)據(jù)對未來值的影響是隨時間間隔的增長而遞減的。所以，更切合實際的方法應(yīng)是對各期觀測值依時間順序進行加權(quán)平均作為預測值。指數(shù)平滑法可滿足這一要求，而且具有簡單的遞推形式。 指數(shù)平滑法根據(jù)平滑次數(shù)的不同，又分為一次指數(shù)平滑法、二次指數(shù)平滑法和三次指數(shù)平滑法等，分別介紹如下。 3.1 一次指數(shù)平滑法 1預測模型 設(shè)時間序列為LL,21 tyyy ， 為加權(quán)系數(shù)， 10 0.00001 Terr=； for j=N+1:m-1 yhat(j)=w*yt(j-1:-1:j-N)； err=yt(j)-yhat(j)； Terr=Terr,abs(err)； w=w+2*k*err*yt(j-1:-1:j-N)； end Terr=max(Terr)； end w, yhat 5.2 kN, 值和初始權(quán)數(shù)的確定 在開始調(diào)整權(quán)數(shù)時，首先要確定權(quán)數(shù)個數(shù) N 和學習常數(shù) k 。一般說來，當時間序列的觀測值呈季節(jié)變動時， N 應(yīng)取季節(jié)性長度值。如序列以一年為周期進行季節(jié)變動時，若數(shù)據(jù)是月度的，則取 12=N ，若季節(jié)是季度的，則取 4=N 。如果時間序列無明顯的周期變動，則可用自相關(guān)系數(shù)法來確定，即取 N 為最高自相關(guān)系數(shù)的滯后時期。 k 的取值一般可定為 N/1 ， 也可以用不同的 k 值來進行計算， 以確定一個能使 S 最小的 k 值。 初始權(quán)數(shù)的確定也很重要，如無其它依據(jù)，也可用 N/1 作為初始權(quán)系數(shù)用，即 ),2,1(1NiNwiL= 自適應(yīng)濾波法有兩個明顯的優(yōu)點：一是技術(shù)比較簡單，可根據(jù)預測意圖來選擇權(quán)數(shù)的個數(shù)和學習常數(shù)，以控制預測。也可以由計算機自動選定。二是它使用了全部歷史數(shù)據(jù)來尋求最佳權(quán)系數(shù)，并隨數(shù)據(jù)軌跡的變化而不斷更新權(quán)數(shù)，從而不斷改進預測。 由于自適應(yīng)濾波法的預測模型簡單，又可以在計算機上對數(shù)據(jù)進行處理，所以這種預測方法應(yīng)用較為廣泛。 6 趨勢外推預測方法 趨勢外推法是根據(jù)事物的歷史和現(xiàn)時資料，尋求事物發(fā)展規(guī)律，從而推測出事物未來狀況的一種比較常用的預測方法。利用趨勢外推法進行預測，主要包括六個階段：（ a）選擇應(yīng)預測的參數(shù)； （ b）收集必要的數(shù)據(jù)； （ c）利用數(shù)據(jù)擬合曲線； （ d）趨勢外推； （ e）預測說明； （ f）研究預測結(jié)果在進行決策中應(yīng)用的可能性。 趨勢外推法常用的典型數(shù)學模型有：指數(shù)曲線、修正指數(shù)曲線、生長曲線、包絡(luò)曲線等。 6.1 指數(shù)曲線法 一般來說，技術(shù)的進步和生產(chǎn)的增長，在其未達飽和之前的新生時期是遵循指數(shù)曲線增長規(guī)律的，因此可以用指數(shù)曲線對發(fā)展中的事物進行預測。 指數(shù)曲線的數(shù)學模型為 Kteyy0= （ 35） 其中系數(shù)0y 和 K 值由歷史數(shù)據(jù)利用回歸方法求得。對式（ 35）取對數(shù)可得 Ktyy +=0lnln （ 36） -297-令 yY ln= ，0ln yA = 則 KtAY += 其中 KA, 可以用最小二乘法求得。 6.2 修正指數(shù)曲線法 利用指數(shù)曲線外推來進行預測時，存在著預測值隨著時間的推移會無限增大的情況。這是不符合客觀規(guī)律的。因為任何事物的發(fā)展都是有一定限度的。例如某種暢銷產(chǎn)品，在其占有市場的初期是呈指數(shù)曲線增長的，但隨著產(chǎn)品銷售量的增加，產(chǎn)品總量接近于社會飽和量時。這時的預測模型應(yīng)改用修正指數(shù)曲線。 ttabKy += （ 37） 在此數(shù)學模型中有三個參數(shù) aK, 和 b 要用歷史數(shù)據(jù)來確定。 修正指數(shù)曲線用于描述這樣一類現(xiàn)象。 （ 1）初期增長迅速，隨后增長率逐漸降低。 （ 2）當 0K ， 0a ， 10 r 。解此微分方程得 rtceLy+=1（ 51） 式中 c 為常數(shù)。 下面我們記 Logistic 曲線的一般形式為 ttabKy+=1， 0K ， 0a ， 10 。假設(shè)),(ARMA110mnA = 模型的殘差ta 之平方和， ),(ARMA221mnA = 模型的殘差ta 之平方和， N 是采集數(shù)據(jù)的數(shù)目，則檢驗準則為： ),(001= NsFNAsAAF ， 其中22mn += ， )(1122mnmns += 。 若這樣得到的 F 值超過由 F 分布查表所得的在 5%置信水平上的 ),( NsF 值，那么由 ),(ARMA11mn 模型改變?yōu)?),(ARMA22mn 時，殘差平方和的改善是顯著的，因而拒絕關(guān)于模型 ),(ARMA11mn 的適用性假設(shè)； F 值低于查表所得之值，就可以認為在該置信水平上這個模型是適用的。 5檢查122,nn 的值是否很小，其置信區(qū)間是否包含零。若不是，則適用的模型就是 )12,2ARMA( nn 。 若122,nn 很小，且其置信區(qū)間包含零，則擬合 )22,12ARMA( nn 。 6利用 F 準則檢驗?zāi)Ｐ?)12,2ARMA( nn 和 )22,12ARMA( nn ，若 F 值不顯著，轉(zhuǎn)入第 7 步；若 F 值顯著，轉(zhuǎn)入第 8 步。 7舍棄小的 MA 參數(shù)，擬合 22 nm 的模型 ),12ARMA( mn ，并用 F 準則進行檢驗。重復這一過程，直到得出具有最小參數(shù)的適用模型為止。 8舍棄小的 MA 參數(shù)，擬合 12 nm 的模型 ),2ARMA( mn ，并用 F 準則進行檢驗。重復這