《利用Mathcad进行时间序列的谱分析Read.docx》由会员分享，可在线阅读，更多相关《利用Mathcad进行时间序列的谱分析Read.docx（14页珍藏版）》请在三一文库上搜索。
1、利用Mathcad进行时间序列谱分析（I）在这一节中我们采用及“利用Excel进行时间序列谱分析（I）” 一节 相同例子，以便计算结果能够相互参照。关于Mathcad基本操作技巧（包 括绘图）己在“回归分析”中备述，本节不再详细说明。【例1】某水文观测站测得一条河流从1979年6月到1980年5月共计12 月份断面平均流量。试借助功率谱分析判断时间序列是否存在周期性特 征，并对周期长度进行估计。第一步，录入或拷入数据可以将在Excel中经过中心化数据直接录入或拷入Mathcad数据表 中，并定义为“data”。data :=010121.141279.1423384.1434117.14453
2、5.1456-62.8667-103.4678-121.9689-135.06910-117.061011-90.461112-5.86图1录入中心化数据需要说明两点：第一，如果是从Excel中拷入数据，一定记住检查最 后一行数据，如果多出一行0数据点，必须将其删除。第二，也可以在Mathcad中对数据中心化。中心化方法非常简单，将原始数据拷入Mathcad数据表并将其定义为“data”（图2）。data :=010119012248图2原始数据表局部显示（数据未经中心化）在被定义过表下定义变量为t = da$然后利用求均值命令mean,键入公式x - mean（x）=回车，立即得到中心化结果
3、（图3）。01021.142179.1422384.1423117.142435.1425-62.8586-103.4587-121.9588-135.0589-117.05810-90.45811-5.858图3中心化数据不过，中心化以后需要重新插入一个数据表并重新定义该表（例如定义为 biao）o将中心化数据拷入biao”中，然后在左方插入单元格（Insert Cells）,再在新插入一列单元格中录入数据序号（图4）。从这个意义上讲,在Mathcad中进行数据处理不如直接将经过Excel预处理即中心化数据拷贝过来。biao :=010121.141279.14图4新建数据表第二步，绘制数据
4、序列变化图假定录入是图1所示中心化数据，将表定义为data以后，在数据表下对变量进行如下定义：t := datax:= datak := 12. 15:= C然后利用Graph中曲线图工具容易绘制曲线图（图5）。 024681c)121I12 T一时序384.142/-135.O58,9Q0图5中心化时间序列变化曲线（1979年6月1980年5月）绘图主要目有两个：其一，可以通过曲线图直观地判断数据变化是否 存在周期，对于本例，当然达不到这种效果；其二，如果拷贝数据多出0 数据行，可以及时将其删除，以免计算出错。第三步，进行快速Fourier变换首先是定义变量。在Mathcad中，系统默认数据
5、排列方式是从0到七 我们有12个月份数据，也就相当于拥有12个样本，即,312,从而A=12-l=llo但是，Fourier变换要求时间序列长度必须是及才个为正 整数），因此，可以考虑在序列末尾加0,将其延长到我2,个。方法是定义 序号4:二12.15,并设照二0,这相当于在Excel加上4个0元素。这种补 充定义已在上文给出。定义快速Fourier变换，命令为fft （注意：在Mathcad中，大写FFT 及小写fft变换结果相差T倍，为了及Excel结果对应，本例采用fft）。 我们是要对x进行fft,令变换结果为匕V将给出对x进行快速Fourier 变换结果；定义为最后一个变换结果序号，
6、则可以通过last自动从对 称点（参见“利用Excel进行时间序列谱分析” 一节图15、图16）截 取前半部结果，加M2=8；定义功率谱为Power；,则有IrJ2Poweg =号注意：如果使用命令FFT,则上式分母应为fo定义频率变化步长为产1/2； 则频率为freq尸力协二力T,并取上0, 1, 2,M全部定义表达式如下：Y:=fft（x）M:=last（Y） M = 8i:=0.MPowerj=（闾w = ；j：=O.M freqj:=w-j定义完成以后，键入“看”，回车，得到aFFT结果；键入Power=或“ |加2=”,回车,得到功率谱密度（图6）。Y.=2.943?0 书2.943
7、?0 T56.3167046.3167042.6497042.649?。41.0377041.03770 45.0370 35.0370 34.62970 34.62970 32.87870 32.87870 35.34970 35.34970 33.1057033.105703i5.425?0 -8232.503+95.3881-74.087+144.9261-75.626+68.2181-67.25-22.525i-50.648-45.43i-0.442-53.645i49.971 -53.401 i55.725图6 Fourier变换结果及功率谱密度第四步，绘制频谱图，识别周期利用Mat
8、hcadGraph工具箱容易绘制频谱图（图7）。在图上点右键， 出现一个选择菜单（图8）；选中Trace,弹出XT Trace对话框（图9）； 用鼠标点击谱线最大值点，则Trace自动给出谱密度极大值户（力陶二63157, 对应频率为=0. 0625；点击Copy X按钮，将0. 0625复制到Mathcad I 作表上，取倒数，立即到周期 片1/十二16。当然，由于时间序列太短，这 个周期是不准确。/ I。4 11116.31610.6 io4 A-00.10.2Pouerjj.1042104一 I s2.94I10 ;lie看”频率f图7例1频谱图.6.316x10；Dteabls Eva
9、luationFormat.Trace.Zoom.Powers.IQ4再匕 cutcop一pas1 *亳B图8右键菜单及Trace位置图9借助Trace显示功率谱密度最大值对应频率【例2】对同一条河流连续两年平均月径流量（1961年6月1963年5月）进行Fourier分析，并识别其周期。步骤及例1相同，不详述。在录入中心化数据以后（图10）,需要重新 定义数据范围。data :=010139.5512168.55图10例2数据表（局部，已中心化）考虑到这次毕24,可以将其延长到伫25=32。于是Fourier变换定义如下:t :=data x:=data lzk:=24.31=0Y :=ff
10、t（x）M :=last（Y） M = 16 i :=0.M2jPower. :=（ |Yj|w:=- j :=0.M freqj := w - j -画出时间序列变化曲线图，从图上可以直观地判断径流量变化具有某种周 期性，但不能确定（图11）。05101520一I时序213.546 300 1000002 1X 由黑空35%0a2524.图11例2时间序列变化图（1961年6月一1963年5月）利用Fourier变换结果（图12）画出频谱图（图13）。借助Trace功 能容易找到最大功率谱密度?（力皿二63157及其对应频率=0.09375 （图 14）；点击Copy X按钮,将0.0937
11、5复制到Mathcad工作表上,取倒数， 可以得到周期尸1/=10.667o在最大值附近存在一个及最大值非常接近功率谱密度值为31107,对应频率为0. 0625o由于这个次最大点同样代表一个突变点一一它及临近值 有很大差距，故须考虑这个数值。根据上例可知，其倒数为16o由此判断, 时间序列周期大约在1016之间。我们知道，河流径流量周期为12个月。可见，对于较短时间序列,功率谱分析会有较大误差。Yi =-3.00570 力31.059 +32.3141154.71 7-65.2011 -23.963+177.338122.68-20.025i17.391 +59.409i -32.591+1
12、2.514i21.468+16.3041-2.775+32.1 73i-3.39-6.952i 11.502+45.611i -17.362+2.902i1.68+37.675i -32.942-6.489i 9.594+28.757i -51.658-3.0681Power.=002.009?032.0097033.12?043.12?043.2027043.202?04915.377915.3773.8327033.832?031.2197031.219703726.679726.6791.043?031.04370 359.81659.8162.2137032.213703309.873
13、09.871.422?031.4227031.127?031.127703918.989918.9892.6787032.678703含摄占图14例2功率谱密度最大值及其对应频率【例3】某海域海面平均高度年际变化功率谱分析。本例时间序列长度为100年，即沪100。将时间序列延长到上2丁128。data :=0101-41.0312-35.03图16海面高度年际变化曲线图15例3数据局部显示（已经中心化）录入数据以后，根据时间序列长度对变量和Fourier变换定义如下:Y :=fft(y) M :=last(Y)7jPower := ( | Yj-w :=-k:=100.127yk：=0M:=6
14、4i:=0.Mj :=0.M freq. := j w J M从原始数据图像可以看到，时间序列可能具有某种内在节律（图16）。108.27446.026-0 . 020406080100年鼎号(M)50，一借助 MathcadTrace户0.09375,于是得到周期约为 61/h10.667 （年）。图18最大功率谱密度值对应频率借助Fourier变换结果，作出频谱图如下（图17）。功能发现功率谱密度最大值对应频率为【例4】在例2中，利用时间序列自相关函数分析其周期特性。在“利用Excel进行时间序列谱分析（I）” 一节，我们谈到可以定义 时间序列周期图为/（力）=7伍,+硝，=1，2，次式中
15、/（）为频率工处强度。以工为横轴，以？（）为纵轴，绘制时间序列 周期图，可以在最大值处找到时间序列周期。实际上，周期图还可以表作N-1/（/；）=讣k）产碘A-1-N这里A为时间序列自协方差函数，即有-1 N-KI&=犷力/州/=1这意味着，周期图是对功率谱直接估计结果。根据上面几个式子及其内在 关系，我们可以利用相关函数估计时间序列周期。我们知道，计算相关系 数公式有两种，本例采用下面公式n-kZ（%-初3-元）/q n 2r-1采用这个公式原因之一是：我们时间序列较短，而这个公式有效性较好； 原因之二是可以直接借助SPSS给出结果。将结果复制到Excel中间，转 换格式以后，再复制到Mat
16、hcad数据表（图19）,然后进行Fourier变换。data :=01010.66120.42230.0334-0.2945-0.5156-0.5767-0.5378-0.3989-0.159100.0910110.311120.412130.4513140.2214150.131516-0.02图19径流量自相关系数进行快速Fourier运算有关定义如下:t:=da别k:= 16.31:=0Y = fft(x)M := last(Y)M = 16i:= 0.MPo阴=（闾）-W ：=-2j:= 0.Mfreq. := w -根据计算结果容易得到频谱图（图20）。利用Trace可知，最大功率谱密 度对应频率为 户0.09375,于是得到周期约为E1/户10.667 （月），及直 接采用中心化数据进行Fourier变换得到结果相同（参见例2）。图20基于径流量自相关系数频谱图前面说过，周期图是对功率谱直接估计结果，而这种估计通常是不准 确。但是，在时间序列较短情况下，利用周期图估计时间序列周期反而更 为可靠。因此建议，在时间序列较短时运用周期图谱一一此时直接采用功 率谱反而不准确；在时间序列较长时运用功率谱一一此时采用周期图计算 工作量较大，而且有误差。可见虽然理论上二者等价，在实用时却各有适 应时间尺度范围。