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

时间序列的定义按照时间的顺序把事件变化发展的过程记录下来就构成了一个时间序列。对时间序列进行观察、研究，找寻它变化发展的规律，预测它将来的走势就是时间序列分析。时间序列例1德国业余天文学家施瓦尔发现太阳黑子的活动具有11年左右的周期时间序列例2上证指数相空间重构如果把一个时间序列看成是由一个确定性的非线性动力系统产生的,要考虑的是以下反问题:如何由时间序列来恢复并刻划原动力系统?由时间序列恢复原系统最常用的方法利用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)可描述原系统。τ称为延迟时间，m称为嵌入维数。由x(n)构造X(n)称为相空间重构。相空间重构法基本思想是:系统中任一分量的演化都是由与之相互作用着的其它分量所决定的,因此这些相关分量的信息就隐含在任一分量的发展过程中。为了重构一个等价的状态空间,只需考察一个分量,并将它在某些固定的时间延迟点上的测量作为新维处理,它们确定了某个多维状态空间中的一点.重复这一过程并测量相对于不同时间的各延迟量,就可以产生出许多这样的点,它可以将原系统的许多性质保存下来,即用系统的一个观察量可以重构出原动力系统模型,可以初步确定原系统的真实信息。相空间重构例Henon映射nnnnnxyyxx3.04.11121该系统虽然有两个状态变量，但如果观测到状态变量Xn的信息，我们可以从Xn建立原系统的模型对状态变量Xn进行相空间重构:Zn=(Xn,Xn-1)由Zn可以重构原来的系统Lorenz系统bzxydtdzyzrxdtdyxydtdx)()(91.37,68.13,34.1538,2810000zyxbr初值，取Lorenz系统的吸引子(x-y-z）-40-20020400204060-20-1001020Lorenz系统的吸引子(x,y相图）-20-15-10-505101520-30-20-100102030Lorenz系统的吸引子(y,z相图）-30-20-10010203005101520253035404550-30-20-10010203005101520253035404550Lorenz系统的吸引子(x,z相图）-20-15-10-50510152005101520253035404550如果只观测到变量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))……-20-1001020-20-1001020-20-1001020-40-20020400204060-20-1001020重构后的相图(x-y-z)原始系统相图(x-y-z)-20-15-10-505101520-20-15-10-505101520-20-15-10-505101520-30-20-100102030重构后的相图(x-y)原系统的相图(x-y)-20-15-10-505101520-20-15-10-505101520重构后的相图(y-z)-30-20-10010203005101520253035404550原系统的相图(y-z)-20-15-10-505101520-20-15-10-505101520重构后的相图(x-z)-20-15-10-50510152005101520253035404550原系统的相图(x-z)如何确定延迟时间和嵌入维数？延迟时间间隔τ的选取主要方法线性自相关函数法平均互信息法线性自相关函数法NnnNnnNnnnxNxxxNxxxxNC11211,)(1))((1)(其中定义自相关函数为选择使得自相关函数C(τ)第一次为零时的τ的值为延迟时间平均互信息法为概率，其中定义平均互信息为对于时间序列),()(,)()(),(log),()(,21nnnnnnnnnNnnxxPxPxPxPxxPxxPIx选择使I(τ)为第一个局部极小的τ为延迟时间间隔。嵌入维数m的选取主要方法虚假邻点法关联积分法奇异值分解法虚假邻点法虚假邻点的定义的虚假邻点。为大很多，认为比如果距离为时当维数增加到距离为的最近邻点，为设当前维数为nnnnnnnnnnnnXXXXXXXXmXXXXm)(m)(1m)(1m)(m)()(,,1,,上面的距离差由于和时间序列数据的大小有关，不太容易确定虚假邻点。实际采用相对度量法之间。和在值的虚假邻点。，其中阈为称若5010/)()(m)(m)(1m)(RXXRXXXXXXnnnnnnnn2)(1/211m)(NnnnnxxNXX标准还需补充以下虚假邻点对于实际的时间序列，虚假邻点法确定嵌入维数对实测时间序列,m从2开始,取R=30,计算虚假最近邻点的比例,然后增加m,直到虚假最近邻点的比例小于5%或虚假最近邻点不再随着m的增加而减少时,可以认为此时的m为合适的嵌入维数。非线性时间序列预测基本思想设时间序列来自确定性系统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)的预测值。基本方法局域预测法局部平均预测法局部线性预测法局部多项式预测法全域预测法神经网络小波网络遗传算法局部平均预测法设时刻T的状态向量为X(T)=(x(T),x(T-τ),…x(T-(m-1)τ))找X(T)的最近邻点X(T1),…X(TK),以X(T1+1),…X(TK+1)的平均值作为X(T+1)的预测值)1(1)1(1^KkkTXKTX这是向量表达式，取第一个分量得)1(1)1(1^KkkTxKTx局部线性预测法局部线性预测模型为为待定系数。为随机误差，其中i210^c))1((...)()())(()1(eemTxcTxcTxccTXgTxm仍设X(T1),…X(TK)为X(T)的K个邻近点，确定Ci的方法：最小二乘法即求Ci使得最小21k21k^)1())(()1()1(KkkKkkTxTXgTxTxyAAACmTxmTxmTxTxTxTxTxTxTxAcccCTxTxTxyTTKKKTmTK12122111021)())1(())1(())1(()()(1)()(1)()(1],,,[)]1(,),1(),1([则记局部多项式预测法为待定系数。为随机误差，其中ij112110201000^c))1(())1(())1(()(...)()()()())1((...)()())(()1(eemTxmTxcmTxTxcTxTxcTxTxcmTxcTxcTxccTXhTxmmmm局部多项式预测模型为（以二次多项式为例）仍以最小二乘法确定系数yAAACmTxmTxmTxTxTxTxmTxmTxmTxTxTxTxTxTxTxAcccCTxTxTxyTTKKKKKTmmTK1)11(1))1(())1(())1(()1()1()1())1(())1(())1(()()(1)()(1)()(11],,,[1)]1(,),1(),1([12222122221212211100021则记预测效果评价为了检验预测的精确性，可以比较预测值与实际观测值之间的差。一次预测可能较好或较差，偶然性较大。为了克服这种偶然性，可以取多个点的预测误差的平均。piiTiTpTTTpTTTyxRMSEyyxxx122121)(p1,,,,y,,,定义均方根误差为的预测值为设如果RMSE比较大，则说明预测效果不好。但是RMSE和观测序列的数值大小有关，为克服这一问题，我们定义正规化均方根误差NRMSEpiTpiTxpxxxRMSENRMSE1i21i21,)(p1,/其中，若NRMSE接近于1，则预测效果不好若NRMSE接近于0，则预测效果较好预测效果的另一个评价标准是相关系数piiTpiiTpipiiTpiiTiTiTyyxxyxyxr1212111)()(相关系数若r接近于1，则预测效果较好上证指数预测数据文件000001.day1990.12.19-2008.06.19共4292个记录每一条记录的长度为40字节:1-4字节为日期5-8字节=开盘指数*10009-12字节=最高指数*100013-16字节=最低指数*100017-20字节=收盘指数*100021-24字节=成交金额(元)/100025-28字节=成交量(手)其余12字节未使用读取数据matlab文件为readdata.m上证指数预测文件为shangzhen1.m相关文件为readdata.mjuli.mdataconstruct1.mreconstruct1.m多变量时间序列以两个变量为例说明多变量情形设给定时间序列x(n),y(n)x(n)的嵌入维数为m1,延迟时间为τ1y(n)的嵌入维数为m2,延迟时间为τ2多变量时间序列的相空间重构X(n)=(x(n),x(n-τ1),…x(n-(m1-1)τ1,,y(n),y(n-τ2),…y(n-(m2-1)τ2)重构后时间序列的维数为m1+m2多变量时间序列预测设时刻T的状态向量为X(T)=(x(T),x(T-τ1),…x(T-(m1-1)τ1,,y(T),y(T-τ2),…y(T-(m2-1)τ2)预测x(T+1)的值以局部多项式（二次）预测为例预测模型为222212221122021021111221111101100)2)1(()2()()()2)1(()()1)1(()1()()()1)1(()())(()1(222111mTycTyTycTycmTycTycmTxcTxTxcTxcmTxcTxccTXgTxmmmmmm设X(T)的K个最近邻点为X(T1),…X(TK)如果系统是确定的，则当X(T)靠近X(Ti)时，X(T+1)应靠近X(Ti+1)以最小二乘估计参数2110022,,mmccc最小使得即求系数21i21100))(()1(,,22KiimmTXgTxcccYAAACmTyTxmTyTxmTyTxATxTxTxYcccCTTKKTKTmm122222222112121100)(,))1(()(1))1(()(1))1(()(1)]1(),1(),1([],,,[22则若记用成交量，收盘指数预测上证指数的文件为shanzhen2.m相关文件为readdata.mjuli.mdataconstruct.mreconstruct.m