利用Mathcad進(jìn)行時間序列的譜分析(I)

1、利用Mathcad進(jìn)行時間序列的譜分析（I）在這一節(jié)中我們采用與“利用Excel進(jìn)行時間序列的譜分析（I）”一節(jié)相同的例子，以便計算結(jié)果能夠相互參照。關(guān)于Mathcad的基本操作技巧（包括繪圖）已在“回歸分析”中備述，本節(jié)不再詳細(xì)說明?！纠?】某水文觀測站測得一條河流從1979年6月到1980年5月共計12月份的斷面平均流量。試借助功率譜分析判斷時間序列是否存在周期性特征，并對周期長度進(jìn)行估計。第一步，錄入或拷入數(shù)據(jù) 可以將在Excel中經(jīng)過中心化的數(shù)據(jù)直接錄入或拷入Mathcad的數(shù)據(jù)表中，并定義為“data”。 圖1 錄入的中心化數(shù)據(jù)需要說明兩點：第一，如果是從Excel中拷入數(shù)據(jù)，一定記

2、住檢查最后一行數(shù)據(jù)，如果多出一行0數(shù)據(jù)點，必須將其刪除。第二，也可以在Mathcad中對數(shù)據(jù)中心化。中心化的方法非常簡單，將原始數(shù)據(jù)拷入Mathcad的數(shù)據(jù)表并將其定義為“data”（圖2）。圖2 原始數(shù)據(jù)表的局部顯示（數(shù)據(jù)未經(jīng)中心化）在被定義過的表下定義變量為然后利用求均值命令mean，鍵入公式回車，立即得到中心化的結(jié)果（圖3）。 圖3 中心化的數(shù)據(jù)不過，中心化以后需要重新插入一個數(shù)據(jù)表并重新定義該表（例如定義為biao）。將中心化的數(shù)據(jù)拷入“biao”中，然后在左方插入單元格（Insert Cells），再在新插入的一列單元格中錄入數(shù)據(jù)序號（圖4）。從這個意義上講，在Mathcad中進(jìn)行數(shù)

3、據(jù)處理不如直接將經(jīng)過Excel預(yù)處理即中心化的數(shù)據(jù)拷貝過來。圖4 新建的數(shù)據(jù)表第二步，繪制數(shù)據(jù)序列的變化圖假定錄入的是圖1所示的中心化的數(shù)據(jù)，將表定義為data以后，在數(shù)據(jù)表下對變量進(jìn)行如下定義：然后利用Graph中的曲線圖工具容易繪制曲線圖（圖5）。 圖5 中心化時間序列的變化曲線（1979年6月1980年5月）繪圖的主要目的有兩個：其一，可以通過曲線圖直觀地判斷數(shù)據(jù)變化是否存在周期，對于本例，當(dāng)然達(dá)不到這種效果；其二，如果拷貝的數(shù)據(jù)多出0數(shù)據(jù)行，可以及時將其刪除，以免計算出錯。第三步，進(jìn)行快速Fourier變換首先是定義變量。在Mathcad中，系統(tǒng)默認(rèn)的數(shù)據(jù)排列方式是從0到K。我們有12

4、個月份的數(shù)據(jù)，也就相當(dāng)于擁有12個樣本，即N=12，從而K=12-1=11。但是，F(xiàn)ourier變換要求時間序列的長度必須是T=2n個（n為正整數(shù)），因此，可以考慮在序列末尾加0，將其延長到T=24個。方法是定義序號k:=12.15，并設(shè)xk=0，這相當(dāng)于在Excel加上4個0元素。這種補(bǔ)充定義已在上文給出。定義快速Fourier變換，命令為fft（注意：在Mathcad中，大寫的FFT與小寫的fft變換結(jié)果相差T倍，為了與Excel的結(jié)果對應(yīng)，本例采用fft）。我們是要對x進(jìn)行fft，令變換的結(jié)果為Y，Y將給出對x進(jìn)行快速Fourier變換的結(jié)果；定義M為最后一個變換結(jié)果的序號，則可以通過l

5、ast自動從對稱點(參見“利用Excel進(jìn)行時間序列的譜分析（I）”一節(jié)的圖15、圖16)截取前半部結(jié)果，M=N/2=8；定義功率譜為Poweri，則有注意：如果使用命令FFT，則上式的分母應(yīng)為T2。定義頻率變化的步長為w=1/2；則頻率為freqj=j/(2M)=j/T，并取j=0,1,2,M。全部定義的表達(dá)式如下：定義完成以后，鍵入“Yi=”，回車，得到x的FFT結(jié)果；鍵入“Poweri=”或“Yi2”，回車，得到功率譜密度（圖6）。 圖6 Fourier變換結(jié)果與功率譜密度第四步，繪制頻譜圖，識別周期利用Mathcad的Graph工具箱容易繪制頻譜圖（圖7）。在圖上點右鍵，出現(xiàn)一個選擇菜

6、單（圖8）；選中Trace，彈出X-Y Trace對話框（圖9）；用鼠標(biāo)點擊譜線的最大值點，則Trace自動給出譜密度的極大值P(f)max=63157，對應(yīng)的頻率為f1=0.0625；點擊Copy X按鈕，將0.0625復(fù)制到Mathcad的工作表上，取倒數(shù)，立即到周期P=1/f1=16。當(dāng)然，由于時間序列太短，這個周期是不準(zhǔn)確的。圖7 例1的頻譜圖圖8 右鍵菜單與Trace的位置 圖9 借助Trace顯示功率譜密度最大值對應(yīng)的頻率【例2】對同一條河流的連續(xù)兩年的平均月徑流量（1961年6月1963年5月）進(jìn)行Fourier分析，并識別其周期。步驟與例1相同，不詳述。在錄入中心化數(shù)據(jù)以后（圖

7、10），需要重新定義數(shù)據(jù)的范圍。圖10 例2的數(shù)據(jù)表（局部，已中心化）考慮到這次N=24，可以將其延長到T=25=32。于是Fourier變換定義如下：畫出時間序列的變化曲線圖，從圖上可以直觀地判斷徑流量變化具有某種周期性，但不能確定（圖11）。圖11 例2的時間序列變化圖（1961年6月1963年5月） 利用Fourier變換結(jié)果（圖12）畫出頻譜圖（圖13）。借助Trace功能容易找到最大功率譜密度P(f)max=63157及其對應(yīng)的頻率f1=0.09375（圖14）；點擊Copy X按鈕，將0.09375復(fù)制到Mathcad的工作表上，取倒數(shù)，可以得到周期P=1/f1=10.667。在最

8、大值附近存在一個與最大值非常接近的功率譜密度值為31107，對應(yīng)的頻率為0.0625。由于這個次最大點同樣代表一個突變點它與臨近值有很大差距，故須考慮這個數(shù)值。根據(jù)上例可知，其倒數(shù)為16。由此判斷，時間序列的周期大約在1016之間。我們知道，河流徑流量的周期為12個月。可見，對于較短的時間序列，功率譜分析會有較大的誤差。 圖12 例2的Fourier變換結(jié)果 圖13 例2的頻譜圖圖14 例2的功率譜密度最大值及其對應(yīng)的頻率【例3】某海域海面平均高度年際變化的功率譜分析。本例時間序列長度為100年，即N=100。將時間序列延長到T=27=128。圖15 例3數(shù)據(jù)的局部顯示（已經(jīng)中心化）錄入數(shù)據(jù)以

9、后，根據(jù)時間序列長度對變量和Fourier變換定義如下：從原始數(shù)據(jù)的圖像可以看到，時間序列可能具有某種內(nèi)在節(jié)律（圖16）。圖16 海面高度的年際變化曲線借助Fourier變換的結(jié)果，作出頻譜圖如下（圖17）。 圖17 海面高度變化的頻譜圖借助Mathcad的Trace功能發(fā)現(xiàn)功率譜密度最大值對應(yīng)的頻率為f=0.09375，于是得到周期約為P=1/f=10.667（年）。圖18 最大功率譜密度值對應(yīng)的頻率【例4】 在例2中，利用時間序列的自相關(guān)函數(shù)分析其周期特性。在“利用Excel進(jìn)行時間序列的譜分析（I）”一節(jié)，我們談到可以定義時間序列的周期圖為，式中I(fi)為頻率fi處的強(qiáng)度。以fi為橫軸

10、，以I(fi)為縱軸，繪制時間序列的周期圖，可以在最大值處找到時間序列的周期。實際上，周期圖還可以表作這里R為時間序列的自協(xié)方差函數(shù)，即有這意味著，周期圖是對功率譜直接估計的結(jié)果。根據(jù)上面幾個式子及其內(nèi)在關(guān)系，我們可以利用相關(guān)函數(shù)估計時間序列的周期。我們知道，計算相關(guān)系數(shù)的公式有兩種，本例采用下面公式采用這個公式的原因之一是：我們的時間序列較短，而這個公式的有效性較好；原因之二是可以直接借助SPSS給出結(jié)果。將結(jié)果復(fù)制到Excel中間，轉(zhuǎn)換格式以后，再復(fù)制到Mathcad的數(shù)據(jù)表（圖19），然后進(jìn)行Fourier變換。圖19 徑流量的自相關(guān)系數(shù)進(jìn)行快速Fourier運算的有關(guān)定義如下：根據(jù)計算結(jié)果容易得到頻譜圖（圖20）。利用Trace可知，最大功率譜密度對應(yīng)的頻率為f=0.09375，于是得到周期約為P=1/f=10.667（月），與直接采用中心化數(shù)據(jù)進(jìn)行Fourier變換得到的結(jié)果相同（參見例2）。 圖20 基于徑流量自相關(guān)系數(shù)的頻譜圖前面說過，