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

动态时间弯曲算法在K线相似度计算中的应用序言在证券交易数据中，股票K线图无疑是一种非常重要的数据．它反映了股票在过去历史中基于开盘价与收盘价的交易价格的变动．古有言，以史为鉴，历史往往存在相似性，对一支股过去历史波动的研究，往往可以对其自身，以及其他股票的未来价格变动作出一些合理预测．而股票价格究其根本是一种时间序列．对于时间序列，是一种以时间为轴，在一些特别规定的时间点上通过采样得到的一系列按照时间顺序排列的，从被观测对象获取到的观测值．通过对时间序列的研究，找到两条时间序列相似程度的的度量方法就被称为时间序列相似性度量，这是时间序列聚类分析中一个不可缺少的步骤，同时也是分类、聚类、规律发现、模式识别等工作的子进程．对于研究股票k线图相似性，是为了对未来进行合理预测，因此度量方法应该考虑其性能对于后期时间序列数据挖掘的效果的的直接影响程度．时间序列的相似程度是由度量距离的大小所决定的．而相似性度量方式的特性又决定了相似性度量的效果．在时间序列相似性度量中，我们最常用的方法就是动态时间弯曲(DynamicTimeWarping,DTW)•这是由Berndt于1994年提出将其应用在时间序列数据挖掘领域中，以此来发现时间序列中的模式.而这刚好适用于股票k线图的相似程度的研究.这是由于动态时间弯曲不仅可以消除欧式距离“点对点”的匹配缺陷，通过弯曲时间来达到时间序列数据点“一对多”的匹配，从而实现不等长时间序列的度量，还对时间序列的偏移，振幅变化等情况具有较强的鲁棒性(鲁棒是Robust的音译，也就是健壮和强壮的意思.鲁棒性指的是遭遇外来干扰，性质保持不变的能力•)•这对于不同股票的k线图在不同的时间跨度而可能形成相同价格形态有着重要意义.一、DTW算法原理动态时间弯曲是一种在语音识别领域得到首次应用的，准确性高并且鲁棒性强的时间序列相似性度量方法.它区别于传统的欧几里的距离，其不同在于动态时间弯曲可以通过弯曲时间序列的时间区域从而对时间序列的数据点进行匹配，这样我们不单单能够得到更好的形态度量的效果，更重要的是我们能够度量两条不等长的时间序列.对于股票K线图的相似性，我们寻求的是价格形态的相似性.例如威廉・欧奈尔提到的一种最普遍的价格形态“带柄茶杯形态”,当我们找到与此形态相似的股票时就要做出准备，这可能是一支带动市场发展的“超级牛股”一如当年的微软与苹果.要想实现股票K线图的这种相似度匹配，依靠欧式距离在度量中讲时间序列进行“一对一“的数据匹配是不够的，尽管它具有高效性，但是它并未能准确的使波峰、波谷匹配起来•而动态时间弯曲则能够通过弯曲时间轴来实现“一对多”的数据匹配•通过这样，动态时间弯曲就能成功将两条不同的股票 K线图的波峰和波谷匹配起来，从而有助于我们度量价格形态的相似程度，体现了动态时间弯曲在形态度量上的优势.1.1动态时间弯曲距离【1】

在介绍动态时间弯曲算法之前，先简单的介绍一下动态时间弯曲距离的定义.定义1给定两条时间序列x={x，x,…，x}和y={y，y,…，y}，计算它们之间的累TOC\o"1-5"\h\z1 2 n "1"2 n积距离D(i-1,j)D(i,j)=d(x冷j)+min{D(i,j-1)1' D(i-1,j- 1)其中d(x.,y.)=IIx.-y.11 (1)ij i j CD为点x.到yj之间的距离，其中i=(12…,n),j=(1,2,…,m),当=2时为欧式距离.得到的累积最小距离就是动态弯曲距离，我们记为 D()在这里我们需要特别注意一点，动态弯曲距离是不符合三角不等式的命题1D()不满足三角不等式warp证明我们可以通过一个反例来证明这个论题，设x?=?0?y?=?1，2? 和？=?1,2,2?z,那我们有：Dwarp(x?,z?)=5>D(x?,y?)+D(y?,z?)warp warp=3+0=3如此命题得证.1.2动态时间弯曲距离的计算计算它的最终累积距离其实可以认为是在距离矩阵 D中寻找一条最优的弯曲路径P，从而使得累积距离达到最小，其中距离矩阵可以表示为以任意两点之间的距离来确立 的距离矩阵D??(???????)???(??????)D=(•• \ 八• ••\ 75D=(???)??(????,??)???(??????),通过寻找弯曲路径P={p1,p2,…pK}(max(n,m)<K<n+m+1)来使得S和Q的累best 丄丄 丄 计距离的值达到最小■其中p表示的是弯曲路径元素在距离矩阵中的位置，k即p=(i,j)表示s.与h之间的匹配关系.则由此可以知道d(p)=d(i,j)k k i ] k•通过观察距离矩阵，我们可以看出一般存在着多条弯曲路径，有效的弯曲路k径P必须符合三个要求：(1)边界性：p=(1,1),p=(n,m)■即路线必须从距离矩阵的第一行第1 K—列出发到达矩阵的第n行第n列(2)单调性：给定pk=(i,j)和pk+1=(x,y)，则x>i,y>j.(3)连续性：给定pk=(i,j)和pk1=(x,y),x<i+1,y<j+1■单调

性和连续性的存在保证了弯曲路径中某个点的下一个点在当前点的上方、右上方、右方.1.3范例给定两条时间序列x={2,5,5,4,6,7} ,y={2,4,6,5,7,7},如下图，计算它们的动态时间弯曲距离:77我们令d（x.訂.）=lx-yl，由此来构造一个距离矩阵，在这里为了方便起见把它）的值.在股票）的值.在股票做成了一个网格形态，每一格中的数字代表距离矩阵相应位置的 d（xi,y由此我们便可以根据动态规划（AP）找到一条最优路径从而得到DOwarp的相似性52231052331丿03001125112/012/0*—023033245：、基于动态时间弯曲的股票k线图相似性计算在股票K线图的相似性计算中，我们通常是将股票进行两两比较.在这里我们将选取其中一支具有特别价格形态的股票作为参考股票，另一支与它进行比较的股票就叫测量股票.2.1将股票数据进行向量化和标准化【6】为了便于后续计算，我们首先要将采集到股票数据进行预先处理我们将参考股票具有研究意义的30天的价格变动以向量方式记录，记为x?（Xi，X2,…,J）同样的我们将测量股票要测量的30天价格变动也用向量方式记录记为天价格变动也用向量方式记录，记为y?=（y^y2,…y）其中n为选取的K线图某时刻的时序标号:n=1为起点时刻，在这里我们选取的是30天，所以n=30为终点时刻.x为第m时刻的K线图的收盘价格.m由于股票价格差别巨大，统一的阙值不好判断，所以要进行规范处理，把所有股票的价格变动都转化为0〜1之间，具体方法如下：xi-mm（x?）Z= imax（x?）-min?（x?）2.2通过构建股票价格的距离矩阵寻找最优路径采用动态时间弯曲算法我们首先要根据参考股票和测试股票的收盘价格来构

建一个距离矩阵D 我们可以用如下办法来快速构建距离矩阵从而找到最优路nxm径.若把测试模板的各个时刻n=1~N在一个二维直角坐标系中的横轴上标出，把参考股票的各时刻m=1~M在纵轴上标出，通过这些表示时刻的整数坐标画出一些纵横线即可形成一个网络，网络中的每一个交叉点（n,m）表示测试股票与参考股票中某一时刻的交汇点.用动态规划（DP）算法来计算并寻找一条通过此网络中若干格点的路径，路径通过的格点即为测试和参考股票中进行计算的时刻.具体如下图：d(x30』1)d（x30』3。）・・・・・・d(x3,y1)叽』1）d('y1)d(\y2)d(x1-y3)-'・…(1(x1，y3。)路径不是随意选择的，首先任何一种股票的价格涨幅都有可能变化，但是其各时刻的先后次序不可能改变，因此所选的路径必定是从左下角出发，在右上角结束为了描述这条路径，假设路径通过的所有格点依次为（「，m】）,……,（nimj）， ，（nN，mM），其中（n]，m]）=（1,1）,（nN，叫）=（N,M）■为了使路径不至于过倾斜，可以约束斜率在0.5~2的范围内，如果路径已经通过了格点（n,m），那么下一个通过的格点（n,m）只可能是下列三种情况之

@（n,m）=（n+1,m+1）{D（n-l，m），D（n-l，m-l），D（n，m-1）}3（n,m）=（n,m+1）用r表示上述三个约束条件.求最佳路径的问题可以归结为满足约束条件r时,求最佳路径，使得沿路径的积累距离达到最小值即D（）warp2.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），这时此路径的积累距离为：D（n- 1,m）D[（n,m）]=d[x,y]+min{???D（n-1,m-1）D（n,m-1）这样可以从（n,m）=（30,30）出发搜索（n,m），对每一个（n,m）都存储相应的距离，这个距离是当前格点的匹配距离与前一个累计距离最小的格点（按照设定的斜率在三个格点中进行比较）•搜索到（1,1）时，只保留一条最佳路径•如果有必要的话，通过逐点向前寻找就可以求得整条路径.DTW算法可以直接按上面描述来实现，即分配两个NXM的矩阵，分别为积累距离矩阵D和时刻匹配距离矩阵d，其中时刻匹配距离矩阵d（i,j）的值为测试股票的第i时刻与参考股票的第j时刻间的距离.D（N,M）即为最佳匹配路径所对应的匹配距离.2・4范例例如我们选取两支股票，一支是中成股份（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)第二步：标准化经过Z经过Zi=xi-min(x?)后,max(x?)-min?(x?)0.820)0.917)X=(0,0.324,0.682,1,0.820)0.917)Y=(0,0.293,0.896,1,第三步：构建距离矩阵并寻找最优路径在这里我们让d（x,y）=|x-y|]j0.3240.8200.5270.0760.3240.8200.5270.0760.1800.097/10.7070.104痢-00.0830.6820.389T0.214/.0.3180.2350.593 打0.031 0.572 0.67600.2930.89610.917最终我们得到一个最优累积距离Dwarp（）=0.046三、DTW在K线图相似度计算机的实际应用【2】在中国A股总共有3470支甚至还有中小企业板，乃至创业板的股票，再涉及到各自的历史数据，这是一个非常庞大的数据库.若要在这个数据可当中找寻于某一支股票相似的其他股票，将是一个计算量非常大的工作.显然不可能人工手算，而是要借助计算机的索引技术.但是如此庞大的数据对索引技术的准确性和高效性有着非常大的挑战.3.1索引技术索引技术是一种快速的文件访问技术，举个例子，好比我们查字典，会根据偏旁部首在索引表中找到这个字在第几页一样，索引技术就是根据事先记录好的能够代表属性的取值将它与物理地址关联，进而在数据库中搜索.之前我们着重提到过一个命题一一D（）不符合三角不等式•这个事实对于warp我们可以用于索引的方法有着重要的意义：任何隐含或明确地符合三角不等式的索引技术都不能避免产生错误的排除.这是一个非常严格的要求：所有的空间访问方法，以及所有使用距离、度量、有利位置树的方法都不能避免错误的排除.唯一保证不会被错误排除的方法是顺序扫描，，对于像股票这种大量长序列集合来说这将是计算量巨大的.为了解决这个问题，人们提出了一种方法，用来避免一小部分的错误排除的同时能够实现的索引的加速.也就是说，人们的目标是提供快速的索引技术，同时尽可能地避免错误的排除.在这里人们引用了两种方法.3.2基于FastMap的访问方法我引用的第一种技术是基于一种名为FastMap的方法【4】它的工作原理如下：给定N个对象和一个距离函数，它将对象映射到k维空间中的N个点，以便保持原始距离，参数k可以由用户给出，或者人们可以在应用中调整以获得更好的系统性能•关键在于是假设对象确实能够表示成n维空间中的点，并且能够仅使用距离函数将这些点投影到相互正交k维空间中（k《n），从而实现降维的目的.如一个三维向量（x,y,z）可以经过函数映射到二维空间（x,y）换用我们如今正在研究的股票来讲，就是我们要将我们包含 30天收盘价格的向量，通过一个距离函数映射到一个二维空间•从而减少数据的处理.在将对象映射到k个点之后，人们可以使用任何空间访问方法来组织它们并搜索范围查询.FastMap在对象的数目N（即序列）上是线性的•而且，将查询序列映射到第k个点需要O（k）时间，也就是说，相对于数据库大小N，时间是恒定的.与其他方法一样（参见命题1），如果不遵守三角不等式，FastMap可能会引入错误的排除•我们观察到，如果我们使用原始距离的平方根即欧几里得距离（很明显它符合三角不等式），我们可以避免更多的错误排除•因此，我们在选择一个距离函数时候使用的就是欧几里得距离.

algorithnfrzm^s^sQarcb.W1w0i/*nasp^nse*//*filtering +/GIysa讥foreacbsecjuenGO玄LnEh七 、ifK巧〔門讪中日哥4《#).thzaddit&/?;/*pGSt^prOCa&filTL.gS-tep*/Forea^hiinR*if(PJ}>*)*TlisB ifr^tnfiiReporr/?;导midlj^xlthiii图来自：【2】算法1描述了如何使用FastMap处理范围查询■如果将FastMap应用于平方根距离，则搜索范围也应该平方根■注意F(s?)表示序列的坐标■在过滤步骤中，根据欧几里得距离而不是时间弯曲距离来比较两个序列，这是由于欧几里得计算的简单性和快捷性，这相当于一步粗筛选•在这一步中不相关的序列被排除在外■—些不符合条件的序列可能包含在内，但在后处理步骤中将删除这些序列.由于两个原因，算法1比单纯的DTW方法更快■首先，它扫描较少的序列■其次，过滤步骤也更快，因为k比序列长度小得多(通常是一些固定常数，比如6)■过滤可能会删除一些合格的序列，导致错误排除，因为我们不能保证在k-d空间中的欧几里得距离越低越好.即使我们使用时间弯曲距离的平方根，情况也是如此，但在实践中错误排除的可能性非常低，我们将在后面看到.3.3下边界技术对于两个给定的序列?=xVX],…,x>和y?=vy],…,y>，让max(?x)和max(y?)表示x?和y?中的最大值.min(x?)和min(y?)的定义相似，但是最小值■—对vmax(x?),min(x?)>定义了序列x?可以控制的范围.warp不失一般性，我们假设max(x?)>max(y?)所提出的方法受以下观察的启发:观察1|max(x?)-max(y?)|<D(x?;?y).warp检查其有效性相当简单■由于max(x?)应该至少匹配y?的一个元素，比如y.,并且我们假设max(x?)>max(y?),|max(?x)-max(y?)|< |max(x?)-y|<D(x?;?y)：i warp这个观察的结果是最大值的绝对差值可以作为一个距离来限制时间弯曲距离■虽然这是真的，但是，它可能不是很有用，因为它可能低于时间弯曲距离太多■因此，我们的目标是找到一个更紧密的下界.我们考虑两个序列的范围的可能排列进行比较.观察2给定两个范围观察2给定两个范围Rx?=vmax(?x),min(?x)>和Ry?=vmax(y?)，min(y?)>(范围的可能排列(见下图))1.overlap:Rx?1.overlap:Rx?和Ry?重叠 (min(x?)Wmax(y?);min(x?)$min(y?)).2.encloses:Rx?包含Ry? (min (?x)vmin(?y))•3.disjoint:Rx?和Ry?是不相交的(min (x?)>max(y?)【2】可以保证范围查询和最近相邻D距离(维度k上Ry■RaKyRy(dIdbjoinc图来自:我们有三种可能的排列iva1.鼻|t|——fimrlyjVi—m«B(f}|n/rJHji—■ 1|jT]— 1讥一 ■【心|ijfhFWKJ"1昭1I舟亠鮎卄7i|rij冊图来自：【2】值得注意的是，当通过扫描一次将序列输入到数据库中时，可以计算序列的最小值和最大值，并且可以将序列存储以供将来使用•而且，通过简单的比较，可以在恒定的时间内确定两个序列的范围的排列•最后，D的定义只需要对每个序lb列进行一次扫描，因此我们可以计算序列长度中线性时间内两个序列之间的D距lb离•这会令速效率有很大的改善，除非D太低估D・lb warp定理1对于任何两个序列X?=VX],・・，x m>和y?=<y n>,y nDib(x?,?y)WD(x?,?y)lb warp证明：见附录A.作为定理1的直接结果，我们得到以下推论.>，如n推论1(没有误诊率)对于任何两个序列X?=<X],・・,X >和>，如nTOC\o"1-5"\h\z& & y果D(x?;?y)W，贝UD(x?;?y)W・warp lb实际距离D与下边界距离d是一种限制条件,2 lb查询不会被错误排除[8]・在过滤步骤中，快速滤除不相关序列，因为可以快速计算lb的线性时间，通常为kW10)•一些不合格的序列可能包含在这个步骤的结果中,因为D小于我们的动态时间弯曲距离D•但是，这些不合格序列在后处理步骤lb warp中被删除.具体算法如下：算法2蛊lgqrlthiq1命pe「一BquacH訥吕蔬a{}；/♦riltaringstep*/Civenry,/o-re^cbsequence>inthedatabase、fif（P»i J►thenaddito{?;/*posl^precessing3t«p*/F&reicLiinR,if（jP.“就讥打IA，）■th餌removeifrom商；fl*portH;endalgorithm图来自：【2】3.3结合两种技术【5】假设有一个股票数据库包含许多任意长度的时间序列并且我们想要找出与某个查询序列相似的所有序列，即时间弯曲距离单位内的所有序列.这种类型的查询是有效范围查询.处理这种查询的直接方式是扫描所有序列并为每个扫描序列计算D （）,以选warp择符合条件的那些序列.虽然非常简单，但速度可能很慢，因为它读取数据库中的每一个序列（因此速度很慢），并且它计算每个股票数据库序列的时间弯曲距离.这个问题的独特之处在于，不仅I/O成本（第一种情况）很高，而且还有计算成本（第二种情况）.因此，我们所需要的技术应该解决这两个问题•为了解决这些问题，人们采用了之前提到的两种技术.1•使用FastMap构建索引结构以加快查询处理速度.这种技术可能会导致一些错误的排除.2•使用低于时间弯曲距离的下边界距离函数. 这种方法保证不会有错误的排除.3•使用这两种技术的组合.由于这两种技术彼此独立，因此可以采用流水线方式进行组合.仔细考虑前面两部分提出的两种技术可以带来更高效的算法.在算法1中，我们仅使用时间弯曲距离筛选过滤的序列.但是，我们可以在计算时间弯曲距离之前使用下边界距离.如果序列确实是合格的序列，这就变成了额外的成本，但是如果不合格就可以节省大量的计算成本. 这样我们就有了一个灵活的多阶段查询处理系统，如图3所示，其中FastMap和D分别作为主要和次要过滤器.lb图来自：【2】它由三个阶段组成，它们以流水线方式连接.第一阶段仅使用FastMap索引排除不相关的序列.在此阶段进行过滤可降低I/O成本和CPU成本.那些通过第一过滤阶段的序列在下一阶段与D()的查询序列进行比较.最后，后处理阶段lb只选择那些真正符合条件的序列.四、计算机程序选取一支股票，例如新钢股份作为参考模板，将其某三十天的收盘价格以向量方式输入.然后另选取三十支股票作为股票池S，股票池中的每一支股票S作为测试模板i与参考模板进行配对，以期待获得两条如下图般相似的股票.将所有股票数据进行单位化xi-min?(xi-min?(x)Z.= imax(x)-min?(x)例如4用D（欧几里得距离）【3】进行第一次筛选2) 2J( )D=x1-y1+?+(xn-yn)A25用D进行第二次筛选lb般票池r >11取那试股孕別输人藝再陨空W 1 戢拥处这T邀行口“"曾址/k进行九皓选二-一一=—丿 输出R中所有的序列s的D (X,s)i warp i6值辅巴5i和凡卄具体程序见附录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,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.0400.017，0.035，0.045，0.052，0.057 ，0.059，0.064，0.0520.061，0.071，0.054，0.0750.096，0.233，0.385，0.552，0.735，0.937，1.000，0.887）Y=(0.036，0.03