一、實驗介紹--點過程基礎及模擬

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 <- 100
params <- 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 過程 ，以及它們各自的模擬方法 ， 下面的一個問題是 ， 當我們拿到一個數據時 ， 如何去擬合出引數 ， 這也是下一節的主要內容。