《利用Mathcad进行时间序列的谱分析-Read》由会员分享，可在线阅读，更多相关《利用Mathcad进行时间序列的谱分析-Read（11页珍藏版）》请在金锄头文库上搜索。
1、利用Mathcad进行时间序列的谱分析（I）在这一节中我们采用与“利用Excel进行时间序列的谱分析（I）”一节相同的例子，以便计算结果能够相互参照。关于Mathcad的基本操作技巧（包括绘图）已在“回归分析”中备述，本节不再详细说明。【例1】某水文观测站测得一条河流从1979年6月到1980年5月共计12月份的断面平均流量。试借助功率谱分析判断时间序列是否存在周期性特征，并对周期长度进行估计。第一步，录入或拷入数据 可以将在Excel中经过中心化的数据直接录入或拷入Mathcad的数据表中，并定义为“data”。 图1 录入的中心化数据需要说明两点：第一，如果是从Excel中拷入数据，一定记
2、住检查最后一行数据，如果多出一行0数据点，必须将其删除。第二，也可以在Mathcad中对数据中心化。中心化的方法非常简单，将原始数据拷入Mathcad的数据表并将其定义为“data”（图2）。图2 原始数据表的局部显示（数据未经中心化）在被定义过的表下定义变量为然后利用求均值命令mean，键入公式回车，立即得到中心化的结果（图3）。 图3 中心化的数据不过，中心化以后需要重新插入一个数据表并重新定义该表（例如定义为biao）。将中心化的数据拷入“biao”中，然后在左方插入单元格（Insert Cells），再在新插入的一列单元格中录入数据序号（图4）。从这个意义上讲，在Mathcad中进行数
3、据处理不如直接将经过Excel预处理即中心化的数据拷贝过来。图4 新建的数据表第二步，绘制数据序列的变化图假定录入的是图1所示的中心化的数据，将表定义为data以后，在数据表下对变量进行如下定义：然后利用Graph中的曲线图工具容易绘制曲线图（图5）。 图5 中心化时间序列的变化曲线（1979年6月1980年5月）绘图的主要目的有两个：其一，可以通过曲线图直观地判断数据变化是否存在周期，对于本例，当然达不到这种效果；其二，如果拷贝的数据多出0数据行，可以及时将其删除，以免计算出错。第三步，进行快速Fourier变换首先是定义变量。在Mathcad中，系统默认的数据排列方式是从0到K。我们有12
4、个月份的数据，也就相当于拥有12个样本，即N=12，从而K=12-1=11。但是，Fourier变换要求时间序列的长度必须是T=2n个（n为正整数），因此，可以考虑在序列末尾加0，将其延长到T=24个。方法是定义序号k:=12.15，并设xk=0，这相当于在Excel加上4个0元素。这种补充定义已在上文给出。定义快速Fourier变换，命令为fft（注意：在Mathcad中，大写的FFT与小写的fft变换结果相差T倍，为了与Excel的结果对应，本例采用fft）。我们是要对x进行fft，令变换的结果为Y，Y将给出对x进行快速Fourier变换的结果；定义M为最后一个变换结果的序号，则可以通过l
5、ast自动从对称点(参见“利用Excel进行时间序列的谱分析（I）”一节的图15、图16)截取前半部结果，M=N/2=8；定义功率谱为Poweri，则有注意：如果使用命令FFT，则上式的分母应为T2。定义频率变化的步长为w=1/2；则频率为freqj=j/(2M)=j/T，并取j=0,1,2,M。全部定义的表达式如下：定义完成以后，键入“Yi=”，回车，得到x的FFT结果；键入“Poweri=”或“Yi2”，回车，得到功率谱密度（图6）。 图6 Fourier变换结果与功率谱密度第四步，绘制频谱图，识别周期利用Mathcad的Graph工具箱容易绘制频谱图（图7）。在图上点右键，出现一个选择菜
6、单（图8）；选中Trace，弹出X-Y Trace对话框（图9）；用鼠标点击谱线的最大值点，则Trace自动给出谱密度的极大值P(f)max=63157，对应的频率为f1=0.0625；点击Copy X按钮，将0.0625复制到Mathcad的工作表上，取倒数，立即到周期P=1/f1=16。当然，由于时间序列太短，这个周期是不准确的。图7 例1的频谱图图8 右键菜单与Trace的位置 图9 借助Trace显示功率谱密度最大值对应的频率【例2】对同一条河流的连续两年的平均月径流量（1961年6月1963年5月）进行Fourier分析，并识别其周期。步骤与例1相同，不详述。在录入中心化数据以后（图
7、10），需要重新定义数据的范围。图10 例2的数据表（局部，已中心化）考虑到这次N=24，可以将其延长到T=25=32。于是Fourier变换定义如下：画出时间序列的变化曲线图，从图上可以直观地判断径流量变化具有某种周期性，但不能确定（图11）。图11 例2的时间序列变化图（1961年6月1963年5月） 利用Fourier变换结果（图12）画出频谱图（图13）。借助Trace功能容易找到最大功率谱密度P(f)max=63157及其对应的频率f1=0.09375（图14）；点击Copy X按钮，将0.09375复制到Mathcad的工作表上，取倒数，可以得到周期P=1/f1=10.667。在最
8、大值附近存在一个与最大值非常接近的功率谱密度值为31107，对应的频率为0.0625。由于这个次最大点同样代表一个突变点它与临近值有很大差距，故须考虑这个数值。根据上例可知，其倒数为16。由此判断，时间序列的周期大约在1016之间。我们知道，河流径流量的周期为12个月。可见，对于较短的时间序列，功率谱分析会有较大的误差。 图12 例2的Fourier变换结果 图13 例2的频谱图图14 例2的功率谱密度最大值及其对应的频率【例3】某海域海面平均高度年际变化的功率谱分析。本例时间序列长度为100年，即N=100。将时间序列延长到T=27=128。图15 例3数据的局部显示（已经中心化）录入数据以
9、后，根据时间序列长度对变量和Fourier变换定义如下：从原始数据的图像可以看到，时间序列可能具有某种内在节律（图16）。图16 海面高度的年际变化曲线借助Fourier变换的结果，作出频谱图如下（图17）。 图17 海面高度变化的频谱图借助Mathcad的Trace功能发现功率谱密度最大值对应的频率为f=0.09375，于是得到周期约为P=1/f=10.667（年）。图18 最大功率谱密度值对应的频率【例4】 在例2中，利用时间序列的自相关函数分析其周期特性。在“利用Excel进行时间序列的谱分析（I）”一节，我们谈到可以定义时间序列的周期图为，式中I(fi)为频率fi处的强度。以fi为横轴
10、，以I(fi)为纵轴，绘制时间序列的周期图，可以在最大值处找到时间序列的周期。实际上，周期图还可以表作这里R为时间序列的自协方差函数，即有这意味着，周期图是对功率谱直接估计的结果。根据上面几个式子及其内在关系，我们可以利用相关函数估计时间序列的周期。我们知道，计算相关系数的公式有两种，本例采用下面公式采用这个公式的原因之一是：我们的时间序列较短，而这个公式的有效性较好；原因之二是可以直接借助SPSS给出结果。将结果复制到Excel中间，转换格式以后，再复制到Mathcad的数据表（图19），然后进行Fourier变换。图19 径流量的自相关系数进行快速Fourier运算的有关定义如下：根据计算结果容易得到频谱图（图20）。利用Trace可知，最大功率谱密度对应的频率为f=0.09375，于是得到周期约为P=1/f=10.667（月），与直接采用中心化数据进行Fourier变换得到的结果相同（参见例2）。 图20 基于径流量自相关系数的频谱图前面说过，周期图是对功率谱直接估计的结果，而这种估计通常是不准确的。但是，在时间序列较短的情况下，利用周期图估计时间序列的周期反而更为可靠。因此建议，在时间序列较短时运用周期图谱此时直接采用功率谱反而不准确；在时间序列较长时运用功率谱此时采用周期图计算工作量较大，而且有误差。可见虽然理论上二者等价，在实用时却各有适应的时间尺度范围。10
