非线性时间序列的例子
1. Logistic模型 xn+1=Rxn(1-xn)
R=1.5时,不管初始状态x0在何处,随时间的演化,系统都将单调地趋向于1/3,
R=2.9时,不管初始状态x0在何处,随时间的演化,系统都将交替地趋向于19/29,
Logistic模型R=3.3时, 不管初始状态x0在何处,随时间的演化,系统都将在0.48和0.82两个状态之间周期性地变化,
R=4时,随时间的演化,系统将出现不规则的振荡,看起来好像是随机的-- ->表明系统对初值具有非常敏感的依赖性,也说明这样的系统只能进行短期预测,要进行较长时间的预测会变得不正确.
2. 弹簧振子受迫振动
3. Lorenz系统
dx/dt=sigma (y-z)
dy/dt=x(r-z)-y
dz/dt=xy-bz
取sigma=10,r=28,b=8/3时,系统是混沌的.
4. 实际问题中的实测时间序列
股票指数时间序列
太阳黑子数时间序列
Chapter 2:单变量非线性时间序列分析
例2.1 Henon映射:
xn+1=1-1.4xn2+yn
yn+1=0.3xn
该系统实际上只与状态变量xn的前两个时刻的状态有关.
相空间重构的基本原理是F.Takens和R.Mane的延迟嵌入定量,它建立了观测信号系统时间波动和动力系统特征之间的桥梁.
基本思想:
通过观测或实验获得单变量时间序列{xn},做以下相空间重构:
xn=(xn,xn-tao,…,xn-(m-1)*tao)
从而形成m维状态空间,在重构的m维状态空间中可以建立数学模型:
xn+1=G(xn)
F.Takens和R.Mane证明了只要适当选取m和tao,原未知数学模型的混沌动力系统的几何特征与重构的m维状态空间的几何特征是等价的,它们具有相同的拓扑结构,这意味着原未知数学模型的混沌动力系统中的任何微分或拓扑不变量可以在重构的状态空间中计算,并且可以通过在重构的m维状态空间中建立数学模型对原未知数学模型的动力系统进行预测,进一步解释/分析/指导原未知数学模型的动力系统.
相空间重构
设动力系统是由非线性差分方程:
zn+1=F(zn)
表示的离散系统,或者是由微分方程
dz(t)/dt = F(z(t))
表示的连续系统,其中zn或z(t)是系统在时刻n或t的状态向量,F(.)是向量值函数.时间序列{xn}是观测到的系统某一维输出,即:
xn=h(zn)+wn,
或:
xn=x(t0+nΔt)=h[z(t0+nΔt)]+wn,
式中,h(.)是多元数量值函数；wn为在观测或者测量过程中由于技术手段不完善或者精度不够引起的测量噪声.
根据F.Takens的定理,当wn=0时,观察到的时间序列{xn}以向量:
xn=(xn,xn-tao,…,xn-(m-1)*tao)
形成m维空间,只要m>=2d+1,动力系统的几何结构可以完全打开,其中d是系统中吸引子的维数,tao是正整数,称为延迟时间间隔.条件m>=2d+1是动力系统重构的充分但不必要条件,获得动力系统重构的整数m叫做嵌入维数.状态空间中xn- ->xn+1的演化反映了未知动力系统zn- ->zn+1或z(t)- ->z(t+1)的演化,并且状态空间Rm中吸引子的几何特征与原动力系统的几何特征等价,这意味着原动力系统中任何微分或拓扑不变量可以重构的状态空间中计算.F.Takens的定理是在无噪声的情况下考虑的,wn=0,后来T.Sauer等把延迟嵌入定理推广到了具有噪声的情形.
对一组长为N的实测时间序列{xn}n=1N,由xn=(xn,xn-tao,…,xn-(m-1)*tao)可构造出m维状态向量:
xn=(xn,xn-tao,…,xn-(m-1)*tao) in Rm, n=N0,N0+1,…,N
其中N0=(m-1)*tao+1,tao是延迟时间间隔.在Rm中在L2或Linfinity范数定义xi到xj的距离.
为了能在重构的Rm空间中刻画原动力系统的性质,需正确地确定延迟时间间隔tao和嵌入维数m.
延迟时间间隔的确定
由延迟嵌入定理可知,在时间序列无限长,无噪声,无限精确的情况相,可以任意选择tao,但实测时间序列是有限长的,且一般都有噪声污染和测量误差,只能根据经验来选择tao.选择tao的基本思想是使xn与xn+tao具有某种程度的独立但又不完全相关,以便它们能在重构的相空间中作为独立的坐标处理.如果tao太大,则xn与xn+tao的值充分靠近,以至于不能区分它们,从实际观点看不能提供两个独立的坐标,导致吸引子重构非常靠近相空间中的”对角线”,重构的相空间可能总是杂乱无规则的；如果tao太大, xn与xn+tao可能会不相关,吸引子轨道会投影在两个完全不相关的方向上,不能反映相空间中轨线的真实演化规则.鉴于此,需要选择一个比较合适的延迟时间间隔tao.目前,可以使用的方法有很多,但从计算复杂性和使用的简便性等角度看,比较常用的方法主要有自相关函数法和平均互信息法.
1. 自相关函数法
基本思想是要考察观测量xn与xn+tao与平均观测量的差之间的线性相关性,即如果假设:
xn+tao – x-bar = CL(tao) (xn – x-bar)
其中:
x-bar=1/N Σn=1N xn (算术平均),
则使:
Σn=1N [xn+tao – x-bar - CL(tao) (xn – x-bar)]2,
最小的CL(tao)为:
CL(tao)=A/B;
A=1/N Σn=1N (xn+tao – x-bar)(xn – x-bar)
B=1/N Σn=1N (xn – x-bar)2;
称这样的CL(tao)为线性自相关函数.取CL(tao)第一次为零时的tao为延迟时间间隔,此时在平均意义下, xn与xn+tao是线性无关的.
自相关函数法是比较简单的寻找时间延迟tao的一种方法,但这种方法只考虑到时间序列中线性关系,至于非线性关系并不清楚,所以并不适合所有情况.特别当自相关函数变化十分缓慢时,选择会非常困难.
2. 平均互信息法
不同于自相关函数法,平均互信息法将非线性关系也考虑在内,这种方法的根据是可从事件bj在B中发生的概率中得到关于aj在集A中发生概率的信息.
I(tao)= 1/N Σn=1N P(xn, xn+tao) log2 [P(xn, xn+tao)/(p(xn)p(xn+tao))]
式中,P(xn),P(xn+tao),P(xn,xn+tao)为概率.
概率P(xn),P(xn+tao)可以通过计算时间序列的直方图获得,联合概率P(xn,xn+tao)可以通过计算时间序列的二维直方图获得,利用计算机可以很方便地计算观测时间序列的平均互信息.
文献[5]建议选择I(tao)的第一个局部最小时的tao为延迟时间间隔,因为此时产生的冗余最小,产生了最大的独立性.与自相关函数法相比,平均互信息法考虑了非线性依赖性,但仍有其局限性,如有时可能无局部最小或对某些例子特别不合适.
除此之外,选择tao的方法还有重构展开法,高阶关联法,通过分析整体和局部混沌吸引子行为获得优化延迟时间的填充因子法等多种方法,这些方法都有各自的特点,但实际应用中用得较多的还是自相关函数法和平均互信息法.
嵌入维数的确定
已从理论上证明m>=2d+1时可获得一个吸引子的嵌入,其中d是吸引子的分形维数,但这只是一个充分条件,对观测时间序列选择m没有帮助.如果仅仅是计算关联维数,已证明了对无噪声,无限长的时间序列,只要取m为大于关联维数d的最小整数即可,但对长度有限且具噪声的时间序列,m要比d大得多.如果m选得太小,则吸引子可能折叠以致在某些地方自相交,这样在相交区域的一个小领域内可能会饮食来自吸引子不同部分的点;如果m选得太大,理论上是可以的,但在实际应用中,随着m的增加会大大增加吸引子几何不变量(如关联维数,Lyapunov指数)的计算工作量,且噪声和舍入误差的影响会大大增加.
1. 试算法
通过逐步增加计算过程中的嵌入维数,观察什么时候某些几何不变量(例如,关联积分/关联维数/Lyapunov指数等)停止变化的方法.从理论上来讲,由于这些几何不变量是吸引子的几何性质,当m大于最小嵌入维数时,几何结构被完全打开,因此这些不变量与嵌入维数无关,取吸引子的几何不变量停止变化时的m为最小维数.这种方法的缺点是对数据要求较高(无噪声),计算量大且比较主观.例如,P.Grassberger通过增加嵌入维数m,计算关联积分CN(r,m,tao),取当关联积分CN(r,m,tao)不再变化时的m为嵌入维数.
2. 虚假邻点法
虚假邻点法建立在以下事实的基础上:选择太小的嵌入维数将导致那些在原相空间中离得比较远的点会在重构的相空间中靠近.其基本思想是当嵌入维数从m变化到m+1时,考察轨线xn邻点中哪些是真实的邻点,哪些是虚假的邻点,当没有虚假邻点时,可以认为几何结构被完全打开.
设xn的最近邻点为xyita(n),当嵌入维数从m增加到m+1时,它们之间的距离从d(m)变为d(m+1),若d(m+1)比d(m)大很多,可以为是由于高维吸引子中两个不相邻的点在投影到低维轨线上时变成相邻的两点造成的,因此这样的邻点是虚假的.阈值Rtol.
观测时间序列通常具有噪声且长度有限,所以仅仅用上面的标准判别虚假最近邻点会不正确,为此M.B.Kennel等提出了增加以下标准,阈值Atol.
让m从1开始增加,计算每个m时的虚假最近邻点的比例,直到虚假最近邻点的比例小于5%或虚假最近邻点不再随着m的增加而减少时,可以认为吸引子几何结构完全打开(???????是为何意???????),此时的m为嵌入维数.
从几何的观点来看,这是一种较好的方法,但在判断虚假邻点时阈值的不同选取会导致不同的结果,因此具有较大的主观性.阈值应当由时间序列的方差确定,因此,不同的时间序列有不同有阈值,这意味洋着要给出一个合适和合理的阈值是非常困难的甚至是完全不可能的.
3. 改进的虚假邻点法,文献[15]
类似虚假最近邻点法的思想,定义:
a(n,m) = d(m+1)infinity / d(m)infinity, 采用Linfinity范数
记所有a(n,m)关于n的均值值为:
E(m)=1/(N-N0+1) Σn=N0N a(n,m),
式子xn=(xn,xn-tao,…,xn-(m-1)*tao) in Rm, n=N0,N0+1,…,N 是向后重构,而文献[15]是向前重构,所以有所不同,但本质上是一样的.
E(m)只依赖于嵌入维数m和延迟时间间隔tao,为了研究嵌入维数从m变为m+1时相空间的变化情况,定义:
E1(m)=E(m+1)/E(m),
如果当m大于某个m0时,E1(m)停止变化,则m0+1就是重构相空间的最小嵌入维数.
对来自于随机系统的时间序列,原则上,随着m的增加,E1(m)将永远不会达到饱和值,但在实际计算中,当m充分大时,很难分辨E1(m)是在缓慢增加还是停止变化.事实上,由于获得的观测时间序列是有限长的,可能会出现虽然是随机时间序列,但E1(m)会在某一m处停止变化.为了解决这一问题,文献[15]定义了量E*(m):
记E2(m)= E*(m+1)/E*(m),对随机时间序列,由于将来的值与过去的值独立,对任何m,E2(m)将等于1,而对确定性系统的观测时间序列,E2(m)显示与m有关,不可能为常数.因此,文献[15]建议可以通过同时计算E1(m)与E2(m)来确定时间序列相空间重构的最小嵌入维数,同时也区分了时间序列是来自于确定性系统还是随机系统.
除上述方法外,确定嵌入维数的方法还有关联积分法,奇异值分解以及上述各种的改进方法.
2.3几何不变量的计算
分形维数的定义有很多种,从时间序列的角度看,关联维数是易于计算的一种分形维数(G-P算法以及Kolmogorov熵公式) .关联维数是系统复杂性程度的一种很好的度量.一般认为,大于关联维数的下一个整数是刻画系统所需的独立变量的个数,这为从时间序列恢复原复杂系统确定了一个框架.Lyapunov指数(轨线法)度量了复杂系统的预测性,定量地刻画了初始靠近的状态空间轨线的指数发散.
在确定性系统中,关联维数就是生成相应复杂系统所必需的独立变量的个数,规则的确定性系统有整数关联维数,而混沌系统有非整数关联维数,大于此关联维数的下一个整数就是系统的独立变量的个数.但是某些关联的随机过程中也有非整数维数[24].文献[25]证明了非整数关联维数不能成为混沌判断的充分条件,因为分形的Browian运动,虽然不是混沌的,但也有非整数关联维数.
关联维数
G-P算法:(经典方法)
关联积分:CN(r,τ,m)
如果在r的某一区间段内,有 CN(r,τ,m)∝rd,则称d是关联维数,这样定义的d就是近似刻产生时间序列的复杂系统复杂程度的某种维数.d就是双对数图ln CN(r,τ,m) – ln r图的线性部分的分辨率,由于连续采样点之间的相关性,会造成”肩峰”效应,为了消除这种效应,文献[27]对关联积分作了修改,忽略了嵌入空间中非常靠近的点对关联积分的贡献,只考虑满足|i-j|≥w的点,其中w>τ(2/N)2/m,特别取w=τ即可.
实际应用中,关联维数的计算是非常耗费时间的,为了改善前面算法的计算效率,文献[35]提出了一种修改算法.在关联积分中,较短的距离起着更有意义的作用,因此选取r的截断距离r0,这样可以把计算时间从O(N2)减到O(Nlog2N).
当关联维数的值较大时,所需的时间序列要求较长,而实际问题中,观测或实验获得的时间序列一般都比较短,为解决这一问题,文献[30]利用极大似然法估计关联维数的方法就有对数据要求相对较少的优点,这一方法的关键在于r0的选取,较小时,可计算的点减少,则估计变得不可靠；较大时,虽然可以产生较稳定的结果,但由于C(r)=(r/r0)d, (0<r≤r0)的距离分布方差使结果产生较大的偏差,文献[30]详细给出了利用χ2测试确定r0的一种方法.
Kolmogorov熵K1和Renyi熵Kq
显然有:limq->∞Kq=K1
Kq是关于q单调减少的,特别q=2时,就是关联积分,因此可由关联积分计算K2.它可以作为K1的一个下界.
Lyapunov指数
混沌系统的基本特点就是系统对初始值的极端敏感性，两个相差无几的初值所产生的轨迹，随着时间的推移按指数方式分离，lyapunov指数就是定量的描述这一现象的.
| 1、这个定义是针对信息技术中的开放性而言的.Lyapunov指数的定义是:取映射F所代表系统在各次迭代点处导数绝对值的对数平均值相空间中的每一维都有其各自的Lyapunov指数 2、Lyapunov指数是指相空间中邻近轨道发散或收敛的平均指数率. 3、Lyapunov指数是指相空间中邻近轨道发散或收敛的平均指数率.考察一个一维映射：Xn＋l二f（n）（l）当初始条件xo出现一个偏差社 4、Lyapunov指数被定义为:λ=1N∑Ni=1log|λmax(Ai)|这里Ai是函数f在点xi处的局部线性化系统的转换矩阵,即f(x)=Aixi+f′(xi)(x-xi),λmax(·)表示矩阵的最大特征值 5、对于一动力学系统x=f(x),其Lyapunov指数定义为:λi=limn→∞1n‖J(xn-1)J(xn-2).J(x0)e0‖(1) 其中J(xi)为系统在xi的Jacobian矩阵,e0为单位向量 6、对于一动力学系统x=f(x),其Lyapunov指数定义为:λi=limn→∞1n‖J(xn-1)J(xn-2).J(x0)e0‖(1) 其中J(xi)为系统在xi的Jacobian矩阵,e0为单位向量 7、在非混沌系统中,相互靠近的轨迹于要么指数地迅速收敛,要么慢于指数迅速地发散,至少在理论上,长期预测的是可能的,这种轨迹收敛或发散的比率,称为Lyapunov指数 8、(5)式中:λ是不依赖于初值x0的数,称为Lyapunov指数.对于非线性映射而言,N为映射迭代的次数,f′(xi)为映射函数在xi点的导数 9、在非混沌系统中,相互靠近的轨迹于要么指数地迅速收敛,要么慢于指数迅速地发散,至少在理论上,长期预测的是可能的,这种轨迹收敛或发散的比率,称为Lyapunov指数 10、v指数考虑一维映射xn+1=f(xn)设它的初值x0的偏差为δx0,则经过n次迭代以后可以得到δxn=|fn(x0+δx0)-fn(x0)|=dfn(x0)dx0δx0=eλnδx0由上式可以得到λ=1nlnδxnδx0=1nlndfn(x0)dxλ称为Lyapunov指数 11、v指数考虑一维映射xn+1=f(xn)(5)设它的初值x0的偏差为δx0,则经过n次迭代以后可以得到δxn=|fn(x0+δx0)-fn(x0)|=dfn(x0)dx0δx0=eλnδx0(6)由上式可以得到λ=1nlnδxnδx0=1nlndfn(x0)dx(7)λ称为Lyapunov指数 12、Lyapunov指数用于量度在相空间中初始条件不同的两条相邻轨迹随时间按指数律吸引或分离的程度,这种轨道收敛或发散的程度称为Lyapunov指数.对于系统xn+1=f(xn),Lyapunov指数可表示为λ(x0)=limn→∞nlndfn(x0)dx(2)式中f为混沌映射 13、Lyapunov指数用于量度在相空间中初始条件不同的两条相邻轨迹随时间按指数律吸引或分离的程度,这种轨道收敛或发散的程度称为Lyapunov指数.对于系统xn+1=f(xn),Lyapunov指数可表示为[4]λ(x0)=limn→∞nlndfn(x0)dx(2)式中f为混沌映射 |
如果最大Lyapunov指数是正的,意味着相邻的轨线按指数发散,即系统是混沌的.
对观测获得的混沌时间序列,最大Lyapunov指数可由A.Wolf等提出的轨线法计算[20],而Lyapunov指数谱可由J.P.Eckmann和R.Brown等提出的矩阵法计算[41,42].
Lyapunov exponentFrom Wikipedia, the free encyclopediaJump to: navigation, search In mathematics the Lyapunov exponent or Lyapunov characteristic exponent of a dynamical system is a quantity that characterizes the rate of separation of infinitesimally close trajectories. Quantitatively, two trajectories in phase space with initial separation where λ is the Lyapunov exponent. The rate of separation can be different for different orientations of initial separation vector. Thus, there is a whole spectrum of Lyapunov exponents—the number of them is equal to the number of dimensions of the phase space. It is common to just refer to the largest one, i.e. to the Maximal Lyapunov exponent (MLE), because it determines the predictability of a dynamical system. A positive MLE is usually taken as an indication that the system is chaotic.
[edit] Definition of the maximal Lyapunov exponentThe maximal Lyapunov exponent can be defined as follows: [edit] Definition of the Lyapunov spectrumFor a dynamical system with evolution equation ft in a n–dimensional phase space, the spectrum of Lyapunov exponents in general, depends on the starting point x0. The Lyapunov exponents describe the behavior of vectors in the tangent space of the phase space and are defined from the Jacobian matrix The Jt matrix describes how a small change at the point x0 propagates to the final point ft(x0). The limit defines a matrix L(x0) (the conditions for the existence of the limit are given by the Oseledec theorem). If Λi(x0) are the eigenvalues of L(x0), then the Lyapunov exponents λi are defined by The set of Lyapunov exponents will be the same for almost all starting points of an ergodic component of the dynamical system. [edit] Basic propertiesIf the system is conservative (i.e. there is no dissipation), a volume element of the phase space will stay the same along a trajectory. Thus the sum of all Lyapunov exponents must be zero. If the system is dissipative, the sum of Lyapunov exponents is negative. If the system is a flow, one exponent is always zero—the Lyapunov exponent corresponding to the eigenvalue of L with an eigenvector in the direction of the flow. [edit] Significance of the Lyapunov spectrumThe Lyapunov spectrum can be used to give an estimate of the rate of entropy production and of the fractal dimension of the considered dynamical system. In particular from the knowledge of the Lyapunov spectrum it is possible to obtain the so-called Kaplan-Yorke dimension DKY, that is defined as follows: where k is the maximum integer such that the sum of the k largest exponents is still non-negative. DKY represents an upper bound for the information dimension of the system [1]. Moreover, the sum of all the positive Lyapunov exponents gives an estimate of the Kolmogorov-Sinai entropy accordingly to Pesin's theorem [2] The inverse of the largest Lyapunov exponent is sometimes referred in literature as Lyapunov time, and defines the characteristic e-folding time. For chaotic orbits, the Lyapunov time will be finite, whereas for regular orbits it will be infinite. [edit] Numerical calculationGenerally the calculation of Lyapunov exponents, as defined above, cannot be carried out analytically, and in most cases one must resort to numerical techniques. The commonly used numerical procedures estimates the L matrix based on averaging several finite time approximations of the limit defining L. One of the most used and effective numerical technique to calculate the Lyapunov spectrum for a smooth dynamical system relies on periodic Gram-Schmidt orthonormalization of the Lyapunov vectors to avoid a misalignement of all the vectors along the direction of maximal expansion [3] [4]. For the calculation of Lyapunov exponents from limited experimental data, various methods have been proposed. [edit] Local Lyapunov exponentWhereas the (global) Lyapunov exponent gives a measure for the total predictability of a system, it is sometimes interesting to estimate the local predictability around a point x0 in phase space. This may be done through the eigenvalues of the Jacobian matrix J 0(x0). These eigenvalues are also called local Lyapunov exponents. The eigenvectors of the Jacobian matrix point in the direction of the stable and unstable manifolds. |
参考资料:
2.4 观测时间序列平稳性的检验
“弱平稳性”的概念
2.5 基于观测时间序列的系统非线性性检验
2.6 基于观测时间序列的系统确定性检验
2.7 观测时间序列噪声处理技术
Chapter 3: 基于观测时间序列的系统非线性性检验
Chapter 4: 多变量时间序列相空间重构
Chapter 5: 多变量非线性时间序列预测方法
Chapter 6: 非线性时间序列分析法在证券市场中的应用
……
没有评论:
发表评论