numerical update rule for discretized hawkes excitation process

Question

So I think I am just misunderstanding some simple notation or something and would appreciate some help. I am trying to replicate this model in an agent based model, but I cannot seem to figure out the right update rule based upon the equations.

I am looking at a paper by Short et. al., entitled "A Statistical Model of Criminal Behavior," where the authors examine a self-exciting process for hotspot formation. So basically there is a regular lattice and at each node point there is a house. Each house has a level of attractiveness for burglary and this level of attractiveness is subject to excitation. But the equations here are throwing me off a bit. The question below defines the level of attractiveness of a single house $s$ at time $t$.

$$ A_s(t) = A_s^0 + B_s(t) $$

So $A_s(t)$ is the attractiveness of house $s$. Then $A_s^0$ is the average level of attractiveness for burglarization of house $s$ (or any other house in the neighborhood). $B_s(t)$ is the stochastic level of attractiveness which can rise for some period of time and then decay back to zero--and the house returns to the average level of attractiveness.

Now, the authors define $B_s(t)$ using the update rule:

$$ B_s(t+\delta t) = B_s(t)(1-\omega \delta t) + \theta E_s(t) $$

So this seems to be a recursively defined function where $B_s(t+\delta t)$ is a function of $B_s(t)$. The $\omega$ represents the rate at which the effect of the excitation diminishes over time as we get farther and farther from some excited level of $B_s(t)$. $E_s(t)$ is just the number of triggering events that increase the excitation, and $\theta$ is a coefficient for the effect of each triggering effect on the excitation component $B_s(t)$.

So here is the problem I am having. I don't understand how to put the update rule into the original equation. I know this is wrong, but I cannot seem to figure out the right way to fix it. I think it is just in the distinction between $t$ and $\delta t$. So if I was to just do something like this:

$$ A_s(t_{n+1}) = A_s^0 + B_s(t_n)(1-\omega t) + \theta E_s(t_n) $$

I mean it seems as if there are 2 time scales happening simultaneously: the $[t: 1,2,...n]$ and then the $[\delta t: 0.1, 0.2, 0.3 . . . 0.9]$ between the $t's$.

So the question is, if this is the wrong update rule, then what is the correct one. Do I have to run an inner loop and outer loop for the $\delta t$ and $t$ respectively? Any suggestions would really be appreciated. Also any references on where to look for stuff like this--as I believe this scheme comes from statistical physics.

Answer 1

I was able to get this to work. So I created the following code.

import numpy as np
from scipy.stats import bernoulli
import matplotlib.pyplot as plt

## set some parameters for the simulation
omega = 0.3
a0 = 0.1
theta = 1.0
p = 0.1
xmax = 50
xsamp = 200

# initialize the simulation
x1 = np.linspace(0, xmax, num = xsamp)
b = np.zeros(len(x1))
a = np.zeros(len(x1))
event = bernoulli.rvs(p,  size = len(x1))

# generate the excitation results
for i in xrange(len(x1)):
    try:
        val = (1-omega)*b[i] + theta*event[i]
        b[i+1] = val
    except IndexError:
        break
a = a0 + b

plt.plot(x1, a)
plt.show