利用Excel進(jìn)行時(shí)間序列的譜分析

1、利用 Excel進(jìn)行時(shí)間序列的譜分析（I）在頻域分析中， 功率譜是揭示時(shí)間序列周期特性的最為有力的工具之一。下面列舉幾個(gè)例子，分別從不同的角度識(shí)別時(shí)間序列的周期。1 時(shí)間序列的周期圖【例 1】某水文觀測(cè)站測(cè)得一條河流從1979 年 6 月到 1980 年 5 月共計(jì) 12 月份的斷面平均流量。試判斷該河流的徑流量變化是否具有周期性，周期長(zhǎng)度大約為多少？分析：假定將時(shí)間序列xt 展開為 Fourier 級(jí)數(shù)，則可表示為kxt( aicos2fi t bi sin 2 f i t )t(1)i1式中 f 為頻率， t 為時(shí)間序號(hào)， k 為周期分量的個(gè)數(shù)即主周期（基波）及其諧波的個(gè)數(shù)，it為標(biāo)準(zhǔn)誤差

2、（白噪聲序列） 。當(dāng)頻率 fi 給定時(shí)，式（ 1）可以視為多元線性回歸模型，可以證明，待定系數(shù) ai、 bi 的最小二乘估計(jì)為a?i2 NtN txt cos2 fi1（ 2）2 N?xt sin 2 fi tbiN t 1這里 N 為觀測(cè)值的個(gè)數(shù)。定義時(shí)間序列的周期圖為I ( fi )N ( ai2bi2 ) ， i1,2, k(3)2式中 I(fi)為頻率 f i 處的強(qiáng)度。以fi 為橫軸，以 I(fi )為縱軸，繪制時(shí)間序列的周期圖，可以在最大值處找到時(shí)間序列的周期。對(duì)于本例，N=12 ，t=1,2, ,N， fi=i/N，下面借助 Excel，利用上述公式，計(jì)算有關(guān)參數(shù)并分析時(shí)間序列的

3、周期特性。第一步，錄入數(shù)據(jù)，并將數(shù)據(jù)標(biāo)準(zhǔn)化或中心化（圖 1）。圖 1 錄入的數(shù)據(jù)及其中心化結(jié)果1中心化與標(biāo)準(zhǔn)化的區(qū)別在于，只需將原始數(shù)據(jù)減去均值，而不必再除以標(biāo)準(zhǔn)差。不難想到，中心化的數(shù)據(jù)均值為0，但方差與原始數(shù)據(jù)相同（未必為1）。第二步，計(jì)算三角函數(shù)值為了借助式（1）計(jì)算參數(shù)ai、 bi，首先需要計(jì)算正弦值和余弦值。取 i 1,2, ,6 ，則頻率為f i i / N 1/ 12,2 / 12, ,6 / 12 （圖 1）。將頻率寫在單元格C3-C14 中（根據(jù)對(duì)稱性，我們只用前6 個(gè)），將中心化的數(shù)據(jù)轉(zhuǎn)置粘貼于第一行的單元格D1-O1中，月份的序號(hào)寫在單元格D2-O2 中（與中心化數(shù)據(jù)對(duì)齊

4、） 。圖 2 計(jì)算余弦值的表格在 D2 單元格中輸入公式“ =COS($B$1*$D$2*C3) ”，回車得到 0.866；按住單元格的右下角右拉至 O3 單元格， 得到 f=1/12=0.083 ，t=1,2, ,12 時(shí)的全部余弦值。 在 D2 單元格中輸入公式“ =COS($B$1*$D$2*C4) ”，回車得到0.5；按住單元格的右下角右拉至O4 單元格，得到 f=2/12=0.167 ，t=1,2, ,12 時(shí)的全部余弦值。 依次類推， 可以算出全部所要的余弦值（在D3-O8 區(qū)域中）。根據(jù)對(duì)稱性，我們的計(jì)算到k=6 為止（ 圖 2）。注意，這里B1 單元格是2=6.28319（圖中

5、未能顯示） 。在上面的計(jì)算中，只要將公式中的“COS”換成“ SIN ”，即可得到正弦值，不過(guò)為了計(jì)算過(guò)程清楚明白，最好在另外一個(gè)區(qū)域給出結(jié)算結(jié)果（在D17-O22 區(qū)域中，參見 圖 3）。圖 3 計(jì)算正弦值的表格第三步，計(jì)算參數(shù) ai、 bi利用中心化的數(shù)據(jù)（仍然表作xt）計(jì)算參數(shù) ai、bi。首先算出 xtcos2fit 和 xt sins2fit。在D9 單元格中輸入公式“ =D1*D3 ”，回車得到18.309；按住單元格的右下角右拉至O9 單元格，得到 f=1/12=0.083 ， t=1,2, ,12 時(shí)的全部xtcos2fit 值；加和得 39.584，再除以6，即得a1=6.5

6、97。在 D10 單元格中輸入公式“ =D1*D4 ”，回車得到 10.571；按住單元格的右下角右拉至 O10 單元格， 得到 f=2/12=0.083 ，t=1,2, ,12 時(shí)的全部 xt cos2fit 值；加和得 -365.25，再除以 6，得到 a2=-60.875 。其余依此類推。將上面公式中的余弦值換成正弦值，即可得到bi 值（見下表） 。上面的計(jì)算過(guò)程相當(dāng)于2采用式（ 2）進(jìn)行逐步計(jì)算。第四步，計(jì)算頻率強(qiáng)度利用式（ 3），非常容易算出I(fi)值。例如I ( f1 )N (a12b12 ) 6 * (6.597 2170.2132 ) 174096 .9142其余依此類推（見

7、圖 4）。圖 4 計(jì)算頻率強(qiáng)度第五步，繪制時(shí)間序列周期圖利用圖 4 中的數(shù)據(jù)，不難畫出周期圖（圖 5）。200000180000160000140000度 120000強(qiáng)100000率頻80000600004000020000000.10.20.30.40.50.6頻率圖 5 某河流徑流量的周期圖（1979 年 6 月 1980年 5月）第六步，周期識(shí)別關(guān)鍵是尋找頻率的極值點(diǎn)或突變點(diǎn)。在本例中，沒有極值點(diǎn)，但在f1 =1/12=0.0833 處，頻率強(qiáng)度突然增加（陡增），而此時(shí) T=1/f1=12 ，故可判斷時(shí)間序列可能存在一個(gè)12 月的周期，即 1 年周期。3【例 2】為了映證上述判斷，我們

8、借助同一條河流的連續(xù)兩年的平均月徑流量（1961 年 6月 1963 年 5 月）。原始數(shù)據(jù)見下圖（圖 6）。圖 6 原始數(shù)據(jù)及部分處理結(jié)果將原始數(shù)據(jù)回車時(shí)間序列變化圖，可以初步估計(jì)具有12 月變化周期，但不能肯定（圖6）。350300250量流200徑均平 150月100500051015202530時(shí)間（月份序號(hào)）圖 6 徑流量的月變化圖（1961 年 6 月 1963 年 5 月）4按照例 1 給出的計(jì)算步驟，計(jì)算參數(shù)數(shù) ai、bi，進(jìn)而計(jì)算頻率強(qiáng)度（結(jié)果將 圖 7）。然后繪制時(shí)間序列的周期圖（ 圖 8）。注意這里， N=24 ，我們?nèi)?k=12 。圖 7 參數(shù)和頻率強(qiáng)度的計(jì)算結(jié)果從圖

9、8 中可以看出，頻率強(qiáng)度的最大值（極值點(diǎn)）對(duì)應(yīng)于頻率f1 =1/12=0.0833 ，故時(shí)間序列的周期判斷為 T=1/f1=12。這與用 12 月的數(shù)據(jù)進(jìn)行估計(jì)的結(jié)果是一致的，但由于例 2 的時(shí)間序列比例1 的時(shí)間序列長(zhǎng)1 倍，故判斷結(jié)果更為可靠。140000120000100000度80000強(qiáng)率頻 60000 40000 20000 000.10.20.30.40.50.6頻率圖 8某河流徑流量的周期圖（1961 年 6 月 1963年 5月）2 時(shí)間序列的頻譜圖【例 3】首先考慮對(duì)例1 的數(shù)據(jù)進(jìn)行功率譜分析。例1 的時(shí)間序列較短，分析的效果不佳，但計(jì)算過(guò)程簡(jiǎn)短。給出這個(gè)例子，主要是幫助大

10、家理解Fourier 變換過(guò)程和方法。為了進(jìn)行 Fourier 分析，需要對(duì)數(shù)據(jù)進(jìn)行預(yù)處理。第一，將數(shù)據(jù)中心化，即用原始數(shù)據(jù)減去其平均值。中心化的數(shù)據(jù)均值為 0，我們對(duì)中心化的數(shù)據(jù)進(jìn)行變換，其周期更為明顯。第二，由于 Fourier 分析通常采用所謂快速 Fourier 變換（ Fast Fourier Transformation ， FFT），而 FFT要求數(shù)據(jù)必須為2n 個(gè)，這里n 為正整數(shù)（ 1,2,3,），而我們的樣本為N=12，它不是 2 的某5個(gè) n 次方。因此，在中心化的數(shù)據(jù)后面加上4 個(gè) 0，這樣新的樣本數(shù)為N =12 4 16 24個(gè)，這才符合 FFT 的需要（ 圖 9）。

11、下面，我們對(duì)延長(zhǎng)后的中心化數(shù)據(jù)進(jìn)行Fourier 變換。圖 9 數(shù)據(jù)的中心化與“延長(zhǎng)”第一步，打開Foureir 分析對(duì)話框沿著主菜單的“工具（ Tools）”“數(shù)據(jù)分析（ Data Analysis ）”路徑打開數(shù)據(jù)分析選項(xiàng)框（圖 10），從中選擇“傅立葉分析（ Fourier Analysis ）”。圖 10 在數(shù)據(jù)分析選項(xiàng)框中選擇Fourier分析第二步，定義變量和輸出區(qū)域確定之后，彈出傅立葉分析對(duì)話框，根據(jù)數(shù)據(jù)在工作表中的分布情況進(jìn)行如下設(shè)置：將光標(biāo)置于“輸入?yún)^(qū)域”對(duì)應(yīng)的空白欄，然后用鼠標(biāo)選中單元格C1-C17 ，這時(shí)空白欄中自動(dòng)以絕對(duì)單元格的形式定義中心化數(shù)據(jù)的區(qū)域范圍（即$C$1-

12、$C$17 ）。選中“標(biāo)志位于第一行”。選中輸出區(qū)域，定義范圍為D2-D17 （圖 11）。注意：如果輸入?yún)^(qū)域的數(shù)據(jù)范圍定義為C2-C17 ，則不要選中“標(biāo)志位于第一行”，這與回歸分析中的原始變量定義是一樣的（ 圖 12）。如果不定義輸出區(qū)域范圍，則變換結(jié)果將會(huì)自動(dòng)給在新的工作表組上。這一點(diǎn)也與回歸分析一樣。6圖 11 選中“標(biāo)志位于第一行”與數(shù)據(jù)輸入范圍的定義圖 12 不選中“標(biāo)志位于第一行”與數(shù)據(jù)輸入范圍的定義第三步，結(jié)果轉(zhuǎn)換定義數(shù)據(jù)輸入輸出區(qū)域完成之后，確定，立即得到Fourier 變換的結(jié)果（圖 13）。圖 13 傅立葉變換的結(jié)果7變換的結(jié)果為一組復(fù)數(shù)，相當(dāng)于將 f(t)變成了 F()

13、，實(shí)際上是將 xt 變成了 XT(f)。我們知道，有了 f(t) 的象函數(shù) F( ) 就可以計(jì)算能量譜密度函數(shù)S()，即有S()F( )F()F (2(4)相應(yīng)地，有了XT(f)也就容易計(jì)算功率譜（密度）2X T ( f )(5)P( f )T式中 f 表示線頻率，與角頻率的轉(zhuǎn)換關(guān)系是 =2/T，這里 T 為數(shù)據(jù)區(qū)間長(zhǎng)度。如果將 XT(f)表作 XT( f)=A+jB （這里 A 為實(shí)部， B 為虛部），則有XT( f) 2XT( f) XT( fi) ( AjB)( AjB)A2B2（ 6）iiiiiiii因此這一步是要分離變換結(jié)果的實(shí)部與虛部。逐個(gè)手動(dòng)提取是非常麻煩而且容易出錯(cuò)的，可以利用

14、Excel 有關(guān)復(fù)數(shù)計(jì)算的命令。提取實(shí)部的Excel 命令是 imreal 。在 H2單元格中輸入命令 “ =IMREAL(D2)”（這里 D2 為變換結(jié)果的第一個(gè)復(fù)數(shù)所在的單元格），回車得到第一個(gè)復(fù)數(shù)的實(shí)部0；點(diǎn)中H2 單元格的右下角，撳住鼠標(biāo)左鍵下拉至H17，得到全部的實(shí)部數(shù)值。提取虛部的命令是imaginary 。在 I2 單元格中輸入公式“=IMAGINARY(D2)”，回車得到第一個(gè)復(fù)數(shù)的虛部0；下拉至I17，得到全部的虛部數(shù)值。根據(jù)式（5）、（ 6），功率譜密度的計(jì)算公式為P( fi )Ai2Bi2（ 7）T考慮到本例中T=N=16，在 J2 單元格中輸入公式“=(H22+I22)

15、/16”，回車得到第一個(gè)功率譜密度 0；下拉至J17，得到全部譜密度數(shù)值（圖 14）?；?FFT 的譜密度分布是對(duì)稱的，可以看出，以J10 所在的3105.28 為對(duì)稱點(diǎn)，上下的數(shù)值完全對(duì)稱。圖 14 功率譜密度的計(jì)算結(jié)果8第四步，繪制頻譜圖線頻率 fi 可以表作f ii / Ti / N ， i0,1,2, N -1顯然 f0 =0/16=0 ， f1 =1/16=0.0625 ， f2=2/16=0.125 ， f15=15/16=0.9375 。在 Excel 中，容易計(jì)算頻率的數(shù)值。將頻率與功率譜對(duì)應(yīng)起來(lái)（圖 15），就可以畫出頻譜圖。如果補(bǔ)上最后一個(gè)頻率數(shù)值f16=1 及其對(duì)應(yīng)的功

16、率0，則可畫出完全對(duì)稱的譜圖（圖 16）。圖 15 功率譜密度與頻率的對(duì)應(yīng)關(guān)系7000060000)50000f(P度40000密譜30000率功20000 10000000.20.40.60.81頻率f圖 16對(duì)稱分布的頻譜圖由于功率譜圖的對(duì)稱性，畫出整個(gè)譜圖實(shí)在沒有必要，因此，在實(shí)際工作中，通常只畫出左半邊（ 圖 17）。9700006000050000)f40000(P率功 30000 20000 10000 000.10.20.30.40.50.6頻率f圖 17 實(shí)用的頻譜圖第五步，功率譜分析頻域分析的主要目標(biāo)之一是判斷時(shí)間序列是否存在周期性。從圖 17 可以看出，功率最大點(diǎn)對(duì)應(yīng)的頻率

17、是f1=0.0625，該頻率對(duì)應(yīng)的周期長(zhǎng)度為16?？梢?，在時(shí)間序列較短的情況下之間用功率譜圖尋找時(shí)間序列的周期不如周期圖準(zhǔn)確。另外可以初步估計(jì)數(shù)據(jù)的性質(zhì)。在圖 17 中，去掉第一個(gè)0 點(diǎn)，剩余的點(diǎn)一般呈冪指數(shù)分布（在雙對(duì)數(shù)坐標(biāo)圖上點(diǎn)列具有直線趨勢(shì)），可以擬合冪指數(shù)函數(shù)如下：P( f ) f（ 8）100000)f(P10000率功P (f ) = 935.68f -1.4952R2 = 0.923310000.010.11頻率f圖 18功率 P(f) 與頻率 f的雙對(duì)數(shù)坐標(biāo)圖結(jié)果得到功率譜指數(shù)=1.4952 1.5。功率譜指數(shù)與時(shí)間序列的Hurst 指數(shù)具有如下關(guān)系2H 1（ 9）據(jù)此估計(jì) H

18、urst 指數(shù)約為0.25。我們知道， Hurst 指數(shù)介于0 1之間，當(dāng) H>0.5 時(shí)，表明時(shí)間序列存在正的自相關(guān)，意味著系統(tǒng)演化具有持久性；當(dāng)H<0時(shí)，表明時(shí)間序列具有負(fù)的10自相關(guān)，意味著系統(tǒng)演化具有反持久性；當(dāng)H=0 時(shí)，表明時(shí)間序列不存在自相關(guān)，過(guò)去與未來(lái)無(wú)關(guān)。對(duì)于這條河流的徑流量而言，H=0.25<0.5，表明時(shí)間序列具有反持久性：過(guò)去的增量意味著今后的減少，過(guò)去的減少意味著未來(lái)的增加。因此，徑流量必然周期性的變化。【例 4】下面對(duì)前述例2 的數(shù)據(jù)進(jìn)行Fourier 變換，方法與例3 相同，但由于N=24，我們?nèi)?T=32=2 5。也就是說(shuō)，對(duì)于中心化的數(shù)據(jù)，要

19、在后面添加8 個(gè) 0 作為補(bǔ)充點(diǎn)數(shù)?；贔FT 的變換結(jié)果如下（見圖 19）。圖 19 例 2 數(shù)據(jù)（經(jīng)中心化處理）的FFT 變換結(jié)果計(jì)算功率譜除例3 講述的方法外，還可以利用Excel 的另外兩個(gè)命令實(shí)現(xiàn)：一是計(jì)算共軛復(fù)數(shù)的命令imconjugate ，首先求出X T ( f ) 的共軛復(fù)數(shù)X T ( f ) ；然后借助復(fù)數(shù)的乘積命令 improduct ，計(jì)算復(fù)數(shù)的X T ( f ) 與 X T ( f ) 的乘積 X T ( f ) 2 ；最后利用式（5）得到功率譜。不過(guò)，此時(shí)的時(shí)間序列長(zhǎng)度視為T=32。例如在 H2 中鍵入公式 “ =IMCONJUGA TE(D2) ”，回車得到第一個(gè)

20、共軛數(shù)， 下拉至 H33，得到全部共軛數(shù)值。在 L2 中鍵入公式“ =IMPRODUCT(D2,H2) ”，回車，得到第一個(gè)復(fù)數(shù)乘22積，下拉至 L33 ，得到全部X T ( f )值。最后用X T ( f )除以 32 得到功率譜密度 （圖 20）。11根據(jù)計(jì)算結(jié)果畫出頻譜圖（圖 21），從圖上可以看出，頻率密度的極值點(diǎn)對(duì)應(yīng)的頻率為0.09375 ，相應(yīng)的周期為T=10.667；在極值點(diǎn)附件存在一個(gè)次最大點(diǎn)，但相對(duì)于其他數(shù)值卻顯然又是突變點(diǎn)，該點(diǎn)對(duì)應(yīng)的頻率為 0.0625，相應(yīng)的周期為 16。故可斷定，該時(shí)間序列的周期比在 10 16 之間。圖 20 功率譜密度值及其相應(yīng)的頻率3500030

21、00025000)f20000(P率功 15000 100005000000.10.20.30.40.50.6頻率f圖 21 例 4的頻譜圖12不過(guò)，由于這里的數(shù)據(jù)包含兩個(gè)周期在內(nèi)， 用它估計(jì)數(shù)據(jù)的自相關(guān)性質(zhì)不太準(zhǔn)確。 最后給一個(gè)較長(zhǎng)時(shí)間序列的分析實(shí)例?！纠?5】某海域測(cè)得多年連續(xù)的海平面年平均高度，發(fā)現(xiàn)大約每隔11 年左右海平面達(dá)到一個(gè)最高值（圖 22），于是科學(xué)家判斷海平面的升降存在一個(gè)11 年周期，與太陽(yáng)黑子（ sunspot）的 11 年周期變化一致，而太陽(yáng)黑子的活動(dòng)正是海平面升降的原因所在。問(wèn)題在于，前述 11 年周期是通過(guò)原始數(shù)據(jù)的變化曲線直觀發(fā)現(xiàn)的，未必可靠。為此，可以進(jìn)行一個(gè)功

22、率譜分析，判斷這種周期是否確實(shí)存在。圖 22 某海域海面年平均高度的時(shí)間序列數(shù)據(jù)首先利用原始數(shù)據(jù)畫出時(shí)間序列變化的曲線圖，觀測(cè)數(shù)據(jù)變化特征，發(fā)現(xiàn)具有周期性波動(dòng)特征，最高峰的時(shí)間間隔大約為10 12 年之間（ 圖 23）。180160140120度高 100面海806040200020406080100年度序號(hào)圖 23 某海域海平面年平均高度的年際變化曲線然后將原始數(shù)據(jù)中心化，取T=128=2 7，這意味著需要將時(shí)間序列延長(zhǎng)到128 位，在計(jì)13算頻率時(shí)用0,1,2,除以 128，在計(jì)算功率譜密度時(shí)，式（ 5）中的序列長(zhǎng)度取T=128。Fourier變換和功率譜密度的計(jì)算步驟與例3、例 4 相同

23、，不贅述。下面直接給出變換結(jié)果（圖 24，圖 25）。圖 24海面高度時(shí)間序列的FFT 變換的部分結(jié)果180001600014000)12000f10000(P率8000功600040002000000.10.20.30.40.50.6頻率f圖 25 海面高度變化的頻譜圖在 Excel 上，將鼠標(biāo)指向功率譜密度的最大值處，立即顯示頻率為f=0.09375，對(duì)應(yīng)的功率為 P(f)=16832.77972 。根據(jù)此處的頻率可以算出周期為T=1/f=1/0.09375=10.667 11。這14表明，海平面高度的變化的確存在一個(gè)為期大約11 年的變化周期，與直觀判斷結(jié)果一致。由于這里采用的時(shí)間序列較

24、長(zhǎng)，故結(jié)果比較可靠。太陽(yáng)黑子的活動(dòng)周期約為11 年稍多一點(diǎn)，二者非常接近。 因此，海面高度變化與太陽(yáng)黑子周期具有對(duì)應(yīng)關(guān)系的猜想，在時(shí)間序列變化的節(jié)律方面，大致是可以接受的。圖 26 功率譜密度最大值對(duì)應(yīng)的頻率顯示結(jié)果【附錄 】我們?cè)谇懊嬲f(shuō)過(guò)，式 (1) 實(shí)則一個(gè)多元線性回歸模型，式 (2)是對(duì)式 (1) 中待定參數(shù)的最小二乘 (OLS) 估計(jì)。下面驗(yàn)證這種判斷。將 圖 3 中計(jì)算出的正弦值 SIN 、余弦值 COS 和中心化的徑流量 Xt 集中在一起， 經(jīng)復(fù)制選擇性粘貼轉(zhuǎn)置， 可將數(shù)據(jù)重新排列如下 （圖 27）。圖 27 重新排列的數(shù)據(jù)若以正弦、余弦值為自變量，以中心化的徑流量為因變量進(jìn)行多元回歸，Excel 拒絕給出結(jié)果，并彈出如下對(duì)話框顯示拒絕計(jì)算的原因。圖 28 Excel拒絕給出回歸結(jié)果15問(wèn)題在于行數(shù)與列數(shù)不能相同，而本例中自變量和樣本數(shù)都是12，這就是問(wèn)題所在?？紤]到最后一個(gè)變量全都是0，不妨去掉最后一個(gè)變量。由于式（1）不含常數(shù)項(xiàng)，在回歸分析選項(xiàng)框中強(qiáng)制“常數(shù)為 0”。確定以后， Excel 給出的多元回歸結(jié)果如下（圖 29）。圖 29 利用圖27 中的前 12 列數(shù)據(jù)進(jìn)行多元回歸的結(jié)果將回歸結(jié)果與利用式（2）直接結(jié)算的結(jié)果進(jìn)行比較，發(fā)現(xiàn)除了SIN6 的系數(shù)沒有給出、COS6 的系數(shù)相差一倍以外，其余的結(jié)果基本一樣（比較圖 29 與下表）。將圖 27