时间序列分析及相空间重构

重复这一过程并测量相对于不同时间的各 延迟量,就可以产生出许多这样的点,它可 以将原系统的许多性质保存下来,即用系 统的一个观察量可以重构出原动力系统 模型,可以初步确定原系统的真实信息。
相空间重构例
Henon 映射
2 xn +1 = 1 − 1.4 xn + yn yn +1 = 0.3 xn
50 45 40 35 30 25 20 15 10 5 0 -20
-15
-10
-5
0
5
10
15
20
如果只观测到变量x的值，利用x作相空间重构 取延迟时间为9，嵌入维数为3 即令(x(1),y(1),z(1))=(x(19),x(10),x(1)) (x(2),y(2),z(2))=(x(20),x(11),x(2)) (x(3),y(3),z(3))=(x(21),x(12),x(3)) ……
N −τ n =1
∑ (x
1 N
N n =1
n +τ N
− x )( xn − x) ,
( xn − x ) 2 ∑
n =1
1 其中x = N
∑x
n
选择使得自相关函数C(τ)第一次为零时的τ的值为延迟时间
平均互信息法
对于时间序列xn , 定义平均互信息为 P ( xn , xn +τ ) I (τ ) = ∑ P ( xn , xn +τ ) log 2 , P ( xn ) P ( xn +τ ) n =1 其中P ( xn )，P ( xn , xn +τ )为概率
预测效果评价
为了检验预测的精确性，可以比较预测值 与实际观测值之间的差。一次预测可能 较好或较差，偶然性较大。为了克服这 种偶然性，可以取多个点的预测误差的 平均。
设xT +1 , xT + 2 , ⋯ , xT + p的预测值为 yT +1 , yT + 2 ,⋯ , yT + p , 定义均方根误差为 1 p 2 RMSE = ∑ ( xT +i − yT +i ) p i =1
m
m
>R
称X η ( n )为X n的虚假邻点。，其中阈值R在10和50之间。
对于实际的时间序列，还需补充以下虚假邻点标准 Xη (n) − X n
m +1
1 / N
( xn − x ) 2 > 2 ∑
n =1
N
虚假邻点法确定嵌入维数
对实测时间序列, m 从2 开始,取R = 30 , 计算虚假最近邻点的比例,然后增加m ,直 到虚假最近邻点的比例小于5 % 或虚假 最近邻点不再随着m 的增加而减少时,可 以认为此时的m 为合适的嵌入维数。
如果RMSE比较大，则说明预测效果不好。但是RMSE 和观测序列的数值大小有关，为克服这一问题，我们定义 正规化均方根误差NRMSE
NRMSE = RMSE / σ , 1 1 2 其中，σ = ∑ ( xT +i − x) , x = ∑ xT +i p i =1 p i =1
2 p p
若NRMSE接近于1，则预测效果不好 若NRMSE接近于0，则预测效果较好
1 X (T + 1) = K
^
∑ X (T
k =1
K
k
+ 1)
这是向量表达式，取第一个分量得
1 x(T + 1) = K
^
∑ x(T
k =1
K
k
+ 1)
局部线性预测法
局部线性预测模型为
x(T + 1) = g ( X (T )) = c0 + c1 x(T ) + c2 x(T − τ ) + ...cm x(T − (m − 1)τ ) + e 其中e为随机误差，ci为待定系数。
0
5
10
重构后的相图(x-y)
原系统的相图(x-y)
50 45 40 35 30 25 20 15 10 5
-15 -10 -5 0 5 10 15
0 -30
20
-20
-10
0
10
20
重构后的相图(y-z)
原系统的相图(y-z)
50
20
45
15 10 5 0 -5
40 35 30 25 20 15
8 取σ = 10，r = 28, b = 3 初值x0 = 15.34, y0 = 13.68, z0 = 37.91
Lorenz系统的吸引子(x-y-z）
20
10
0
-10
-20 60 40 20 -20 0 -40 0 40 20
Lorenz系统的吸引子(x,y相图）
30 20
10
0
-10
τ称为延迟时间，m称为嵌入维数。由x(n)构造X(n) 称 为相空间重构。
相空间重构法基本思想是: 系统中任一分量的演化都是由与之相互作 用着的其它分量所决定的,因此这些相关 分量的信息就隐含在任一分量的发展过 程中。
为了重构一个等价的状态空间,只需考察一 个分量,并将它在某些固定的时间延迟点 上的测量作为新维处理,它们确定了某个 多维状态空间中的一点.
非线性时间序列预测
基本思想 设时间序列来自确定性系统 X(n)=F(X(n-1)),F(.)为连续函数。 若 X(n)和X(j)距离很小，则F(X(n))和F(X(j))距 离也应很小，即X(n+1)和X(j+1)间的距离很 小，从而 可以用X(j+1)作为X(n+1)的预测值。
基本方法 局域预测法
预测效果的另一个评价标准是相关系数
相关系数r =
∑x
i =1 p i =1
p
T +i
yT +i − ∑ xT +i ∑ yT +i
i =1 i =1
p
p
( xT +i − x) 2 ∑
( yT +i − y ) 2 ∑
i =1
p
若r接近于1，则预测效果较好
上证指数预测
数据文件000001.day 1990.12.19-2008.06.19 共4292个记录 每一条记录的长度为40字节: 1-4字节为日期 5-8字节=开盘指数*1000 9-12字节=最高指数*1000 13-16字节=最低指数*1000 17-20字节=收盘指数*1000 21-24字节=成交金额(元)/1000 25-28字节=成交量(手) 其余12字节未使用 读取数据matlab文件为readdata.m
局部平均预测法 局部线性预测法 局部多项式预测法 全域预测法 神经网络 小波网络 遗传算法
局部平均预测法
设时刻T的状态向量为
X(T)=(x(T),x(T-τ),…x(T-(m-1)τ)) X(T)=(x(T),x(Tx(T-(m-1)
找X(T)的最近邻点X(T1),…X(TK), 以X(T1+1),…X(TK+1)的平均值作为X(T+1)的预测值
m m +1 m m +1
, 如果 X η ( n ) − X n
比
X η ( n ) − X n 大很多，认为X η ( n )为X n的虚假邻点。
上面的距离差由于和时间序列数据的大 小有关，不太容易确定虚假邻点。实际采用 相对度量法
若( X η ( n ) − X n
m +1
− Xη (n) − X n ) / X η (n) − X n
^
仍设X(T1),…X(TK)为X(T)的K个邻近点， 确定Ci的方法：最小二乘法 即 求Ci使得
∑
K
k =1
x (T k + 1) − x (T k + 1) =
^
2
∑
K
2
k =1
g ( X (T k )) − x (T k + 1) 最小
记y = [ x(T1 + 1), x(T2 + 1),⋯ , x(TK + 1)]T C = [c0 , c1 , ⋯ , cm ]T 1 x(T1 ) x(T1 − τ ) ⋯ x(T1 − (m − 1)τ ) 1 x(T ) x(T − τ ) ⋯ x(T − (m − 1)τ ) 2 2 2 A= ⋮ ⋮ ⋮ ⋮ ⋮ 1 x(TK ) x(TK − τ ) ⋯ x(TK − (m − 1)τ ) 则C = ( AT A) −1 AT y
20
20
10
10
0
0
-10
-10
-20 20 10 0 -10 -20 -20 0 -10 10 20
-20 60 40 20 -20 0 0 -40 20
重构后的相图(x-y-z)
原始系统相图(x-y-z)
30
20
10
0
-10
-20
-15
-10
-5
0
5
10
15
20
-30 -20
-15
-10
-5
^
仍以最小二乘法确定系数
记y = [ x(T1 + 1), x(T2 + 1),⋯ , x(TK + 1)]T C1 = [c00 , c10 ,⋯ , cmm ]T 1 x(T1 ) x(T1 − τ ) ⋯ x(T1 − (m − 1)τ ) x(T1 + 1) 2 ⋯ x(T1 − (m − 1)τ ) 2 1 x(T2 ) x(T2 − τ ) ⋯ x(T2 − (m − 1)τ ) x(T2 + 1) 2 ⋯ x(T2 − (m − 1)τ ) 2 A1 = ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ 1 x(TK ) x(TK − τ ) ⋯ x(TK − (m − 1)τ ) x(TK + 1) 2 ⋯ x(TK − (m − 1)τ ) 2 则C1 = ( A1T A1) −1 A1T y
上证指数预测文件为shangzhen1.m 相关文件为readdata.m juli.m dataconstruct1.m reconstruct1.m
多变量时间序列
以两个变量为例说明多变量情形 设给定时间序列x(n),y(n) x(n)的嵌入维数为m1,延迟时间为τ1 y(n)的嵌入维数为m2,延迟时间为τ2
-10
10
-15 -20 -20
5 0 15 -20
-15
-10
-5
0
5
10
20 -15
-10
-5
0
5
10
原系统的相图(x-z) 重构后的相图(x-z)
如何确定延迟时间和嵌入维数？
延迟时间间隔τ的选取
主要方法 线性自相关函数法 平均互信息法
线性自相关函数法
定义自相关函数为
C (τ ) =
1 N
由时间序列恢复原系统最常用的方法利用Takens 的 延迟嵌入定理 对于一个非线性系统,通过观测,可以得到一组测量值 x ( n)，n=1,2,…N 利用此测量值可以构造一组m 维向量 X( n) = ( x ( n) , x ( n -τ) , ⋯,x ( n -( m - 1)τ) ) n= ( m - 1)τ) +1,…N 如果参数τ, m 选择恰当,则X( n) 可描述原系统。
时间序列的定义
按照时间的顺序把事件变化发展的过程 记录下来就构成了一个时间序列。对时 间序列进行观察、研究，找寻它变化发 展的规律，预测它将来的走势就是时间 序列分析。
时间序列例1
德国业余天文学家施瓦尔发现太阳黑子的活动具有11年左右的周期
时间序列例2
上证指数
相空间重构
如果把一个时间序列看成是由一个确定性 的非线性动力系统产生的,要考虑的是以 下反问题:如何由时间序列来恢复并刻划 原动力系统?
局部多项式预测法
局部多项式预测模型为（以二次多项式为例）
x(T + 1) = h( X (T )) = c00 + c10 x(T ) + c20 x(T − τ ) + ...cm 0 x(T − (m − 1)τ ) + c11 x(T ) x(T ) + c12 x(T ) x(T − τ ) + ...c1m x(T ) x(T − (m − 1)τ ) + ⋯ + cmm x(T − (m − 1)τ ) x(T − (m − 1)τ ) + e 其中e为随机误差，cij为待定系数。
-20
来自百度文库
-30 -20
-15
-10
-5
0
5
10
15
20
Lorenz系统的吸引子(y,z相图）
50 50
45 45
40 35
40 35 30
30
25
25 20
15 20 10 15
10 5
5 0 -30
-20
-10
0
10
20
30
0 -30
-20
-10
0
10
20
30
Lorenz系统的吸引子(x,z相图）
该系统虽然有两个状态变量，但如果观测到状态变量 Xn的信息，我们可以从Xn建立原系统的模型 对状态变量Xn进行相空间重构:Zn=(Xn,Xn-1) 由Zn 可以重构原来的系统
Lorenz系统
dx dt = σ ( y − x) dy = x(r − z ) − y dt dz dt = xy − bz
选择使I(τ) 为第一个局部极小的τ为延迟时间间隔。
N
嵌入维数m的选取
主要方法 虚假邻点法 关联积分法 奇异值分解法
虚假邻点法
虚假邻点的定义
设当前维数为m, X η ( n )为X n的最近邻点， 距离为 X η ( n ) − X n , 当维数增加到m + 1时, 距离为 X η ( n ) − X n