当前位置:文档之家› 利用Excel进行时间序列的谱分析

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

利用Excel 进行时间序列的谱分析(I )在频域分析中,功率谱是揭示时间序列周期特性的最为有力的工具之一。

下面列举几个例子,分别从不同的角度识别时间序列的周期。

1 时间序列的周期图【例1】某水文观测站测得一条河流从1979年6月到1980年5月共计12月份的断面平均流量。

试判断该河流的径流量变化是否具有周期性,周期长度大约为多少?分析:假定将时间序列x t 展开为Fourier 级数,则可表示为∑=++=ki 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 的最小二乘估计为∑∑====Nt i tiNt i ti tf xNb tf xN a112sin 2ˆ2cos 2ˆππ (2)这里N 为观测值的个数。

定义时间序列的周期图为)(2)(22i 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)(2221211=+=+=b a N f I 其余依此类推(见图4)。

图4 计算频率强度第五步,绘制时间序列周期图利用图4中的数据,不难画出周期图(图5)。

02000040000600008000010000012000014000016000018000020000000.10.20.30.40.50.6频率频率强度图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)。

050100150200250300350051015202530时间(月份序号)月平均径流量图6 径流量的月变化图(1961年6月-1963年5月)按照例1给出的计算步骤,计算参数数a i 、b i ,进而计算频率强度(结果将图7)。

然后绘制时间序列的周期图(图8)。

注意这里,N =24,我们取k =12。

图7 参数和频率强度的计算结果从图8中可以看出,频率强度的最大值(极值点)对应于频率f 1=1/12=0.0833,故时间序列的周期判断为T =1/f 1=12。

这与用12月的数据进行估计的结果是一致的,但由于例2的时间序列比例1的时间序列长1倍,故判断结果更为可靠。

02000040000600008000010000012000014000000.10.20.30.40.50.6频率频率强度图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的某个n次方。

因此,在中心化的数据后面加上4个0,这样新的样本数为N′=12+4=16=24个,这才符合FFT的需要(图9)。

下面,我们对延长后的中心化数据进行Fourier变换。

图9 数据的中心化与“延长”第一步,打开Foureir分析对话框沿着主菜单的“工具(Tools)”→“数据分析(Data Analysis)”路径打开数据分析选项框(图10),从中选择“傅立叶分析(Fourier Analysis)”。

图10 在数据分析选项框中选择Fourier分析第二步,定义变量和输出区域确定之后,弹出傅立叶分析对话框,根据数据在工作表中的分布情况进行如下设置:将光标置于“输入区域”对应的空白栏,然后用鼠标选中单元格C1-C17,这时空白栏中自动以绝对单元格的形式定义中心化数据的区域范围(即$C$1-$C$17)。

选中“标志位于第一行”。

选中输出区域,定义范围为D2-D17(图11)。

注意:如果输入区域的数据范围定义为C2-C17,则不要选中“标志位于第一行”,这与回归分析中的原始变量定义是一样的(图12)。

如果不定义输出区域范围,则变换结果将会自动给在新的工作表组上。

这一点也与回归分析一样。

图11 选中“标志位于第一行”与数据输入范围的定义图12 不选中“标志位于第一行”与数据输入范围的定义第三步,结果转换定义数据输入-输出区域完成之后,确定,立即得到Fourier变换的结果(图13)。

图13 傅立叶变换的结果变换的结果为一组复数,相当于将f (t )变成了F (ω),实际上是将x t 变成了X T (f )。

我们知道,有了f (t )的象函数F (ω)就可以计算能量谱密度函数S (ω),即有2)()()()(ωωωωF F F S == (4)相应地,有了X T (f )也就容易计算功率谱(密度)Tf X f P T 2)()(=(5)式中f 表示线频率,与角频率ω的转换关系是ω=2π/T ,这里T 为数据区间长度。

如果将X T (f )表作X T (f )=A +jB (这里A 为实部,B 为虚部),则有222))(()()()(i i i i i i i T i T i T B A jB A jB A f X f X f X +=-+== (6)因此这一步是要分离变换结果的实部与虚部。

逐个手动提取是非常麻烦而且容易出错的,可以利用Excel 有关复数计算的命令。

提取实部的Excel 命令是imreal 。

在H2单元格中输入命令“=IMREAL(D2)”(这里D2为变换结果的第一个复数所在的单元格),回车得到第一个复数的实部0;点中H2单元格的右下角,揿住鼠标左键下拉至H17,得到全部的实部数值。

提取虚部的命令是imaginary 。

在I2单元格中输入公式“=IMAGINARY(D2)”,回车得到第一个复数的虚部0;下拉至I17,得到全部的虚部数值。

根据式(5)、(6),功率谱密度的计算公式为TB A f P i i i 22)(+= (7)考虑到本例中T =N =16,在J2单元格中输入公式“=(H2^2+I2^2)/16”,回车得到第一个功率谱密度0;下拉至J17,得到全部谱密度数值(图14)。

基于FFT 的谱密度分布是对称的,可以看出,以J10所在的3105.28为对称点,上下的数值完全对称。

图14 功率谱密度的计算结果第四步,绘制频谱图线频率f i 可以表作N i T i f i //==,N i ,,2,1,0 =-1显然f 0=0/16=0,f 1=1/16=0.0625,f 2=2/16=0.125,…,f 15=15/16=0.9375。

在Excel 中,容易计算频率的数值。

将频率与功率谱对应起来(图15),就可以画出频谱图。

相关主题