动态时间弯曲算法在K线相似度计算中的应用

1、动态时间弯曲算法在K线相似度计算中的应用序言在证券交易数据中，股票K线图无疑是一种非常重要的数据它反映了股票在过去历史中基于开盘价与收盘价的交易价格的变动古有言，以史为鉴，历史往往存在相似性，对一支股过去历史波动的研究，往往可以对其自身，以及其他股票的未来价格变动作出一些合理预测而股票价格究其根本是一种时间序列对于时间序列，是一种以时间为轴，在一些特别规定的时间点上通过采样得到的一系列按照时间顺序排列的，从被观测对象获取到的观测值通过对时间序列的研究，找到两条时间序列相似程度的的度量方法就被称为时间序列相似性度量，这是时间序列聚类分析中一个不可缺少的步骤，同时也是分类、聚类、规律发现、模式识别

2、等工作的子进程对于研究股票k线图相似性，是为了对未来进行合理预测，因此度量方法应该考虑其性能对于后期时间序列数据挖掘的效果的的直接影响程度时间序列的相似程度是由度量距离的大小所决定的而相似性度量方式的特性又决定了 相似性度量的效果在时间序列相似性度量中，我们最常用的方法就是动态时间弯曲（Dynamic Time Warping，DTW）这是由Berndt于1994年提出将其应用在时间序列数据挖掘领域中，以此来发现时间序列中的模式而这刚好适用于股票k线图的相似程度的研究这是由于动态时间弯曲不仅可以消除欧式距离“点对点”的匹配缺陷，通过弯曲时间来达到时间序列数据点“一对多”的匹配，从而实现不等长时

3、间序列的度量，还对时间序列的偏移，振幅变化等情况具有较强的鲁棒性（鲁棒是Robust的音译，也就是健壮和强壮的意思鲁棒性指的是遭遇外来干扰，性质保持不变的能力）这对于不同股票的k线图在不同的时间跨度而可能形成相同价格形态有着重要意义一、 DTW算法原理动态时间弯曲是一种在语音识别领域得到首次应用的，准确性高并且鲁棒性强的时间序列相似性度量方法它区别于传统的欧几里的距离，其不同在于动态时间弯曲可以通过弯曲时间序列的时间区域从而对时间序列的数据点进行匹配，这样我们不单单能够得到更好的形态度量的效果，更重要的是我们能够度量两条不等长的时间序列对于股票K线图的相似性，我们寻求的是价格形态的相似性例如威

4、廉·欧奈尔提到的一种最普遍的价格形态“带柄茶杯形态”,当我们找到与此形态相似的股票时就要做出准备，这可能是一支带动市场发展的“超级牛股”一如当年的微软与苹果要想实现股票K线图的这种相似度匹配，依靠欧式距离在度量中讲时间序列进行“一对一“的数据匹配是不够的，尽管它具有高效性，但是它并未能准确的使波峰、波谷匹配起来而动态时间弯曲则能够通过弯曲时间轴来实现“一对多”的数据匹配通过这样，动态时间弯曲就能成功将两条不同的股票K线图的波峰和波谷匹配起来，从而有助于我们度量价格形态的相似程度，体现了动态时间弯曲在形态度量上的优势1.1动态时间弯曲距离【1】在介绍动态时间弯曲算法之前，先简单的介绍一

5、下动态时间弯曲距离的定义定义1 给定两条时间序列x=x1,x2,xn和y=y1,y2,yn,计算它们之间的累积距离D(i,j)=d(xi.yj)+minD(i1,j)Di,j1D(i1,j1)其中 d(xi,yj)= | xi-y j| (1)为点xi到yj之间的距离，其中i=（1,2,n）,j=(1,2,m),当=2时为欧式距离得到的累积最小距离就是动态弯曲距离，我们记为Dwarp（）在这里我们需要特别注意一点，动态弯曲距离是不符合三角不等式的命题1 Dwarp()不满足三角不等式证明我们可以通过一个反例来证明这个论题，设x=0 y=1,2 和z=1,2,2， 那我们有： Dwarp（x，z

6、）=5> Dwarp（x，y）+Dwarp（y，z）=3+0=3如此命题得证1.2动态时间弯曲距离的计算 计算它的最终累积距离其实可以认为是在距离矩阵D中寻找一条最优的弯曲路径P，从而使得累积距离达到最小，其中距离矩阵可以表示为以任意两点之间的距离来确立n×m的距离矩阵D D= d（x1,yn)d(xn,yn)d(x1,y1)d(xn,y1) 通过寻找弯曲路径Pbest=p1,p2,pK (max(n,m) Kn+m+1)来使得S和Q的累计距离的值达到最小其中pk表示的是弯曲路径元素在距离矩阵中的位置，即pk=（i，j）k表示si与hj之间的匹配关系则由此可以知道d（pk）=d

7、（i，j）k通过观察距离矩阵，我们可以看出一般存在着多条弯曲路径，有效的弯曲路径P必须符合三个要求：（1） 边界性：p1=（1,1），pK=（n，m）即路线必须从距离矩阵的第一行第一列出发到达矩阵的第n行第n列（2） 单调性：给定pk=（i，j）和pk+1=（x，y），则xi，yj（3） 连续性: 给定pk=（i，j）和pk+1=（x，y），xi+1，yj+1 单调性和连续性的存在保证了弯曲路径中某个点的下一个点在当前点的上方、右上方、右方1.3范例给定两条时间序列x=2,5,5,4,6,7，y=2,4,6,5,7,7,如下图，计算它们的动态时间弯曲距离：我们令d(xi,yj)=|x-y|，由

8、此来构造一个距离矩阵，在这里为了方便起见把它做成了一个网格形态，每一格中的数字代表距离矩阵相应位置的d(xi,yj)的值 由此我们便可以根据动态规划（AP）找到一条最优路径从而得到Dwarp（）在股票的相似性522310523310300112511201201023033245二、基于动态时间弯曲的股票k线图相似性计算在股票K线图的相似性计算中，我们通常是将股票进行两两比较在这里我们将选取其中一支具有特别价格形态的股票作为参考股票，另一支与它进行比较的股票就叫测量股票2.1将股票数据进行向量化和标准化【6】为了便于后续计算，我们首先要将采集到股票数据进行预先处理我们将参考股票具有研究意义的3

9、0天的价格变动以向量方式记录，记为x=（x1,x2,xn）同样的我们将测量股票要测量的30天价格变动也用向量方式记录记为天价格变动也用向量方式记录，记为y=（y1,y2,yn）其中n为选取的K线图某时刻的时序标号，n=1为起点时刻，在这里我们选取的是30天，所以n=30为终点时刻xm为第m时刻的K线图的收盘价格由于股票价格差别巨大，统一的阙值不好判断，所以要进行规范处理，把所有股票的价格变动都转化为01之间，具体方法如下：Zi=ximin(x)maxxmin(x)2.2通过构建股票价格的距离矩阵寻找最优路径采用动态时间弯曲算法我们首先要根据参考股票和测试股票的收盘价格来构建一个距离矩阵Dn&#

10、215;m我们可以用如下办法来快速构建距离矩阵从而找到最优路径若把测试模板的各个时刻n=1N在一个二维直角坐标系中的横轴上标出，把参考股票的各时刻m=1M在纵轴上标出，通过这些表示时刻的整数坐标画出一些纵横线即可形成一个网络，网络中的每一个交叉点（n，m）表示测试股票与参考股票中某一时刻的交汇点用动态规划（DP）算法来计算并寻找一条通过此网络中若干格点的路径，路径通过的格点即为测试和参考股票中进行计算的时刻具体如下图：d(x30,y1)d(x30,y30)d(x3,y1)d(x2,y1)d（x1,y1）d（x1,y2）d(x1.y3)d(x1,y30)路径不是随意选择的，首先任何一种股票的价格

11、涨幅都有可能变化，但是其各时刻的先后次序不可能改变，因此所选的路径必定是从左下角出发，在右上角结束为了描述这条路径，假设路径通过的所有格点依次为（n1 ，m1 ），（ni ，mj ），（nN ，mM ），其中（n1，m1 ）=（1，1），（nN ，mM ）=（N，M）为了使路径不至于过倾斜，可以约束斜率在0.52的范围内，如果路径已经通过了格点（n ，m ），那么下一个通过的格点（n ，m ）只可能是下列三种情况之一：（n ，m ）=（n +1，m）（n ，m ）=（n +1，m +1）D(n-1,m),

12、D(n-1,m-1),D(n,m-1)（n ，m ）=（n ，m+1 ）用r表示上述三个约束条件求最佳路径的问题可以归结为满足约束条件r时，求最佳路径，使得沿路径的积累距离达到最小值，即：Dwarp（）2.3搜索该路径的方法搜索从（N, M）点出发，可以展开若干条满足r的路径，假设可计算每条路径达到（n, m）点时的总的积累距离，具有最小累积距离者即为最佳路径易于证明，限定范围的任一格点（n, m）只可能有一条搜索路径通过对于（n, m），其可达到该格点的前一个格点只可能是（n-1, m）、（n-1, m -1）和（n, m-1），那么（n, m）一定选择这3个距离之路径延伸而通过（n, m）

13、，这时此路径的积累距离为：D(n, m)=dxn,ym+minD(n1,m) D(n1,m1)D(n,m1)这样可以从（n ，m ）=（30，30）出发搜索（n ，m ），对每一个（n ，m ）都存储相应的距离，这个距离是当前格点的匹配距离与前一个累计距离最小的格点（按照设定的斜率在三个格点中进行比较）搜索到（1 ，1）时，只保留一条最佳路径如果有必要的话，通过逐点向前寻找就可以求得整条路径DTW算法可以直接按上面描述来实现，即分配两个N×M的矩阵，分别为积累距离矩阵D和时刻匹配距离矩阵d，其中时刻匹配距离矩阵d（i，j）的值为测试股票的第i时刻与参考股票的第j时刻间的距离D（N，M

14、）即为最佳匹配路径所对应的匹配距离24范例例如我们选取两支股票，一支是中成股份（000151）作为参考股份x，另一支是翼凯股份（002691）作为测试股份y，为了计算方便我们只取五天的数据第一步：提取数据X=（10.53，11.58，12.74，13.77，13.19）Y=（27.99，29.4，32.3，32.8，32.4）第二步：标准化经过Zi=ximin(x)maxxmin(x)后，X=( 0，0.324，0.682，1，0.820)Y=( 0，0.293，0.896，1，0.917)第三步：构建距离矩阵并寻找最优路径在这里我们让d(xi,yj)=|x-y|0.8200.5270.076

15、0.1800.09710.7070.10400.0830.6820.3890.2140.3180.2350.3240.0310.5720.6760.59300.2930.89610.917最终我们得到一个最优累积距离Dwarp（）=0.046三、 DTW在K线图相似度计算机的实际应用【2】在中国A股总共有3470支甚至还有中小企业板，乃至创业板的股票，再涉及到各自的历史数据，这是一个非常庞大的数据库若要在这个数据可当中找寻于某一支股票相似的其他股票，将是一个计算量非常大的工作显然不可能人工手算，而是要借助计算机的索引技术但是如此庞大的数据对索引技术的准确性和高效性有着非常大的挑战 3.1索引技

16、术索引技术是一种快速的文件访问技术，举个例子，好比我们查字典，会根据偏旁部首在索引表中找到这个字在第几页一样，索引技术就是根据事先记录好的能够代表属性的取值将它与物理地址关联，进而在数据库中搜索之前我们着重提到过一个命题Dwarp（）不符合三角不等式这个事实对于我们可以用于索引的方法有着重要的意义：任何隐含或明确地符合三角不等式的索引技术都不能避免产生错误的排除 这是一个非常严格的要求：所有的空间访问方法，以及所有使用距离、度量、有利位置树的方法都不能避免错误的排除 唯一保证不会被错误排除的方法是顺序扫描，对于像股票这种大量长序列集合来说这将是计算量巨大的为了解决这个问题，人们提出了一种方法，

17、用来避免一小部分的错误排除的同时能够实现的索引的加速 也就是说，人们的目标是提供快速的索引技术，同时尽可能地避免错误的排除 在这里人们引用了两种方法3.2基于FastMap的访问方法我引用的第一种技术是基于一种名为FastMap的方法【4】它的工作原理如下： 给定N个对象和一个距离函数，它将对象映射到k维空间中的N个点，以便保持原始距离，参数k可以由用户给出，或者人们可以在应用中调整以获得更好的系统性能关键在于是假设对象确实能够表示成n维空间中的点，并且能够仅使用距离函数将这些点投影到相互正交k维空间中（kn），从而实现降维的目的如一个三维向量（x,y,z）可以经过函数映射到二维空间（x，y）

18、换用我们如今正在研究的股票来讲，就是我们要将我们包含30天收盘价格的向量，通过一个距离函数映射到一个二维空间从而减少数据的处理在将对象映射到k个点之后，人们可以使用任何空间访问方法来组织它们并搜索范围查询 FastMap在对象的数目N（即序列）上是线性的而且，将查询序列映射到第k个点需要O（k）时间，也就是说，相对于数据库大小N，时间是恒定的与其他方法一样（参见命题1），如果不遵守三角不等式，FastMap可能会引入错误的排除我们观察到，如果我们使用原始距离的平方根即欧几里得距离（很明显它符合三角不等式），我们可以避免更多的错误排除因此，我们在选择一个距离函数时候使用的就是欧几里得距离算法图来

19、自：【2】算法1描述了如何使用FastMap处理范围查询如果将FastMap应用于平方根距离，则搜索范围也应该平方根注意F（s）表示序列的坐标在过滤步骤中，根据欧几里得距离而不是时间弯曲距离来比较两个序列，这是由于欧几里得计算的简单性和快捷性，这相当于一步粗筛选在这一步中不相关的序列被排除在外一些不符合条件的序列可能包含在内，但在后处理步骤中将删除这些序列由于两个原因，算法1比单纯的DTW方法更快首先，它扫描较少的序列其次，过滤步骤也更快，因为k比序列长度小得多（通常是一些固定常数，比如6）过滤可能会删除一些合格的序列，导致错误排除，因为我们不能保证在k -d空间中的欧几里得距离越低越好即使我

20、们使用时间弯曲距离的平方根，情况也是如此，但在实践中错误排除的可能性非常低，我们将在后面看到3.3下边界技术对于两个给定的序列x= <x1, xm>和y= <y1, yn>，让max（x）和max（y）表示x和y中的最大值 min（x）和min（y）的定义相似，但是最小值一对<max（x）,min（x）>定义了序列x可以控制的范围不失一般性，我们假设max（x）max（y）所提出的方法受以下观察的启发：观察1 max（x）-max（y）Dwarp（x; y）检查其有效性相当简单由于max（x）应该至少匹配y的一个元素，比如yi，并且我们假设max（x）max

21、（y），max（x）-max（y） max（x）-yi Dwarp（x; y）：这个观察的结果是最大值的绝对差值可以作为一个距离来限制时间弯曲距离虽然这是真的，但是，它可能不是很有用，因为它可能低于时间弯曲距离太多因此，我们的目标是找到一个更紧密的下界我们考虑两个序列的范围的可能排列进行比较观察2给定两个范围，Rx = <max（x）, min（x）>和Ry = <max(y),min(y)>（范围的可能排列（见下图）图来自：【2】我们有三种可能的排列1.overlap :Rx和Ry重叠 （min（x）max（y）; min（x）min（y）2.encloses :Rx

22、包含Ry (min（x）<min（y）3.disjoint :Rx和Ry是不相交的 (min（x）> max（y）图来自：【2】值得注意的是，当通过扫描一次将序列输入到数据库中时，可以计算序列的最小值和最大值，并且可以将序列存储以供将来使用而且，通过简单的比较，可以在恒定的时间内确定两个序列的范围的排列最后，Dlb的定义只需要对每个序列进行一次扫描，因此我们可以计算序列长度中线性时间内两个序列之间的Dlb距离这会令速效率有很大的改善，除非Dlb太低估Dwarp定理1对于任何两个序列x = <x1,., xm>和y= <y1,., yn>，Dlb（x， y）D

23、warp（x， y）证明：见附录A.作为定理1的直接结果，我们得到以下推论推论1（没有误诊率）对于任何两个序列x = <x1,., xm>和y= <y1,., yn>，如果Dwarp（x; y），则Dlb（x; y）实际距离D2与下边界距离dlb是一种限制条件，可以保证范围查询和最近相邻查询不会被错误排除8在过滤步骤中，快速滤除不相关序列，因为可以快速计算Dlb距离（维度k上的线性时间，通常为k 10）一些不合格的序列可能包含在这个步骤的结果中，因为Dlb小于我们的动态时间弯曲距离Dwarp但是，这些不合格序列在后处理步骤中被删除具体算法如下：算法2图来自：【2】3.3

24、结合两种技术【5】假设有一个股票数据库包含许多任意长度的时间序列并且我们想要找出与某个查询序列相似的所有序列，即时间弯曲距离单位内的所有序列 这种类型的查询是有效范围查询处理这种查询的直接方式是扫描所有序列并为每个扫描序列计算Dwarp（），以选择符合条件的那些序列虽然非常简单，但速度可能很慢，因为它读取数据库中的每一个序列（因此速度很慢），并且它计算每个股票数据库序列的时间弯曲距离这个问题的独特之处在于，不仅I / O成本（第一种情况）很高，而且还有计算成本（第二种情况） 因此，我们所需要的技术应该解决这两个问题为了解决这些问题，人们采用了之前提到的两种技术1.使用FastMap构建索引结构

25、以加快查询处理速度 这种技术可能会导致一些错误的排除 2.使用低于时间弯曲距离的下边界距离函数 这种方法保证不会有错误的排除 3.使用这两种技术的组合 由于这两种技术彼此独立，因此可以采用流水线方式进行组合仔细考虑前面两部分提出的两种技术可以带来更高效的算法 在算法1中，我们仅使用时间弯曲距离筛选过滤的序列 但是，我们可以在计算时间弯曲距离之前使用下边界距离 如果序列确实是合格的序列，这就变成了额外的成本，但是如果不合格就可以节省大量的计算成本 这样我们就有了一个灵活的多阶段查询处理系统，如图3所示，其中FastMap和Dlb分别作为主要和次要过滤器图来自：【2】它由三个阶段组成，它们以流水线

26、方式连接第一阶段仅使用FastMap索引排除不相关的序列 在此阶段进行过滤可降低I / O成本和CPU成本 那些通过第一过滤阶段的序列在下一阶段与Dlb（）的查询序列进行比较 最后，后处理阶段只选择那些真正符合条件的序列 四、计算机程序1. 选取一支股票，例如新钢股份作为参考模板，将其某三十天的收盘价格以向量方式输入2. 然后另选取三十支股票作为股票池S，股票池中的每一支股票Si作为测试模板与参考模板进行配对，以期待获得两条如下图般相似的股票 3. 将所有股票数据进行单位化Zi=ximin（x）maxxmin(x)例如4用D2（欧几里得距离）【3】进行第一次筛选D2=2x1y12+(xnyn)

27、25用Dlb进行第二次筛选6输出R中所有的序列si的Dwarp（X，si）值具体程序见附录B5实验我们仍然选取以中成股份（000151）为参考股票，以翼凯股份（002691）作为测试股份第一步：读取数据X=（7.71，7.73，7.76，7.81，7.85，7.61，7.36，7.43，7.48，7.59，7.46，7.56，7.62，7.66，7.69，7.7，7.73，7.66，7.71，7.77，7.67，7.79，7.91，8.7，9.57，10.53，11.58，12.74，13.1，12.45）Y=（23.23,23.26,23.26,23.29,23.38,23.56,23.22

28、,23.12,23.05,23.06,23,22.95,22.96,22.93,24.12,24.45,23.3,23.27,23.2,23.28,23.14,22.87,23.89,25.4,26.46,27.99,29.4,32.3,32.8,32.4）第二步：标准化X=（0.061 ，0.064 ，0.070 ，0.078 ，0.085 ，0.044 ，0.000 ，0.012 ，0.021 ，0.040 0.017 ，0.035， 0.045 ，0.052 ，0.057 ，0.059 ，0.064 ，0.052 ，0.061 ，0.071 ，0.054 ，0.075 ，0.096 ，0.233 ，0.385 ，0.552 ，0.735 ，0.937 ，1.000 ，0.887）Y=（0.036 ，0.039 ，0.039 ，0.042 ，0.051 ，0.069 0.035 ，0.025 ，0.018 ，0.019 ，0.013 ，0.008 ，0.009 ，0.006 ，0.126 ，0.159 ，0.043 ，0.040 ，0.