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

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

cos2f

tbi

sin2f

t)i

(1)式中fi

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

xcos2ftt i

〔2〕b? i

2NN

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

2 i

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

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

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

)2XT

(f)Xi

(f)(Ai

jBi

)(Ai

jBi

i

B2i

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

i iT

〔7〕考慮到本例中N=1，在J2=(H2^2+I2^2)/1譜密度；下拉至J1，得到全部譜密度數(shù)值圖1?；贔FT可以看出，以J103105.28為對(duì)稱點(diǎn)，上下的數(shù)值完全對(duì)稱。14功率譜密度的計(jì)算結(jié)果第四步，繪制頻譜圖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計(jì)算頻率的數(shù)值。將頻率與功率譜對(duì)應(yīng)起來圖1，就可以畫出頻譜圖。假設(shè)補(bǔ)上最終一個(gè)頻率數(shù)值f16=1及其對(duì)應(yīng)的功率，則可畫出完全對(duì)稱的譜圖圖1。15功率譜密度與頻率的對(duì)應(yīng)關(guān)系700007000060000)50000f(P度密譜率功40000300002023010000000.20.40.60.81頻率f16對(duì)稱分布的頻譜圖出左半邊圖1。70000700006000050000)f(40000P率功300002023010000000.10.20.30.40.50.6頻率f17有用的頻譜圖第五步，功率譜分析頻域分析的主要目標(biāo)之一是推斷時(shí)間序列是否存在周期性。從圖17可以看出，功率最大點(diǎn)對(duì)應(yīng)的頻率是f1=0.0625，該頻率對(duì)應(yīng)的周期長度為16?？梢姡跁r(shí)間序列較短的狀況下之間用功率譜圖查找時(shí)間序列的周期不如周期圖準(zhǔn)確。另外可以初步估量數(shù)據(jù)的性質(zhì)。在170點(diǎn)，剩余的點(diǎn)一般呈冪指數(shù)分布〔在雙對(duì)數(shù)坐標(biāo)圖上點(diǎn)列具有直線趨勢(shì)，可以擬合冪指數(shù)函數(shù)如下：P(f) f 〔8〕100000100000)f(P率功10000P(f)=935.68f-1.4952R2=0.923310000.010.11頻率f18功率P(ff結(jié)果得到功率譜指數(shù)=1.495≈1.Hurst指數(shù)具有如下關(guān)系2H1 〔9〕據(jù)此估量Hurst0.25。我們知道，Hurst0～1H>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í)間序列具有反長期性：過去的增量意味著今后的削減，過去的削減意味著將來的增加。因此，徑流量必定周期性的變化。42的數(shù)據(jù)進(jìn)展Fourier3N=24，我們?nèi)=32=258個(gè)0作為補(bǔ)充點(diǎn)數(shù)?；贔FT的變換結(jié)果如下〔圖19。19例2〔經(jīng)中心化處理〕FFT3表達(dá)的方法外，還可以利用Excel的另外兩個(gè)命令實(shí)現(xiàn)：一是計(jì)算共軛復(fù)數(shù)的命令imconjugateXT

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

(f)2XT

(f)2除以32圖20。依據(jù)計(jì)算結(jié)果畫出頻譜圖圖210.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)的頻率35000350003000025000)f(20230P率功15000100005000000.10.20.30.40.50.6頻率f214給一個(gè)較長時(shí)間序列的分析實(shí)例?！纠?11年左右海平面到達(dá)一個(gè)最高值〔圖22，于是科學(xué)家推斷海平面的升降存在一個(gè)11年周期，與太陽黑子〔sunspot〕的11年周期變化全都，而太陽黑子的活動(dòng)正是海平面升降的緣由所在。問題在于，前述11年周期是通過原始數(shù)據(jù)的變化曲線直觀覺察的，未必牢靠。為此，可以進(jìn)展一個(gè)功率譜分析，推斷這種周期是否確實(shí)存在。22某海疆海面年平均高度的時(shí)間序列數(shù)據(jù)180160140120度高面海100806040200020406080180160140120度高面海100806040200020406080100年度序號(hào)23某海疆海平面年平均高度的年際變化曲線然后將原始數(shù)據(jù)中心化，取T=128=27128位，在計(jì)0,1,2,128〔5〕T=128Fourier34一樣，不贅述。下面直接給出變換結(jié)果〔圖24，圖25。24海面高度時(shí)間序列的FFT1800018000160001400012023)f(10000P率功8000600040002023000.10.20.30.40.50.6頻率f25海面高度變化的頻譜圖在Excelf=0.09375，對(duì)應(yīng)的功率為P(f)=16832.77972T=1/f=1/0.09375=10.667≈1111年的變化周期，與直觀推斷結(jié)果全都。11年稍多一點(diǎn)，的節(jié)律方面，大致是可以承受的。26功率譜密度最大值對(duì)應(yīng)的頻率顯示結(jié)果【附錄】我們?cè)谇懊嬲f過，式(1)實(shí)則一個(gè)多元線性回歸模型，式(2)是對(duì)式(1)中待定參數(shù)的最小二乘(OLS)3中計(jì)算出的正弦值SIN、余弦值COS和中心化的徑流量Xt〔圖227重排列的數(shù)據(jù)假設(shè)以正弦、余弦值為自變量，以中心化的徑流量為因變量進(jìn)展多元回歸，Excel拒絕給出結(jié)果，并彈出如下對(duì)話框顯示拒絕計(jì)算的緣由。28Excel問題在于行數(shù)與列數(shù)不能一樣，而本例中自變量和樣本數(shù)都是12，這就是問題所在。考慮到最終一個(gè)變量全都是0，不妨去掉最終一個(gè)變量。由于式〔1〕不含常數(shù)項(xiàng)，在回歸常數(shù)為Excel給出的多元回歸結(jié)果如下圖2。29利用圖2712將回歸結(jié)果與利用式〔2〕直接結(jié)算的結(jié)果進(jìn)展比較，覺