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

利用Excel 进行时间序列的谱分析（I ）

在频域分析中，功率谱是揭示时间序列周期特性的最为有力的工具之一。下面列举几个例子，分别从不同的角度识别时间序列的周期。

1 时间序列的周期图

【例1】某水文观测站测得一条河流从1979年6月到1980年5月共计12月份的断面平均流量。试判断该河流的径流量变化是否具有周期性，周期长度大约为多少？

分析：假定将时间序列x t 展开为Fourier 级数，则可表示为

∑=++=

k

i t i i i i

t t f b t f a

x 1

)2sin 2cos (εππ (1)

式中f i 为频率，t 为时间序号，k 为周期分量的个数即主周期（基波）及其谐波的个数，εt 为标准误差（白噪声序列）。当频率f i 给定时，式（1）可以视为多元线性回归模型，可以证明，待定系数a i 、b i 的最小二乘估计为

∑∑====N

t i t

i

N

t i t

i t

f x

N

b t

f x

N a 112sin 2ˆ2cos 2ˆππ （2）

这里N 为观测值的个数。定义时间序列的周期图为

)(2

)(2

2

i i i b a N f I +=

，k i ,,2,1 = (3)

式中I (f i )为频率f i 处的强度。以f i 为横轴，以I (f i )为纵轴，绘制时间序列的周期图，可以在

最大值处找到时间序列的周期。对于本例，N =12，t =1,2,…,N ，f i =i /N ，下面借助Excel ，利用上述公式，计算有关参数并分析时间序列的周期特性。

第一步，录入数据，并将数据标准化或中心化（图1）。

图1 录入的数据及其中心化结果

中心化与标准化的区别在于，只需将原始数据减去均值，而不必再除以标准差。不难想到，中心化的数据均值为0，但方差与原始数据相同（未必为1）。 第二步，计算三角函数值

为了借助式（1）计算参数a i 、b i ，首先需要计算正弦值和余弦值。

取6,,2,1 =i ，则频率为12/6,,12/2,12/1/ ==N i f i （图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时的全部余弦值。依次类推，可以算出全部所要的余弦值（在D3-O8区域中）。根据对称性，我们的计算到k =6为止（图2）。注意，这里B1单元格是2π=6.28319（图中未能显示）。

在上面的计算中，只要将公式中的“COS ”换成“SIN ”，即可得到正弦值，不过为了计算过程清楚明白，最好在另外一个区域给出结算结果（在D17-O22区域中，参见图3）。

图3 计算正弦值的表格

第三步，计算参数a i 、b i

利用中心化的数据（仍然表作x t ）计算参数a i 、b i 。首先算出x t cos2πf i t 和x t sins2πf i t 。在D9单元格中输入公式“=D1*D3”，回车得到18.309；按住单元格的右下角右拉至O9单元格，得到f =1/12=0.083，t =1,2,…,12时的全部x t cos2πf i t 值；加和得39.584，再除以6，即得a 1=6.597。在D10单元格中输入公式“=D1*D4”，回车得到10.571；按住单元格的右下角右拉至O10单元格，得到f =2/12=0.083，t =1,2,…,12时的全部x t cos2πf i t 值；加和得-365.25，再除以6，得到a 2=-60.875。其余依此类推。

将上面公式中的余弦值换成正弦值，即可得到b i 值（见下表）。上面的计算过程相当于

采用式（2）进行逐步计算。

第四步，计算频率强度

利用式（3），非常容易算出I (f i )值。例如

914.174096)213.170597

.6(*6)(2

)(2

2

2

1211=+=+=b a N f I

其余依此类推（见图4）。

图4 计算频率强度

第五步，绘制时间序列周期图

利用图4中的数据，不难画出周期图（图5）。

图5 某河流径流量的周期图（1979年6月－1980年5月）

第六步，周期识别

关键是寻找频率的极值点或突变点。在本例中，没有极值点，但在f 1=1/12=0.0833处，频率强度突然增加（陡增），而此时T =1/f 1=12，故可判断时间序列可能存在一个12月的周期，即1年周期。

【例2】为了映证上述判断，我们借助同一条河流的连续两年的平均月径流量（1961年6月－1963年5月）。原始数据见下图（图6）。

图6 原始数据及部分处理结果

将原始数据回车时间序列变化图，可以初步估计具有12月变化周期，但不能肯定（图6）。

图6 径流量的月变化图（1961年6月－1963年5月）

按照例1给出的计算步骤，计算参数数a i、b i，进而计算频率强度（结果将图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的时间序列较短，分析的效果不佳，但计算过程简短。给出这个例子，主要是帮助大家理解Fourier变换过程和方法。为了进行Fourier分析，需要对数据进行预处理。第一，将数据中心化，即用原始数据减去其平均值。中心化的数据均值为0，我们对中心化的数据进行变换，其周期更为明显。第二，由于Fourier分析通常采用所谓快速Fourier变换（Fast Fourier Transformation，FFT），而FFT 要求数据必须为2n个，这里n为正整数（1,2,3,…），而我们的样本为N=12，它不是2的某