# A Stochastic Model for the Early Stages of Highly Contagious Epidemics by
using a State-Dependent Point Process

Jonathan Chávez-Casillas Department of Mathematics and Applied Mathematical
Sciences, University of Rhode Island, USA (jchavezc@uri.edu).

###### Abstract

The recent COVID-19 pandemic has shown that when the reproduction number is
high and there are no proper measurements in place, the number of infected
people can increase dramatically in a short time, producing a phenomenon that
many stochastic SIR-like models cannot describe: overdispersion of the number
of infected people (i.e., the variance of the number of infected people during
any interval is very high compared to the average). To address this issue, in
this paper we explore the possibility of modeling the total number of
infections as a state dependent self-exciting point process. In this way,
infections are not independent among themselves, but any infection will
increase the likelihood of a new infection while also the number of currently
infected and recovered individuals are included into determining the
likelihood of new infections, Since long term simulation is extremely
computationally intensive, exact expressions for the moments of the processes
determining the number of infected and recovered individuals are computed,
while also simulation algorithms for these state-dependent processes are
provided.

Keywords: SIRQ-model, overdispersion, Self-exciting process, Quarantine
distribution, stochastic-intensity, state-dependent.

##  1 Introduction

As the, now seemingly endemic, disease of COVID-19 is being incorporated into
our daily lives, we have been reminded that epidemics are a phenomena that can
expand quickly and have devastating economic and social consequences. Highly
contagious diseases can produce huge damage, and so the development of
accurate and effective models for epidemics is important. However, as noted in
Adiga et al., (2020), different kinds of models are appropriate at different
stages, and for addressing different kinds of questions. For example, some
statistical methods based on machine learning techniques are very useful in
predicting the behaviour on the short-term. However, they are not very
effective for the long term predictions nor for describing the evolving
circumstances. Simple compartmental-type models, and their extensions, that
is, structured meta-population models, are useful for several population-level
questions. However, once the outbreak has spread, and complex societal
interactions are at play, stochastic agent-based models could be seen as a
more robust and effective tool, since they allow for a more systematic
representation of complex social interactions, individual and collective
behavioral adaptation, and public policies. With this idea of trying to be
more precise in describing different stages of a pandemic the current paper
tries to shed some light into how to incorporate certain salient features
observed in the data to model future outbreaks of epidemics with a high
control reproduction number.

Stochastic modeling of epidemics becomes relevant when the environmental or
demographic variability -such as transmission, recovery, births, deaths, or
environmental impacts- are too large to account for in a deterministic model
(see Allen, (2017)).

The ongoing pandemic has returned epidemic modeling at the forefront of
worldwide public policy-making Bertozzi et al., (2020). Many scientist from
different fields, particularly medicine, biology, mathematics, physics and
chemistry (for a few, non-exhaustive, list of examples, see Cheng et al.,
(2020),Chiang et al., (2022),dos Santos Gomes and de Oliveira Serra,
(2021),Franco, (2020), Hazarika and Gupta, (2020), Liu et al., (2022), Van den
Driessche, (2017))

Due to their versatility and capabilities to model properties that regular
stochastic SIR models can’t, point process models have been recently explored
as a viable option. Models using these processes are data driven and are
flexible enough to allow for parametric or nonparametric estimation of the
reproduction number and the time scale at which the contagions occur (see
Bertozzi et al., (2020). Furthermore, at some level, they can also be viewed
as stochastic versions of popular compartmental models used in epidemiology.
Since these processes have just been recently introduced, before stating the
model assumptions and results, in remainder of this introduction, an
introduction to point processes is provided as well as how they have recently
been used to model the novel COVID-19 epidemics.

###  1.1 An introduction to Point processes

Point processes can be understood roughly as a random set of points in a space
$\mathcal{X}$. The topology and properties of the space can be very general,
making them a very versatile tool for modeling different phenomena. In this
setting, these points, may represent anything, but the most common usage of
point processes are occurrences of events in time, location of objects in the
space or a mixture of the previous two. Some of the most recognized and
applied types of point processes are:

  * •

_Spatial Point process:_ Here, the space $\mathcal{X}$ represents some fixed
“geographical location”, generally, $\mathcal{X}\subseteq\mathbb{R}^{n}$. In
this instance the number of events occurring and their location are random.

  * •

_Temporal point process:_ In this case, the space $\mathcal{X}$ is a totally
ordered set. Generally, $\mathcal{X}=[0,T]\subset\mathbb{R}$ and each point
represents the times at which certain events occurred. In these kind of
processes the number of events and the times at which those events occurred
are random.

  * •

_Spatio-Temporal point process:_ Here, the space is a combination of the
previous two cases and each point represents a place and time for an event and
usually $\mathcal{X}\subseteq\mathbb{R}_{+}\times\mathbb{R}^{n}$.

It also important to mention that the previous cases are sometimes called
“unmarked point processes”. However, sometimes, we want to measure or model an
additional characteristic of the point (other than time or location). In this
instance, we can attach a “mark” to each random point to classify them as
belonging to a certain class. This flexibility comes almost at no cost, since
“marked point processes” are merely a generalization where each point belongs
to a space $\mathcal{Y}=\mathcal{X}\times\mathcal{K}$, where $\mathcal{X}$ is
one of the spaces described above and $\mathcal{K}$ is the so-called mark
space.

Generally speaking, when talking about temporal point processes there are
mainly two -not disjoint- ways to describe them. The first one was born in the
neurophysiologist community where point processes are seen as random variables
in a complicated functional space and are characterized and analyzed via their
moments (such as the average, variance, skewness, kurtosis, etc.). The
advantage of this point of view is that these moments can be associated to
statistical quantities which can then be estimated, fitted and calibrated
using data. The second way to characterize point processes is through its
stochastic intensity, which provides a summary, at any given time $t$, of the
likelihood that some new future event arrives in the time interval $[t,t+h)$
given its history, i.e., given the times of all past events up until time $t$.
This notion of intensity has been used extensively because of its analytical
tractability. The resulting family of processes that possess an intensity
function is called (stochastic) intensity-based point processes, and actually,
contains almost all the point processes which have some practical interest.
Exhaustive treatments containing all the important probabilistic and
statistical features of Point processes can be found in, e.g. Daley and Vere-
Jones, (2003, 2008); Brémaud, (2020).

A realization of a point process over $[0,\infty)$ is a sequence
$\\{T_{0}\\}_{n\geq 1}$ in $[0,\infty]$ such that

$T_{0}=0,\qquad\qquad T_{n}<\infty\qquad\text{and}\qquad T_{n}<T_{n+1}.$

and the point process is considered to be non-explosive if
$\lim\limits_{n\to\infty}T_{n}=\infty$. Further, for each realization of the
point process, there is a counting function, or counting process, $N_{t}$
defined as

$N_{t}=\left\\{\begin{array}[]{rcl}n&\text{ if }&t\in[T_{n},T_{n+1}),\ n\geq
0\\\ \infty&\text{ if }&t\geq\lim_{n\to\infty}T_{n}\end{array}\right.$

The above definition is not that intuitive, but basically, the counting
process $N_{t}$ tells us how many events have occurred up to time $t$. To see
this, we can also define the inter-arrival times between events
$\tau_{n}:=T_{n}-T_{n-1}$ for $n\geq 1$ and define $\tau_{0}:=0$. Thus,
$N_{t}$ can be written as

$N_{t}=\max\limits_{n\geq 0}\\{\tau_{0}+\tau_{1}+\ldots+\tau_{n}\leq t\\}.$

Many of the point processes that are found in the literature are simple, which
means that only one event arrives each time. That is,

$\lim\limits_{h\to 0}\mathbb{P}\left[\frac{N(t+h)-N(t)>1}{h}\right]=0$ (1)

A point process is called _regular_ if it possess an intensity function
$\lambda(t)$, which is defined, when the limit exist, as

$\lambda(t)=\lim\limits_{h\to
0}\mathbb{E}\left.\left[\frac{N(t+h)-N(t)}{h}\right|\sigma(N_{s})_{0\leq s\leq
t}\right],$

where $\sigma(N_{s})_{0\leq s\leq t}$ denotes the past history of the process
$(N_{t})_{t\geq 0}$. Further, whenever the point process is simple, the above
definition is equivalent to,

$\lambda(t)=\lim\limits_{h\to
0}\mathbb{P}\left.\left[\frac{N(t+h)-N(t)>0}{h}\right|\sigma(N_{s})_{0\leq
s\leq t}\right]$

Notice that by definition, as $\lambda(t)$ is a conditional expectation over a
sigma algebra, it is itself a random process, but sometimes, it can be a
deterministic function of time as in the case of the Poisson process. Since
probabilities are always non-negative, the intensity function is always non-
negative and thus the cumulative intensity function,
$\Lambda(t)=\int_{0}^{t}\lambda(s)ds$ is a non-decreasing function. The
cumulative intensity function always exists (even if the point process is not
regular). In fact, it can be shown that the cumulative intensity function is
the compensator of the point process $N(t)$.

A technical remark might be in place here: since a temporal point process is
always a càdlàg (i.e., right continuous with left limits) submartingale, it
can always be then decomposed as $N(t)=M(t)+A(t)$, where $M$ is a local
martingale and $A$ is an increasing predictable process called the
compensator. As can be inferred by the name, it turns out that
$\Lambda(t)=A(t)$. Moreover, whenever the measure generated by $\Lambda(t)$ is
absolutely continuous with respect to the Lebesgue measure, the “normal”
intensity function $\lambda(t)$ exists and so the limit above exists.

To explain the above definitions and before proceeding to the canonical
example of interest we can look at the simplest and most common point process:
the Poisson process. In this case, the intensity is completely deterministic,
i.e. $\lambda(t):[0,\infty)\to[0,\infty)$ and

$\mathbb{P}[N(t)-N(s)=k]=e^{-\int_{s}^{t}\lambda(u)du}\frac{1}{k!}\left[\int_{s}^{t}\lambda(u)du\right]^{k}.$

Since the intensity is deterministic we always know the infinitesimal
probability of a new event. That is, the probability that we observe a new
event in the interval $[t,t+h)$ is roughly $\approx h\lambda(t)$ and in this
case this probability is actually independent of the process’ history. In this
case, it can be shown that $\mathbb{E}[N_{t}]=\int_{0}^{t}\lambda(s)ds$.

The Poisson process and many others point processes are good to model events
where the rate at which new events appear is independent of how many or when
past events occurred, but sometimes the phenomena that wants to be analyzed
exhibits a “clustering effect”, where the appearance of one new event will
trigger the occurrence of more new events. An important example that possess
this property, and a process that will be used thoroughly during this
proposal, is the family of Hawkes processes, which was first studied
thoroughly by Ogata (see Ogata, (1983)) to model the occurrence of earthquakes
in Japan. His motivation to use this type of processes was that usually after
a big earthquake, many replicas will follow. In the Hawkes Processes family,
the (stochastic) intensity function, $\Lambda(t)$, “feeds” itself from the
point process $N_{t}$. That is, the intensity increases when a point arrives,
which in turn will trigger more points to arrive. In the most classical
example, the intensity can be characterized by

$\Lambda(t)=\lambda_{0}+\int_{0}^{t}\phi(t-s)dN_{s}=\lambda_{0}+\sum\limits_{T_{i}<t}\phi(t-T_{i}),$
(2)

where $N_{t}=\\{T_{1},T_{2},T_{3},\ldots,T_{N_{t}}\\}$ is itself the random
point process and $T_{i}$ is the time of the $i-$th occurrence. As mentioned
before, the intensity is also stochastic (random) but the source of its
randomness is not exogenous but comes from the point process itself. Moreover,
if $\phi(t)=Qe^{-rt}$, then the process $(N_{t},\lambda_{t})$ is a Markov
Process (a process whose evolution depends solely on the current state the
system and not in the history of how it arrived to such current state) and the
intensity function decays exponentially between events. As before, we can also
compute the expected number of events in the interval $[0,T]$ which turns out
to be $\mathbb{E}[N_{t}]=\frac{\lambda_{0}}{1-\int_{0}^{\infty}\phi(s)ds}$.

Figure 1: Six realizations of point process with $T=20$ and
$\mathbb{E}[N_{T}]=25$. On the left, three different inhomogeneous Poisson
processes with intensities proportional to a constant, $x^{2}$ and $\sin(2x)$
are shown, whereas on the right three realizations of the same Hawkes process
are displayed.

Figure 1, tries to illustrate the concepts and processes described above. In
all the plots the intensity function and a realization of the point process
are shown. At the same time, all the processes below were “standardized” by
setting the expected number of points to be the same in all cases. The left
column shows 3 different inhomogeneous Poisson processes with different
intensities. Since the intensity is deterministic, it does not matter how many
simulations of the point process are performed, the intensity will the the
same. As it can be seen, when the intensity is higher, the likelihood of more
points increases. In the right column, 3 different realizations of the same
Hawkes process are plotted. In this case since the intensity is stochastic, it
will change between realizations, but more importantly the clustering
phenomena can be observed. Whenever a new event arrives, the intensity spikes
and thus the likelihood of more events happening increases and to balance the
process, in between the arrival of events the intensity decreases
exponentially reducing the likelihood of more events coming. However, there is
always a baseline intensity for which the intensity process cannot go below
guaranteeing that there will be a new event at some point. To illustrate this,
we present next a heatmap of the Hawkes process (i.e. the process counting the
number of events up to time $t$) and one of its intensity.

(a) The Hawkes process

(b) The intensity process

Figure 2: Heatmap of the Hawkes process and its intensity. The frequency is
given in a logarithmic scale.

###  1.2 Linking point processes to Epidemiology

Epidemic models have been trying to model, understand and predict different
features and generalities of the pandemics that humanity has lived with.
However, it was not until around 1930 that the first stochastic epidemic model
was created. Before that, all the models were deterministic were based on
creating different systems of ODEs, but some desirable random effects were
missing. Nonetheless, those deterministic models provided the skeleton upon
which the behavior of corresponding stochastic systems are built. Furthermore,
many times, when populations are large, the limiting behavior of such
stochastic systems converge to a corresponding system of ODEs. The main models
discussed in this proposal will be variations of the famous stochastic
susceptible-infected-recovered (SIR) models and belong to the class of the so-
called compartmental models, where each individual is placed into one of the
compartments and the individuals move between the compartments under
predetermined rules. For a survey on the description of many stochastic
models, see Greenwood and Gordillo, (2009).

The most common stochastic epidemic model, and the building block for many of
the more detailed and precise models produced nowadays, is the so-called
general stochastic model,

  * •

There are three compartments where an individual can be: Susceptible, Infected
and Recovered. Moreover, each individual belongs to one and only one
compartment.

  * •

The number of susceptible, infected and recovered individuals at time $t$ are
denoted by $S_{t},I_{t},R_{t}$. Also, the population is constant at a level
$N$ and thus, $N=S_{t}+I_{t}+R_{t}$.

  * •

A recovered individual cannot contract the disease again.

  * •

During the interval $[t,t+\Delta t]$, one of the following things has to
happen:

    1. 1.

A susceptible individual gets infected with probability $\beta\frac{I_{t}}{N}$
because $I_{t}/N$ is the likelihood of contacting one of the current number of
infected individuals, $I_{t}$, multiplied by the likelihood $\beta$ of
contracting the virus once you got in contact with a person. Since there are
$S_{t}$ susceptible people at that moment,

$\mathbb{P}\Big{[}(S_{t+\Delta t},I_{t+\Delta
t})-(S_{t},I_{t})=(-1,1)\Big{]}=\beta S_{t}\frac{I_{t}}{N}\Delta_{t}+o(\Delta
t)$ (3)

    2. 2.

An infected individual recovers from the virus at any time with likelihood
$\gamma$. Thus, since there are $I_{t}$ infected individuals with the same
likelihood to recover,

$\mathbb{P}\Big{[}(S_{t+\Delta t},I_{t+\Delta
t})-(S_{t},I_{t})=(0,-1)\Big{]}=\gamma I_{t}\Delta_{t}+o(\Delta t)$ (4)

    3. 3.

Nothing happens and the system remains the same. This is the complementary
event to the union of the two actions above,

$\mathbb{P}\Big{[}(S_{t+\Delta t},I_{t+\Delta
t})-(S_{t},I_{t})=(0,0)\Big{]}=1-\left(\beta\frac{S_{t}}{N}+\gamma\right)I_{t}\Delta_{t}+o(\Delta
t)$ (5)

Many of the compartmental models that are available nowadays, to allow for an
easy simulation and mathematical tractability, remain within the framework of
the Markovian world. As such, those models have deterministic limits and
diffusion approximations. The main results and techniques to obtain such
limits are detailed in Kurtz, (1981); Meyn and Tweedie, (2012). These limits
can be used to estimate and calibrate the models quickly and constantly which
can be of great help for public policy makers, since many public health
policies might be influenced by predictions of how large an epidemic might be.

In a more abstract form, the model (3-5) above can be written in terms of the
difference of some Poisson processes and in fact, at their core, many
stochastic SIR-type models can be established as functions of Poisson
processes or as limits of them, which leads to an immediate question of
whether the point process leading the epidemic model can be generalized to
include more features observed empirically in the data. Besides, there have
been few models that try to explore the dynamics of the frequency at which
individuals get infected and its effects in the likelihood of future people
getting infected. That is, many of the current models either try to elaborate
on the number of compartments or in the finesse of the conditions required for
a certain individual to transition from one compartment to another but there
are few models that try to get a more realistic description of when such
transitions occur and how the frequency of such transitions affects the future
evolution of the system.

Indeed, in most of the current compartmental models such as (3-5), each
susceptible individual is equally likely to get infected and that probability
is solely dependent on two factors: the proportion of infected individuals and
the number of susceptible individuals at time $t$. To understand better
paradigm assume the following two scenarios

  * •

Assume a population of a 100 individuals and at time 0, there are 5 infected
people and 95 susceptible but at time $t=20$, it is observed that
$(S_{20},I_{20},R_{20})=(50,30,20)$.

  * •

Assume the same population as above with the same initial conditions of
$(S_{0},I_{0},R_{0})=(95,5,0)$ but at time $t=5$, it is observed that
$(S_{5},I_{5},R_{5})=(50,30,20)$.

As it can be seen, in both cases you start with the same number of infected
and susceptible people and at at some point in time you also have the state
$(50,30,20)$. The difference is that in the first case it took longer to reach
that point, whereas in the second, such state was reached much faster. When
this happens it might be better to assume that the driving process counting
the number of infected people can cluster due to the fact that in periods
where there are many infected people, the likelihood of more people to get
infected is higher than in periods where few people have been recently
infected. In this spirit, various empirical studies have revealed that when an
arrival processes shows display over-dispersion, i.e. the variance of the
number of arrivals in a given interval exceeds the corresponding expected
value, the standard assumption of having a Poisson process driving the number
of infections is not valid and new models are required. As Figure 3 suggests,
COVID-19 has proven to be a very contagious virus where a high level of over-
dispersion could be observed. This implies that rather than having a
stochastic epidemic model featuring Poisson processes, a Hawkes point process
might be a much more suitable candidate to model it.

Figure 3: Plots showing the over-dispersion in a 2,7 and 30 days window of the
number of newly infected people during COVID-19 in 3 US states: New York,
California and Rhode Island. The $x$ axis shows the time in days, starting
from February of 2020, while the $y$ axis shows the quotient between the
variance and the average of the new cases within that time period. A quotient
greatly larger than 1, indicates high over-dispersion in the “arrival process”
of new infect people.

Our first step will be the creation of a model that generalizes the Poissonian
flow. In contrast to deterministic models, rather than modeling the number of
infected people, $I_{t}$, directly the contagion process $C_{t}$ (counting the
number of infected people up to time $t$) is the process that needs to be
modeled. Reason being that $I_{t}$ can increase (with a new infection) and
decrease (with a recovery) but the counting process associated to a point
process must be non-decreasing.

One of the first attempts of incorporating Hawkes processes into an SIR model
is provided in Rizoiu et al., (2018). In there, the authors explain a
relationship between the general stochastic SIR model and a finite population
Hawkes-SIR model (they account for the fact that the population is capped at
$N$ and therefore $C_{t}\leq N$), but this model has some limitations. In
particular they only explore a model where the conditional expectation of the
intensity of $C_{t}$ conditioned on all the times at which there is a recovery
coincides with a finite Hawkes process with a baseline intensity
$\lambda_{0}(t)=0$. These restrictions cripples the model in the sense that a
true Hawkes-SIR model would require the intensity of $C_{t}$ to be of the form
(2) with $\lambda_{0}(t)\neq 0$ and also not be the result of a conditional
expectation on past and future unknown times.

Other models that have considered an epidemic and relating the total number of
infected people with a Hawkes process are presented in Escobar, (2020);
Garetto et al., (2021); Chiang et al., (2020) among others. Nonetheless, there
is missing a model that becomes a natural extension of an stochastic SIR
model. That is, there is no model in the literature where people are divided
into the $S,I,R$ compartments and people start arriving to the “infected”
compartment according to a Hawkes process. Recovered individuals should arrive
via a Poisson or renewal process, since contrary to infections, the recovery
of a person does not affect how other people recover. One of the aims of this
proposal is to create and analyze such model as well as find conditions under
which a disease introduced into a community will develop into a large
outbreak, and if it does, conditions under which the disease may become
endemic. This condition is linked to the so-called basic reproductive number,
$R_{0}$, defined as the expected number of secondary infective cases per
primary cases in a susceptible population.

As it can be noted, many of the models that have started to incorporate Hawkes
processes are very recent and inspired by the highly infectious rate of
COVID-19, making this a novel area to propose, analyze and create models that
can shed some light into the question of how fast this disease propagates,
hopefully contributing to the decision process on policy-making by providing
more accurate estimates when incorporating some of the salient feature of
these types of models.

Scope and Limitations of the current model: In its current form, the model
presented in this article provides a way to model the overdispersion and high
variability observed in the data. However, due to the limiting ability of
incorporating a stronger dependency on the current state of the epidemic and
the different rates at which different variants have shown to be transmitted,
this model is more suitable for keeping track of the pandemic more at the
beginning stage when not so many social and political reactions have occurred.
To account for this, the author is currently investigating a regime switching
alternative. Further, as it is, these type of models should only be used to
forecast a short to medium term prediction.

Organization of this article: In section 2, a brief description of the model
is presented while in Section 3 we compute the moments of the process which
can be used to predict the medium term dynamics of the current number of
infected people given the probabilistic distribution of the quarantined people
(which in terms can be used for public policy). Section 4 provides a brief
survey and ideas of how to perform efficient estimation on a state-dependent
self-exciting process. However, since estimation is out of the scope of this
paper, it will be treated in subsequent research. In Section 5, we provide
different numerical examples to explore how the different parameters of the
model play a role in it and their sensitivity. The conclusions and elements of
further research are presented in Section 6. Finally two appendices are also
provided. In Appendix A all the proofs to the technical lemmas and theorems of
Section 3 are provided while in Appendix B the pseudo-code of the different
algorithms used to simulate the different processes are provided.

##  2 A Brief Description of The Model

In the model presented in this article, an important assumption will be that
the the population is rather large (technically, it is assumed that the
possible number of susceptible people is infinite to facilitate the analysis
of the moments of the number of infected people $I_{t}$). This assumption
provides a tractable background for the problem and is not extremely
unrealistic given that the model is specifically designed to simulate a
pandemic in its early stages. Considering that up May 1st, 2022 only 6.6% of
the population has been reportedly infected111Information taken from
https://covid19.who.int/ (although it is of general consensus that this number
is very underestimated, specially now with the development of at-home test
kits), we can think of the number of susceptible individuals as a rather large
number and the model will still be accurate. In fact, a precise finite model
can be created, but there is little to gain and the closed formulas obtained
in Section 3 become very cumbersome. This compartmental model follows a
similar SIR dynamics than the models in Greenwood and Gordillo, (2009), where
$S_{t},R_{t}$ and $I_{t}$ represent the number of susceptible, recovered and
infected people at time $t$. Obviously, since, as discussed above, the number
of susceptible people is assumed to be infinite, we cannot have a classic
relationship such as $S_{t}+I_{t}+R_{t}=N$, but we will rather only focus on
the number of infected individuals $I_{t}$ and the number of recovered
individuals $R_{t}$.

The main modeling assumption in this paper will be that the (historical)
number of infected people up to time $t$, denoted by $C_{t}$, will be driven
by a (Hawkes-like) counting process with stochastic intensity

$\Lambda_{t}=\lambda^{\infty}(I_{t},R_{t})+\int_{0}^{t}Q_{s}\phi(s,t)dC_{s},$
(6)

where $\\{Q_{n}\\}_{n=1}^{\infty}$ is a sequence of i.i.d. random variables
independent among themselves and from every other process. As described below,
these random variables $Q$ will be interpreted as the level of quarantine that
the $n-$th infected person will have, modulating the probability that this
individual will produce future infections.

In this paper we will assume that the baseline intensity
$\lambda^{\infty}(i,\rho)$ of the Hawkes-like process $C_{t}$ given in
Equation (6) reacts to the state of the epidemics at time $t$. Indeed, we will
assume that if $I_{t}=i$ and $R_{t}=\rho$. Then,

$\lambda^{\infty}(i,\rho)=\lambda_{0}+i\log(\alpha)+\rho\log(\beta),$ (A1)

where $\alpha>1$ and $\beta,\lambda_{0}$ are positive numbers.

The form of the baseline intensity given in assumption A1 provides great
flexibility to consider multiple situations. We describe some of them next.

  * $\bullet$

A pure Hawkes approach. In this setting, the intensity takes the form of a
classical Hawkes process. This could be a basic model to use for explaining
the over-dispersion observed in the data. This can be achieved by setting
$\alpha=\beta=1$. In this sense, the baseline intensity will be independent of
the number of infected and recovered and the general intensity of the point
process will only depend on the cumulative number of infected people and times
of infection of these. In this case, the intensity becomes

$\Lambda_{t}=\lambda_{0}+\int_{0}^{t}Q_{s}\phi(s,t)dC_{s}.$ (7)

This case is considered in Lesage, (2020) where the author analyzes, simulates
and fits a classical Hawkes process for the total number of infected people in
France but to signify a change in public policy modifies the intensity to
mimic a lock-down as in China. As another example, in Escobar, (2020) where
the author discretizes the intensity to analyze and fit the epidemic model to
the gathered data of several countries, including Mexico. This model is very
popular because there are well established algorithms to simulate and to
calibrate the model to existing data.

  * •

A proportional-to-infections approach. In this case, the baseline intensity
will be proportional to the current number of infected people people. Here,
not only can over-dispersion be captured but the intensity, and thus the
likelihood, of observing a new infection will increase or decrease according
to current number of infected people. This can be achieved by setting
$\lambda_{0}=0$ and $\beta=1$ and the intensity becomes

$\Lambda_{t}=\lambda_{0}+i\log(\alpha)+\int_{0}^{t}Q_{s}\phi(s,t)dC_{s},$ (8)

  * •

A state dependent Hawkes baseline intensity. This is the more general case and
the one that will be considered in this paper. In this case, assumption A1 can
be used to focus on different aspects of the model depending on whether
$0<\beta<1$, $\beta=1$ or if $\beta>1$. Indeed, if $0<\beta<1$, then the
quantity $\rho\log(\beta)$ is negative and in this case we can think of a
situation where the recovered individuals suppress the epidemics. This is more
in line with the classical assumptions of an SIR model where by increasing the
number of recovered people the infectivity of the virus decays and new cases
are less likely. However, this behaviour hinges on a critical assumption:
recovered people cannot contribute more to the epidemics and that they become
isolated from it. This assumption is not completely true for the COVID-19
epidemics where many recovered people become reinfected and can contribute to
the spread of the disease. One should be careful here though, if the number of
recovered individuals become very large and $\beta<1$, the baseline intensity
$\lambda^{\infty}(i,\rho)$ might become negative and this is not permitted.
Several solutions could be considered at this point, the two more relevant
being to stop the process before $\lambda^{\infty}(i,\rho)<0$ or change the
form Assumption A1 to
$\lambda^{\infty}(i,\rho)=\bigg{(}\lambda_{0}+i\log(\alpha)+\rho\log(\beta)\bigg{)}_{+}$,
where $f(x)=(x)_{+}$ represents the positive part of $x$, being this latter
option a very interesting case left out for further research. If $\beta=1$,
then the number of recovered people does not affect the model and the process
$R_{t}$ becomes a homogenous Poisson process independent of $C_{t}$. Finally,
whenever $\beta>1$ the recoveries can also contribute to the infection of new
individuals (rather than suppressing it). As discussed above, in the current
pandemic many individuals can become reinfected and keep transmitting the
virus to to other people. This feature might comes in handy specifically at
the current situation with COVID-19 as universities and other workplaces treat
recovered individuals as non-infectious but indeed they might become
infectious again but at a different rate than that of susceptible or currently
infected individuals.

It is also important to determine the parameters of the model:

  * •

The baseline infective rate $\lambda_{0}$, which is associated to the
likelihood that someone new gets infected due to exogenous factors to the
model such as migrations, population dynamics, etc. This is not affected by
the amount of infected or recovered people.

  * •

The infective rate $\alpha>1$. This parameter (or more precisely, the
logarithm of it) measures how the likelihood at time $t$ of a new infection
will increase due to the number of infected people at that particular instant
regardless of the total number of people that has been infected by the virus.
The higher $\alpha$, the more likely a new infection will occur at the level
of current infections $I_{t}$.

  * •

The situational rate $\beta>1$. This parameter (or more precisely, the
logarithm of it) measures how the likelihood at time $t$ of a new infection
will increase or decrease due to the number of recovered people up to that
particular instant. As $\beta\to 0$ the likelihood that a new infection will
occur at the level of current of $R_{t}$ decreases while if $\beta\to\infty$
the likelihood will increase.

  * •

The mean recovery time $\mu$. This parameter is the expected recovery time
period from the infection. There are several studies about this quantity and
it depends on several factors, but according to Chowdhury et al., (2021) a
good approximation could be around 9 days.

  * •

The historical influence parameter $r$ (also called the reversion coefficient
$r$). This parameter will measure the influence of previous (historical)
infections in the arrival of a new infection and indicates the speed at which
the likelihood of a new infection decays to the current level of the baseline
intensity on the absence of a new infection. In fact, for out model, the
increase in the likelihood of a new infection due to a previous infection at
time $s<t$ is proportional to $e^{-r(t-s)}$. This means that more recent an
infection, the more it will contribute to the likelihood of a further
infection. This behavior is the key difference with respect to the classical
epidemiological models and what allows for the clustering and over-disperssion
observed in the data for the current pandemic.

  * •

The probability distribution of the “quarantine effect” $Q$. This random
variable specifies the level of “quarantine” each individual will have.
Indeed, our model specifies some stochastic dynamics where the probability of
a new infection at time $t$ is driven by the amount of people that has been
infected up to that time and how recent their infections have been. This
quarantine factor will basically determine the proportion at which infected
individual person will contribute to a new infection. From the modeling
perspective, this random variable can be thought as a measure of quarantining.
The lower the (random variable) $Q$, the lower the contribution of such
infected individual to a new infection. This random variable can be discrete
or continuous, but its support has to be over the positive numbers. In other
words, an infected individual must have a positive contribution to the general
likelihood of generating a new infection even if it is small That is, there
cannot be a “perfect quarantine”, which is consistent because people have to
go to the groceries or buy basic services and even interact with delivery
services by receiving goods at home. For this work we will impose the mild
restriction that the Moment Generating Function (MGF) of $Q$ exists on a
neighbourhood of 0.

######  Remark 1.

It is also known (see Chapter 3.3 in Laub et al., (2022)) than when the random
variable $Q$ is constant, and the intensity of the point process is the
classical Hawkes process in Equation 2 with $\phi(t)=Qe^{-rt}$, then the
branching ratio of the process would be given by

$n=\frac{Q}{r}$

where as above, $r$ is the reversion coefficient. In the SIR process where the
number total (historical) number of people that has been infected by the virus
is driven by a Hawkes process the branching ration $n$ has the interpretation
that when $0<n<1$, it becomes the ratio of the number of people that one
individual will infect relative to the entire population; that is, it is
related to the epidemiological basic reproduction number $\mathcal{R}_{0}$.

######  Remark 2.

Assumption (A1) specifies the form of the so-called baseline intensity. This
quantity will not depend on the past number of infections but solely on the
present number of infected and recovered individuals. In fact, this baseline
intensity remains constant between events, that is, between new infections or
recoveries, the probability of a new infection is exponentially distributed
with rate $\lambda^{\infty}(i,\rho)$ and thus the process dictating the
arrival of a new infection is equal in distribution to a homogeneous Poisson
process with rate $\lambda^{\infty}(i,\rho)$.

Also, notice that $\alpha$ and $\beta$ represent the factors that drives the
baseline intensity according to the state of the system. The more infected
people there is at the moment, the more likely a new infection will occur, and
the more recovered individuals there are in the present moment, the less
likely a new infection will occur.

Further, since we don’t want a “degenerate” Hawkes process, we will assume
that $Q$ does not have an atom at 0. That is,

$\mathbb{P}[Q\leq 0]=0$ (A2)

##  3 Derivation of the Moments of the Number of Infected individuals $I_{t}$

As mentioned in the previous section, one of the features that the epidemic
model presented her has is that every person is assumed to have a random level
of quarantine, whose law is given by $Q$. As an easy example, assume that
$supp(Q)=\\{0.25,0.75,1\\}$. In this case, $Q=0.25$ will imply a higher level
of quarantine (contributing less to new infections) while $Q=1$ would mean a
low level of quarantine (so that this person will increase the likelihood of a
new infection happening). Obviously, $Q=0.5$ would be an intermediate case.

To compare how this different levels of quarantine and other parameters affect
the model, it might be worth to look at the average behaviour of $I_{t}$ and
$R_{t}$. This is because depending on Var$[I_{t}]$, some comparisons might not
be depicted accurately by plotting some trajectories of the process.

Recall from Equation (6) that the total number of infected people up to time
$t$, denoted by $C_{t}$, is a determined by a counting process with stochastic
intensity given by

$\Lambda_{t}=\lambda^{\infty}(I_{t},R_{t})+\int_{0}^{t}Q_{s}\phi(s,t)dC_{s},$
(9)

where $\\{Q_{n}\\}_{n=1}^{\infty}$ is a sequence of i.i.d. random variables
independent among themselves and from every other process denoting the level
of quarantine that the $n-$th infected person will have, modulating the
probability that this individual will produce future infections.

We are interested in computing the generating function of the triplet
$(I_{t},R_{t},\lambda(t))$

$\mathbb{E}\left[z^{I_{t}}w^{R_{t}}e^{-s\lambda_{t}}\right]$

Here, we will assume that all the stochastic processes are Markovian, which
can be attained if

  * •

The recovery time is exponentially distributed with rate parameter $\mu$. That
is, if $\tau$ is the recovery time for an infected individual, then

$\mathbb{P}[\tau\leq t]=1-e^{-\mu t}.$ (A3)

This implies in particular that if $I_{t}=k$, then

$\mathbb{P}[R_{t+\Delta
t}-R_{t}=1]=\mathbb{P}[\min\\{\tau_{1},\ldots,\tau_{k}\\}\leq\Delta
t]=1-e^{k\mu\Delta t}=k\mu\Delta t+o(\Delta t)$ (10)

  * •

The self-exciting kernel
$\phi(\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ )$ is
an exponential function. That is,

$\phi(s,t)=e^{-r(t-s)}$ (A4)

These assumptions can be relaxed, but the analytical tractability will be lost
and different techniques would have to be employed. The analysis presented
here is inspired by the one presented in Koops et al., (2018) but here is
generalized to allow a state dependent baseline intensity.

In order to compute the joint distribution of $(I_{t},R_{t},\Lambda(t))$, we
need to solve a system of differential equations presented next.

######  Theorem 3.

Let
$F(t,i,\rho,\lambda)=\mathbb{P}[I_{t}=i,R_{t}=\rho,\Lambda_{t}\leq\lambda]$
and $f(t,i,\rho,\lambda)=\dfrac{\partial
F}{\partial\lambda}(t,i,\rho,\lambda)$. Then, under assumptions (A3)-(A4),
$f(t,i,\rho,\lambda)$ satisfy the Partial difference differential equation
(PDDE):

$\displaystyle\dfrac{\partial}{\partial
t}f(t,i,\rho,\lambda)-\frac{\partial}{\partial\lambda}(r\lambda$
$\displaystyle
f(t,i,\rho,\lambda))+r\lambda^{\infty}(i,\rho)\frac{\partial}{\partial\lambda}f(t,i,\rho,\lambda)=$
(11)
$\displaystyle\int_{0}^{\lambda-\Delta\lambda^{\infty,1}_{i,\rho}}yf(t,i-1,\rho,y)\frac{\partial}{\partial\lambda}\mathbb{P}[Q\leq\lambda-\Delta\lambda^{\infty,1}_{i,\rho}-y]dy+$
$\displaystyle(i+1)\mu
f(t,i+1,\rho-1,\lambda-\Delta\lambda^{\infty,2}_{i,\rho})-(i\mu+\lambda)f(t,i,\rho,\lambda)$

where

$\Delta\lambda^{\infty,1}_{i,\rho}=\lambda^{\infty}(i,\rho)-\lambda^{\infty}(i-1,\rho)\qquad\qquad\text{and}\qquad\qquad\Delta\lambda^{\infty,2}_{i,\rho}=\lambda^{\infty}(i,\rho)-\lambda^{\infty}(i+1,\rho-1)$

The next objective is to compute the generating function of the triple
$(I_{t},R_{t},\Lambda(t))$, which can then be used to compute the moments and
other quantities of interest. However, since $I_{t}$ and $R_{t}$ are discrete
and $\Lambda(t)$ is continuous, we will have to compute the transformations of
those random variables separately.

######  Theorem 4.

Let
$\varphi(t,z,w,s)=\mathbb{E}\left[z^{I_{t}}w^{R_{t}}e^{-s\Lambda(t)}\right]$
for $|z|\leq 1$ and $|w|\leq 1$. Then, $\varphi:=\varphi(t,z,w,s)$ satisfies
the PDE

$\displaystyle\dfrac{\partial}{\partial
t}\varphi+\left[rs-1+\left(\frac{1}{\alpha}\right)^{s}M_{Q}(s)z\right]\frac{\partial}{\partial
s}\varphi$
$\displaystyle+\left[\left(\mu+rs\log\left(\alpha\right)\right)z-\mu\left(\frac{\alpha}{\beta}\right)^{s}w\right]\frac{\partial}{\partial
z}\varphi$ (12)
$\displaystyle+\bigg{[}\log(\beta)rsw\bigg{]}\frac{\partial}{\partial
w}\varphi=-\lambda_{0}rs\varphi$

with initial condition

$\varphi(0,z,w,s)=z^{i_{0}}e^{-s(\lambda_{0}+i_{0}\log(\alpha))}$

Since the resulting PDE is linear of first order, we can apply the method of
characteristics to simplify the problem to a system of 2 by 2 ODEs. However,
the solution can be quite messy and numerical methods are very likely to be
needed to solve such system which we just point out for the sake of
completeness.

######  Corollary 5.

Let $(z,w,s)$ be a fixed point in $[0,\infty)^{3}$ and let
$\Xi(t):=\varphi(t,z,w,s)=\mathbb{E}\left[z^{I_{t}}w^{R_{t}}e^{-s\Lambda(t)}\right]$.
Then, $\Xi(t)$ is the solution of the system of ODEs

$\displaystyle\left\\{\begin{array}[]{ccll}\dfrac{d\vartheta}{dt}&=&-1+r\vartheta(t)+\left(\frac{1}{\alpha}\right)^{\vartheta(t)}M_{Q}(\vartheta(t))\bigg{[}&z\exp\left(\mu
t+r\log\left(\alpha\right)\displaystyle\int_{0}^{t}\vartheta(u)du\right)\\\\[11.38092pt]
&&&-w\mu\displaystyle\int_{0}^{t}\left(\frac{\alpha}{\beta}\right)^{\vartheta(u)}\exp\left(-\mu
u+r\log(\beta)\displaystyle\int_{0}^{u}\vartheta(y)dy\right)du\bigg{]}\\\\[11.38092pt]
\dfrac{d\Xi}{dt}&=&-\lambda_{0}r\vartheta(t)\hfil\lx<EMAIL_ADDRESS>\vartheta(0)&=&s&\\\\[14.22636pt]
\Xi(0)&=&z^{i_{0}}\exp\left(-\vartheta(\lambda_{0}+i_{0}\log(\alpha))\right)&\end{array}\right.$
(13)

Our objective is to be able to compute the moments of the random variables
$I_{t}$ and $R_{t}$, however the system of ODEs (13) might not prove that
useful for this task, therefore we turn our attention to the PDE (12) again.
In order to obtain the joint $(i,j,k)-$th moment of $I_{t}$, $R_{t}$ and
$\Lambda(t)$ we take the $i-$th derivative of the PDE (12) with respect to
$z$, the $j-$th derivative with respect to $w$ and the $k-$th derivative with
respect to $s$ plug in the values $z=w=1$ and $s=0$. However, for simplicity,
we will proceed to give an explicit formula for the first and second moments
of the processes $I_{t},R_{t},\Lambda(t)$ since this are the ones that might
be used the most.

Unfortunately, the moments of $I_{t}$ and $R_{t}$ are not independent of those
of $\Lambda(t)$, therefore to compute the first moment, we will need to solve
a system of 3 linear differential equations. Each equation will be obtained by
differentiating the PDE (12) with respect to one parameter $z,w$ or $s$ and
plugging the values mentioned before. We show the procedure in the next two
lemmas.

######  Lemma 6.

For any time $t>0$, let
$\overline{\mathbf{X}}_{t}=\big{[}\mathbb{E}[\Lambda(t)],\mathbb{E}[I_{t}],\mathbb{E}[R_{t}]\big{]}^{T}$.
Then, $\overline{\mathbf{X}}_{t}$ is the solution to the system of
differential equations

$\displaystyle\dfrac{d}{dt}\overline{\mathbf{X}}_{t}$
$\displaystyle=\mathbf{A}_{1}\overline{\mathbf{X}}_{t}+\mathbf{b}_{1}$ (14)
$\displaystyle\overline{\mathbf{X}}_{0}$
$\displaystyle=\left[\begin{array}[]{c}\lambda_{0}+i_{0}\log(\alpha)\\\
i_{0}\\\ 0\end{array}\right]$ (18)

where

$\mathbf{A}_{1}=\left[\begin{array}[]{ccc}\log\left(\alpha\right)+\mathbb{E}[Q]-r&(r-\mu)\log\left(\alpha\right)+\mu\log\left(\beta\right)&r\log(\beta)\\\
1&-\mu&0\\\
0&\mu&0\end{array}\right]\qquad\text{and}\qquad\mathbf{b}_{1}=\left[\begin{array}[]{c}\lambda_{0}r\\\
0\\\ 0\end{array}\right]$

That is,

$\overline{\mathbf{X}}_{t}=\overline{\mathbf{X}}_{0}e^{\mathbf{A}_{1}t}+\int_{0}^{t}e^{\mathbf{A}_{1}(t-s)}\mathbf{b}_{1}ds$
(19)

We can use the same approach as in Lemma 6 to compute the second moments.
Besides, we will need to consider all 6 possible double products of the random
variables $I_{t},R_{t}$ and $\Lambda(t)$.

For the following Lemma, it will be useful to simplify the notation so that
the matrix used to compute second moments is displayed nicely. Indeed, let

$\displaystyle C_{1}$ $\displaystyle:=\log\left(\alpha\right)+\mathbb{E}[Q]-r$
$\displaystyle C_{2}$
$\displaystyle:=(r-\mu)\log\left(\alpha\right)+\mu\log(\beta)$ $\displaystyle
C_{3}$ $\displaystyle:=r\log(\beta)$ $\displaystyle C_{4}$
$\displaystyle:=r\log(\alpha)$

An important remark is that without too much work we can transition from the
models presented in Section 2 by just modifying the parameters of the original
model. For example, if a pure-Hawkes (pH) model with intensity given by
Equation (7) wants to be considered; then by setting $\alpha=\beta=1$ in the
previous result it is obtained that the following:

$\displaystyle\dfrac{d}{dt}\overline{\mathbf{X}}_{t}^{pH}$
$\displaystyle=\mathbf{A}_{1}^{pH}\overline{\mathbf{X}}_{t}+\mathbf{b}_{1}^{pH}$
(20) $\displaystyle\overline{\mathbf{X}}_{0}^{pH}$
$\displaystyle=\left[\begin{array}[]{c}\lambda_{0}\\\ i_{0}\\\
0\end{array}\right]$ (24)

where

$\mathbf{A}_{1}^{pH}=\left[\begin{array}[]{ccc}\mathbb{E}[Q]-r&0&0\\\
1&-\mu&0\\\
0&\mu&0\end{array}\right]\qquad\text{and}\qquad\mathbf{b}_{1}=\left[\begin{array}[]{c}\lambda_{0}r\\\
0\\\ 0\end{array}\right]$

which agrees with the result provided in Koops et al., (2018). Further, if a
model whose infection rate increases or decreases proportional to the number
of infected people in its baseline intensity (pI) as given by Equation (8) is
sought, then it should be set $\beta=1$ and the resulting system of ODEs is

$\displaystyle\dfrac{d}{dt}\overline{\mathbf{X}}_{t}^{pI}$
$\displaystyle=\mathbf{A}_{1}\overline{\mathbf{X}}_{t}^{pI}+\mathbf{b}_{1}^{pI}$
(25) $\displaystyle\overline{\mathbf{X}}_{0}^{pI}$
$\displaystyle=\left[\begin{array}[]{c}\lambda_{0}+i_{0}\log(\alpha)\\\
i_{0}\\\ 0\end{array}\right]$ (29)

where

$\mathbf{A}_{1}^{pI}=\left[\begin{array}[]{ccc}\log\left(\alpha\right)+\mathbb{E}[Q]-r&(r-\mu)\log\left(\alpha\right)&0\\\
1&-\mu&0\\\
0&\mu&0\end{array}\right]\qquad\text{and}\qquad\mathbf{b}_{1}=\left[\begin{array}[]{c}\lambda_{0}r\\\
0\\\ 0\end{array}\right]$

######  Remark 7.

Notice that in the two particular cases above given by systems of ODEs in
Equations (20)-(29) since $\beta=1$, then the intensity becomes independent of
$R_{t}$ and thus the entire system becomes independent of $R_{t}$ which is
reflected by the third column of the matrices $\mathbf{A}_{1}^{pH}$ and
$\mathbf{A}_{1}^{pI}$ being 0.

Finally, a third important differentiation from the base case $\alpha,\beta>1$
is when the number of recovered individuals actually _decreases_ the
likelihood of a new infection, such as in a classical SIR model and which is
also described in Section 2. In this case, and as mentioned earlier, provided
that the process is stopped whenever the intensity becomes negative, the
corresponding system of differential equations for the first moments can be
obtained from Lemma 6 by setting $\beta=1/\tilde{\beta}$ with
$\tilde{\beta}>1$ so that, in this case, $\beta<1$ and the baseline intensity
given by Assumption (A1) becomes

$\lambda^{\infty}(i,\rho)=\lambda_{0}+i\log(\alpha)+\rho\log(\beta)=\lambda_{0}+i\log(\alpha)-\rho\log(\tilde{\beta}).$
(30)

In this case, the corresponding system of differential equations governing the
first moments of this system is given by

$\displaystyle\dfrac{d}{dt}\overline{\mathbf{X}}_{t}^{SIR}$
$\displaystyle=\mathbf{A}_{1}\overline{\mathbf{X}}_{t}^{SIR}+\mathbf{b}_{1}^{SIR}$
(31) $\displaystyle\overline{\mathbf{X}}_{0}^{SIR}$
$\displaystyle=\left[\begin{array}[]{c}\lambda_{0}+i_{0}\log(\alpha)\\\
i_{0}\\\ 0\end{array}\right]$ (35)

where

$\mathbf{A}_{1}^{SIR}=\left[\begin{array}[]{ccc}\log\left(\alpha\right)+\mathbb{E}[Q]-r&(r-\mu)\log\left(\alpha\right)-\mu\log(\tilde{\beta})&-r\log(\tilde{\beta})\\\
1&-\mu&0\\\
0&\mu&0\end{array}\right]\qquad\text{and}\qquad\mathbf{b}_{1}=\left[\begin{array}[]{c}\lambda_{0}r\\\
0\\\ 0\end{array}\right]$

Next, a characterization of the second moments is provided. The techniques and
methods to obtain them follow from the ones used in Lemma 6 but to solve this
system it is necessary to obtain the solution of the system (14)-(35), since
its solution its dependent (as expected) on the first moments of
$(\Lambda(t),I_{t},R_{t})$.

######  Lemma 8.

For any time $t>0$, let
$\overline{\mathbf{Y}}_{t}=\Big{[}\mathbb{E}[\Lambda^{2}(t),\mathbb{E}[I_{t}(I_{t}-1)],\mathbb{E}[R_{t}(R_{t}-1)],\mathbb{E}[\Lambda(t)I_{t}],\mathbb{E}[\Lambda(t)R_{t}],\mathbb{E}[I_{t}R_{t}]\Big{]}^{T}$.
Then, $\overline{\mathbf{Y}}_{t}$ is the solution to the system of equations

$\displaystyle\dfrac{d}{dt}\overline{\mathbf{Y}}_{t}$
$\displaystyle=\mathbf{A}_{2}\overline{\mathbf{Y}}_{t}+\mathbf{b}_{2}$ (36)
$\displaystyle\overline{\mathbf{Y}}_{0}$
$\displaystyle=\left[(\lambda_{0}+i_{0}\log(\alpha)))^{2},i_{0}^{2},0,i_{0}(\lambda_{0}+i_{0}\log(\alpha)),0,0\right]^{T}$
(37)

where,

$\mathbf{A}_{2}=\left[\begin{array}[]{cccccc}2C_{1}&0&0&2C_{2}&2C_{3}&0\\\
0&-2\mu&0&2&0&0\\\ 0&0&0&0&0&2\mu\\\ 1&C_{2}&0&C_{1}&0&C_{3}\\\
0&0&C_{3}&-\mu&C_{1}&C_{2}\\\ 0&\mu&0&0&1&-\mu\\\ \end{array}\right]$

and

$\mathbf{b}_{2}=\left[\begin{array}[]{c}\left[\log^{2}\left(\alpha\right)+2\log(\alpha)+\mathbb{E}[Q^{2}]+2\lambda_{0}r\right]\mathbb{E}[\Lambda(t)]+\mu\log^{2}\left(\frac{\alpha}{\beta}\right)\mathbb{E}[I_{t}]\\\
0\\\ 0\\\
(C_{1}+r)\mathbb{E}[\Lambda(t)]+r\lambda_{0}\mathbb{E}[I_{t}]+C_{4}\mathbb{E}[R_{t}]\\\
(r\lambda_{0}+C_{3})\mathbb{E}[R_{t}]+(C_{2}-C_{4})\mathbb{E}[I_{t}]\\\
0\end{array}\right]$

That is,

$\overline{\mathbf{Y}}_{t}=\overline{\mathbf{Y}}_{0}e^{\mathbf{A}_{2}t}+\int_{0}^{t}e^{\mathbf{A}_{2}(t-s)}\mathbf{b}_{2}ds$
(38)

##  4 Some considerations regarding the Estimation of the State Dependent
Hawkes Process

As with many many point processes, there are mainly three types of estimation
procedures: MLE, (generalized)MoM and LSE-type methods. Below, a discussion on
each of these three process, their advantages and disadvantages is presented.
It is important to emphasize that there is not a definite answer into which
method provides an advantage when doing inference over the parameters of the
model and this intricate topic is left out for future research.

  * •

Maximum Likelihood Estimation (MLE). This method is one of the most used and
regarded within the academic literature. This is due in part because there
there are many theoretical results that guarantee that the MLE methods will
yield an optimal solution. Unfortunately, for many Markovian point processes
the evaluation of the log-likelihood function is of order $O(N_{t}^{2})$.
Several algorithms such as EM are used to improve this (see Lewis et al.,
(2012)), but in general, the biggest problem of MLE is that the likelihood
curve is very flat, with many local maxima and without a clear way to decide
which is the global maximum. An important step in improving this is discussed
in Lewis et al., (2012) and Veen and Schoenberg, (2008). Many of these methods
apply to the exponential kernel form used in this paper (see Assumption (A4)),
but also have a constant baseline intensity. It is still an ongoing research
topic of this author to generalize such methods to a state dependent
intensity.

  * •

(Generalized) Method of Moments (gMoM). These methods are derived under the
assumption that the process is in its limiting stationary state, and as any
Method of Moments, the idea is to form a system of equations of $n\times n$,
where $n$ is the number of parameters to be determined and $n$ linear
independent equations relating the moments are used. However, for Hawkes
processes, usually there are not enough linear independent conditions on the
moments and thus the Autocovariance function is introduced to generate other
equations. A prime example is given in the work by Foschi et al., (2020).
Further, these methods are inefficient when the number of observations is
“small” and they usually does not work properly on higher dimensions. However,
it is important to mention that in our case, these methods might not be the
best since we are trying to model the early stages of the epidemic where the
processes are far away from the stationary state.

  * •

Least Square Estimation (LSE) type methods. These methods have not been
explored until recently, mainly because the order of these methods is similar
to the ones of MLE. However, the recent work by Cartea et al., (2021) show
that unlike MLE methods, LSE methods can possess certain algebraic properties
that help with the stochastic approximation of the kernels to then maximize or
minimize the LSE functional.

The nature of the process we are dealing with requires special care in the
estimation of parameters and as such, this paper will not try to just follow
the MLE method presented in Daley Vere Jones, Section , but rather this
sensitive topic is left as an object of further research. Specially if
Assumption (A4) is relaxed and we allow for a more general kernel and also a
state-dependent self-exciting kernel.

##  5 Empirical and Numerical Examples

The objective of this section is to provide numerical evidence of the
behaviour that the model has under the different scenarios proposed as well as
the implications of increasing or decreasing the level of quarantine provided
by the random variable $Q$ and other change in the parameters.

For all the experiments there will be four levels of Quarantine: high
($Q_{H}$), medium-high ($Q_{M_{H}}$), medium-low ($Q_{M_{L}}$) and low
($Q_{L}$). To completely specify the levels of quarantine provided by these
random variables, their distribution is presented next.

$Q_{H}=\begin{cases}0.05&\text{w.p. }0.92\\\ 0.45&\text{w.p. }0.03\\\
0.60&\text{w.p. }0.02\\\ 0.95&\text{w.p. }0.03\end{cases},\quad
Q_{M_{H}}=\begin{cases}0.05&\text{w.p. }0.11\\\ 0.45&\text{w.p. }0.18\\\
0.60&\text{w.p. }0.46\\\ 0.95&\text{w.p. }0.25\end{cases},\quad
Q_{M_{L}}=\begin{cases}0.05&\text{w.p. }0.12\\\ 0.45&\text{w.p. }0.46\\\
0.60&\text{w.p. }0.32\\\ 0.95&\text{w.p. }0.10\end{cases},\quad
Q_{L}=\begin{cases}0.05&\text{w.p. }0.02\\\ 0.45&\text{w.p. }0.02\\\
0.60&\text{w.p. }0.02\\\ 0.95&\text{w.p. }0.91\end{cases}$

As it can be seen, the random variables take the same values in all cases but
their probabilities change. They are specified so that
$\mathbb{E}[Q_{L}]=0.9$, $\mathbb{E}[Q_{M_{L}}]=0.6$,
$\mathbb{E}[Q_{M_{H}}]=0.5$ and $\mathbb{E}[Q_{H}]=0.1$.

Next, we provide several figures with various simulations of the different
cases mentioned in Section 3. In particular, it is important to notice how the
difference on the parameters affect the speed at which the number of infected
individuals grow. Also, it is important to remark that although the simulation
algorithm 2 is a variation of the thinning algorithm by Ogata described in
Appendix B, to the author’s knowledge, an explicit algorithm to simulate a
Hawkes process with a state dependent intensity is not readily available.
Thus, as part of this work, a detailed pseudo-code for simulating these kind
of processes is provided in the Appendix B. Furthermore, due to the immense
amount of simulations and computations required to simulate these processes
-at the end the thinning algorithm is a variant of an acceptance-rejection
method and as such, many simulations are rejected- the process cannot be
easily simulated for large time intervals.

To understand the behaviour of the process counting the number of infected
people some different parameters of the model are changed, in Figures 4-7
below, the number of infected people is simulated for the small time interval
$[0,4]$ together with the corresponding solution of the differential equation
under different scenarios, where either $\alpha$, $\beta$ or the Quarantine
distribution changes. This exercise has two purposes: verify that the
simulation algorithm and the differential equations yield similar results
(performing a cross-verification) and showing how the change in different
parameters yield logical conclusions as well as exploring how such changes
affect the speed at which the number of infected people grows.

To understand the behaviour of the process counting the number of infected
people under a general purely Hawkes process (i.e. we set $\alpha=\beta=1$),
in Figure 4, $I_{t}$ is plotted under 12 different scenarios and each scenario
under the 4 different Quarantine scenarios. In all the scenarios, the
parameters were set to be the same except for the decaying parameter $r$ (see
Equation (A4)) and the constant baseline intensity $\lambda_{0}$. As expected,
the number of infected people decreases as $r$ gets larger, meaning that an
infectious person will contribute to a new infection significantly only during
a short time span.

Two interesting cases arrive in the general state dependent Hawkes baseline
intensity case. First, as can be seen in Figure 5, when $\alpha>\beta$ (in
this case $\alpha=1.4$ while $\beta=1.2$, we see that there are more
infections as compared to the purely Hawkes case reflected in Figure 4.
However, when $\beta>\alpha>1$, we see that the infections also grow as
compared within the baseline case of the purely Hawkes model but they are
actually even higher than in the previous case where $\alpha>\beta$. The
reason for this is that the number of infected people becomes, as time
progresses, comparative smaller to the number of people that have recovered
and continue propagating the virus. This is a classical behaviour of viruses
that do not create immunity as is the case of COVID-19. Finally, in Figure 7
we plot the case where $\alpha>1$ but $\beta<1$. This is to illustrate how a
$\beta$ less than 1 can mitigate greatly the size of the epidemics. Indeed, as
a good comparison, note that, with all the parameters save $\alpha$ and
$\beta$ kept constant, the maximum range up to time 4 for the number of
infected individuals in the Pure Hawkes case is of roughly 400, while on the
state-dependent Hawkes cases with $\alpha>\beta>1$ is of 600; when
$\beta>\alpha>1$ is of 1200 while on the case where $\beta<1<\alpha$ it is of
200.

Finally, in Figure 8, the long-run behavior of the state-dependent Hawkes
process with $\beta<1$ is displayed. However, as mentioned above the
simulation of the process for larger time windows is not feasible and as such
only the solution of the differential equation is provided. The purpose of
this figure is to show that under this case, we can observe the typical
behaviour of an SIR model. Since this was computed in the long run, we only
plot $t$ vs $\mathbb{E}[I_{t}]$.

Figure 4: Figure showing 12 simulations and their behaviour for different
parameters under a purely Hawkes approach. Plots in the same row share the
same value of $r$ while plots in the same column share the same value of
$\lambda_{0}$. The rest of the parameters, as well as the scale, where kept
the same for comparison purposes. The Quarantine level is also represented,
with High Quarantine ($Q_{H}$) in purple, Medium-High Quarantine ($Q_{M_{H}}$)
in blue, Medium-Low Quarantine ($Q_{M_{L}}$) in green and Low Quarantine
($Q_{L}$) in red. Furthermore, the mean (deterministic) behaviour of the
system is plotted with the same color as the simulated paths. Figure 5: Figure
showing 12 simulations and their behaviour for different parameters under a
general state dependent Hawkes intensity with $\alpha=1.4$ and $\beta=1.2$.
Plots in the same row share the same value of $r$ while plots in the same
column share the same value of $\lambda_{0}$. The rest of the parameters, as
well as the scale, where kept the same for comparison purposes. Furthermore,
the Quarantine level is also represented, with High Quarantine ($Q_{H}$) in
purple, Medium-High Quarantine ($Q_{M_{H}}$) in blue, Medium-Low Quarantine
($Q_{M_{L}}$) in green and Low Quarantine ($Q_{L}$) in red. Figure 6: Figure
showing 12 simulations and their behaviour for different parameters under a
general state dependent Hawkes intensity with $\alpha=1.2$ and $\beta=1.4$.
Plots in the same row share the same value of $r$ while plots in the same
column share the same value of $\lambda_{0}$. The rest of the parameters, as
well as the scale, where kept the same for comparison purposes. Furthermore,
the Quarantine level is also represented, with High Quarantine ($Q_{H}$) in
purple, Medium-High Quarantine ($Q_{M_{H}}$) in blue, Medium-Low Quarantine
($Q_{M_{L}}$) in green and Low Quarantine ($Q_{L}$) in red. Figure 7: Figure
showing 12 simulations and their behaviour for different parameters under a
general state dependent Hawkes intensity with $\alpha=1.2$ and $\beta=0.8$.
Plots in the same row share the same value of $r$ while plots in the same
column share the same value of $\lambda_{0}$. The rest of the parameters, as
well as the scale, where kept the same for comparison purposes. Furthermore,
the Quarantine level is also represented, with High Quarantine ($Q_{H}$) in
purple, Medium-High Quarantine ($Q_{M_{H}}$) in blue, Medium-Low Quarantine
($Q_{M_{L}}$) in green and Low Quarantine ($Q_{L}$) in red. Figure 8: Figure
showing the long term behaviour of the mean for the number of infected people
under a general state dependent Hawkes intensity.

##  6 Conclusions and further Research Direction

A phenomena that has been observed in the current COVID-19 epidemics is the
overdispersion of process describing the number of contagions $C_{t}$. That
is, that for different time windows, the variance of the number of new
contagions within that time interval is great larger than their average. I.e.,
that for any real number $\tau>0$,

$\text{Var}[C_{t+\tau}-C_{t}]\gg\mathbb{E}[C_{t+\tau}-C_{t}]$

This work tries to construct a model that takes into account the
overdispersion observed in epidemics with a high control reproduction number
that the classic stochastic-SIR models fail to capture due to their Poissonian
nature (where the expectation and variance are roughly the same for any time
interval).

One possible solution to model the overdispersion is through the usage of
regular point processes, in particular self-exciting point processes. Not only
can they capture this phenomena but the rationale behind them makes sense.
Every new infection will increase the likelihood of a new infection occurring.
However, one of the main difficulties, and a very active are of research
currently, is the efficient estimation of its parameters. While it is part of
the ultimate goal of this work to research efficient estimation methods, it is
not the primary intent of this paper. The basic estimation method is in terms
of the likelihood function, but there is no clear way to determine whether the
method found a local maxima or a global one. Furthermore, the likelihood curve
is very flat and although theoretical convergence is guaranteed, in practice
it is hard to achieve.

To continue this research, there are two interesting directions to pursue.
First to explore a regime-switching model where the parameters of the model
change or switch randomly to different scenarios considering possible social,
economical and political factors. Particularly of interest is the level of
quarantining for the general population, where an stochastic optimization
problem can be devised for minimizing the economic impact of a disease like
COVID-19. Another direction of research is creating a full compartmental model
where the transitions to different compartments are guided by point processes.
This approach to classical stochastic compartmental models is of interest due
to the flexibility and the statistical properties of point processes.

##  Appendix A Proofs of Section 3

###### Proof of Theorem 3.

Let $\Delta t>0$ denote an infinitesimal time step. Due to the dynamics of the
SIR model, from time $t$ to time $t+\Delta t$, and since we assume that the
point processes defining the arrivals of events is simple (see Equation
1),only one of following four (disjoint) events can happen:

  * $\blacksquare$

Event of type (I): There is a new infection but no recoveries. That is,

$C_{t+\Delta t}-C_{t}=1\qquad\qquad\text{and}\qquad\qquad R_{t+\Delta
t}-R_{t}=0$

  * $\blacksquare$

Event of type (II): There is no new infection but a recovery. Then,

$C_{t+\Delta t}-C_{t}=0\qquad\qquad\text{and}\qquad\qquad R_{t+\Delta
t}-R_{t}=1$

  * $\blacksquare$

Event of type (III): There is a new infection and a recovery. Thus,

$C_{t+\Delta t}-C_{t}=1\qquad\qquad\text{and}\qquad\qquad R_{t+\Delta
t}-R_{t}=1$

  * $\blacksquare$

Event of type (IV): There are no new infections neither recoveries. Therefore,

$C_{t+\Delta t}-C_{t}=0\qquad\qquad\text{and}\qquad\qquad R_{t+\Delta
t}-R_{t}=0$

Recall that $\\{t_{i}\\}_{i=1}^{\infty}$ are the jump times of the counting
process $(C_{t})_{t\geq 0}$. Then, define the function

$J(t)=\int_{0}^{t}Q_{s}e^{-r(t-s)}dC_{s}=\sum\limits_{t_{i}\leq
t}Q_{i}e^{-r(t-t_{i})}.$

Fix the terminal state at time $t+\Delta t$, and computing

$\displaystyle F\Big{(}t+\Delta t,$ $\displaystyle
i,\rho,\lambda-r(\lambda-\lambda^{\infty}(i,\rho))\Delta
t\Big{)}=\mathbb{P}\Big{[}I_{t+\Delta t}=i,R_{t+\Delta
t}=\rho,\Lambda(t+\Delta
t)\leq\lambda-r(\lambda-\lambda^{\infty}(i,\rho))\Delta t\Big{]}$
$\displaystyle=\mathbb{P}\left[I_{t+\Delta t}=i,R_{t+\Delta
t}=\rho,\lambda^{\infty}(I_{t+\Delta t},R_{t+\Delta t})+\sum\limits_{t_{i}\leq
t+\Delta t}Q_{i}e^{-r(t+\Delta
t-t_{i})}\leq\lambda-r(\lambda-\lambda^{\infty}(i,\rho))\Delta t\right]$
$\displaystyle=\mathbb{P}\Big{[}I_{t+\Delta t}=i,R_{t+\Delta
t}=\rho,\lambda^{\infty}(i,\rho)+J(t+\Delta t)\leq\lambda(1-r\Delta t)+r\Delta
t\lambda^{\infty}(i,\rho)\Big{]}$ $\displaystyle=\mathbb{P}\Big{[}I_{t+\Delta
t}=i,R_{t+\Delta t}=\rho,(1-r\Delta t)\lambda^{\infty}(i,\rho)+J(t+\Delta
t)\leq\lambda(1-r\Delta t)\Big{]}$ $\displaystyle=\mathbb{P}\Big{[}I_{t+\Delta
t}=i,R_{t+\Delta t}=\rho,(1-r\Delta
t)\big{[}\lambda^{\infty}(I_{t},R_{t})+J(t)\big{]}+(1-r\Delta
t)\big{[}\lambda^{\infty}(i,\rho)-\lambda^{\infty}(I_{t},R_{t})\big{]}$
$\displaystyle\qquad\qquad\qquad+\big{(}J(t+\Delta t)-(1-r\Delta
t)J(t)\big{)}\leq\lambda(1-r\Delta t)\Big{]}$
$\displaystyle=\mathbb{P}\Big{[}I_{t+\Delta t}=i,R_{t+\Delta
t}=\rho,(1-r\Delta t)\Lambda(t)+(1-r\Delta
t)\big{[}\lambda^{\infty}(i,\rho)-\lambda^{\infty}(I_{t},R_{t})\big{]}$
$\displaystyle\qquad\qquad\qquad+\Big{(}\sum\limits_{t_{i}\leq t+\Delta
t}\left\\{(1-r\Delta t)Q_{i}e^{-r(t-t_{i})}+o(\Delta t)\right\\}-(1-r\Delta
t)J(t)\Big{)}\leq\lambda(1-r\Delta t)\Big{]}$
$\displaystyle=\mathbb{P}\left[I_{t+\Delta t}=i,R_{t+\Delta
t}=\rho,\Lambda(t)+\big{[}\lambda^{\infty}(i,\rho)-\lambda^{\infty}(I_{t},R_{t})\big{]}+\sum\limits_{t<t_{i}\leq
t+\Delta t}Q_{i}e^{-r(t-t_{i})}+o(\Delta t)\leq\lambda\right]$
$\displaystyle=\mathbb{P}\left[I_{t+\Delta t}=i,R_{t+\Delta
t}=\rho,\Lambda(t)\leq\lambda-\big{[}\lambda^{\infty}(i,\rho)-\lambda^{\infty}(I_{t},R_{t})\big{]}-\sum\limits_{t<t_{i}\leq
t+\Delta t}Q_{i}+o(\Delta t)\right]$

Since we are later going to take a difference quotient, divide by $\Delta t$
and take the limit as $dt\to 0$, for ease of notation, we will disregard all
the terms of lower order after displaying them once and we will introduce the
following notation

$\displaystyle\Delta\lambda^{\infty,1}_{i,\rho}$
$\displaystyle=\lambda^{\infty}(i,\rho)-\lambda^{\infty}(i-1,\rho)$
$\displaystyle\Delta\lambda^{\infty,2}_{i,\rho}$
$\displaystyle=\lambda^{\infty}(i,\rho)-\lambda^{\infty}(i+1,\rho-1)$
$\displaystyle\Delta\lambda^{\infty,3}_{i,\rho}$
$\displaystyle=\lambda^{\infty}(i,\rho)-\lambda^{\infty}(i,\rho-1)$
$\displaystyle\Delta\lambda^{\infty,4}_{i,\rho}$
$\displaystyle=\lambda^{\infty}(i,\rho)-\lambda^{\infty}(i,\rho)=0$

Define

$E(t+\Delta t;t,i,\rho,\lambda)=\left\\{I_{t+\Delta t}=i,R_{t+\Delta
t}=\rho,\Lambda(t)\leq\lambda-\big{[}\lambda^{\infty}(i,\rho)-\lambda^{\infty}(I_{t},R_{t})\big{]}-\sum\limits_{t<t_{i}\leq
t+\Delta t}Q_{i}\right\\}$

Then, by recalling that the point processes are simple and that
$\\{Q,Q_{i}\\}_{i=1}^{\infty}$ is a sequence of i.i.d. random variables, we
have that

$E(t+\Delta t;t,i,\rho,\lambda)=\left\\{I_{t+\Delta t}=i,R_{t+\Delta
t}=\rho,\Lambda(t)\leq\lambda-\big{[}\lambda^{\infty}(i,\rho)-\lambda^{\infty}(I_{t},R_{t})\big{]}-Q\right\\}$

and thus we can defined the events

$\displaystyle E_{I}(t,i,\rho,\lambda)$
$\displaystyle=\left\\{I_{t}=i-1,R_{t}=\rho,\Lambda(t)\leq\lambda-\Delta\lambda^{\infty}_{1}(i,\rho)-Q,\;C_{t+\Delta
t}-C_{t}=1,R_{t+\Delta t}-R_{t}=0\right\\}$ $\displaystyle
E_{II}(t,i,\rho,\lambda)$
$\displaystyle=\bigg{\\{}I_{t}=i+1,R_{t}=\rho-1,\Lambda(t)\leq\lambda-\Delta\lambda^{\infty}_{2}(i,\rho),\;C_{t+\Delta
t}-C_{t}=0,R_{t+\Delta t}-R_{t}=1\bigg{\\}}$ $\displaystyle
E_{III}(t,i,\rho,\lambda)$
$\displaystyle=\left\\{I_{t}=i,R_{t}=\rho-1,\Lambda(t)\leq\lambda-\Delta\lambda^{\infty}_{3}(i,\rho)-Q,\;C_{t+\Delta
t}-C_{t}=1,R_{t+\Delta t}-R_{t}=1\right\\}$ $\displaystyle
E_{IV}(t,i,\rho,\lambda)$
$\displaystyle=\bigg{\\{}I_{t}=i,R_{t}=\rho,\Lambda(t)\leq\lambda-\Delta\lambda^{\infty}_{4}(i,\rho),\;C_{t+\Delta
t}-C_{t}=0,R_{t+\Delta t}-R_{t}=0\big{]}\bigg{\\}}$

By considering the likelihood of the 4 types of events happening and the
dynamics of the SIR model, it is clear that

$\mathbb{P}\Big{[}E(t+\Delta
t;t,i,\rho,\lambda)\Big{]}=\mathbb{P}\Big{[}E_{I}(t,i,\rho,\lambda)\Big{]}+\mathbb{P}\Big{[}E_{II}(t,i,\rho,\lambda)\Big{]}+\mathbb{P}\Big{[}E_{III}(t,i,\rho,\lambda)\Big{]}+\mathbb{P}\Big{[}E_{IV}(t,i,\rho,\lambda)\Big{]}$
(39)

In the following, for exposition purposes, we will compute each of the fours
summands separately, and we will only consider operations up to the order
$\Delta t$. Further, Then, it follows that

$\displaystyle\mathbb{P}\Big{[}E_{I}(t,i,\rho,\lambda)\Big{]}$
$\displaystyle=\int_{0}^{\lambda-\Delta\lambda^{\infty,1}_{i,\rho}}f(t,i-1,\rho,y)\mathbb{P}[C_{t+\Delta
t}-C_{t}=1]\mathbb{P}[R_{t+\Delta
t}-R_{t}=0]\mathbb{P}\Big{[}Q\leq\lambda-\Delta\lambda^{\infty,1}_{i,\rho}-y\Big{]}dy$
$\displaystyle=\int_{0}^{\lambda-\Delta\lambda^{\infty,1}_{i,\rho}}f(t,i-1,\rho,y)\Big{[}y\Delta
t+o(\Delta t)\Big{]}e^{-(i-1)\mu\Delta
t}\mathbb{P}\Big{[}Q\leq\lambda-\Delta\lambda^{\infty,1}_{i,\rho}-y\Big{]}dy$
$\displaystyle=\int_{0}^{\lambda-\Delta\lambda^{\infty,1}_{i,\rho}}f(t,i-1,\rho,y)\Big{[}y\Delta
t+o(\Delta t)\Big{]}\Big{[}1-(i-1)\mu\Delta t+o(\Delta
t)\Big{]}\mathbb{P}\Big{[}Q\leq\lambda-\Delta\lambda^{\infty,1}_{i,\rho}-y\Big{]}dy$
$\displaystyle=\int_{0}^{\lambda-\Delta\lambda^{\infty,1}_{i,\rho}}y\Delta
tf(t,i-1,\rho,y)\mathbb{P}\Big{[}Q\leq\lambda-\Delta\lambda^{\infty,1}_{i,\rho}-y\Big{]}dy+o(\Delta
t),$ $\displaystyle\mathbb{P}\Big{[}E_{II}(t,i,\rho,\lambda)\Big{]}$
$\displaystyle=\int_{0}^{\lambda-\Delta\lambda^{\infty,2}_{i,\rho}}f(t,i+1,\rho-1,y)\mathbb{P}[C_{t+\Delta
t}-C_{t}=0]\mathbb{P}[R_{t+\Delta t}-R_{t}=1]dy$
$\displaystyle=\int_{0}^{\lambda-\Delta\lambda^{\infty,2}_{i,\rho}}f(t,i+1,\rho-1,y)\Big{[}1-y\Delta
t+o(\Delta t)\Big{]}\left[1-e^{(i+1)\mu\Delta t}\right]dy$
$\displaystyle=\int_{0}^{\lambda-\Delta\lambda^{\infty,2}_{i,\rho}}f(t,i+1,\rho-1,y)\Big{[}1-y\Delta
t+o(\Delta t)\Big{]}\left[1-(1-(i+1)\mu\Delta t+o(\Delta t))\right]dy$
$\displaystyle=\int_{0}^{\lambda-\Delta\lambda^{\infty,2}_{i,\rho}}f(t,i+1,\rho-1,y)\Big{[}1-y\Delta
t+o(\Delta t)\Big{]}(i+1)\mu\Delta tdy+o(\Delta t)$
$\displaystyle=\int_{0}^{\lambda-\Delta\lambda^{\infty,2}_{i,\rho}}(i+1)\mu\Delta
tf(t,i+1,\rho-1,y)dy+o(\Delta t),$

$\displaystyle\mathbb{P}\Big{[}E_{III}(t,i,\rho,\lambda)\Big{]}$
$\displaystyle=\int_{0}^{\lambda-\Delta\lambda^{\infty,3}_{i,\rho}}f(t,i,\rho-1,y)\mathbb{P}[C_{t+\Delta
t}-C_{t}=1]\mathbb{P}[R_{t+\Delta
t}-R_{t}=1]\mathbb{P}\Big{[}Q\leq\lambda-\Delta\lambda^{\infty,3}_{i,\rho}-y\Big{]}dy$
$\displaystyle=\int_{0}^{\lambda-\Delta\lambda^{\infty,3}_{i,\rho}}f(t,i,\rho-1,y)\Big{[}y\Delta
t+o(\Delta t)\Big{]}\Big{[}1-e^{-i\mu\Delta
t}\Big{]}\mathbb{P}\Big{[}Q\leq\lambda-\Delta\lambda^{\infty,3}_{i,\rho}-y\Big{]}dy$
$\displaystyle=\int_{0}^{\lambda-\Delta\lambda^{\infty,3}_{i,\rho}}f(t,i,\rho-1,y)\Big{[}y\Delta
t+o(\Delta t)\Big{]}\Big{[}i\mu\Delta t+o(\Delta
t)\Big{]}\mathbb{P}\Big{[}Q\leq\lambda-\Delta\lambda^{\infty,3}_{i,\rho}-y\Big{]}dy$
$\displaystyle=\int_{0}^{\lambda-\Delta\lambda^{\infty,3}_{i,\rho}}yi\mu(\Delta
t)^{2}f(t,i,\rho-1,y)\mathbb{P}\Big{[}Q\leq\lambda-\Delta\lambda^{\infty,3}_{i,\rho}-y\Big{]}dy+o(\Delta
t)$ $\displaystyle=o(\Delta t)$

and

$\displaystyle\mathbb{P}\Big{[}E_{IV}(t,i,\rho,\lambda)\Big{]}$
$\displaystyle=\int_{0}^{\lambda-\Delta\lambda^{\infty,4}_{i,\rho}}f(t,i,\rho,y)\mathbb{P}[C_{t+\Delta
t}-C_{t}=0]\mathbb{P}[R_{t+\Delta t}-R_{t}=0]dy$
$\displaystyle=\int_{0}^{\lambda}f(t,i,\rho,y)\Big{[}1-y\Delta t+o(\Delta
t)\Big{]}e^{i\mu\Delta t}dy$
$\displaystyle=\int_{0}^{\lambda}f(t,i+1,\rho-1,y)\Big{[}1-y\Delta t+o(\Delta
t)\Big{]}\Big{[}1-i\mu\Delta t+o(\Delta t)\Big{]}dy$
$\displaystyle=(1-i\mu\Delta
t)\int_{0}^{\lambda}f(t,i+1,\rho-1,y)\Big{[}1-y\Delta t+o(\Delta
t)\Big{]}dy+o(\Delta t)$ $\displaystyle=(1-i\mu\Delta
t)F(t,i,\rho,y)-\int_{0}^{\lambda}y\Delta tf(t,i,\rho,y)dy+o(\Delta t)$

Therefore, by using Equation (39) and recalling that $F\Big{(}t+\Delta
t,i,\rho,\lambda-r(\lambda-\lambda^{\infty}(i,\rho))\Delta
t\Big{)}=\mathbb{P}\Big{[}E(t+\Delta t;t,i,\rho,\lambda)\Big{]}$, we have that

$\displaystyle F\Big{(}t+\Delta
t,i,\rho,\lambda-r(\lambda-\lambda^{\infty}(i,\rho))\Delta t\Big{)}$
$\displaystyle=\int_{0}^{\lambda-\Delta\lambda^{\infty,1}_{i,\rho}}y\Delta
tf(t,i-1,\rho,y)\mathbb{P}\Big{[}Q\leq\lambda-\Delta\lambda^{\infty,1}_{i,\rho}-y\Big{]}dy$
$\displaystyle+\int_{0}^{\lambda-\Delta\lambda^{\infty,2}_{i,\rho}}(i+1)\mu\Delta
tf(t,i+1,\rho-1,y)dy$ $\displaystyle+(1-i\mu\Delta
t)F(t,i,\rho,y)-\int_{0}^{\lambda}y\Delta tf(t,i,\rho,y)dy+o(\Delta t)$

Dividing by $\Delta t$ and letting $\Delta t\to 0$,

$\displaystyle\frac{\partial}{\partial
t}F\Big{(}t,i,\rho,\lambda\Big{)}-r(\lambda-\lambda^{\infty}(i,\rho))\frac{\partial}{\partial\lambda}F(t,i,\rho,\lambda)$
$\displaystyle=\int_{0}^{\lambda-\Delta\lambda^{\infty,1}_{i,\rho}}yf(t,i-1,\rho,y)\mathbb{P}\Big{[}Q\leq\lambda-\Delta\lambda^{\infty,1}_{i,\rho}-y\Big{]}dy$
$\displaystyle+\int_{0}^{\lambda-\Delta\lambda^{\infty,2}_{i,\rho}}(i+1)\mu
f(t,i+1,\rho-1,y)dy$ $\displaystyle-i\mu
F(t,i,\rho,y)-\int_{0}^{\lambda}yf(t,i,\rho,y)dy$

Differentiating with respect to $\lambda$; applying Leibniz Integral Rule
where appropriate; and using Assumption (A2), we get the result. ∎

###### Proof of Theorem 4.

Define

$\displaystyle\psi(t,i,\rho,s)$
$\displaystyle:=\int_{0}^{\infty}e^{-s\lambda}f(t,i,r,\lambda)d\lambda$ (40)
$\displaystyle\varphi(t,z,w,s)$
$\displaystyle:=\mathbb{E}\left[z^{I_{t}}w^{R_{t}}e^{-s\Lambda(t)}\right].$
(41)

By definition,

$\displaystyle\varphi(t,z,w,s)$
$\displaystyle=\mathbb{E}\left[z^{I_{t}}w^{R_{t}}e^{-s\Lambda(t)}\right]$
$\displaystyle=\sum\limits_{i=0}^{\infty}\sum\limits_{\rho=0}^{\infty}\int_{0}^{\infty}z^{i}w^{\rho}e^{-s\lambda}f(t,i,\rho,\lambda)d\lambda$
$\displaystyle=\sum\limits_{i=0}^{\infty}\sum\limits_{\rho=0}^{\infty}z^{i}w^{\rho}\int_{0}^{\infty}e^{-s\lambda}f(t,i,\rho,\lambda)d\lambda$
$\displaystyle=\sum\limits_{i=0}^{\infty}\sum\limits_{\rho=0}^{\infty}z^{i}w^{\rho}\psi(t,i,\rho,s)$

Thus, the strategy is to first use Theorem 3 to find a PDDE equation that
characterizes the function
$\psi(\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ )$ and
then transform that PDDE as the corresponding one for
$\varphi(\,\begin{picture}(-1.0,1.0)(-1.0,-2.0)\circle*{2.0}\end{picture}\ )$.

Multiplying equation 11 by $e^{-s\lambda}$ and integrating over $(0,\infty)$
we get:

$\displaystyle\dfrac{\partial}{\partial
t}\int_{0}^{\infty}f(t,i,\rho,\lambda)e^{s\lambda}d\lambda$
$\displaystyle-\int_{0}^{\infty}\frac{\partial}{\partial\lambda}(r\lambda
f(t,i,\rho,\lambda))e^{-s\lambda}d\lambda+r\lambda^{\infty}(i,\rho)\int_{0}^{\infty}e^{-s\lambda}\frac{\partial}{\partial\lambda}f(t,i,\rho,\lambda)d\lambda=$
$\displaystyle\int_{0}^{\infty}\int_{0}^{\lambda-\Delta\lambda^{\infty}_{1}(i,\rho)}yf(t,i-1,\rho,y)d\mathbb{P}[Q\leq\lambda-\Delta\lambda^{\infty}_{1}(i,\rho)-y]dye^{-s\lambda}d\lambda+$
$\displaystyle(i+1)\mu\int_{0}^{\infty}f(t,i+1,\rho-1,\lambda-\Delta\lambda^{\infty}_{2}(i,\rho))e^{s\lambda}d\lambda-$
$\displaystyle\int_{0}^{\infty}(i\mu+\lambda)f(t,i,\rho,\lambda)e^{-s\lambda}d\lambda$

Using integration by parts and Fubini’s theorem, and relation (40)

$\displaystyle\dfrac{\partial}{\partial t}\psi(t,i,\rho,s)$
$\displaystyle-\int_{0}^{\infty}rs\lambda
f(t,i,\rho,\lambda)e^{-s\lambda}d\lambda+r\lambda^{\infty}(i,\rho)\int_{0}^{\infty}sf(t,i,\rho,\lambda)e^{-s\lambda}d\lambda=$
$\displaystyle\int_{0}^{\infty}\int_{y+\Delta\lambda^{\infty}_{1}(i,\rho)}^{\infty}yf(t,i-1,\rho,y)d\mathbb{P}[Q\leq\lambda-\Delta\lambda^{\infty}_{1}(i,\rho)-y]e^{-s\lambda}d\lambda
dy+$ $\displaystyle(i+1)\mu
e^{-s\Delta\lambda^{\infty}_{2}(i,\rho)}\int_{-\Delta\lambda^{\infty}_{2}(i,\rho)}^{\infty}f(t,i+1,\rho-1,\lambda)e^{s\lambda}d\lambda-$
$\displaystyle i\mu\psi(t,i,\rho,s)-$ $\displaystyle\int_{0}^{\infty}\lambda
f(t,i,\rho,\lambda)e^{-s\lambda}d\lambda$

Let $M_{Q}(s)=\mathbb{E}[e^{-sQ}]$ denote the moment generating function of
$Q$, which by assumption it exists on a neighbourhood of 0. Then, by using the
fact that $f(t,i,\rho,\lambda)=0$ if $\lambda<0$, definition (40) and Fubini’s
Theorem we can further simplify the above equation as

$\displaystyle\dfrac{\partial}{\partial t}\psi(t,i,\rho,s)$
$\displaystyle+rs\frac{\partial}{\partial
s}\psi(t,i,\rho,s)+rs\lambda^{\infty}(i,\rho)\psi(t,i,\rho,s)=$ $\displaystyle
e^{-s\Delta\lambda^{\infty}_{1}(i,\rho)}M_{Q}(s)\frac{\partial}{\partial
s}\psi(t,i-1,\rho,s)+(i+1)\mu
e^{-s\Delta\lambda^{\infty}_{2}(i,\rho)}\psi(t,i+1,\rho-1,s)$
$\displaystyle-i\mu\psi(t,i,\rho,s)+\frac{\partial}{\partial
s}\psi(t,i,\rho,s)$

Finally, rearranging terms we have

$\displaystyle\dfrac{\partial}{\partial t}\psi(t,i,\rho,s)$
$\displaystyle+(rs-1)\frac{\partial}{\partial
s}\psi(t,i,\rho,s)+rs\lambda^{\infty}(i,\rho)\psi(t,i,\rho,s)+e^{-s\Delta\lambda^{\infty}_{1}(i,\rho)}M_{Q}(s)\frac{\partial}{\partial
s}\psi(t,i-1,\rho,s)=$ (42) $\displaystyle(i+1)\mu
e^{-s\Delta\lambda^{\infty}_{2}(i,\rho)}\psi(t,i+1,\rho-1,s)-i\mu\psi(t,i,\rho,s)$

At this point, we use the particular form of the baseline intensity of the
Hawkes process $C_{t}$, $\lambda^{\infty}(i,\rho)$. By using Assumption (A1),
we have that

$\lambda^{\infty}(i,\rho)=\lambda_{0}+i\log(\alpha)+\rho\log(\beta)$

and substituting this into the PDDE (42) we have

$\displaystyle\dfrac{\partial}{\partial t}\psi(t,i,\rho,s)$
$\displaystyle+(rs-1)\frac{\partial}{\partial
s}\psi(t,i,\rho,s)+rs\lambda^{\infty}(i,\rho)\psi(t,i,\rho,s)+\left(\frac{1}{\alpha}\right)^{s}M_{Q}(s)\frac{\partial}{\partial
s}\psi(t,i-1,\rho,s)=$ (43)
$\displaystyle(i+1)\mu\left(\alpha\beta\right)^{s}\psi(t,i+1,\rho-1,s)-i\mu\psi(t,i,\rho,s)$

By doing some algebra and using the fact $\psi(t,i,\rho,s)=0$ if $i<0$ or
$\rho<0$, we have the following relations

$\displaystyle\sum\limits_{i=0}^{\infty}\sum\limits_{\rho=0}^{\infty}iz^{i}w^{\rho}\psi(t,i,\rho,s)$
$\displaystyle=z\frac{\partial}{\partial z}\varphi(t,z,w,s)$ (44)
$\displaystyle\sum\limits_{i=0}^{\infty}\sum\limits_{\rho=0}^{\infty}\rho
z^{i}w^{\rho}\psi(t,i,\rho,s)$ $\displaystyle=w\frac{\partial}{\partial
w}\varphi(t,z,w,s)$ (45)
$\displaystyle\sum\limits_{i=0}^{\infty}\sum\limits_{\rho=0}^{\infty}z^{i}w^{\rho}\psi(t,i-1,\rho,s)$
$\displaystyle=z\varphi(t,z,w,s)$ (46)
$\displaystyle\sum\limits_{i=0}^{\infty}\sum\limits_{\rho=0}^{\infty}(i+1)z^{i}w^{\rho}\psi(t,i+1,\rho-1,s)$
$\displaystyle=w\frac{\partial}{\partial z}\varphi(t,z,w,s)$ (47)

Multiplying both sides by $z^{i}w^{\rho}$; adding those terms as a series;
using Assumption (A1); and using relation (41), we obtain the result.

∎

###### Proof of Lemma 6.

Taking equation (12), substituting the value function $\varphi(t,z,w,s)$ and
computing some of the derivatives within,

$\displaystyle\dfrac{\partial}{\partial
t}\mathbb{E}\left[z^{I_{t}}w^{R_{t}}e^{-s\Lambda(t)}\right]$
$\displaystyle+\left[rs-1+\left(\frac{1}{\alpha}\right)^{s}M_{Q}(s)z\right]\mathbb{E}\left[-\Lambda(t)z^{I_{t}}w^{R_{t}}e^{-s\Lambda(t)}\right]$
(48)
$\displaystyle+\bigg{[}\left(\mu+rs\log\left(\alpha\right)\right)z-\mu\left(\alpha\beta\right)^{s}w\bigg{]}\mathbb{E}\left[I_{t}z^{I_{t}-1}w^{R_{t}}e^{-s\Lambda(t)}\right]$
$\displaystyle-\bigg{[}rs\log(\beta)w\bigg{]}\mathbb{E}\left[R_{t}z^{I_{t}}w^{R_{t}-1}e^{-s\Lambda(t)}\right]=-rs\lambda_{0}\mathbb{E}\left[z^{I_{t}}w^{R_{t}}e^{-s\Lambda(t)}\right]$

Taking partial derivatives with respect to $s$ from equation (48) and plugging
in $(z,w,s)=(1,1,0)$,

$-\dfrac{\partial}{\partial
t}\mathbb{E}\left[\Lambda(t)\right]-\bigg{[}r-\log\left(\alpha\right)-\mathbb{E}[Q]\bigg{]}\mathbb{E}\left[\Lambda(t)\right]+\bigg{[}r\log\left(\alpha\right)-\mu\log\left(\alpha\beta\right)\bigg{]}\mathbb{E}[I_{t}]-\bigg{[}r\log(\beta)\bigg{]}\mathbb{E}\left[R_{t}\right]=-r\lambda_{0}$
(49)

where we have used that $M_{Q}(0)=1$ and
$\left.\frac{d}{ds}M_{Q}(s)\right|_{s=0}=-\mathbb{E}[Q]$. Similarly, taking
partial derivatives with respect to $z$ from equation (48) and plugging in
$(z,w,s)=(1,1,0)$,

$\dfrac{\partial}{\partial
t}\mathbb{E}\left[I_{t}\right]-\mathbb{E}\left[\Lambda(t)\right]+\mu\mathbb{E}\left[I_{t}\right]=0$
(50)

Finally, taking partial derivatives with respect to $w$ from equation (48) and
plugging in $(z,w,s)=(1,1,0)$,

$\dfrac{\partial}{\partial
t}\mathbb{E}\left[R_{t}\right]-\mu\mathbb{E}\left[I_{t}\right]=0$ (51)

Arranging equations (49)-(51) in matrix form, we get the result. ∎

##  Appendix B Simulation Algorithms

When the baseline intensity $\lambda^{\infty}(i,\rho)$ is constant, several
simulation algorithms are available. However, it is rare to find simulation
algorithms for more general intensities.

However, by using the classical thinning algorithm by Ogata, we can simulate
our process, but some care and considerations need to be taken.

Before showing the main algorithm, we need to be able to evaluate the
intensity at any time $t$ given the history. In this case, the history is
provided in an array tIRarray which is a $3\times n$ array with Row 1 having
all the times at which $I_{t}$ or $R_{t}$ changed. Rows 2 and 3 have the value
of $I_{t}$ and $R_{t}$, respectively, at the times on row 1.

Result: Compute
$\Lambda(t)=\lambda^{\infty}(i,\rho)+\int_{0}^{t}Q_{s}e^{-r(t-s)}dN_{s}$, for
$t>0$;

1 Inputs: $C_{t}$,$Q_{t}$, tIRarray,$t$,$\ldots$;

2 Output: Intensity=function();

3 Initialization;

4

5$n\;\leftarrow\;$length$(C_{t})=||C_{t}||$;

6 $aux\;\leftarrow\;sum(C_{t}<t)=\int_{0}^{t}dC_{s}$;

7 Computation of the intensity function;

8 if  _$(n=0\ ||\ aux=0)$_ then

9 return $\lambda^{\infty}(i_{0},\rho_{0})$;

10

11else

12 $Qs\;\leftarrow\;Q_{t}[1:aux]$;

13 $Cs\;\leftarrow\;C_{t}[1:aux]$;

14 $exps\;\leftarrow\;\exp(-r(t-Cs))$;

15 $integ\;\leftarrow\;sum(Qs*exps)$;

16 $lam\;\leftarrow\;\lambda^{\infty}(I_{t},R_{t})$;

17 return $lam+integ$;

18

19 end if

Algorithm 1 Computation of the Intensity function at time $t$.

Next, we will show the algorithm to simulate the state dependent Hawkes
process. Notice how the intensity function needs to be updated between new
infections as recoveries affect the intensity function also.

Result: Simulation of State-Dependent Hawkes Process with intensity function
$\Lambda(t)=\lambda^{\infty}(i,\rho)+\int_{0}^{t}Q_{s}e^{-r(t-s)}dN_{s}$

1 Inputs: $T$, Qvals, Qprobs, $i_{0},\rho_{0},\ldots$;

2 Output: SDHP=function();

3 Initialization;

4 tIRarray $\;\leftarrow\;[0,i_{0},\rho_{0}]^{T}$;

5 $t\;\leftarrow\;0$;

6 while  _$t\leq T$_ do

7 $M\;\leftarrow\;$Intensity $(t,\ldots)$;

8 if  _$M\leq 0$_ then

9 break;

10

11 end if

12 $E\sim\text{Exp}(M)$;

13 $t\;\leftarrow\;t+E$;

14 $U\sim\text{Unif}(0,M)$;

15

16 end while

17## Since $\lambda^{\infty}(i,\rho)$ can change between increments of $C_{t}$
(via only a recovery), we need to update our intensity function between 0 and
the next time candidate for a jump of $C_{t}$ which is $E$;

18 $\text{Raux}\;\leftarrow\;0$;

19 while  _$\text{Raux} <E$_ do

20 $\text{PossibleR}\sim\text{Exp}(\;\mu\cdot\text{Icurrent}\;)$;

21 if  _$(\text{PossibleR+Raux}) >E$_ then

22 break

23 end if

24 $\text{Raux}\;\leftarrow\;\text{Raux}+\text{PossibleR}$;

25 Update tIRarray;

26

27 end while

28if  _$t <E$_ then

29 ## Perform a thinning routine to accept new Infections;

30
$\Lambda_{N}\;\leftarrow\;{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}{\tt
Intensity}}(t,\ldots)$;

31 if  _$U <\Lambda_{N}$_ then

32 Accept $t$ as time of Infection;

33 Update all parameters;

34

35 end if

36

37 end if

Algorithm 2 Simulation of a trajectory of the State-dependent Hawkes SIR model

## References

  * Adiga et al., (2020) Adiga, A., Chen, J., Marathe, M., Mortveit, H., Venkatramanan, S., and Vullikanti, A. (2020).  Data-driven modeling for different stages of pandemic response.  Journal of the Indian Institute of Science, 100(4):901–915.
  * Allen, (2017) Allen, L. J. (2017).  A primer on stochastic epidemic models: Formulation, numerical simulation, and analysis.  Infectious Disease Modelling, 2(2):128–142.
  * Bertozzi et al., (2020) Bertozzi, A. L., Franco, E., Mohler, G., Short, M. B., and Sledge, D. (2020).  The challenges of modeling and forecasting the spread of covid-19.  Proceedings of the National Academy of Sciences, 117(29):16732–16738.
  * Brémaud, (2020) Brémaud, P. (2020).  Point Process Calculus in Time and Space: An Introduction with Applications, volume 98.  Springer Nature.
  * Cartea et al., (2021) Cartea, Á., Cohen, S. N., and Labyad, S. (2021).  Gradient-based estimation of linear hawkes processes with general kernels.  arXiv preprint arXiv:2111.10637.
  * Cheng et al., (2020) Cheng, X., Han, Z., Abba, B., and Wang, H. (2020).  Regional infectious risk prediction of covid-19 based on geo-spatial data.  PeerJ, 8:e10139.
  * Chiang et al., (2020) Chiang, W.-H., Liu, X., and Mohler, G. (2020).  Hawkes process modeling of covid-19 with mobility leading indicators and spatial covariates.  medRxiv.
  * Chiang et al., (2022) Chiang, W.-H., Liu, X., and Mohler, G. (2022).  Hawkes process modeling of covid-19 with mobility leading indicators and spatial covariates.  International journal of forecasting, 38(2):505–520.
  * Chowdhury et al., (2021) Chowdhury, A., Shahbaz, M., Karim, R., Islam, J., Dan, G., and Shuixiang, H. (2021).  A comparative study on ivermectin-doxycycline and hydroxychloroquine-azithromycin therapy on covid-19 patients.  Eurasian J. Med. Oncol, 5:63–70.
  * Daley and Vere-Jones, (2003) Daley, D. J. and Vere-Jones, D. (2003).  An introduction to the theory of point processes: volume I: elementary theory and methods.  Springer.
  * Daley and Vere-Jones, (2008) Daley, D. J. and Vere-Jones, D. (2008).  An introduction to the theory of point processes: volume II: general theory and structure.  Springer New York.
  * dos Santos Gomes and de Oliveira Serra, (2021) dos Santos Gomes, D. C. and de Oliveira Serra, G. L. (2021).  Machine learning model for computational tracking and forecasting the covid-19 dynamic propagation.  IEEE Journal of Biomedical and Health Informatics, 25(3):615–622.
  * Escobar, (2020) Escobar, J. V. (2020).  A hawkes process model for the propagation of covid-19: Simple analytical results.  EPL (Europhysics Letters), 131(6):68005.
  * Foschi et al., (2020) Foschi, R., Lilla, F., and Mancini, C. (2020).  Warnings about future jumps: properties of the exponential hawkes model.  Available at SSRN 3639050.
  * Franco, (2020) Franco, E. (2020).  A feedback sir (fsir) model highlights advantages and limitations of infection-dependent mitigation strategies.  arXiv preprint arXiv:2004.13216.
  * Garetto et al., (2021) Garetto, M., Leonardi, E., and Torrisi, G. L. (2021).  A time-modulated hawkes process to model the spread of covid-19 and the impact of countermeasures.  arXiv preprint arXiv:2101.00405.
  * Greenwood and Gordillo, (2009) Greenwood, P. E. and Gordillo, L. F. (2009).  Stochastic epidemic modeling.  In Mathematical and statistical estimation approaches in epidemiology, pages 31–52. Springer.
  * Hazarika and Gupta, (2020) Hazarika, B. B. and Gupta, D. (2020).  Modelling and forecasting of covid-19 spread using wavelet-coupled random vector functional link networks.  Applied Soft Computing, 96:106626.
  * Koops et al., (2018) Koops, D. T., Saxena, M., Boxma, O. J., and Mandjes, M. (2018).  Infinite-server queues with hawkes input.  Journal of Applied Probability, 55(3):920–943.
  * Kurtz, (1981) Kurtz, T. G. (1981).  Approximation of population processes.  SIAM.
  * Laub et al., (2022) Laub, P. J., Lee, Y., and Taimre, T. (2022).  The Elements of Hawkes Processes.  Springer.
  * Lesage, (2020) Lesage, L. (2020).  A Hawkes process to make aware people of the severity of COVID-19 outbreak: application to cases in France.  PhD thesis, Université de Lorraine; University of Luxembourg.
  * Lewis et al., (2012) Lewis, E., Mohler, G., Brantingham, P. J., and Bertozzi, A. L. (2012).  Self-exciting point process models of civilian deaths in iraq.  Security Journal, 25(3):244–264.
  * Liu et al., (2022) Liu, Q.-H., Zhang, J., Peng, C., Litvinova, M., Huang, S., Poletti, P., Trentini, F., Guzzetta, G., Marziano, V., Zhou, T., et al. (2022).  Model-based evaluation of alternative reactive class closure strategies against covid-19.  Nature communications, 13(1):1–10.
  * Meyn and Tweedie, (2012) Meyn, S. P. and Tweedie, R. L. (2012).  Markov chains and stochastic stability.  Springer Science & Business Media.
  * Ogata, (1983) Ogata, Y. (1983).  Estimation of the parameters in the modified omori formula for aftershock frequencies by the maximum likelihood procedure.  Journal of Physics of the Earth, 31(2):115–124.
  * Rizoiu et al., (2018) Rizoiu, M.-A., Mishra, S., Kong, Q., Carman, M., and Xie, L. (2018).  Sir-hawkes: Linking epidemic models and hawkes processes to model diffusions in finite populations.  In Proceedings of the 2018 world wide web conference, pages 419–428.
  * Van den Driessche, (2017) Van den Driessche, P. (2017).  Reproduction numbers of infectious disease models.  Infectious Disease Modelling, 2(3):288–303.
  * Veen and Schoenberg, (2008) Veen, A. and Schoenberg, F. P. (2008).  Estimation of space–time branching process models in seismology using an em–type algorithm.  Journal of the American Statistical Association, 103(482):614–624.