利用Ecel進行時間序列的譜分析

10利用Excel進展時間序列的譜分析〔I〕例子，分別從不同的角度識別時間序列的周期。時間序列的周期圖【1197961980512月份的斷面平均流量。試推斷該河流的徑流量變化是否具有周期性，周期長度大約為多少？x開放為Fourier級數(shù)，則可表示為tx k (at i1

cos2f

tbi

sin2f

t)i

(1)式中fi

為頻率，t為時間序號，k為周期重量的個數(shù)即主周期〔基波〕及其諧波的個數(shù)，εt為標準誤差〔白噪聲序列f給定時，式〔1〕可以視為多元線性回歸模型，可以證ia、b的最小二乘估量為i ia? iNN

xcos2ftt i

〔2〕b? i

2NN

xsin2ftt iN為觀測值的個數(shù)。定義時間序列的周期圖為NI(f)i

2 i

b2)，i1,2, ,k (3)iI(f)ffI(f)為縱軸，繪制時間序列的周期圖，可以在i i i iN=1，=1,2,，f=N，下面借助Exce，利i用上述公式，計算有關參數(shù)并分析時間序列的周期特性。第一步，錄入數(shù)據(jù)，并將數(shù)據(jù)標準化或中心化圖1。1錄入的數(shù)據(jù)及其中心化結果到，中心化的數(shù)據(jù)均值為，但方差與原始數(shù)據(jù)一樣〔未必為。其次步，計算三角函數(shù)值為了借助式〔1〕a、b，首先需要計算正弦值和余弦值。i i取i,2,,6，則頻率為fii/N1/12,2/12,,6/12〔圖。將頻率寫在單元格C3-C14中〔依據(jù)對稱性，我們只用前6個，將中心化的數(shù)據(jù)轉置粘貼于第一行的單元格D1-O1中，月份的序號寫在單元格D2-O2中〔與中心化數(shù)據(jù)對齊。2計算余弦值的表格在D2=COS($B$1*$D$2*C30.86下角右拉至O3f=1/12=0.083，t=1,2,?,12D2單元格中輸=COS($B$1*$D$2*C40.；按住單元格的右下角右拉至O4單元格，得到f=2/12=0.167，t=1,2,?,12〔在D3-O8區(qū)域中。依據(jù)對稱性，我們的計算到=6為止圖。留意，這里B1單元格是2π=6.28319〔圖中未能顯示。在上面的計算中，只要將公式中的“COS”換成“SI在D17-O22圖3計算正弦值的表格a、bi i利用中心化的數(shù)據(jù)〔照舊表作xt〕計算參數(shù)ai、bi。首先算出xtcos2πfit和xtsins2πfit。在D9=D1*D18.30；按住單元格的右下角右拉至O9單元f=1/12=0.083，t=1,2,?,12xtcos2πfit39.5846，即得a1=6.59。在D10=D1*D10.57；按住單元格的右下角右拉至O10f=2/12=0.083，t=1,2,?,12xtcos2πfit值；加和得-365.25，再6a2=-60.875。其余依此類推。將上面公式中的余弦值換成正弦值，即可得到bi值〔見下表。上面的計算過程相當于承受式〔2〕進展逐步計算。第四步，計算頻率強度利用式3，格外簡潔算出f值。例如iNI(f) (a2b2)6*(6.5972170.2132)174096.9141 2 1 1其余依此類推〔圖。4計算頻率強度第五步，繪制時間序列周期圖強率頻10000080000600004000020230000.10.20.3強率頻10000080000600004000020230000.10.20.3頻率0.40.50.6202300180000160000140000度1202305某河流徑流量的周期圖〔1979619805〕第六步，周期識別關鍵是查找頻率的極值點或突變點。在本例中，沒有極值點，但在f1=1/12=0.0833處，頻率強度突然增加〔陡增，而此時T=11=12，故可推斷時間序列可能存在一個12月的周1年周期。【2】為了映證上述推斷，我們借助同一條河流的連續(xù)兩年的平均月徑流量〔196161963年5月。原始數(shù)據(jù)見以以以下圖圖6。6原始數(shù)據(jù)及局部處理結果將原始數(shù)據(jù)回車時間序列變化圖，可以初步估量具有12月變化周期，但不能確定〔圖6。350350300250量流徑均平月200150100500051015202530時間〔月份序號〕6徑流量的月變化圖〔1961619635〕依據(jù)例1給出的計算步驟，計算參數(shù)數(shù)ab，進而計算頻率強度〔圖。然后i i繪制時間序列的周期圖圖8N=2，我們取=1。7參數(shù)和頻率強度的計算結果f8中可以看出，頻率強度的最大值〔極值點〕對應于頻率=1/12=0.0833，故時間f11T=1/f=12122的111倍，故推斷結果更為牢靠。140000140000120230100000度強率頻80000600004000020230000.10.20.30.40.50.6頻率8某河流徑流量的周期圖〔1961619635〕2時間序列的頻譜圖【3】首先考慮對例1的數(shù)據(jù)進展功率譜分析。例1的時間序列較短，分析的效果不佳，但計算過程簡短。給出這個例子，主要是幫助大家理解Fourier變換過程和方法。為了進展Fourier分析，需要對數(shù)據(jù)進展預處理。第一，將數(shù)據(jù)中心化，即用原始數(shù)據(jù)減去其平均值。中心化的數(shù)據(jù)均值為0，我們對中心化的數(shù)據(jù)進展變換，其周期更為明顯。其次，由于Fourier分析通常承受所謂快速Fourier變換FastFourierTransformation，F(xiàn)F，而FFT要求數(shù)據(jù)必需為2n個，這里n為正整數(shù)1,2,3?，而我們的樣本為N=1，它不是2的某n次方。因此，在中心化的數(shù)據(jù)后面加上40N′=12＋4＝16＝24個，這才符合FFT的需要圖。下面，我們對延長后的中心化數(shù)據(jù)進展Fourier變換。9數(shù)據(jù)的中心化與“延長”第一步，翻開Foureir分析對話框沿著主菜單的“工具Tool”→“數(shù)據(jù)分析DataAnalysi框〔圖10，從中選擇“傅立葉分析FourierAnalysi。10在數(shù)據(jù)分析選項框中選擇Fourier其次步，定義變量和輸出區(qū)域確定之后，彈出傅立葉分析對話框，依據(jù)數(shù)據(jù)在工作表中的分布狀況進展如下設置：將光標置于“輸入?yún)^(qū)域”對應的空白欄，然后用鼠標選中單元格C1-C17，這時空白欄中自動以確定單元格的形式定義中心化數(shù)據(jù)的區(qū)域范圍〔$C$1-$C$1選中“標志位于第一行選中輸出區(qū)域，定義范圍為D2-D1圖1。留意：假設輸入?yún)^(qū)域的數(shù)據(jù)范圍定義為C2-C1，則不要選中“標志位于第一行圖1自動給在的工作表組上。這一點也與回歸分析一樣。11選中“標志位于第一行”與數(shù)據(jù)輸入范圍的定義12不選中“標志位于第一行”與數(shù)據(jù)輸入范圍的定義第三步，結果轉換定義數(shù)據(jù)輸入－輸出區(qū)域完成之后，確定，馬上得到Fourier變換的結果〔圖13。13傅立葉變換的結果f(t)F(ω)xX(f)。我們知t Tf(t)F(ω)S(ω)，即有S()F()F()F()2X(f)也就簡潔計算功率譜〔密度〕T

(4)X (f)2P(X (f)2T

(5)式中f表示線頻率，與角頻率ω的轉換關系是=T，這里T假設將X表作X)=+j〔這里AB為虛部，則有T TX (fT i

)2XT

(f)Xi

(f)(Ai

jBi

)(Ai

jBi

i

B2i

〔6〕因此這一步是要分別變換結果的實部與虛部。逐個手動提取是格外麻煩而且簡潔出錯Excel有關復數(shù)計算的命令。提取實部的ExcelimrealH2單元格=IMREAL(D2〔這里D2為變換結果的第一個復數(shù)所在的單元格0；點中H2單元格的右下角，撳住鼠標左鍵下拉至H17，得到全部的實部數(shù)值。提取虛部的命令是imaginar。在I2=IMAGINARY(D2；下拉至I1，得到全部的虛部數(shù)值。依據(jù)式5〔6密度的計算公式為P(f)i

i iT

〔7〕考慮到本例中N=1，在J2=(H2^2+I2^2)/1譜密度；下拉至J1，得到全部譜密度數(shù)值圖1?；贔FT可以看出，以J103105.28為對稱點，上下的數(shù)值完全對稱。14功率譜密度的計算結果第四步，繪制頻譜圖f可以表作i

f i/Ti/N，i0,1,2, ,N-1if f f 明顯=0/16=0，=1/16=0.0625，=2/16=0.125，?，=15/16=0.9375Excelf f f 0 1 2 15計算頻率的數(shù)值。將頻率與功率譜對應起來圖1，就可以畫出頻譜圖。假設補上最終一個頻率數(shù)值f16=1及其對應的功率，則可畫出完全對稱的譜圖圖1。15功率譜密度與頻率的對應關系700007000060000)50000f(P度密譜率功40000300002023010000000.20.40.60.81頻率f16對稱分布的頻譜圖出左半邊圖1。70000700006000050000)f(40000P率功300002023010000000.10.20.30.40.50.6頻率f17有用的頻譜圖第五步，功率譜分析頻域分析的主要目標之一是推斷時間序列是否存在周期性。從圖17可以看出，功率最大點對應的頻率是f1=0.0625，該頻率對應的周期長度為16?？梢姡跁r間序列較短的狀況下之間用功率譜圖查找時間序列的周期不如周期圖準確。另外可以初步估量數(shù)據(jù)的性質。在170點，剩余的點一般呈冪指數(shù)分布〔在雙對數(shù)坐標圖上點列具有直線趨勢，可以擬合冪指數(shù)函數(shù)如下：P(f) f 〔8〕100000100000)f(P率功10000P(f)=935.68f-1.4952R2=0.923310000.010.11頻率f18功率P(ff結果得到功率譜指數(shù)=1.495≈1.Hurst指數(shù)具有如下關系2H1 〔9〕據(jù)此估量Hurst0.25。我們知道，Hurst0～1H>0.5時，說明時間序列存在正的自相關，意味著系統(tǒng)演化具有長期性；當H<0時，說明時間序列具有負的自相關，意味著系統(tǒng)演化具有反長期性；當H=0時，說明時間序列不存在自相關，過去與將來無關。對于這條河流的徑流量而言，H=0.25<0.5，說明時間序列具有反長期性：過去的增量意味著今后的削減，過去的削減意味著將來的增加。因此，徑流量必定周期性的變化。42的數(shù)據(jù)進展Fourier3N=24，我們取T=32=258個0作為補充點數(shù)?；贔FT的變換結果如下〔圖19。19例2〔經(jīng)中心化處理〕FFT3表達的方法外，還可以利用Excel的另外兩個命令實現(xiàn)：一是計算共軛復數(shù)的命令imconjugateXT

(f)的共軛復數(shù)X (f)；然后借助復數(shù)的乘積命T令improduct，計算復數(shù)的X (f)與X (f)的乘積X (f)2；最終利用式〔5〕得到功率T T T譜。不過，此時的時間序列長度視為T=32。例如在H2=IMCONJUGTE(D2得到全部共軛數(shù)值。在L2=IMPRODUCT(D2,H2L33XT

(f)2XT

(f)2除以32圖20。依據(jù)計算結果畫出頻譜圖圖210.09375，相應的周期為T=10.667；在極值點附件存在一個次最大點，但相對于其他數(shù)值卻明顯又是突變點，該點對應的頻率為0.0625，相應的周期為16。故可斷定，該時間序列的10～16之間。20功率譜密度值及其相應的頻率35000350003000025000)f(20230P率功15000100005000000.10.20.30.40.50.6頻率f214給一個較長時間序列的分析實例。【例511年左右海平面到達一個最高值〔圖22，于是科學家推斷海平面的升降存在一個11年周期，與太陽黑子〔sunspot〕的11年周期變化全都，而太陽黑子的活動正是海平面升降的緣由所在。問題在于，前述11年周期是通過原始數(shù)據(jù)的變化曲線直觀覺察的，未必牢靠。為此，可以進展一個功率譜分析，推斷這種周期是否確實存在。22某海疆海面年平均高度的時間序列數(shù)據(jù)180160140120度高面海100806040200020406080180160140120度高面海100806040200020406080100年度序號23某海疆海平面年平均高度的年際變化曲線然后將原始數(shù)據(jù)中心化，取T=128=27128位，在計0,1,2,128〔5〕T=128Fourier34一樣，不贅述。下面直接給出變換結果〔圖24，圖25。24海面高度時間序列的FFT1800018000160001400012023)f(10000P率功8000600040002023000.10.20.30.40.50.6頻率f25海面高度變化的頻譜圖在Excelf=0.09375，對應的功率為P(f)=16832.77972T=1/f=1/0.09375=10.667≈1111年的變化周期，與直觀推斷結果全都。11年稍多一點，的節(jié)律方面，大致是可以承受的。26功率譜密度最大值對應的頻率顯示結果【附錄】我們在前面說過，式(1)實則一個多元線性回歸模型，式(2)是對式(1)中待定參數(shù)的最小二乘(OLS)3中計算出的正弦值SIN、余弦值COS和中心化的徑流量Xt〔圖227重排列的數(shù)據(jù)假設以正弦、余弦值為自變量，以中心化的徑流量為因變量進展多元回歸，Excel拒絕給出結果，并彈出如下對話框顯示拒絕計算的緣由。28Excel問題在于行數(shù)與列數(shù)不能一樣，而本例中自變量和樣本數(shù)都是12，這就是問題所在?？紤]到最終一個變量全都是0，不妨去掉最終一個變量。由于式〔1〕不含常數(shù)項，在回歸常數(shù)為Excel給出的多元回歸結果如下圖2。29利用圖2712將回歸結果與利用式〔2〕直接結算的結果進展比較，覺