利用Excel进行时间序列的谱分析(I).pdf

《利用Excel进行时间序列的谱分析(I).pdf》由会员分享，可在线阅读，更多相关《利用Excel进行时间序列的谱分析(I).pdf（17页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

1、 1利用利用 Excel 进行时间序列的谱分析（进行时间序列的谱分析（I）在频域分析中，功率谱是揭示时间序列周期特性的最为有力的工具之一。下面列举几个例子，分别从不同的角度识别时间序列的周期。1 时间序列的周期图时间序列的周期图【例 1】某水文观测站测得一条河流从 1979 年 6 月到 1980 年 5 月共计 12 月份的断面平均流量。试判断该河流的径流量变化是否具有周期性，周期长度大约为多少？分析：假定将时间序列 xt展开为 Fourier 级数，则可表示为=+=kitiiiittfbtfax1)2sin2cos(1)式中 fi为频率，t 为时间序号，k 为周期分量的个数即主周期（基波）

2、及其谐波的个数，t为标准误差（白噪声序列）。当频率 fi给定时，式（1）可以视为多元线性回归模型，可以证明，待定系数 ai、bi的最小二乘估计为=NtitiNtititfxNbtfxNa112sin22cos2 （2）这里 N 为观测值的个数。定义时间序列的周期图为)(2)(22iiibaNfI+=，ki,2,1=(3)式中 I(fi)为频率 fi处的强度。以 fi为横轴，以 I(fi)为纵轴，绘制时间序列的周期图，可以在最大值处找到时间序列的周期。对于本例，N=12，t=1,2,N，fi=i/N，下面借助 Excel，利用上述公式，计算有关参数并分析时间序列的周期特性。第一步，录入数据，并将

3、数据标准化或中心化（图 1）。图 1 录入的数据及其中心化结果 图 1 录入的数据及其中心化结果 2中心化与标准化的区别在于，只需将原始数据减去均值，而不必再除以标准差。不难想到，中心化的数据均值为 0，但方差与原始数据相同（未必为 1）。第二步，计算三角函数值 为了借助式（1）计算参数 ai、bi，首先需要计算正弦值和余弦值。取6,2,1=i，则频率为12/6,12/2,12/1/=Nifi（图 1）。将频率写在单元格 C3-C14 中（根据对称性，我们只用前 6 个），将中心化的数据转置粘贴于第一行的单元格 D1-O1 中，月份的序号写在单元格 D2-O2 中（与中心化数据对齐）。图 2

4、计算余弦值的表格 图 2 计算余弦值的表格 在 D2 单元格中输入公式“=COS($B$1*$D$2*C3)”，回车得到 0.866；按住单元格的右下角右拉至 O3 单元格，得到 f=1/12=0.083，t=1,2,12 时的全部余弦值。在 D2 单元格中输入公式“=COS($B$1*$D$2*C4)”，回车得到 0.5；按住单元格的右下角右拉至 O4 单元格，得到 f=2/12=0.167，t=1,2,12 时的全部余弦值。依次类推，可以算出全部所要的余弦值（在D3-O8 区域中）。根据对称性，我们的计算到 k=6 为止（图 2）。注意，这里 B1 单元格是2=6.28319（图中未能显示

5、）。在上面的计算中，只要将公式中的“COS”换成“SIN”，即可得到正弦值，不过为了计算过程清楚明白，最好在另外一个区域给出结算结果（在 D17-O22 区域中，参见图 3）。图 3 计算正弦值的表格 图 3 计算正弦值的表格 第三步，计算参数 ai、bi 利用中心化的数据（仍然表作 xt）计算参数 ai、bi。首先算出 xtcos2fit 和 xtsins2fit。在D9 单元格中输入公式“=D1*D3”，回车得到 18.309；按住单元格的右下角右拉至 O9 单元格，得到 f=1/12=0.083，t=1,2,12 时的全部 xtcos2fit 值；加和得 39.584，再除以 6，即得a

6、1=6.597。在 D10 单元格中输入公式“=D1*D4”，回车得到 10.571；按住单元格的右下角右拉至 O10 单元格，得到 f=2/12=0.083，t=1,2,12 时的全部 xtcos2fit 值；加和得-365.25，再除以 6，得到 a2=-60.875。其余依此类推。将上面公式中的余弦值换成正弦值，即可得到 bi值（见下表）。上面的计算过程相当于 3采用式（2）进行逐步计算。第四步，计算频率强度 利用式（3），非常容易算出 I(fi)值。例如 914.174096)213.170597.6(*6)(2)(2221211=+=+=baNfI 其余依此类推（见图 4）。图 4

7、计算频率强度 图 4 计算频率强度 第五步，绘制时间序列周期图 利用图 4 中的数据，不难画出周期图（图 5）。02000040000600008000010000012000014000016000018000020000000.10.20.30.40.50.6频率频率强度 图 5 某河流径流量的周期图（1979 年 6 月1980 年 5 月）图 5 某河流径流量的周期图（1979 年 6 月1980 年 5 月）第六步，周期识别 关键是寻找频率的极值点或突变点。在本例中，没有极值点，但在 f1=1/12=0.0833 处，频率强度突然增加（陡增），而此时 T=1/f1=12，故可判断时间

8、序列可能存在一个 12 月的周期，即 1 年周期。4【例 2】为了映证上述判断，我们借助同一条河流的连续两年的平均月径流量（1961 年 6月1963 年 5 月）。原始数据见下图（图 6）。图 6 原始数据及部分处理结果 图 6 原始数据及部分处理结果 将原始数据回车时间序列变化图，可以初步估计具有 12 月变化周期，但不能肯定（图6）。050100150200250300350051015202530时间（月份序号）月平均径流量 图 6 径流量的月变化图（1961 年 6 月1963 年 5 月）图 6 径流量的月变化图（1961 年 6 月1963 年 5 月）5按照例 1 给出的计算步

9、骤，计算参数数 ai、bi，进而计算频率强度（结果将图 7）。然后绘制时间序列的周期图（图 8）。注意这里，N=24，我们取 k=12。图 7 参数和频率强度的计算结果 图 7 参数和频率强度的计算结果 从图 8 中可以看出，频率强度的最大值（极值点）对应于频率 f1=1/12=0.0833，故时间序列的周期判断为 T=1/f1=12。这与用 12 月的数据进行估计的结果是一致的，但由于例 2 的时间序列比例 1 的时间序列长 1 倍，故判断结果更为可靠。02000040000600008000010000012000014000000.10.20.30.40.50.6频率频率强度 图 8 某

10、河流径流量的周期图（1961 年 6 月1963 年 5 月）图 8 某河流径流量的周期图（1961 年 6 月1963 年 5 月）2 时间序列的频谱图时间序列的频谱图【例 3】首先考虑对例 1 的数据进行功率谱分析。例 1 的时间序列较短，分析的效果不佳，但计算过程简短。给出这个例子，主要是帮助大家理解 Fourier 变换过程和方法。为了进行 Fourier 分析，需要对数据进行预处理。第一，将数据中心化，即用原始数据减去其平均值。中心化的数据均值为 0，我们对中心化的数据进行变换，其周期更为明显。第二，由于 Fourier 分析通常采用所谓快速 Fourier 变换（Fast Four

11、ier Transformation，FFT），而 FFT要求数据必须为 2n个，这里 n 为正整数（1,2,3,），而我们的样本为 N=12，它不是 2 的某 6个 n 次方。因此，在中心化的数据后面加上 4 个 0，这样新的样本数为 N=1241624个，这才符合 FFT 的需要（图 9）。下面，我们对延长后的中心化数据进行 Fourier 变换。图 9 数据的中心化与“延长”图 9 数据的中心化与“延长”第一步，打开 Foureir 分析对话框 沿着主菜单的“工具（Tools）”“数据分析（Data Analysis）”路径打开数据分析选项框（图 10），从中选择“傅立叶分析（Fouri

12、er Analysis）”。图 10 在数据分析选项框中选择 Fourier 分析 图 10 在数据分析选项框中选择 Fourier 分析 第二步，定义变量和输出区域 确定之后，弹出傅立叶分析对话框，根据数据在工作表中的分布情况进行如下设置：将光标置于“输入区域”对应的空白栏，然后用鼠标选中单元格 C1-C17，这时空白栏中自动以绝对单元格的形式定义中心化数据的区域范围（即$C$1-$C$17）。选中“标志位于第一行”。选中输出区域，定义范围为 D2-D17（图 11）。注意：如果输入区域的数据范围定义为 C2-C17，则不要选中“标志位于第一行”，这与回归分析中的原始变量定义是一样的（图 1

13、2）。如果不定义输出区域范围，则变换结果将会自动给在新的工作表组上。这一点也与回归分析一样。7 图 11 选中“标志位于第一行”与数据输入范围的定义 图 11 选中“标志位于第一行”与数据输入范围的定义 图 12 不选中“标志位于第一行”与数据输入范围的定义 图 12 不选中“标志位于第一行”与数据输入范围的定义 第三步，结果转换 定义数据输入输出区域完成之后，确定，立即得到 Fourier 变换的结果（图 13）。图 13 傅立叶变换的结果 图 13 傅立叶变换的结果 8变换的结果为一组复数，相当于将 f(t)变成了 F()，实际上是将 xt变成了 XT(f)。我们知道，有了 f(t)的象函

14、数 F()就可以计算能量谱密度函数 S()，即有 2)()()()(FFFS=(4)相应地，有了 XT(f)也就容易计算功率谱（密度）TfXfPT2)()(=(5)式中 f 表示线频率，与角频率 的转换关系是=2/T，这里 T 为数据区间长度。如果将 XT(f)表作 XT(f)=A+jB（这里 A 为实部，B 为虚部），则有 222)()()()(iiiiiiiTiTiTBAjBAjBAfXfXfX+=+=（6）因此这一步是要分离变换结果的实部与虚部。逐个手动提取是非常麻烦而且容易出错的，可以利用 Excel 有关复数计算的命令。提取实部的 Excel 命令是 imreal。在 H2 单元格中

15、输入命令“=IMREAL(D2)”（这里 D2 为变换结果的第一个复数所在的单元格），回车得到第一个复数的实部 0；点中 H2 单元格的右下角，揿住鼠标左键下拉至 H17，得到全部的实部数值。提取虚部的命令是 imaginary。在 I2 单元格中输入公式“=IMAGINARY(D2)”，回车得到第一个复数的虚部 0；下拉至 I17，得到全部的虚部数值。根据式（5）、（6），功率谱密度的计算公式为 TBAfPiii22)(+=（7）考虑到本例中 T=N=16，在 J2 单元格中输入公式“=(H22+I22)/16”，回车得到第一个功率谱密度 0；下拉至 J17，得到全部谱密度数值（图 14）。

16、基于 FFT 的谱密度分布是对称的，可以看出，以 J10 所在的 3105.28 为对称点，上下的数值完全对称。图 14 功率谱密度的计算结果 图 14 功率谱密度的计算结果 9第四步，绘制频谱图 线频率 fi可以表作 NiTifi/=，Ni,2,1,0=-1 显然 f0=0/16=0，f1=1/16=0.0625，f2=2/16=0.125，f15=15/16=0.9375。在 Excel 中，容易计算频率的数值。将频率与功率谱对应起来（图 15），就可以画出频谱图。如果补上最后一个频率数值 f16=1 及其对应的功率 0，则可画出完全对称的谱图（图 16）。图 15 功率谱密度与频率的对应

17、关系 图 15 功率谱密度与频率的对应关系 01000020000300004000050000600007000000.20.40.60.81频率f功率谱密度P(f)图 16 对称分布的频谱图 图 16 对称分布的频谱图 由于功率谱图的对称性，画出整个谱图实在没有必要，因此，在实际工作中，通常只画出左半边（图 17）。1001000020000300004000050000600007000000.10.20.30.40.50.6频率f功率P(f)图 17 实用的频谱图 图 17 实用的频谱图 第五步，功率谱分析 频域分析的主要目标之一是判断时间序列是否存在周期性。从图 17 可以看出，功率

18、最大点对应的频率是 f1=0.0625，该频率对应的周期长度为 16。可见，在时间序列较短的情况下之间用功率谱图寻找时间序列的周期不如周期图准确。另外可以初步估计数据的性质。在图 17 中，去掉第一个 0 点，剩余的点一般呈幂指数分布（在双对数坐标图上点列具有直线趋势），可以拟合幂指数函数如下：ffP)(（8）P(f)=935.68f-1.4952R2=0.92331000100001000000.010.11频率f功率P(f)图 18 功率 P(f)与频率 f 的双对数坐标图 图 18 功率 P(f)与频率 f 的双对数坐标图 结果得到功率谱指数=1.49521.5。功率谱指数与时间序列的

19、Hurst 指数具有如下关系 12+=H （9）据此估计 Hurst 指数约为 0.25。我们知道，Hurst 指数介于 01 之间，当 H0.5 时，表明时间序列存在正的自相关，意味着系统演化具有持久性；当 H0 时，表明时间序列具有负的 11自相关，意味着系统演化具有反持久性；当 H=0 时，表明时间序列不存在自相关，过去与未来无关。对于这条河流的径流量而言，H=0.250.5，表明时间序列具有反持久性：过去的增量意味着今后的减少，过去的减少意味着未来的增加。因此，径流量必然周期性的变化。【例 4】下面对前述例 2 的数据进行 Fourier 变换，方法与例 3 相同，但由于 N=24，我

20、们取 T=32=25。也就是说，对于中心化的数据，要在后面添加 8 个 0 作为补充点数。基于FFT 的变换结果如下（见图 19）。图 19 例 2 数据（经中心化处理）的 FFT 变换结果 图 19 例 2 数据（经中心化处理）的 FFT 变换结果 计算功率谱除例 3 讲述的方法外，还可以利用 Excel 的另外两个命令实现：一是计算共轭复数的命令 imconjugate，首先求出)(fXT的共轭复数)(fXT；然后借助复数的乘积命令 improduct，计算复数的)(fXT与)(fXT的乘积2)(fXT；最后利用式（5）得到功率谱。不过，此时的时间序列长度视为 T=32。例如在 H2 中键

21、入公式“=IMCONJUGATE(D2)”，回车得到第一个共轭数，下拉至 H33，得到全部共轭数值。在 L2 中键入公式“=IMPRODUCT(D2,H2)”，回车，得到第一个复数乘积，下拉至 L33，得到全部2)(fXT值。最后用2)(fXT除以 32 得到功率谱密度（图 20）。12根据计算结果画出频谱图（图 21），从图上可以看出，频率密度的极值点对应的频率为0.09375，相应的周期为 T=10.667；在极值点附件存在一个次最大点，但相对于其他数值却显然又是突变点，该点对应的频率为 0.0625，相应的周期为 16。故可断定，该时间序列的周期比在 1016 之间。图 20 功率谱密度

22、值及其相应的频率 图 20 功率谱密度值及其相应的频率 0500010000150002000025000300003500000.10.20.30.40.50.6频率f功率P(f)图 21 例 4 的频谱图 图 21 例 4 的频谱图 13不过，由于这里的数据包含两个周期在内，用它估计数据的自相关性质不太准确。最后给一个较长时间序列的分析实例。【例 5】某海域测得多年连续的海平面年平均高度，发现大约每隔 11 年左右海平面达到一个最高值（图 22），于是科学家判断海平面的升降存在一个 11 年周期，与太阳黑子（sunspot）的 11 年周期变化一致，而太阳黑子的活动正是海平面升降的原因所在

23、。问题在于，前述 11 年周期是通过原始数据的变化曲线直观发现的，未必可靠。为此，可以进行一个功率谱分析，判断这种周期是否确实存在。图 22 某海域海面年平均高度的时间序列数据 图 22 某海域海面年平均高度的时间序列数据 首先利用原始数据画出时间序列变化的曲线图，观测数据变化特征，发现具有周期性波动特征，最高峰的时间间隔大约为 1012 年之间（图 23）。020406080100120140160180020406080100年度序号海面高度 图 23 某海域海平面年平均高度的年际变化曲线 图 23 某海域海平面年平均高度的年际变化曲线 然后将原始数据中心化，取 T=128=27，这意味着

24、需要将时间序列延长到 128 位，在计 14算频率时用 0,1,2,除以 128，在计算功率谱密度时，式（5）中的序列长度取 T=128。Fourier变换和功率谱密度的计算步骤与例 3、例 4 相同，不赘述。下面直接给出变换结果（图 24，图 25）。图 24 海面高度时间序列的 FFT 变换的部分结果图 24 海面高度时间序列的 FFT 变换的部分结果 02000400060008000100001200014000160001800000.10.20.30.40.50.6频率f功率P(f)图 25 海面高度变化的频谱图 图 25 海面高度变化的频谱图 在 Excel 上，将鼠标指向功率谱

25、密度的最大值处，立即显示频率为 f=0.09375，对应的功率为 P(f)=16832.77972。根据此处的频率可以算出周期为 T=1/f=1/0.09375=10.66711。这 15表明，海平面高度的变化的确存在一个为期大约 11 年的变化周期，与直观判断结果一致。由于这里采用的时间序列较长，故结果比较可靠。太阳黑子的活动周期约为 11 年稍多一点，二者非常接近。因此，海面高度变化与太阳黑子周期具有对应关系的猜想，在时间序列变化的节律方面，大致是可以接受的。图 26 功率谱密度最大值对应的频率显示结果 图 26 功率谱密度最大值对应的频率显示结果 【附录】我们在前面说过，式(1)实则一个

26、多元线性回归模型，式(2)是对式(1)中待定参数的最小二乘(OLS)估计。下面验证这种判断。将图 3 中计算出的正弦值 SIN、余弦值 COS 和中心化的径流量 Xt集中在一起，经复制选择性粘贴转置，可将数据重新排列如下（图 27）。图 27 重新排列的数据 图 27 重新排列的数据 若以正弦、余弦值为自变量，以中心化的径流量为因变量进行多元回归，Excel 拒绝给出结果，并弹出如下对话框显示拒绝计算的原因。图 28 Excel 拒绝给出回归结果 图 28 Excel 拒绝给出回归结果 16问题在于行数与列数不能相同，而本例中自变量和样本数都是 12，这就是问题所在。考虑到最后一个变量全都是

27、0，不妨去掉最后一个变量。由于式（1）不含常数项，在回归分析选项框中强制“常数为常数为 0”。确定以后，Excel 给出的多元回归结果如下（图 29）。图 29 利用图 27 中的前 12 列数据进行多元回归的结果 图 29 利用图 27 中的前 12 列数据进行多元回归的结果 将回归结果与利用式（2）直接结算的结果进行比较，发现除了 SIN6 的系数没有给出、COS6 的系数相差一倍以外，其余的结果基本一样（比较图 29 与下表）。将图 27 中的数据复制到 SPSS 的工作表中，利用 SPSS 进行回归。由于式(1)中没有截距（intercept），即常数项为零，故在回归分析的选项设置中必

28、须强制回归线通过坐标原点，即回归分析选项（Linear Regression:Options）不选“包括常数项（Include constant in equation）”(图 30)，这与 Excel 中的设置“常数为 0”是一样道理。图 29 将数据复制到 SPSS 中 图 29 将数据复制到 SPSS 中 从回归结果可以看出，COS6 的系数依然不对，SIN6 的系数有极大偏差。可见强行进行某种违反规则的回归必然出错。我们进行上述对比，旨在说明，式（1）本质上是一个多元线性回归模型，在一定条件下可以利用多元回归进行参数估计。但是，在绝大多数情况下，我们无法直接进行多元回归，17即便回归也

29、会引起某些误会。因此之故，功率谱分析必须沿着功率谱计算的套路进行。有关技术将会在不断锻炼的过程中逐步掌握。图图 30 强制回归线通过原点强制回归线通过原点 Coefficientsa,b6.597.000.033.-60.875.000-.305.15.017.000.075.45.092.000.226.6.886.000.035.-18.575.000-.132.170.213.000.853.-10.089.000-.051.-44.833.000-.225.2.295.000.012.44.553.000.223.-426.169.000.000.COS1COS2COS3COS4COS5COS6SIN1SIN2SIN3SIN4SIN5SIN6Model1BStd.ErrorUnstandardizedCoefficientsBetaStandardizedCoefficientstSig.Dependent Variable:XTa.Linear Regression through the Originb.SPSS给出的回归结果 SPSS给出的回归结果