-
利用
Excel
进行时间序列的谱分析（
I
）
在频域分析中
，
功率谱是揭示时间序列周期特性的最为有力的工具之一。
下面
列举几个
例子，分别从不同的角度识别时间序列的周期。
1
时间序列的周期图
【
例
1
】某水文观测站测得
一条河流从
1979
年
6
月到
1980
年
5
月共计
12
月份的断面平均
流量。试判断该河流的径流量变化是否具有周期性，周期长度大约为多少？
分析：假定将时间序列
x
t
< br>展开为
Fourier
级数，则可表示为
x
t
?
?
(
a
i
cos
2
?
f
i
t
?
b
i<
/p>
sin
2
?
f<
/p>
i
t
)
?
?
t
(1)
i
?
1
k
式中
f
i
为频率，
t
为时间序号，
k
为周期分量的个数即
主周期（基波）及其谐波的个数，
ε
t
为标准误差（白噪声序列）
。当频率
f
i
给定时，式（
1
）可以视为多元线性
回归模型，可以证
明，待定系数
a
i<
/p>
、
b
i
的最小二
乘估计为
2
N
?
i
?
?
x
t
cos
2
?
f
i
t
a
N
t
?
1
2
N
?
b
i
?
?
x
< br>t
sin
2
?
< br>f
i
t
N
t
?
1
这里
N
为观测值的个数。定义时间序列的周期图为
（
2
）
I
(
f
i
)
?
N
2
(
a
i
?
< br>b
i
2
)
，
i
?
1
,
2
,
?
,
k
(3)
2
式中
I
(
f
i
)
为频率
f
i
处的强度。以
f
i
为横轴，以
I
(
f
i
)
为纵轴，绘制时间序列的周期图，可以在
最大值处找到时间序列的周期。对于本例，
N
=
12
，
t
=1,2,
< br>…
,
N
，
f
i
=
i
/
N
，下面借助
Excel
，利
用上述公式，计算有关参数并分析时间序列的周期特性。
第一步，录入数据，并将数据标准化或中心化（<
/p>
图
1
）
。
图
1
录入的数据及其中心化结果
1
中心化与标准化的区别在于，
只需
将原始数据减去均值，
而不必再除以标准差。
不难想
到，中心化的数据均值为
0
，但方差与原始数据相
同（未必为
1
）
。
第二步，计算三角函数值
为了
借助式（
1
）计算参数
a
i
、
b
i
< br>，首先需要计算正弦值和余弦值。
取
< br>i
?
1
,
2
,
?
,
6
，则频率为
f
i
?
i
/
N
?
1
/
12
,<
/p>
2
/
12
,
?
,
6
/
12
（
图
1
）
。
将频率写在单元格
C3-C14
中（根据对称性，我们只用前
6
个）
，将中心化的数据转置
粘贴
于第一行的单元格
D1-O1
中，月份的序号写在单元格
D2-O2
中（与中心化数据对齐）
。
图
2
计算余弦值的表格
在
D2
单元格中输入公式“
=COS($$B$$1*$$D$$
2*C3)
”
，回车得到
0.866<
/p>
；按住单元格的右
下角右拉至
O3
单元格，
得到
f
=1
/12=0.083
，
t
=1,2,<
/p>
…
,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<
/p>
”
，即可得到正弦值，不过为了
计算过程
清楚明白，最好在另外一个区域给出结算结果（在
D17-O22
区域中，参见
图
3
）
。
图
3
计算正弦值的表格
第三步，计算参数
a
i
、
b
i
利用中心化的数据（仍然表作
x
t
）计算参数
a<
/p>
i
、
b
i
。首先算出
x
t
co
s2
π
f
i
t
和
x
t
sin
s2
π
f
i
t
。在
D9
单元格中输入公式“
=D1*D3
”
，回车得到
18.309
；按住单元格的右下角右拉至
O9
单元
格，得到
f
=1/1
2=0.083
，
t
=1,2,
…
,12
时的全部
x
t
cos2
π
f
i
t
值；加和得
39.584
，再除以
6
，即得<
/p>
a
1
=6.597
。在
D10
单元格中输入公式“
=D
1*D4
”
，回车得到
10.571<
/p>
；按住单元格的右下角右
拉至
O10
单元格，
得到
f
=
2/12=0.083
，
t
=1,2,
…
,12
时的全部
x
t
cos2
π
< br>f
i
t
值；
加和得
-365.25
，
再
除以
6
，得到
a
2
=-60.875
。其余依此类推。
将上面公式中的余弦值换成正弦值，即可得到
b
i
值（见下表）
。上面的计算过
程相当于
2
采用式（
2
）进行逐步计算。
第四步，计算频率强度
利用式（
3
）
，非常容易算出
I
(
f
i
)
值。例如
I
(
f
1
)
?<
/p>
N
2
(
a
1
?
b
1
2
)
?
6
*
(
6
.
597
2
?
170
.
213
2
)
?
174096
.
914
2
其余依此类推（见
图
4
）
。
图
4
计算频率强度
第五步，绘制时间序列周期图
利用图
4
中的数据，不难画出周期图（
图
5
）
。
200000
180000
160000
140000
频
率
强
度
120000
100000
80000
60000
40000
20000
0
0
0.1
0.2
0.3
0.4
0.5
0.6
频率
图
5
某河流径流量的周期图（
1979
年
6
月－<
/p>
1980
年
5
月
）
第六步，周期识别
关键是寻找频率的极值点或突变点。在本例中，没有极值点，但在
f
1
=1/12=0.0833
处，
频率强度突然增加（陡增）
，而此时
T
=1/
f
1
=12
，故可判断时间序列可能存在一个
12
月的周<
/p>
期，即
1
年周期。
3
【
例
2
】为了映证上述判断，我们借助
同一条河流的连续两年的平均月径流量（
1961
年
6
月－
1963
年
5
月）
。原始数据见下图（
< br>图
6
）
。
图
6
原始数据及部分处理结果
将原始数据回车时间序列变化图，可以初步估计具有
12
月变化周期，但不能肯定（
图
6
）
。
350
300
月
平
均
径
流
量
250
200
150
100
50
0
0
5
10
15
20
25
30
时间（月份序号）
图
6
径流量的月变化图（
1961
年
6
月－
1963
年
5
月）<
/p>
4
按照例
1
给出的计算步骤，计算参数数
a
i
、
b
i
，进而计算频率强度（结果将
图
7
）
。然后
绘制时间序列的周期图（
< br>图
8
）
。注意这里，
N
=24
，我们取
k
=12
。
图
7
参数和频率强度的计算结果
从图
8
中可以看出，频率强度的最大值（极值
点）对应于频率
f
1
=1/12=0.
0833
，故时间
序列的周期判断为
T
=1/
f
1
=
12
。
这与用
12
月的数据进行估计的结果是一致的，
但由于例
2
的
时间序列比例
1
的时间
序列长
1
倍，故判断结果更为可靠。
140000
120000
10000
0
频
率
强
度<
/p>
80000
60000
40000
20000
0
0
0.
1
0.2
0.3
0.4
0.5
0.6
频率
图
8
某河流径流量的周期图（
1961
年
6
月－<
/p>
1963
年
5
月
）
2
时间序列的频谱图
【
例
3
】首先考虑对例
1
的数据进行功率谱分析。例
1
的时间序列较短
，分析的效果不
佳，但计算过程简短。给出这个例子，主要是帮助大家理解
Fourier
变换过程和方法。为了
进行
Fourier
分析，需要对数据进行预处理。第一，将数据中心化，即
用原始数据减去其平
均值。中心化的数据均值为
0
，我们对中心化的数据进行变换，其周期更为明显。第二，由
于
Fourier
分析通常采用所谓快速
Fourie
r
变换（
Fast Fourier Transformat
ion
，
FFT
）
，而
FFT
n
要求数据必须为
2
个，这里
n
为正
整数（
1,2,3,
…）
，而我们的样
本为
N
=12
，它不是
2
的某
5
个
n
次方。因此，在中心化的数据后面加上
4
个
0
，这样新的样本数
为
N
′
=12
＋
4
＝
16
＝
2
4
个，这才符合
FFT
的需要（
图
9
）
。下面，我们对延长后的中心化数据进行
Four
ier
变换。
图
9
数据的中心化与“延长”
第一步，打
开
Foureir
分析对话框
沿着主菜单的“工具（
Tools
）
”→“数据分析（
Data Analysis
）
”路径打开数据分析选项
框（图
10
）
，从中选择“傅立叶分析（
Four
ier Analysis
）
”
。
图
10
在数据分析选项框中选择
Fourier
分
析
第二步，定义变量和输出区域
确定之
后，弹出傅立叶分析对话框，根据数据在工作表中的分布情况进行如下设置：
将光标置于“输入区域”对应的空白栏，然后用鼠标选中单元格
C1-C
17
，这时空白栏
中自动以绝对单元格的形式定义中心化数据的
区域范围（即
$$C$$1-$$C$$17
）
。
选中“标志位于第一行”
。
选中输出区域，定义范围为
D2-D17
（
图
11
）
。
注意：如果输入区域的数据范围定义为
C2
-C17
，则不要选中“标志位于第一行”
，这与
回归分析中的原始变量定义是一样的（
图
12
）
。如果不定义输出区域范围，则变换结果将会
自动给在新的工作表组上。这一点也与回归分析一样。
6
图
11
选中“标志位于第一行”与数据输入范围的定义
图
12
不选中“标志位于第一行”与数据输入范围的定义
第三步，结果转换
定义数据输入－输出区域完成之后，确定，立即得到
Fourier
变换的结果（
图
13
）<
/p>
。
图
13
傅立叶变换的结果
7
变换的结果为一组复数，
相当于将
f
(
t
)
变成了
F
(
ω
)
，
实际上是将
x
t
变成了
X
T
(
f
)
。
我们知
道，有了
f
(
t
)
的象函数
F
(
ω
)
就
可以计算能量谱密度函数
S
(
ω
)
，即有
S
(
?
)
?
F
(
?
)
F
(
?
)
?
F
(
?
)
(4)
相应地，有了
X
T
(
< br>f
)
也就容易计算功率谱（密度）
P
(
f
)
?
2
X
T
(
f
)
T
2
(5)
式中
f
表示线频率，与角频率
ω
的转换关系是
ω
=2
π
/
T
，这里
T
为数据区间长度。
如果将
X
T
(
f
)
表作
X
T
(
f
)=
A
+
jB
（这里
A
为实部，
B
为虚部）
，则有
X<
/p>
T
(
f
i
)
?
X
T
(
f
i
)
X
T
(
f
i
)
?
(
A
i
?
jB
i
)(
A
i
?<
/p>
jB
i
)
?
A
i
2
?
B
i
2
（
6
）
因此这一步是要分离变换结果的实部与虚部。逐个手动提取是非常麻烦而且容易出错
的，可以利用
Excel
有关复数计算的命令。
提取实部的
Excel
命令是
imre
al
。在
H2
单元格
< br>中输入命令
“
=IMREAL(D2)
< br>”
（这里
D2
为变换结果的第一
个复数所在的单元格）
，
回车得到
第一
个复数的实部
0
；点中
H2
单元格的右下角，揿住鼠标左键下拉至
H17
，得
到全部的实
部数值。提取虚部的命令是
imaginary
。在
I2
单元格中输入公式“
=IMAGINARY(D2)
”
，回
车得到第一个复数的虚部
0
；下拉至
I17
，得到全部的虚部数值。根据式（
5
）
、
（
6
）
，功率谱
密度的计算公式为
2
A
i
2
?
B
i
2
P
(
f
i
)
?
（
7
）
T
考虑到本例中
T
=
N
=16
，在
J2
单元格中输入公式“
=(H2^2+I2^2)/16<
/p>
”
，回车得到第一个功率
谱密度
0
；下拉至
J17
，得
到全部谱密度数值（
图
14
）
。基于
FFT
的谱密度分布是对称的，
可以看出，以
J10
所在的
3105.28
为对称点，上下的数值完全对称。
图
14
功率谱密度的计算结果
8
-
-
-
-
-
-
-
-
-
上一篇：初二英语上册短语大全
下一篇：英语必修四短语