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

利用Excel进行时间序列的谱分析（I）在频域分析中，功率谱是揭示时间序列周期特性的最为有力的工具之一下面列举几个例子，分别从不同的角度识别时间序列的周期1时间序列的周期图【例1】某水文观测站测得一条河流从1979年6月到1980年5月共计12月份的断面平均流量试判断该河流的径流量变化是否具有周期性，周期长度大约为多少？分析：假定将时间序列xt展开为Fourier级数，则可表示为kxt二'（aicos2fitbisin2二£t）亠-t⑴i=4式中fi为频率，t为时间序号，k为周期分量的个数即主周期（基波）及其谐波的个数，为标准误差（白噪声序列）当频率fi给定时，式（1）可以视为多元线性回归模型，可以证明，待定系数ai、bi的最小二乘估计为?=一送xtcos2兀£tNtm2N?xtsin2二ftNy这里N为观测值的个数定义时间序列的周期图为N22I（fiH2（aibi）,i=1,2,,k⑶式中I（fi）为频率fi处的强度以fi为横轴，以1（町为纵轴，绘制时间序列的周期图，可以在最大值处找到时间序列的周期对于本例，N=12,t=1,2,…,N,fi=i/N,下面借助Excel，利用上述公式，计算有关参数并分析时间序列的周期特性。

第一步，录入数据，并将数据标准化或中心化（图1）ABCDIE■F1月份月份序号#径流量航中心化周期序号r频率£26119021.1421CL083333|37224879.14220.1666674_83553384.14230.25594286117.14240.333333610520435.14250.4166677116106-62.85860.5812765.4-103.45870.583333918迢9-121.958S0.666667102933.8-135.05890.751131051.8-117.053100.833333121'4：1178.4-90.45B110.91666713512163-5.858121图1录入的数据及其中心化结果中心化与标准化的区别在于，只需将原始数据减去均值，而不必再除以标准差不难想到，中心化的数据均值为0，但方差与原始数据相同（未必为1）第二步，计算三角函数值为了借助式（1）计算参数ai、bi，首先需要计算正弦值和余弦值将频率写在单元格粘贴于第一行的单元格取i=1,2,,6，则频率为fj=i/N=1/12,2/12,,6/12(图1)。

C3-C14中（根据对称性，我们只用前6个），将中心化的数据转置D1-O1中，月份的序号写在单元格D2-O2中（与中心化数据对齐）1CDEFGHIJELNN0中心牝藝卅21.1416779,14167384-1417117.141735.1^167-62.8S33-103,46-121,558-135,05S-117.058-90.4583-S.85832頻率f1234567891011123山083333H86&65006000-0.5007866-1.000■686&-0.500Oi00065000..8661■1D,166^670.500-o.500-：.0000.5001..0000.500-0.500-1.000-0.5000.500150.250.000-1..0000.0001.0000.000-1.00060001.0000.000-1.0000..0001—&0.333333-0-5001-000-0.500*5001..COO-0.500-0.500

在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单元格是2n=6.28319（图中未能显示）在上面的计算中，只要将公式中的“COS”换成“SIN”，即可得到正弦值，不过为了计算过程清楚明白，最好在另外一个区域给出结算结果（在D17-O22区域中，参见图3）16170,0833331B0.166667190s25200.333333210.4166670.522123450”500O',@66L0000,B6SCL5000.8660.8660.000-0.866~~^6.8661.0006DOC-1ooo0,0001.0000.E&678660.ooo0.866-0-B665500-0i866L000-Q-866CL5000.000Oiooo0.0000-oooCL00010.57168,539384.142101,44317.571IE.68.5390.000-101.448-30.43421.1^20,0(W-384,1420,ooo35.14218.309-BBL5390..ooo10L44£-30.43410.571-68,5$9384.142-10L44817.5710.Q000.ooo0.oooCL00060000.5S33330.66666T(X1&0.3333330.915667678910'11120,000-a500-a邂-L000-a3&e-0.5000.00050000.8660.8660.0002^866-0.8660.0000,000-1,ooo0.000LOOO6000-L0000,000OuOOO0.866-0.8&60.0000.866-0.8660.000notw-0.5000.8&6-1.0000„8&5-0..5000.00Q0-0000.0000..000此000(1OOOm0000.000aooci51.725105.619135,05810L375452290,0000-000-89.598-105.619IX000101.37578.3390.000aooo103.4580..000-135,053HOOO90.4580P000也000-89.598105.619山000-101.37578.3390..ooo04OCI51,729-105.619m-10L3754&r2290.0000.0000.0000.000山0000.ooo0.ODO0.000图3计算正弦值的表格第三步，计算参数ai、bi利用中心化的数据（仍然表作xj计算参数ai、bi。

首先算出Xtcos2nt和xtsins2nt在D9单元格中输入公式“=D1*D3”，回车得到18.309;按住单元格的右下角右拉至O9单元格，得到f=1/12=0.083,t=1,2,…,12时的全部Xtcos2nt值；加和得39.584，再除以6,即得a1=6.597在D10单元格中输入公式“=D1*D4”，回车得到10.571;按住单元格的右下角右拉至O10单元格，得到f=2/12=0.083,t=1,2,…,12时的全部Xtcos2nt值；加和得-365.25,再除以6,得到a2=-60.875其余依此类推bi值（见下表）上面的计算过程相当于将上面公式中的余弦值换成正弦值，即可得到采用式（2）进行逐步计算££hIP39.5846-5971021.281170-213-365.250-60.875~~-60.535-1008990.10015.017-269.000-44.833270.55045-09213.7702.29541,3166-8S&267,31944.553-222.900-37.150Q.0000.000第四步，计算频率强度利用式（3）,非常容易算出l（fi）。