利用Excel進行時間序列的譜分析(I)

1、利用Excel進行時間序列的譜分析（I）在頻域分析中，功率譜是揭示時間序列周期特性的最為有力的工具之一。下面列舉幾個例子，分別從不同的角度識別時間序列的周期。1 時間序列的周期圖【例1】某水文觀測站測得一條河流從1979年6月到1980年5月共計12月份的斷面平均流量。試判斷該河流的徑流量變化是否具有周期性，周期長度大約為多少？分析：假定將時間序列xt展開為Fourier級數(shù)，則可表示為 (1)式中fi為頻率，t為時間序號，k為周期分量的個數(shù)即主周期（基波）及其諧波的個數(shù)，t為標準誤差（白噪聲序列）。當頻率fi給定時，式（1）可以視為多元線性回歸模型，可以證明，待定系數(shù)ai、bi的最小二乘估計

2、為 （2）這里N為觀測值的個數(shù)。定義時間序列的周期圖為， (3)式中I(fi)為頻率fi處的強度。以fi為橫軸，以I(fi)為縱軸，繪制時間序列的周期圖，可以在最大值處找到時間序列的周期。對于本例，N=12，t=1,2,N，fi=i/N，下面借助Excel，利用上述公式，計算有關(guān)參數(shù)并分析時間序列的周期特性。第一步，錄入數(shù)據(jù)，并將數(shù)據(jù)標準化或中心化（圖1）。圖1 錄入的數(shù)據(jù)及其中心化結(jié)果中心化與標準化的區(qū)別在于，只需將原始數(shù)據(jù)減去均值，而不必再除以標準差。不難想到，中心化的數(shù)據(jù)均值為0，但方差與原始數(shù)據(jù)相同（未必為1）。第二步，計算三角函數(shù)值為了借助式（1）計算參數(shù)ai、bi，首先需要計算正弦

3、值和余弦值。取，則頻率為（圖1）。將頻率寫在單元格C3-C14中（根據(jù)對稱性，我們只用前6個），將中心化的數(shù)據(jù)轉(zhuǎn)置粘貼于第一行的單元格D1-O1中，月份的序號寫在單元格D2-O2中（與中心化數(shù)據(jù)對齊）。圖2 計算余弦值的表格在D2單元格中輸入公式“=COS($B$1*$D$2*C3)”，回車得到0.866；按住單元格的右下角右拉至O3單元格，得到f=1/12=0.083，t=1,2,12時的全部余弦值。在D2單元格中輸入公式“=COS($B$1*$D$2*C4)”，回車得到0.5；按住單元格的右下角右拉至O4單元格，得到f=2/12=0.167，t=1,2,12時的全部余弦值。依次類推，可以算

4、出全部所要的余弦值（在D3-O8區(qū)域中）。根據(jù)對稱性，我們的計算到k=6為止（圖2）。注意，這里B1單元格是2=6.28319（圖中未能顯示）。在上面的計算中，只要將公式中的“COS”換成“SIN”，即可得到正弦值，不過為了計算過程清楚明白，最好在另外一個區(qū)域給出結(jié)算結(jié)果（在D17-O22區(qū)域中，參見圖3）。圖3 計算正弦值的表格第三步，計算參數(shù)ai、bi利用中心化的數(shù)據(jù)（仍然表作xt）計算參數(shù)ai、bi。首先算出xtcos2fit和xtsins2fit。在D9單元格中輸入公式“=D1*D3”，回車得到18.309；按住單元格的右下角右拉至O9單元格，得到f=1/12=0.083，t=1,2,

5、12時的全部xtcos2fit值；加和得39.584，再除以6，即得a1=6.597。在D10單元格中輸入公式“=D1*D4”，回車得到10.571；按住單元格的右下角右拉至O10單元格，得到f=2/12=0.083，t=1,2,12時的全部xtcos2fit值；加和得-365.25，再除以6，得到a2=-60.875。其余依此類推。將上面公式中的余弦值換成正弦值，即可得到bi值（見下表）。上面的計算過程相當于采用式（2）進行逐步計算。第四步，計算頻率強度利用式（3），非常容易算出I(fi)值。例如其余依此類推（見圖4）。圖4 計算頻率強度第五步，繪制時間序列周期圖利用圖4中的數(shù)據(jù)，不難畫出周

6、期圖（圖5）。圖5 某河流徑流量的周期圖（1979年6月1980年5月）第六步，周期識別關(guān)鍵是尋找頻率的極值點或突變點。在本例中，沒有極值點，但在f1=1/12=0.0833處，頻率強度突然增加（陡增），而此時T=1/f1=12，故可判斷時間序列可能存在一個12月的周期，即1年周期。【例2】為了映證上述判斷，我們借助同一條河流的連續(xù)兩年的平均月徑流量（1961年6月1963年5月）。原始數(shù)據(jù)見下圖（圖6）。圖6 原始數(shù)據(jù)及部分處理結(jié)果將原始數(shù)據(jù)回車時間序列變化圖，可以初步估計具有12月變化周期，但不能肯定（圖6）。圖6 徑流量的月變化圖（1961年6月1963年5月）按照例1給出的計算步驟，計

7、算參數(shù)數(shù)ai、bi，進而計算頻率強度（結(jié)果將圖7）。然后繪制時間序列的周期圖（圖8）。注意這里，N=24，我們?nèi)=12。圖7 參數(shù)和頻率強度的計算結(jié)果從圖8中可以看出，頻率強度的最大值（極值點）對應于頻率f1=1/12=0.0833，故時間序列的周期判斷為T=1/f1=12。這與用12月的數(shù)據(jù)進行估計的結(jié)果是一致的，但由于例2的時間序列比例1的時間序列長1倍，故判斷結(jié)果更為可靠。圖8 某河流徑流量的周期圖（1961年6月1963年5月）2 時間序列的頻譜圖【例3】首先考慮對例1的數(shù)據(jù)進行功率譜分析。例1的時間序列較短，分析的效果不佳，但計算過程簡短。給出這個例子，主要是幫助大家理解Fouri

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

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

10、入?yún)^(qū)域的數(shù)據(jù)范圍定義為C2-C17，則不要選中“標志位于第一行”，這與回歸分析中的原始變量定義是一樣的（圖12）。如果不定義輸出區(qū)域范圍，則變換結(jié)果將會自動給在新的工作表組上。這一點也與回歸分析一樣。圖11 選中“標志位于第一行”與數(shù)據(jù)輸入范圍的定義圖12 不選中“標志位于第一行”與數(shù)據(jù)輸入范圍的定義第三步，結(jié)果轉(zhuǎn)換定義數(shù)據(jù)輸入輸出區(qū)域完成之后，確定，立即得到Fourier變換的結(jié)果（圖13）。 圖13 傅立葉變換的結(jié)果變換的結(jié)果為一組復數(shù)，相當于將f(t)變成了F()，實際上是將xt變成了XT(f)。我們知道，有了f(t)的象函數(shù)F()就可以計算能量譜密度函數(shù)S()，即有 (4)相應地，有了

11、XT(f)也就容易計算功率譜（密度） (5)式中f表示線頻率，與角頻率的轉(zhuǎn)換關(guān)系是=2/T，這里T為數(shù)據(jù)區(qū)間長度。如果將XT(f)表作XT(f)=A+jB（這里A為實部，B為虛部），則有 （6）因此這一步是要分離變換結(jié)果的實部與虛部。逐個手動提取是非常麻煩而且容易出錯的，可以利用Excel有關(guān)復數(shù)計算的命令。提取實部的Excel命令是imreal。在H2單元格中輸入命令“=IMREAL(D2)”（這里D2為變換結(jié)果的第一個復數(shù)所在的單元格），回車得到第一個復數(shù)的實部0；點中H2單元格的右下角，撳住鼠標左鍵下拉至H17，得到全部的實部數(shù)值。提取虛部的命令是imaginary。在I2單元格中輸入公

12、式“=IMAGINARY(D2)”，回車得到第一個復數(shù)的虛部0；下拉至I17，得到全部的虛部數(shù)值。根據(jù)式（5）、（6），功率譜密度的計算公式為 （7）考慮到本例中T=N=16，在J2單元格中輸入公式“=(H22+I22)/16”，回車得到第一個功率譜密度0；下拉至J17，得到全部譜密度數(shù)值（圖14）?；贔FT的譜密度分布是對稱的，可以看出，以J10所在的3105.28為對稱點，上下的數(shù)值完全對稱。圖14 功率譜密度的計算結(jié)果第四步，繪制頻譜圖線頻率fi可以表作，-1顯然f0=0/16=0，f1=1/16=0.0625，f2=2/16=0.125，f15=15/16=0.9375。在Excel

13、中，容易計算頻率的數(shù)值。將頻率與功率譜對應起來（圖15），就可以畫出頻譜圖。如果補上最后一個頻率數(shù)值f16=1及其對應的功率0，則可畫出完全對稱的譜圖（圖16）。圖15 功率譜密度與頻率的對應關(guān)系圖16 對稱分布的頻譜圖由于功率譜圖的對稱性，畫出整個譜圖實在沒有必要，因此，在實際工作中，通常只畫出左半邊（圖17）。圖17 實用的頻譜圖第五步，功率譜分析頻域分析的主要目標之一是判斷時間序列是否存在周期性。從圖17可以看出，功率最大點對應的頻率是f1=0.0625，該頻率對應的周期長度為16?？梢姡跁r間序列較短的情況下之間用功率譜圖尋找時間序列的周期不如周期圖準確。另外可以初步估計數(shù)據(jù)的性質(zhì)。在

14、圖17中，去掉第一個0點，剩余的點一般呈冪指數(shù)分布（在雙對數(shù)坐標圖上點列具有直線趨勢），可以擬合冪指數(shù)函數(shù)如下： （8）圖18 功率P(f)與頻率f的雙對數(shù)坐標圖結(jié)果得到功率譜指數(shù)=1.49521.5。功率譜指數(shù)與時間序列的Hurst指數(shù)具有如下關(guān)系 （9）據(jù)此估計Hurst指數(shù)約為0.25。我們知道，Hurst指數(shù)介于01之間，當H>0.5時，表明時間序列存在正的自相關(guān)，意味著系統(tǒng)演化具有持久性；當H<0時，表明時間序列具有負的自相關(guān)，意味著系統(tǒng)演化具有反持久性；當H=0時，表明時間序列不存在自相關(guān)，過去與未來無關(guān)。對于這條河流的徑流量而言，H=0.25<0.5，表明時間序

15、列具有反持久性：過去的增量意味著今后的減少，過去的減少意味著未來的增加。因此，徑流量必然周期性的變化。【例4】下面對前述例2的數(shù)據(jù)進行Fourier變換，方法與例3相同，但由于N=24，我們?nèi)=32=25。也就是說，對于中心化的數(shù)據(jù)，要在后面添加8個0作為補充點數(shù)。基于FFT的變換結(jié)果如下（見圖19）。圖19 例2數(shù)據(jù)（經(jīng)中心化處理）的FFT變換結(jié)果計算功率譜除例3講述的方法外，還可以利用Excel的另外兩個命令實現(xiàn)：一是計算共軛復數(shù)的命令imconjugate，首先求出的共軛復數(shù)；然后借助復數(shù)的乘積命令improduct，計算復數(shù)的與的乘積；最后利用式（5）得到功率譜。不過，此時的時間序列

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

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

18、動特征，最高峰的時間間隔大約為1012年之間（圖23）。圖23 某海域海平面年平均高度的年際變化曲線 然后將原始數(shù)據(jù)中心化，取T=128=27，這意味著需要將時間序列延長到128位，在計算頻率時用0,1,2,除以128，在計算功率譜密度時，式（5）中的序列長度取T=128。Fourier變換和功率譜密度的計算步驟與例3、例4相同，不贅述。下面直接給出變換結(jié)果（圖24，圖25）。圖24 海面高度時間序列的FFT變換的部分結(jié)果圖25 海面高度變化的頻譜圖在Excel上，將鼠標指向功率譜密度的最大值處，立即顯示頻率為f=0.09375，對應的功率為P(f)=16832.77972。根據(jù)此處的頻率可以

19、算出周期為T=1/f=1/0.09375=10.66711。這表明，海平面高度的變化的確存在一個為期大約11年的變化周期，與直觀判斷結(jié)果一致。由于這里采用的時間序列較長，故結(jié)果比較可靠。太陽黑子的活動周期約為11年稍多一點，二者非常接近。因此，海面高度變化與太陽黑子周期具有對應關(guān)系的猜想，在時間序列變化的節(jié)律方面，大致是可以接受的。圖26 功率譜密度最大值對應的頻率顯示結(jié)果【附錄】我們在前面說過，式(1)實則一個多元線性回歸模型，式(2)是對式(1)中待定參數(shù)的最小二乘(OLS)估計。下面驗證這種判斷。將圖3中計算出的正弦值SIN、余弦值COS和中心化的徑流量Xt集中在一起，經(jīng)復制選擇性粘貼轉(zhuǎn)置，可將數(shù)據(jù)重新排列如下（圖27）。圖27 重新排列的數(shù)據(jù)若以正弦、余弦值為自變量，以中心化的徑流量為因變量進行多元回歸，Excel拒絕給出結(jié)果，并彈出如下對話框顯示拒絕計算的原因。圖28 Excel拒絕給出回