精品資料（2021-2022年收藏）利用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ù)，則可表示為 (1)式中fi為頻率，t為時(shí)間序號(hào)，k為周期分量的個(gè)數(shù)即主周期（基波）及其諧波的個(gè)數(shù)，t為標(biāo)準(zhǔn)誤差（白噪聲序列）。當(dāng)頻率fi給定時(shí)，式（1）可以視為多元線性回歸模型，可以證明，待定系數(shù)ai、bi的最小二乘估計(jì)

2、為 （2）這里N為觀測(cè)值的個(gè)數(shù)。定義時(shí)間序列的周期圖為， (3)式中I(fi)為頻率fi處的強(qiáng)度。以fi為橫軸，以I(fi)為縱軸，繪制時(shí)間序列的周期圖，可以在最大值處找到時(shí)間序列的周期。對(duì)于本例，N=12，t=1,2,N，fi=i/N，下面借助Excel，利用上述公式，計(jì)算有關(guān)參數(shù)并分析時(shí)間序列的周期特性。第一步，錄入數(shù)據(jù)，并將數(shù)據(jù)標(biāo)準(zhǔn)化或中心化（圖1）。圖1 錄入的數(shù)據(jù)及其中心化結(jié)果中心化與標(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ì)算正弦

3、值和余弦值。取，則頻率為（圖1）。將頻率寫在單元格C3-C14中（根據(jù)對(duì)稱性，我們只用前6個(gè)），將中心化的數(shù)據(jù)轉(zhuǎn)置粘貼于第一行的單元格D1-O1中，月份的序號(hào)寫在單元格D2-O2中（與中心化數(shù)據(jù)對(duì)齊）。圖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í)的全部余弦值。依次類推，可以算

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

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

6、期圖（圖5）。圖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年周期。【例2】為了映證上述判斷，我們借助同一條河流的連續(xù)兩年的平均月徑流量（1961年6月1963年5月）。原始數(shù)據(jù)見下圖（圖6）。圖6 原始數(shù)據(jù)及部分處理結(jié)果將原始數(shù)據(jù)回車時(shí)間序列變化圖，可以初步估計(jì)具有12月變化周期，但不能肯定（圖6）。圖6 徑流量的月變化圖（1961年6月1963年5月）按照例1給出的計(jì)算步驟，計(jì)

7、算參數(shù)數(shù)ai、bi，進(jìn)而計(jì)算頻率強(qiáng)度（結(jié)果將圖7）。然后繪制時(shí)間序列的周期圖（圖8）。注意這里，N=24，我們?nèi)=12。圖7 參數(shù)和頻率強(qiáng)度的計(jì)算結(jié)果從圖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é)果更為可靠。圖8 某河流徑流量的周期圖（1961年6月1963年5月）2 時(shí)間序列的頻譜圖【例3】首先考慮對(duì)例1的數(shù)據(jù)進(jìn)行功率譜分析。例1的時(shí)間序列較短，分析的效果不佳，但計(jì)算過程簡(jiǎn)短。給出這個(gè)例子，主要是幫助大家理解Fouri

8、er變換過程和方法。為了進(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，F(xiàn)FT），而FFT要求數(shù)據(jù)必須為2n個(gè)，這里n為正整數(shù)（1,2,3,），而我們的樣本為N=12，它不是2的某個(gè)n次方。因此，在中心化的數(shù)據(jù)后面加上4個(gè)0，這樣新的樣本數(shù)為N=1241624個(gè)，這才符合FFT的需要（圖9）。下面，我們對(duì)延長(zhǎng)后的中心化數(shù)據(jù)進(jìn)行Fourier變換。圖9 數(shù)據(jù)的中心化與“

9、延長(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-$C$17）。選中“標(biāo)志位于第一行”。選中輸出區(qū)域，定義范圍為D2-D17（圖11）。注意：如果輸

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

11、XT(f)也就容易計(jì)算功率譜（密度） (5)式中f表示線頻率，與角頻率的轉(zhuǎn)換關(guān)系是=2/T，這里T為數(shù)據(jù)區(qū)間長(zhǎng)度。如果將XT(f)表作XT(f)=A+jB（這里A為實(shí)部，B為虛部），則有 （6）因此這一步是要分離變換結(jié)果的實(shí)部與虛部。逐個(gè)手動(dòng)提取是非常麻煩而且容易出錯(cuò)的，可以利用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單元格中輸入公

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

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

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

15、列具有反持久性：過去的增量意味著今后的減少，過去的減少意味著未來的增加。因此，徑流量必然周期性的變化?！纠?】下面對(duì)前述例2的數(shù)據(jù)進(jìn)行Fourier變換，方法與例3相同，但由于N=24，我們?nèi)=32=25。也就是說，對(duì)于中心化的數(shù)據(jù)，要在后面添加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，首先求出的共軛復(fù)數(shù)；然后借助復(fù)數(shù)的乘積命令improduct，計(jì)算復(fù)數(shù)的與的乘積；最后利用式（5）得到功率譜。不過，此時(shí)的時(shí)間序列

16、長(zhǎng)度視為T=32。例如在H2中鍵入公式“=IMCONJUGATE(D2)”，回車得到第一個(gè)共軛數(shù)，下拉至H33，得到全部共軛數(shù)值。在L2中鍵入公式“=IMPRODUCT(D2,H2)”，回車，得到第一個(gè)復(fù)數(shù)乘積，下拉至L33，得到全部值。最后用除以32得到功率譜密度（圖20）。根據(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í)間序列的周期比在1016之間。圖20 功率譜密度值及其相應(yīng)的頻率圖21

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

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

19、算出周期為T=1/f=1/0.09375=10.66711。這表明，海平面高度的變化的確存在一個(gè)為期大約11年的變化周期，與直觀判斷結(jié)果一致。由于這里采用的時(shí)間序列較長(zhǎng)，故結(jié)果比較可靠。太陽黑子的活動(dòng)周期約為11年稍多一點(diǎn)，二者非常接近。因此，海面高度變化與太陽黑子周期具有對(duì)應(yīng)關(guān)系的猜想，在時(shí)間序列變化的節(jié)律方面，大致是可以接受的。圖26 功率譜密度最大值對(duì)應(yīng)的頻率顯示結(jié)果【附錄】我們?cè)谇懊嬲f過，式(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é)