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

1、利用Excel进行时间序列的谱分析（I）在频域分析中，功率谱是揭示时间序列周期特性的最为有力的工具之一。下面列举几个例子，分别从不同的角度识别时间序列的周期。1 时间序列的周期图【例1】某水文观测站测得一条河流从1979年6月到1980年5月共计12月份的断面平均流量。试判断该河流的径流量变化是否具有周期性，周期长度大约为多少？分析：假定将时间序列xt展开为Fourier级数，则可表示为 (1)式中fi为频率，t为时间序号，k为周期分量的个数即主周期（基波）及其谐波的个数，t为标准误差（白噪声序列）。当频率fi给定时，式（1）可以视为多元线性回归模型，可以证明，待定系数ai、bi的最小二乘估计

2、为 （2）这里N为观测值的个数。定义时间序列的周期图为， (3)式中I(fi)为频率fi处的强度。以fi为横轴，以I(fi)为纵轴，绘制时间序列的周期图，可以在最大值处找到时间序列的周期。对于本例，N=12，t=1,2,N，fi=i/N，下面借助Excel，利用上述公式，计算有关参数并分析时间序列的周期特性。第一步，录入数据，并将数据标准化或中心化（图1）。图1 录入的数据及其中心化结果中心化与标准化的区别在于，只需将原始数据减去均值，而不必再除以标准差。不难想到，中心化的数据均值为0，但方差与原始数据相同（未必为1）。第二步，计算三角函数值为了借助式（1）计算参数ai、bi，首先需要计算正弦

3、值和余弦值。取，则频率为（图1）。将频率写在单元格C3-C14中（根据对称性，我们只用前6个），将中心化的数据转置粘贴于第一行的单元格D1-O1中，月份的序号写在单元格D2-O2中（与中心化数据对齐）。图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时的全部余弦值。依次类推，可以算

4、出全部所要的余弦值（在D3-O8区域中）。根据对称性，我们的计算到k=6为止（图2）。注意，这里B1单元格是2=6.28319（图中未能显示）。在上面的计算中，只要将公式中的“COS”换成“SIN”，即可得到正弦值，不过为了计算过程清楚明白，最好在另外一个区域给出结算结果（在D17-O22区域中，参见图3）。图3 计算正弦值的表格第三步，计算参数ai、bi利用中心化的数据（仍然表作xt）计算参数ai、bi。首先算出xtcos2fit和xtsins2fit。在D9单元格中输入公式“=D1*D3”，回车得到18.309；按住单元格的右下角右拉至O9单元格，得到f=1/12=0.083，t=1,2,

5、12时的全部xtcos2fit值；加和得39.584，再除以6，即得a1=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值（见下表）。上面的计算过程相当于采用式（2）进行逐步计算。第四步，计算频率强度利用式（3），非常容易算出I(fi)值。例如其余依此类推（见图4）。图4 计算频率强度第五步，绘制时间序列周期图利用图4中的数据，不难画出周

6、期图（图5）。图5 某河流径流量的周期图（1979年6月1980年5月）第六步，周期识别关键是寻找频率的极值点或突变点。在本例中，没有极值点，但在f1=1/12=0.0833处，频率强度突然增加（陡增），而此时T=1/f1=12，故可判断时间序列可能存在一个12月的周期，即1年周期。【例2】为了映证上述判断，我们借助同一条河流的连续两年的平均月径流量（1961年6月1963年5月）。原始数据见下图（图6）。图6 原始数据及部分处理结果将原始数据回车时间序列变化图，可以初步估计具有12月变化周期，但不能肯定（图6）。图6 径流量的月变化图（1961年6月1963年5月）按照例1给出的计算步骤，计

7、算参数数ai、bi，进而计算频率强度（结果将图7）。然后绘制时间序列的周期图（图8）。注意这里，N=24，我们取k=12。图7 参数和频率强度的计算结果从图8中可以看出，频率强度的最大值（极值点）对应于频率f1=1/12=0.0833，故时间序列的周期判断为T=1/f1=12。这与用12月的数据进行估计的结果是一致的，但由于例2的时间序列比例1的时间序列长1倍，故判断结果更为可靠。图8 某河流径流量的周期图（1961年6月1963年5月）2 时间序列的频谱图【例3】首先考虑对例1的数据进行功率谱分析。例1的时间序列较短，分析的效果不佳，但计算过程简短。给出这个例子，主要是帮助大家理解Fouri

8、er变换过程和方法。为了进行Fourier分析，需要对数据进行预处理。第一，将数据中心化，即用原始数据减去其平均值。中心化的数据均值为0，我们对中心化的数据进行变换，其周期更为明显。第二，由于Fourier分析通常采用所谓快速Fourier变换（Fast Fourier Transformation，FFT），而FFT要求数据必须为2n个，这里n为正整数（1,2,3,），而我们的样本为N=12，它不是2的某个n次方。因此，在中心化的数据后面加上4个0，这样新的样本数为N=1241624个，这才符合FFT的需要（图9）。下面，我们对延长后的中心化数据进行Fourier变换。图9 数据的中心化与“

9、延长”第一步，打开Foureir分析对话框沿着主菜单的“工具（Tools）”“数据分析（Data Analysis）”路径打开数据分析选项框（图10），从中选择“傅立叶分析（Fourier Analysis）”。图10 在数据分析选项框中选择Fourier分析第二步，定义变量和输出区域确定之后，弹出傅立叶分析对话框，根据数据在工作表中的分布情况进行如下设置：将光标置于“输入区域”对应的空白栏，然后用鼠标选中单元格C1-C17，这时空白栏中自动以绝对单元格的形式定义中心化数据的区域范围（即$C$1-$C$17）。选中“标志位于第一行”。选中输出区域，定义范围为D2-D17（图11）。注意：如果输

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

11、XT(f)也就容易计算功率谱（密度） (5)式中f表示线频率，与角频率的转换关系是=2/T，这里T为数据区间长度。如果将XT(f)表作XT(f)=A+jB（这里A为实部，B为虚部），则有 （6）因此这一步是要分离变换结果的实部与虚部。逐个手动提取是非常麻烦而且容易出错的，可以利用Excel有关复数计算的命令。提取实部的Excel命令是imreal。在H2单元格中输入命令“=IMREAL(D2)”（这里D2为变换结果的第一个复数所在的单元格），回车得到第一个复数的实部0；点中H2单元格的右下角，揿住鼠标左键下拉至H17，得到全部的实部数值。提取虚部的命令是imaginary。在I2单元格中输入公

12、式“=IMAGINARY(D2)”，回车得到第一个复数的虚部0；下拉至I17，得到全部的虚部数值。根据式（5）、（6），功率谱密度的计算公式为 （7）考虑到本例中T=N=16，在J2单元格中输入公式“=(H22+I22)/16”，回车得到第一个功率谱密度0；下拉至J17，得到全部谱密度数值（图14）。基于FFT的谱密度分布是对称的，可以看出，以J10所在的3105.28为对称点，上下的数值完全对称。图14 功率谱密度的计算结果第四步，绘制频谱图线频率fi可以表作，-1显然f0=0/16=0，f1=1/16=0.0625，f2=2/16=0.125，f15=15/16=0.9375。在Excel

13、中，容易计算频率的数值。将频率与功率谱对应起来（图15），就可以画出频谱图。如果补上最后一个频率数值f16=1及其对应的功率0，则可画出完全对称的谱图（图16）。图15 功率谱密度与频率的对应关系图16 对称分布的频谱图由于功率谱图的对称性，画出整个谱图实在没有必要，因此，在实际工作中，通常只画出左半边（图17）。图17 实用的频谱图第五步，功率谱分析频域分析的主要目标之一是判断时间序列是否存在周期性。从图17可以看出，功率最大点对应的频率是f1=0.0625，该频率对应的周期长度为16。可见，在时间序列较短的情况下之间用功率谱图寻找时间序列的周期不如周期图准确。另外可以初步估计数据的性质。在

14、图17中，去掉第一个0点，剩余的点一般呈幂指数分布（在双对数坐标图上点列具有直线趋势），可以拟合幂指数函数如下： （8）图18 功率P(f)与频率f的双对数坐标图结果得到功率谱指数=1.49521.5。功率谱指数与时间序列的Hurst指数具有如下关系 （9）据此估计Hurst指数约为0.25。我们知道，Hurst指数介于01之间，当H>0.5时，表明时间序列存在正的自相关，意味着系统演化具有持久性；当H<0时，表明时间序列具有负的自相关，意味着系统演化具有反持久性；当H=0时，表明时间序列不存在自相关，过去与未来无关。对于这条河流的径流量而言，H=0.25<0.5，表明时间序

15、列具有反持久性：过去的增量意味着今后的减少，过去的减少意味着未来的增加。因此，径流量必然周期性的变化。【例4】下面对前述例2的数据进行Fourier变换，方法与例3相同，但由于N=24，我们取T=32=25。也就是说，对于中心化的数据，要在后面添加8个0作为补充点数。基于FFT的变换结果如下（见图19）。图19 例2数据（经中心化处理）的FFT变换结果计算功率谱除例3讲述的方法外，还可以利用Excel的另外两个命令实现：一是计算共轭复数的命令imconjugate，首先求出的共轭复数；然后借助复数的乘积命令improduct，计算复数的与的乘积；最后利用式（5）得到功率谱。不过，此时的时间序列

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

17、例4的频谱图不过，由于这里的数据包含两个周期在内，用它估计数据的自相关性质不太准确。最后给一个较长时间序列的分析实例。【例5】某海域测得多年连续的海平面年平均高度，发现大约每隔11年左右海平面达到一个最高值（图22），于是科学家判断海平面的升降存在一个11年周期，与太阳黑子（sunspot）的11年周期变化一致，而太阳黑子的活动正是海平面升降的原因所在。问题在于，前述11年周期是通过原始数据的变化曲线直观发现的，未必可靠。为此，可以进行一个功率谱分析，判断这种周期是否确实存在。图22 某海域海面年平均高度的时间序列数据首先利用原始数据画出时间序列变化的曲线图，观测数据变化特征，发现具有周期性波

18、动特征，最高峰的时间间隔大约为1012年之间（图23）。图23 某海域海平面年平均高度的年际变化曲线 然后将原始数据中心化，取T=128=27，这意味着需要将时间序列延长到128位，在计算频率时用0,1,2,除以128，在计算功率谱密度时，式（5）中的序列长度取T=128。Fourier变换和功率谱密度的计算步骤与例3、例4相同，不赘述。下面直接给出变换结果（图24，图25）。图24 海面高度时间序列的FFT变换的部分结果图25 海面高度变化的频谱图在Excel上，将鼠标指向功率谱密度的最大值处，立即显示频率为f=0.09375，对应的功率为P(f)=16832.77972。根据此处的频率可以

19、算出周期为T=1/f=1/0.09375=10.66711。这表明，海平面高度的变化的确存在一个为期大约11年的变化周期，与直观判断结果一致。由于这里采用的时间序列较长，故结果比较可靠。太阳黑子的活动周期约为11年稍多一点，二者非常接近。因此，海面高度变化与太阳黑子周期具有对应关系的猜想，在时间序列变化的节律方面，大致是可以接受的。图26 功率谱密度最大值对应的频率显示结果【附录】我们在前面说过，式(1)实则一个多元线性回归模型，式(2)是对式(1)中待定参数的最小二乘(OLS)估计。下面验证这种判断。将图3中计算出的正弦值SIN、余弦值COS和中心化的径流量Xt集中在一起，经复制选择性粘贴转置，可将数据重新排列如下（图27）。图27 重新排列的数据若以正弦、余弦值为自变量，以中心化的径流量为因变量进行多元回归，Excel拒绝给出结果，并弹出如下对话框显示拒绝计算的原因。图28 Excel拒绝给出回