R语言对高频交易订单流进行建模分析 1

一、实验介绍--点过程基础及模拟

1.1 实验知识点

• 泊松过程及其模拟
• Hawkes 过程及其模拟

1.2 实验环境

• R 3.4.1
• Rstudio

二、点过程基础

假设你蹲在一个交通站台后面，看着人来人往。你觉得乘客的到达似乎存在某种数学规律， 于是你把每个人到达的时刻记录了下来。

有什么办法可以对这些人到达的时刻进行建模？你渐渐进入了沉思状态。

也许提炼这些点形成的集合所具有的特征是一个好办法。

你想到乘客到来的速率肯定是一个重要特征，如果是在一个偏僻的小公交站，可能半天也看不到一个人 ； 到了市中心的大车站，人潮涌动可能让你难以计数。

不同人到来的间隔是另外一个有意思的特征，乘客不是从工厂出来的产品，肯定不会乖乖地等间隔的到来，那么不同乘客到来的间隔有什么规律呢？

要回答这些问题，必须要借助概率的语言，更确切地说，是点过程。下面我们就从泊松过程开始我们的旅程。

2.1 泊松过程

泊松过程有以下几个性质：

• 不相交的时间段上到来的数量是相互独立的；
• 两个点几乎肯定不会同时到达；
• 在某个给定的时间段到达的数量服从泊松分布，分布均值正比于时间段的长度。

我们从数学层面来描述具有这样性质的过程，首先我们从第二个性质开始 ， 我们用 N(a,b] 表示 a < t <= b 这段时间发生的事件数。

对于一个趋近于0的 Δt , 我们声明 , 对任意t

P(N(t, t + Δt]=1)=λΔt

由于 λ 的含义是单位时间内事件的数量 ， 所以可以定义为事件发生的强度。

由于两个点不会同时到达 , 在小段时间里发生两次的概率约等于0

P(N(t, t + Δt >= 2) → 0

那么对于任意时间段 (a,b] , 我们可以先将其划分为多个小时间段 ， 然后由不同时间段的 独立性， 用二项分布来计算概率分布 ， 再用泊松分布近似：

P( N(a , b] = k ) = (n!/(k!(n-k)!) ) (λΔt)^k (1-λΔt)^{((b-a)/Δt) - k)} → (λ(b-a))^ke^(-λ(b-a)) / k!

可以看到 N(a,b] 近似服从参数为 λ(b-a) 的泊松分布。

现在我们来看一看两个事件的间隔服从什么分布 , 间隔为 $ t_0 $ 意味着这$ t_0 $时间段没有事件发生， 那么可以很容易的进行计算：

P(interval > t₀)= P(N(t, t + t₀]=0)= e^−λt₀

我们可以计算出 interval 的密度函数：

f(interval = t₀)=λe^−λt₀

可以发现这正是指数分布。

泊松分布还有很多有趣的性质 ， 在这里我们就不一一列举了 ， 下面来看一下如何对泊松分布进行模拟 。

我们可以利用间隔服从指数分布的性质，模拟服从泊松过程的事件,第k个事件的时刻就是第 k-1 个事件时刻加参数为λ 的指数分布的随机变量：

举个例子， 我们可以模拟一个 λ 为 0.5 的泊松过程 ， 总共模拟 50 个事件， 可以画出事件与时间的关系：

x<-cumsum(rexp(50,rate=0.5))y <- rep(0,50)plot(x,y,main="events arrival")

还可以画出累积事件与时间的关系 ， 按照我们的估算 ， 发生 50 个 λ 为 0.5 的事件大约要用 100 的时间，我们可以从图中进行验证。

y_cum<-cumsum(c(0,rep(1,50))) plot(stepfun(x,y_cum),do.points = F,main="cumsum of events")

2.2 Hawkes过程

在泊松过程中，强度保持恒定，事件的发生遵循“无记忆性”的原则，在现实世界中，很多情况 都不符合这样的假定,例如犯罪行为往往具有空间上的聚集性，这是由于罪犯在得手后倾向于在附近继续作案；而在高频交易中趋势交易者会跟踪大订单,使得市场在短时间涌入大量订单。

在这些系统中，事件发生的速率都是不均匀的。如何描述这种空间上的聚集性，或者说是正反馈的机制呢？

我们需要对模型进行扩展，不再把 λ 固定为一个确定的值，而是让他成为一个关于时间的函数，即 λ(t) 。

比较精确的定义是当 Δt 趋近为 0 时：

P(N(t, t + Δt]=1)=λ(t)Δt

其余的假设相似 ， 在小间隔内发生2次或以上事件的概率趋近为0 。

而 λ(t) 的定义则为：

这个式子中的 λ₀(t) 代表的是背景的强度 ， 而 v(t - t_i) 则代表发生在 t 时刻之前的事件对时刻 t 产生的正向影响 ， v 函数就是核函数 ， 简单来说 ， t_i离t 越近 ， 对t时刻造成的影响就应更大。

我们先使用一个比较简单的核函数：指数函数来看一看 Hawkes 过程究竟有什么特性：

可以定义

v(t − s)=αe^{−β(t − s)}

我们把背景强度设为恒定的 μ ， 那么强度随时间变化的函数可以表示为

λ(t)=μ + ∑_ti < tαe^{−β(t − s)}

更详尽的介绍可以参考slides

理论框架搭设完毕，我们尝试来模拟一下 Hawkes 过程，这里我们使用 Ogata(1981) 的方法进行模拟，他的方法的大意是先用一个 λ 的上界作为参数喂到指数分布里给出一个间隔 ，进而算出新的事件发生的时间。 但是实际情况在间隔里 λ(t) 是应该不断变小的（核函数的距离不断拉长），所以我们只能以一定概率来接受这个新的点，如果被拒绝应该模拟新的点。换句话说，有一些模拟出来的点需要被删掉， 使得间距和真实的更接近。更详细的解释可参考原论文 [1]：

假设μ = 1.2 , α = 0.6 , β = 0.8。

set.seed(728)n <- 100params <- c(0.5 , 0.3 , 1.2)mu <- params[1]alpha <- params[2]beta <- params[3]lambda_intensity <- function(event , t , params) {  mu <- params[1]  alpha <- params[2]  beta <- params[3]  partial <- event[event < t]  lambda_int <- mu + sum(alpha * exp(-beta *(t - partial)) )  lambda_int}event_time <- numeric(n) for (i in 1:n) {  if (i==1) {    lambda_star <- mu    event_time[1] <- rexp(1 , rate=lambda_star)  }  else {    lambda_star <- lambda_intensity(event_time , event_time[i-1] , params) + alpha    event_time[i] <- event_time[i-1] + rexp(1 , rate=lambda_star)    while ( runif(1) > lambda_intensity(event_time , event_time[i] , params)/lambda_star) {      # thinning procedure , delete this point and make new       lambda_star <- lambda_intensity(event_time , event_time[i] , params)      event_time[i] <- event_time[i] + rexp(1 , rate=lambda_star)    }  }}

我们模拟了 100 个事件，可以把事件和对应的强度画在一张图上：

time <- seq(min(event_time) , max(event_time) , 0.001)intensity <- unlist( lapply(time , function(x) lambda_intensity(event_time , x , params)) )plot( event_time , rep(0,n) , col="blue" , ylim=c(-1,5) , main="simulated hawkes process" , ylab="lambda intensity")lines(time , intensity , type='l' ,col="red")

可以通过数值验算一下我们的结果是否合理：

强度的实际平均值为 N_t / t大约为 100/150 ,

而我们可以推导理论的平均值：

几个需要注意的点

• 图片里的 λ₀ 是我们定义的 μ ， 而图片里的 μ 代表我们要算来和实际对比的强度的平均值 E[λ(t)]

• 注意从第一排到第二排的变换是根据 λ(t) 的定义。

所以可以计算出理论值

E[λ(t)] = μ / (1 - (α/β)) = 2/3

与实际的较为接近,说明模拟是正确的。

三.总结

在上面的课程中 ， 我们从理论层面介绍了泊松过程和 hawkes 过程 ，以及它们各自的模拟方法 ， 下面的一个问题是 ， 当我们拿到一个数据时 ， 如何去拟合出参数 ， 这也是下一节的主要内容。