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

《利用Excel进行时间序列的谱分析.docx》由会员分享，可在线阅读，更多相关《利用Excel进行时间序列的谱分析.docx（18页珍藏版）》请在三一办公上搜索。

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

2、多元线性回归模型，可以证 明，待定系数a.、b.的最小二乘估计为2 y=一 x cos 2写 t（2）N 乙 sin2Kftt=1t=1这里N为观测值的个数。定义时间序列的周期图为I（f ） = N （a2 + b2），i = 1,2, ,k（3）i 211式中If）为频率f.处的强度。以f,为横轴，以If）为纵轴，绘制时间序列的周期图，可以在 最大值处找到时间序列的周期。对于本例，N=12，t=1,2,N，f=i/N，下面借助Excel，利 用上述公式，计算有关参数并分析时间序列的周期特性。第一步，录入数据，并将数据标准化或中心化（图1）。ABCDEF1月份月份序号丈径流量也中心化周期序号f

3、频率26119021. 14：210. 0S333337224879. 14220. 1666674S3553304. 14230. 25594286117. 14=240. 333333610520435. 14250. 4166677116106-62. 85860. 5312765. 4-103. i5S70. 58333391S46.9-121.958S0. 666667102933. 8-135.05890. 751131051. 0-117.050100. 833333124117S. 4-90. 4：58110. 91666713512163-5. 858121图1录入的数据及其

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

5、7-62.游83-103.花-L21.958-135. 05S-117- 05STO, 4够3-5. S5E32频率f1234&6I891112:.: 83333JU：:.J.1.-1.:.Ji: :-:Ji.: ,i.：.L ：1a0.1666670- 500P. 500-1-000-0. 500Cl 5001. 000LL 500-0.500-1. D00-CL 5000. E00150- 250. 0001. . 7.1. :L 0001.: JI.-1.:CL 0001.0. 000Z. 7JI.7.1. : .116CL 333333-0. 500-o. 5如LODO-0, 500-

6、0, 51X1L DOO-0, 500-Db5D0L DOO-CL 51)(1-o.50D171. Ill-1 . .r :.r.l.:l.-l. 190- 251. ii0. 000II0. ooo1. 0i：!UD. OCidj-LOOO0. ooo1-0000. U00-1. 0LU0. ooo20J.0 8667. Ji: :1. 7 JI.7. Ji: :7. Ji: :1.J.J. . . JJ. *：.L ：J210-4166670.500Y, 866L ODOY, 8660. 500OS-0.5000. 866-LOODIX 866-0. 5000. ooo220,50. DO

7、OAWDDOOD0. 000o. ooo0. ooo0. ooo23J. :*:10. &716S.眺30.142101.i4S1?.571tk ooo&L i29IOS.619135.涮101L3?5g 22&0. ooo0.66666718.3096S. 5龄(L ODO-101.44S-30. 434AWD-89.59B -105. 6L9WL3T578. 3於0. ooo250.7521. 142U. JOO-38 A 142u. JOO35.1420. OLI JLOS. 4580. oao -135- 05S：L OODgo.蛟0. ooo260.83333319.309-68.

8、539CL 000101.4Q3-30. 434o. ooo-89.598105. 619侦Mi-101.37578. 3390. ooo270B 91666710B 571-68. 5393E3.142-101.448IT. 571儿ooo51.729 -105. 619135.058-101.37545, 2290. ooo2B10. ooo0. 0000. ooo0. oooU. ILIOO0. OLI J0. 000o. oaoLL 0000. OOD0. am0. ooo图3计算正弦值的表格第二步，计算参数a、 b.利用中心化的数据（仍然表作气）计算参数a、b。首先算出气cos2f

9、和气sins2f。在 D9单元格中输入公式“ =D1*D3”，回车得到18.309；按住单元格的右下角右拉至O9单元 格，得到f=1/12=0.083，t=1,2,“,12时的全部气cos2f值；加和得39.584，再除以6，即得 a1=6.597。在D10单元格中输入公式“=D1*D4”，回车得到10.571；按住单元格的右下角右 拉至O10单元格，得到f=2/12=0.083, t=1,2,“,12时的全部气cos2f 值；加和得-365.25,再 除以6，得到a2=-60.875。其余依此类推。将上面公式中的余弦值换成正弦值，即可得到b.值（见下表）。上面的计算过程相当于采用式（2）进行

10、逐步计算。39. 5346. 5971021.231170_ 213-365. 250-60. 875-60.53510- 08990. 10015. 017-269.000-44.833270. 55045. 09213.7702- 29541. 3166. 886267.31944-553-222. 900-37. 1500.0000- 000第四步，计算频率强度利用式（3），非常容易算出/（/；）值。例如I（f ） = （a2 + b2） = 6*（6.5972 +170.2132） = 174096.9141211其余依此类推（见图4）。1 QR|_ ST频率胃必0. 0836. 59

11、7170. 213174096.9140. 167-60. 875-10. 08922845.3450. 25015.017-44. 83313413.1680. 33345.0922. 29512231.1520. 4176. S8644. 55312194.4230. 500-37.1500. 0008280.735图4计算频率强度第五步，绘制时间序列周期图利用图4中的数据，不难画出周期图（图5）。图5某河流径流量的周期图（1979年6月一 1980年5月）第六步，周期识别关键是寻找频率的极值点或突变点。在本例中，没有极值点，但在九=1/12=0.0833处， 频率强度突然增加（陡增），而

12、此时T=1/f=12，故可判断时间序列可能存在一个12月的周 期，即1年周期。【例2】为了映证上述判断，我们借助同一条河流的连续两年的平均月径流量1961年6月一1963年5月）。原始数据见下图（图6）。ABCDEFGH16. 28319中心化致据径流量39. 545S3168.545858. 5458366. 545832月份径流量中心化数据频率f12343114939. 545830. 0416670. 9660. 8660. 7070. 50042278168.54580. 0833330. 8660.5000.000-0.5005316858. 545830.1250. 7070. 0

13、00-0. 707-1. 0006417666. 545830.1666670. 500-0. 500-1.000-0. 5007512212. 545830. 2083330. 259-0. 866-0. 7070. 500S672-37. 45420. 250. 000-1.0000. 0001. 0009750-59. 45420. 291667-0. 259-0. 8660. 7070. 50010843. 6-65. 85420. 333333-0. 500-0. 5001. 000-0.50011928-81. 45420. 375-0. 7070. 0000. 707-1. 00

14、0121041. 4-68. 05420. 416667-0. 8660. 5000. 000-0.500131156-53. 45420. 458333-0. 9660.866-0. 7070. 500141265-44. 45420. 5-1. 0001. 000-1. 0001. 000151315242. 545330. 54166738.198145.96541.398S3. 2731614229119. 54580. 58333334. 24884. 2730. 000-33. 2731715323213.54580. 62527. 9630.000-41.398-66. 5461

15、81614434. 545830. 66666719. 773-84. 273-58. 546-33. 273191719181. 545830. 70833310. 235-145. 965-41. 39833. 273201884-25. 45420. 750. 000-168. 5460. 00066. 546211935. 8-73. 65420. 791667-10. 235-145. 96541. 39833. 273222029.4-80. 05420. 833333-19. 773-84. 27358. 546-33. 273232125.1-84. 35420. 875-27

16、. 9630. 00041. 398-66. 546242225.1-84. 35420. 916667-34. 24884.2730.000-33. 273252350-59. 45420. 958333-38.198145. 965-41. 39833. 273262489. 5-19. 95421-39. 546168. 546-58. 54666. 546图6原始数据及部分处理结果将原始数据回车时间序列变化图，可以初步估计具有12月变化周期，但不能肯定（图 6）。051015202530时间（月份序号）350300量流径均平月250200150100500图6径流量的月变化图（1961

17、年6月一 1963年5月）按照例1给出的计算步骤，计算参数数a、b,进而计算频率强度（结果将图7）。然后 绘制时间序列的周期图（图8）。注意这里，N=24，我们取k=12。| AE匚AF1 AG频率f可了与少）10. 0417-6.6800-S.25391352. 9950. 08337. 9876103.694129808. 3420. 125014. 23500. SSS02441.0060. 1667-26.162511.640S9839. 8230. 20332. 83666. 40EE5S9. 4190. 2500-15. 1667-1.30832780. 8740. 2917-1.

18、98979.63551161. 6190. 3333-2.6292-16.2S853266. 7310. 3750-13. 1183-11.99S33792. 5990. 41676. 80404.4673795. 0220. 458316. 9665-18.36407501.1780. 5000-6.07500.0000442. 868图7参数和频率强度的计算结果从图8中可以看出，频率强度的最大值（极值点）对应于频率4=1/12=0.0833，故时间 序列的周期判断为T=1/f=12.这与用12月的数据进行估计的结果是一致的，但由于例2的 时间序列比例1的时间序列长1倍，故判断结果更为可靠。

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

20、3,），而我们的样本为N=12，它不是2的某个n次方。因此，在中心化的数据后面加上4个0,这样新的样本数为N =12+4=16=24 个，这才符合FFT的需要（图9）。下面，我们对延长后的中心化数据进行Fourier变换。ABC1月份径流量中心化2119021.1423224879. 14243553384. 14254286117. 1426S20435. 14276106-62. 858E765. 4-103. 4589846. 9-121.95310933. 3-135. 058111051. 8-117. 058121178. 4-90. 4581312163-5. 858141301

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

22、出区域，定义范围为D2-D17 （图11）。注意：如果输入区域的数据范围定义为C2-C17，则不要选中“标志位于第一行”这与 回归分析中的原始变量定义是一样的（图12）。如果不定义输出区域范围，则变换结果将会 自动给在新的工作表组上。这一点也与回归分析一样。图11选中“标志位于第一行”与数据输入范围的定义图12不选中“标志位于第一行”与数据输入范围的定义第三步，结果转换定义数据输入一输出区域完成之后，确定，立即得到Fourier变换的结果（图13）。ABC|_ D _1 E|_ FG1月份径流量中心化Foim.改换结果2119021. 142013224379. 14293土 01384461

23、9619-381. 55316&656499i43迎384. 142-296. 348675854492-579. 70485245&Q13i54286117. 142-302. 505754499164-272. 8702594TS912i6520435. 142-269+90.li76106-62. S5S-202.5关180343623+181. 7199S5749509iST&5. 4-103. 45S-1.76799081217553+214. 5784SaST732igS46. 9-121.涮199.884090223169+213. 60374523359110933. -135

24、. 05S222.111051. 8-117. 058199. S84090223168-213. 603刊5%8589i1211T8. 4-90.如R-1. 76799081217561-214. 5784SaST732i1312163-5. S5S-202. 592180343624-181. 7199S5749509i1413026990.li15140-302.505754499164+272. 8702594TS913i16150-296. 348675854491+579. 70485245601411716a030.013844619619+381. 55316665649Ti图

25、13傅立叶变换的结果变换的结果为一组复数，相当于将人。变成了 F(),实际上是将孔变成了 Xf 我们知 道，有了人。的象函数F(时就可以计算能量谱密度函数S(/)，即有S(3) = F ()F ()=|F ()|2相应地，有了乂朋也就容易计算功率谱(密度)P( f)=Xt (f )|2T式中f表示线频率，与角频率刃的转换关系是刃=2兄T，这里T为数据区间长度。如果将XTf)表作XT(f)=A+jB (这里A为实部，B为虚部)，则有X (f )|2= X (f )Xf = (A + jB )(A - jB ) = A2 + B2(6)T iT i T ii i i i i i因此这一步是要分离变

26、换结果的实部与虚部。逐个手动提取是非常麻烦而且容易出错 的，可以利用Excel有关复数计算的命令。提取实部的Excel命令是imreal。在H2单元格 中输入命令“=IMREAL(D2)”(这里D2为变换结果的第一个复数所在的单元格)，回车得到 第一个复数的实部0；点中H2单元格的右下角，揿住鼠标左键下拉至H17，得到全部的实 部数值。提取虚部的命令是imaginary。在I2单元格中输入公式=IMAGINARY(D2)”，回 车得到第一个复数的虚部0；下拉至I17,得到全部的虚部数值。根据式(5)、(6)，功率谱 密度的计算公式为(7)考虑到本例中T=N=16，在J2单元格中输入公式“=(H

27、2A2+I2A2)/16，回车得到第一个功率 谱密度0；下拉至J17，得到全部谱密度数值(图14)。基于FFT的谱密度分布是对称的， 可以看出，以J10所在的3105.28为对称点，上下的数值完全对称。DEFGH I IJ I1Founer换结果imrealimaginaryW200003930.0138446161-301. 5531666564991930. 013S-381.55363156. 794-296.348675S54492-579. 704852456013i-296. 349-579. 70526492. 525-302.505754499164-272. 87025947

28、8912i-302. 506-272. 8710372. 996-269+90.11-26990. 15029. 9337-202.5921S0343623+1S1. 7199857495091-202. 592181. 724629. 1093-1.76799031217553+214. 57848087732i-1. 76799214. 57352377. 941g199.884010223169+213.603745238591199. 8841213. 60375348. 76310222.9222. 903105. 285348. 76311199.88409022316S-213.

29、 603745238589i199. 8841-213. 60412-1.767990S1217561-214. 57848087732i-1. 76799-214. 5782877. 94113-202.5921S0343624-1S1. 7199857495091-202. 592-LSI. 724629. 10914-269-90.H-269-90. 15029. 93315-302.5057544164+272. S70259478913i-302. 506272. 370310372. 9916-296.34867F8F4491+F79. 7048524560141-296. 349

30、579. 704926492. 5217930.0138446161+301. 553166656497i930. 013S381.553263156. 79图14功率谱密度的计算结果第四步，绘制频谱图线频率f可以表作f = i / T = i / N， i = 0,1,2, , N-1密谱率功AAi / inL1L /1 jLiiiij 170000600005000040000300002000010000000.20.40.60.81频率f显然 f =0/16=0, f=1/16=0.0625, f2=2/16=0.125，5=15/16=0.9375。在 Excel 中，容易 计算频率

31、的数值。将频率与功率谱对应起来（图15），就可以画出频谱图。如果补上最后 个频率数值f16=1及其对应的功率0，则可画出完全对称的谱图（图16）。KH L M 11序号I频率/功率谱密度分布户（壬）2000310. 062563156. 78564420. 12526492. 51585530. 187510372. 99438640. 255029. 938125750. 31254629.109047860. 3752S77.94064970. 43755348. 76309410K0. 53105. 275&251190. 562553犯 76309412100. 6252877.940

32、6413110. 6S754629. 10904714120.花5029. 93812515130. 312510372.9943S16140. 87526492. 5158517150. 937563156. 78564图15功率谱密度与频率的对应关系图16对称分布的频谱图由于功率谱图的对称性，画出整个谱图实在没有必要，因此，在实际工作中，通常只画 出左半边（图17）。00.10.20.30.4频率f图17实用的频谱图700006000050000400003000020000100000.50.6第五步，功率谱分析频域分析的主要目标之一是判断时间序列是否存在周期性。从图17可以看出，功率最

33、 大点对应的频率是f=0.0625，该频率对应的周期长度为16。可见，在时间序列较短的情况 下之间用功率谱图寻找时间序列的周期不如周期图准确。另外可以初步估计数据的性质。在图17中，去掉第一个0点，剩余的点一般呈幂指数 分布（在双对数坐标图上点列具有直线趋势），可以拟合幂指数函数如下：（8）P (f) = 935.6fR2 = 0.92330.10.01频率f图18功率P（f）与频率f的双对数坐标图结果得到功率谱指数俭1.49521.5。功率谱指数与时间序列的Hurst指数具有如下关系(9)据此估计Hurst指数约为0.25。我们知道，Hurst指数介于01之间，当H0.5时，表明时 间序列存

34、在正的自相关，意味着系统演化具有持久性；当H0时，表明时间序列具有负的自相关，意味着系统演化具有反持久性；当H=0时，表明时间序列不存在自相关，过去与 未来无关。对于这条河流的径流量而言，H=0.250.5，表明时间序列具有反持久性：过去的 增量意味着今后的减少，过去的减少意味着未来的增加。因此，径流量必然周期性的变化。【例4】下面对前述例2的数据进行Fourier变换，方法与例3相同，但由于N=24，我 们取7=32=25。也就是说，对于中心化的数据，要在后面添加8个0作为补充点数。基于 FFT的变换结果如下(见图19)。ACDEI FI G1月份径流量中心化Fourier换结果211493

35、9. 54583032278168. 545S175.697209544&61-182. 79376905547114316S58. 54583875.211791457S11+481. 96764957647215417666. 54583-135. 552842511169-1003. 1725467126516512212. 5458312S.296969619&7+113. 2782612482317672-37.454298.3772036360729-336. 06927332013718750-59.4542-184. 362200004434-70. 78851551733H9843.6-65.8542121.43913435022-92. 2294676618315110928-81.4542-15. 7000000000001-1821111041. 4-68.0542-19.1752816836816+39. 32466038537991121156-53.454265.0634777431574-258. 0121759909451131265-44.4542-93. 215776868303-16. 416749686007H141315242. 545839.50303038032981-213. 1217387517711514