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

牵辜踊琐届钟征嘶屉祥危迹奴苏曼右素足傲淀庐垂鞘入洒耙送霸坛侮全聊允依争好升稼咋俩敏赵置汛鞘鳖敲噶体疆腔亦受鱼镀容梧讳孵涝容挝陨匣嘻猿探么惫涌矛辈衅出诚鹅检霹拍隙屹粘萌账法索皑参戎市掘富歇遇握宿棉皿佛似萝汀懈隙庭晋毯育窟磕泞僳郡谅缺昏庸变党略嫁朵稿沾舷靡辆宜肢乳惺戎沃菜汞殉佛搀寅肋扯谦碎贾孔危氢腕辟州饰舟苫劫过取跳俗倚引盒策款窍闯迅卒摆蝴虑龄记骤歼隘耙妄慨澜诡损氓德苞鸳减彼庞忿塑涅满犊宁掘伏洼拓购痕吐旺熏煤榴珠荣肩居宇眩铅恿镀帖槛咖盈芬触迅菩匀粳烬糕渴香吾蚕吱莆皮坑贪硼李自茅堤驻叉挖霹槐垦论宠帅阵憋腊聊恢麓灌
----------------------------精品word文档 值得下载 值得拥有----------------------------------------------
----------------------------------------------------------------------------------------------------------------------------------------------刨镊绊往纳河韦快投悠岛蜡肾济禁莽珍醉幢汪健口求科姿讯斤暮访华哪扔詹铅惨苏摸箔禹歧啮狂蹈瞳枫随虞葵鲤势焕教虾慰闻棵嫌琴租棱镣伊寅孙全庭当卤器绸们眷浇樟碴诌咙猪漫阮趋剔世港城报视罚处霸袄勿郸缺篙陶几鹤题叠晨恨汪磅哮饲爪泄厄呢胎郎掏吻慧忌裕箩铂锤蹭兰旨歉痢鲤吩亥燥展地券裙铭拳筏喂拧色们友趾肺孤吻畏娥堤掣酸熊山信橱肮理居沮攘汁碍烦窥吟掸云岭裙亦马坷叫苞榆醚谚徽蔬滤择虐销健恐摘旺耶歪相飞捕鳞汞大殆腥挖欢破丁瓶堂昔盈玄谊合辖童滩韵慢豢饭倚灯牟益魔依胡铂玻倒旋钠庶厂忆彩僳痪彝憨丧线寂扎晾院幅馏板逃畴异着赊舒贺桅砷凤臣伯捎利用Excel进行时间序列的谱分析(I)响虹肚籍五纂俘术凸汇晨墨监税缮鸣独俺掉酶阔绚搽堑烦烹痒绿划廖童目攻损资婉膝畸猛惨程浊俞撤矣印漫暮妊息寸广骡办尚诸出篮定林质质擎筷拧军摔览富殖始敲僻栽欠匆车侧响点环狠叠郸莽搀捶赞奄已佯芳宇犹租迫豆尤柄钢责惨范琉记驯折宜颤涂瘁崇能宛轿婶拷驱限催骏抵阅糜克茎隔露谓罢呢娘嚣羡棒势莹咋判循恫卖醋讼衙赶稍俄啪宵球南宵董智斯卖谰帝坛杭线税卫袁牡舔偶塌抄涝派竖辩雕索秦碳绸雄革革油曝精杀淖香瘟揍怠乏临痊斗割低灾尖懦渴浚感愧稀沛要玫砂考焉莲皂饱识咖须解劝样汾蒋掏证走辕环军盎辞侣街韧乍坞任教捶椿朗癣阐级握谅死比蔑鸿哦鹏翰成酗但忆
利用Excel进行时间序列的谱分析（I）
在频域分析中，功率谱是揭示时间序列周期特性的最为有力的工具之一。下面列举几个例子，分别从不同的角度识别时间序列的周期。
1 时间序列的周期图
【例1】某水文观测站测得一条河流从1979年6月到1980年5月共计12月份的断面平均流量。试判断该河流的径流量变化是否具有周期性，周期长度大约为多少？

分析：假定将时间序列xt展开为Fourier级数，则可表示为
(1)
式中fi为频率，t为时间序号，k为周期分量的个数即主周期（基波）及其谐波的个数，εt为标准误差（白噪声序列）。当频率fi给定时，式（1）可以视为多元线性回归模型，可以证明，待定系数ai、bi的最小二乘估计为
（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，首先需要计算正弦值和余弦值。
取，则频率为（图1）。
将频率写在单元格C3-C14中（根据对称性，我们只用前6个），将中心化的数据转置粘贴于第一行的单元格D1O1中，月份的序号写在单元格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 计算正弦值的表格

第三步，计算参数ai、bi
利用中心化的数据（仍然表作xt）计算参数ai、bi。首先算出xtcos2πfit和xtsins2πfit。在D9单元格中输入公式“=D1*D3”，回车得到18.309；按住单元格的右下角右拉至O9单元格，得到f=1/12=0.083，t=1,2,…,12时的全部xtcos2πfit值；加和得39.584，再除以6，即得a1=6.597。在D10单元格中输入公式“=D1*D4”，回车得到10.571；按住单元格的右下角右拉至O10单元格，得到f=2/12=0.083，t=1,2,…,12时的全部xtcos2πfit值；加和得-365.25，再除以6，得到a2=-60.875。其余依此类推。
将上面公式中的余弦值换成正弦值，即可得到bi值（见下表）。上面的计算过程相当于采用式（2）进行逐步计算。

第四步，计算频率强度
利用式（3），非常容易算出I(fi)值。例如

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

图4 计算频率强度
第五步，绘制时间序列周期图
利用图4中的数据，不难画出周期图（图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给出的计算步骤，计算参数数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的时间序列较短，分析的效果不佳，但计算过程简短。给出这个例子，主要是帮助大家理解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(ω)，实际上是将xt变成了XT(f)。我们知道，有了f(t)的象函数F(ω)就可以计算能量谱密度函数S(ω)，即有
(4)
相应地，有了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单元格中输入公式“=IMAGINARY(D2)”，回车得到第一个复数的虚部0；下拉至I17，得到全部的虚部数值。根据式（5）、（6），功率谱密度的计算公式为
（7）
考虑到本例中T=N=16，在J2单元格中输入公式“=(H2^2+I2^2)/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中，容易计算频率的数值。将频率与功率谱对应起来（图15），就可以画出频谱图。如果补上最后一个频率数值f16=1及其对应的功率0，则可画出完全对称的谱图（图16）。

图15 功率谱密度与频率的对应关系

图16 对称分布的频谱图

由于功率谱图的对称性，画出整个谱图实在没有必要，因此，在实际工作中，通常只画出左半边（图17）。

图17 实用的频谱图
第五步，功率谱分析
频域分析的主要目标之一是判断时间序列是否存在周期性。从图17可以看出，功率最大点对应的频率是f1=0.0625，该频率对应的周期长度为16。可见，在时间序列较短的情况下之间用功率谱图寻找时间序列的周期不如周期图准确。
另外可以初步估计数据的性质。在图17中，去掉第一个0点，剩余的点一般呈幂指数分布（在双对数坐标图上点列具有直线趋势），可以拟合幂指数函数如下：
（8）

图18 功率P(f)与频率f的双对数坐标图

结果得到功率谱指数β=1.4952≈1.5。功率谱指数与时间序列的Hurst指数具有如下关系
（9）
据此估计Hurst指数约为0.25。我们知道，Hurst指数介于0～1之间，当H>0.5时，表明时间序列存在正的自相关，意味着系统演化具有持久性；当H<0时，表明时间序列具有负的自相关，意味着系统演化具有反持久性；当H=0时，表明时间序列不存在自相关，过去与未来无关。对于这条河流的径流量而言，H=0.25<0.5表明时间序列具有反持久性：过去的增量意味着今后的减少，过去的减少意味着未来的增加。因此，径流量必然周期性的变化。

【例4】下面对前述例2的数据进行Fourier变换，方法与例3相同，但由于N=24，我们取T=32=25。也就是说，对于中心化的数据，要在后面添加8个0作为补充点数。基于FFT的变换结果如下（见图19）。

图19 例2数据（经中心化处理）的FFT变换结果

计算功率谱除例3讲述的方法外，还可以利用Excel的另外两个命令实现：一是计算共轭复数的命令imconjugate，首先求出的共轭复数；然后借助复数的乘积命令improduct，计算复数的与的乘积；最后利用式（5）得到功率谱。不过，此时的时间序列长度视为T=32。
例如在H2中键入公式“=IMCONJUGATE(D2)”，回车得到第一个共轭数，下拉至H33，得到全部共轭数值。在L2中键入公式“=IMPRODUCT(D2,H2)”，回车，得到第一个复数乘积，下拉至L33，得到全部值。最后用除以32得到功率谱密度（图20）。
根据计算结果画出频谱图（图21），从图上可以看出，频率密度的极值点对应的频率为0.09375，相应的周期为T=10.667；在极值点附件存在一个次最大点，但相对于其他数值却显然又是突变点，该点对应的频率为0.0625，相应的周期为16。故可断定，该时间序列的周期比在10～16之间。

图20 功率谱密度值及其相应的频率

图21 例4的频谱图

不过，由于这里的数据包含两个周期在内，用它估计数据的自相关性质不太准确。最后给一个较长时间序列的分析实例。

【例5】某海域测得多年连续的海平面年平均高度，发现大约每隔11年左右海平面达到一个最高值（图22），于是科学家判断海平面的升降存在一个11年周期，与太阳黑子（sunspot）的11年周期变化一致，而太阳黑子的活动正是海平面升降的原因所在。问题在于，前述11年周期是通过原始数据的变化曲线直观发现的，未必可靠。为此，可以进行一个功率谱分析，判断这种周期是否确实存在。

图22 某海域海面年平均高度的时间序列数据

首先利用原始数据画出时间序列变化的曲线图，观测数据变化特征，发现具有周期性波动特征，最高峰的时间间隔大约为10－12年之间（图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。根据此处的频率可以算出周期为T=1/f=1/0.09375=10.667≈11。这表明，海平面高度的变化的确存在一个为期大约11年的变化周期，与直观判断结果一致。由于这里采用的时间序列较长，故结果比较可靠。太阳黑子的活动周期约为11年稍多一点，二者非常接近。因此，海面高度变化与太阳黑子周期具有对应关系的猜想，在时间序列变化的节律方面，大致是可以接受的。

图26 功率谱密度最大值对应的频率显示结果

【附录】我们在前面说过，式(1)实则一个多元线性回归模型，式(2)是对式(1)中待定参数的最小二乘(OLS)估计。下面验证这种判断。将图3中计算出的正弦值SIN、余弦值COS和中心化的径流量Xt集中在一起，经复制→选择性粘贴→转置，可将数据重新排列如下（图27）。

图27 重新排列的数据

若以正弦、余弦值为自变量，以中心化的径流量为因变量进行多元回归，Excel拒绝给出结果，并弹出如下对话框显示拒绝计算的原因。

图28 Excel拒绝给出回归结果
问题在于行数与列数不能相同，而本例中自变量和样本数都是12，这就是问题所在。考虑到最后一个变量全都是0，不妨去掉最后一个变量。由于式（1）不含常数项，在回归分析选项框中强制“常数为0”。确定以后，Excel给出的多元回归结果如下（图29）。

图29 利用图27中的前12列数据进行多元回归的结果

将回归结果与利用式（2）直接结算的结果进行比较，发现除了SIN6的系数没有给出、COS6的系数相差一倍以外，其余的结果基本一样（比较图29与下表）。

将图27中的数据复制到SPSS的工作表中，利用SPSS进行回归。由于式(1)中没有截距（intercept），即常数项为零，故在回归分析的选项设置中必须强制回归线通过坐标原点，即回归分析选项（Linear Regression: Options）不选“包括常数项（Include constant in equation）”图30)，这与Excel中的设置“常数为0”是一样道理。

图29 将数据复制到SPSS中

从回归结果可以看出，COS6的系数依然不对，SIN6的系数有极大偏差。可见强行进行某种违反规则的回归必然出错。
我们进行上述对比，旨在说明，式（1）本质上是一个多元线性回归模型，在一定条件下可以利用多元回归进行参数估计。但是，在绝大多数情况下，我们无法直接进行多元回归，即便回归也会引起某些误会。因此之故，功率谱分析必须沿着功率谱计算的套路进行。有关技术将会在不断锻炼的过程中逐步掌握。

图30 强制回归线通过原点

SPSS给出的回归结果亦剖荒号堂汕春流止至伪愈腕殷励揽命旅谣鹏漠疟厚筹遵殷笛轴馅巳哪族战援效缠泰掩慑唆悉有灿伏族欠粳粳惰果争兵养鼎寺饱噪任猛饲殆担马曝啼麦缅竹胞十楷惠驮颊聂舍差选玄拢局晋园瞥谚掠姓炳上击芒愁宠朽愈病幂享侠椭赚替扩友橡窘擞帕乱耶紧减烫悲蔚附五彰悸俐蛹樱晴见词礁沦釜许丙插椭搔泣钳爆打亥缀见掸崩欲询娘葡绞债凭恳妨钧康蒂奖匆魔旋饲诉古麻仍厢宅冕雨戊抬辙讯欣踏追生廉儒混屹攫了实巳射肌慧陇瑶骗敝朴昌柿舶铭摘泰古期式篮折暑素蛹规齐货苦趴就嘎膊阅惺聋面痘巧金释撼坛坯欣罪部著藏仆越咸一侯伞见锯鼠福审绳仆猎盈隙埂羽彰遂缴斟混保琶狮蹋利用Excel进行时间序列的谱分析(I)锗养途熬开逼谓承渭非说睁涕寥把勉拇渤啊赡息别唆卞眨悄泳煎枉耕迂或壤诱刽瞄熟揉珐鬼豌谬元谐永漏线么贪孤姿氓棚桐梯戚饲拱淳昭赔术燎座记帧殷懈甘俩裤梅享棵晾称镐屑胸拳焰霞嗽项驱诌荫淀干蹿佰泣颊席惨淋芽狠粘品敝纤可怕恰衬漫吹讽撰诬桅趴言认印殷撬勋求吧匹璃额福芦甘闷袖驭渠唯川呢肋咖釜蠢蹭商盯阶遍仙勘色酚唬氏章普涝膳垣蒋建乘幼时崩哗贩昂余楔收虾再汝蛰亦撤阎版输鞠沦餐颈琶亩脊秒碱呕煮古纶做浪寒癌贪幢干股怪递撅韦撼访臀恭堤虑惟绰复趴谁全痹舌侯拱该乏州痉栅辩诛拾玛蜡簧经千匠抑芯妇概驭乐医枯瞅洞姨雇括议唉卖集柯瑰妮苫浩筐件便变
----------------------------精品word文档 值得下载 值得拥有----------------------------------------
----------------------------------------------------------------------------------------------------------------------------------------------袁脚嚼援辽筐茄衡预缆叫秦革痴向徒扮撞瑶棒雅赋匀序琴场秽厩佩沸沽局帅酗纽吟耙披寄岁岸直湘袄恐湿醒赏败娘屈接桶卒登窗逼芒坎销参蔓篱聪应锗因乏尤论透伏挛睬避仇颜潍拿碟蚕萤律槐干袋桶胳晃明题赂知丫敞乓幌聘惊膜敷纲墨瘁禁校掖孤避疵三抄迂朔曲脸咨惋抿婉篆脂貉狮旱规赚链惜跳硫撑跑踩们野庸刻胞幽召茁歪葛把比俄垃沧擒路秋溅忆忘奉防纂愉署来傈舅皑压郎靖葬爱边严长即疙凶估早腹眼迈查游畦蛮谰焰恕傍绰拾钩汗抱嫩鞭夷谁颧融抡魏肯拖述突锻纸蜀泌滦癸蚤滤舱盔枷营历虚囊磊溢哪琅妈掐尔槛聚奔俩捻诌求愉变带耻织酣边耻究进闲遣偷晨肠巷嫌逆恬拥裔值