Bootstrap Inference for Hawkes and General Point Processes

\begingroup\leftskip=0.8cm\rightskip=0.8cm

Inference and testing in general point process models such as the Hawkes model is predominantly based on asymptotic approximations for likelihood-based estimators and tests. As an alternative, and to improve finite sample performance, this paper considers bootstrap-based inference for interval estimation and testing. Specifically, for a wide class of point process models we consider a novel bootstrap scheme labeled `fixed intensity bootstrap' (FIB), where the conditional intensity is kept fixed across bootstrap repetitions. The FIB, which is very simple to implement and fast in practice, extends previous ideas from the bootstrap literature on time series in discrete time, where the so-called `fixed design' and `fixed volatility' bootstrap schemes have shown to be particularly useful and effective. We compare the FIB with the classic recursive bootstrap, which is here labeled `recursive intensity bootstrap' (RIB). In RIB algorithms, the intensity is stochastic in the bootstrap world and implementation of the bootstrap is more involved, due to its sequential structure. For both bootstrap schemes, we provide new bootstrap (asymptotic) theory which allows to assess bootstrap validity, and propose a `non-parametric' approach based on resampling time-changed transformations of the original waiting times. We also establish the link between the proposed bootstraps for point process models and the related autoregressive conditional duration (ACD) models. Lastly, we show effectiveness of the different bootstrap schemes in finite samples through a set of detailed Monte Carlo experiments, and provide applications to both financial data and social media data to illustrate the proposed methodology.

Keywords: Self-exciting point processes; conditional intensity; bootstrap inference; Hawkes process; autoregressive conditional duration models.

JEL Classification: C32.

\numberwithin{equation}{section} \endgroup

Introduction

Point processes are well-known to be useful tools to characterize dynamics of event occurrence times. This includes the homogeneous Poisson process where the intensity process is constant over time, the inhomogeneous Poisson process, where the intensity is a deterministic (or strictly exogenous) time-varying function, as well as the class of `self-exciting' point processes, such as the well-known and much applied Hawkes process. In particular, for the Hawkes process, the conditional intensity process depends on all past history of the events and thereby allows for (exponential or fractional) memory features, similar to autoregressive or fractional time-series processes in discrete time series econometrics. The self-exciting class of models, which are the focus of this paper, were originally proposed for modelling earthquake sequences; see O1988 and the references therein. Later, they have been put to use in a wide range of applications such as financial transactions B2007, BH2009, financial contagion ASCDL2015, monetary policy DMD2002, criminal fights and relations MSBST2011, forecasting electricity price spikes CHH2015 and the rich literature on social network information diffusion RLMX2017, among others. Hawkes processes are also closely related to the class of autoregressive conditional duration [ACD] models of ER1998, which are well known and much used in financial economics; see Sections (ref) and (ref) below for the relation between the two classes of processes.

Inference for self-exciting point process models is generally performed through classic, likelihood-based asymptotic inference and testing\footnote{As an alternative, the general methods of moments (GMM) has also been used, see e.g. ASCDL2015.}, as originally discussed in O1978. However, see e.g. R2018 and WSJ2010, the finite sample performance of asymptotic inference is not always satisfactory. This is in general the case because the finite sample distributions of the estimators are often very skewed and far from the Gaussian asymptotic distribution.

In this framework, a key motivation for the results presented in the paper is to provide a simple to implement, and theoretically well-grounded, bootstrap approach to inference in self-exciting point process models. We do this by providing six main contributions.

The first contribution is to propose a novel (non-)parametric bootstrap scheme for such point process models, which we label as `fixed intensity bootstrap' (FIB). The FIB is simple and fast to implement in practice -- particularly so when compared to existing (recursive) applications of the bootstrap. The key difference between the new and the classic bootstrap schemes is how to generate the sequence of waiting times in the bootstrap world. Specifically, while for standard, recursive bootstrap schemes, the bootstrap event times are generated recursively through the past bootstrap events, for the novel bootstrap scheme the bootstrap event times are generated using a `fixed' conditional intensity function, which entirely depends on the event times in the original world. Therefore, the FIB contrasts with existing implementations of the bootstrap, see e.g. ELL2011 and SNCWWSB2011, which utilize a (possibly highly complex and time consuming) sequential update of the bootstrap conditional intensities.

The second contribution is to provide bootstrap (first-order) asymptotic theory, including establishing bootstrap validity for inference and testing in point process models for both the novel (FIB) bootstrap and for the classic recursive bootstrap (for which no theory exists in the literature). We show that the bootstrap based on the FIB is valid under regularity conditions which are milder than those required for validity of recursive bootstrap schemes (hereafter, RIB).

The third contribution is to introduce novel `non-parametric' implementations of the FIB and RIB schemes, which are based on resampling time-changed transformations of the original waiting times, rather than generating the (transformed) waiting times through a parametric model (usually the exponential distribution), as is in the literature. These implementations are likely to be robust to model misspecifications which generate non-exponential transformed waiting times. We show how to scale the original time-changed waiting times properly and to resample them; we also show that, in the homogeneous case, validity of the implied bootstrap follows from a time-change functional central limit theorem derived by B1968 which, as far as we are aware, has never been applied to the bootstrap of self-exciting point process models.

The fourth contribution of the paper is a detailed Monte Carlo simulation study on the performance of the bootstrap for self-exciting Hawkes processes. Possibly due to the high computational costs involved in the implementation of a simulation study for the bootstrap in this framework, to the best of our knowledge studies like ours have not been attempted in the literature. We show that for Hawkes processes with exponential kernels, the coverage probabilities of confidence intervals based on the Gaussian asymptotic approximation may be well below the nominal level. In contrast, the bootstrap is able to correct this and, in particular, FIB implementations are particularly well performing in terms of coverage probabilities.

The fifth contribution is to provide two real data examples where we illustrate the key differences between asymptotic and the various bootstrap inference methods in applications. The first refers to the problem of modeling and predicting extreme financial results, see ELL2011. We use this example, including the data sample considered in ELL2011, to compare the outcome of the four different bootstrap schemes discussed in the paper. The second example is based on social media data and considers the flows of tweets and re-tweets proceeding and following a political announcement. Specifically, using recent tweets related to the COVID-19 pandemic in Denmark, we show how bootstrap-based inference is able to detect structural breaks (in the mean intensity as well as in the decay rate of intensity) induced by the announcement, which may not be detected based on asymptotic inference.

Sixth, we discuss the link between the proposed bootstraps based on the point process representation and bootstrap inference for autoregressive conditional duration (ACD) models. We establish the relation between our proposed bootstrap schemes and bootstrap algorithms based on the ACD representation. Specifically, we show that our recursive bootstrap corresponds to a residual-based bootstrap in the ACD world PS2021, with the crucial difference that the number of events generated through our scheme is random, rather than being fixed. This is a key improvement, as our bootstrap ensures that the sum of the bootstrap waiting times always cover the original time interval. We also show that residual-based implementation of the bootstrap in the ACD world corresponds to our proposed non-parametric bootstrap with re-sampling based on the (estimated) transformed waiting times. Finally, we discuss the relation between our proposed fixed intensity bootstrap and a bootstrap in the ACD framework, which is novel in the literature, where the conditional duration in the bootstrap world is fixed to the estimated conditional duration in the original sample.

Structure of the paper

The paper is organized as follows. In Section (ref) likelihood-based analysis for point processes inference is presented, and in Section (ref) the novel fixed intensity, as well as the recursive intensity, bootstraps are discussed, with theory and validity results in Section (ref). Non-parametric bootstrap is discussed in Section (ref), and the relation between our bootstraps and the bootstraps for ACD models is discussed in Section (ref). Section (ref) provides a Monte Carlo study of the different schemes. Section (ref) contains two empirical illustrations, and Section (ref) concludes. All proofs are contained in the Appendix.

Notation

We use the counting process $N(t)$ to characterize the total number of events occurring before and including time $t$, with $N(s,t]$ and $N[s,t)$ the numbers of events in the interval $(s,t]$ and $[s,t)$, respectively, for $s<t$. For a right-continuous natural filtration $(\mathcal{F}_{t} )_{t\in\mathbb{R}}$ of a continuous time stochastic process, we denote by $\mathcal{F}_{t-}$ the left limit of $\mathcal{F}_{t}$, which contains all the information before but not including time $t$. We use $\mathbb{I}(\cdot)$ to denote the indicator function, and define $\mathbb{R}^{+}:=(0,\infty)$ and $\mathbb{R}_{+}:=[0,\infty)$. For $x\in\mathbb{R}$, $\lfloor x\rfloor :=\max_{z\in\mathbb{Z}}\{z\leq x\}$. For the bootstrap, as is standard, we denote by $P^{\ast}$ the probability measure induced by the bootstrap; expectation and variance computed under $P^{\ast}$ are denoted by $E^{\ast}$ and $V^{\ast}$, respectively. For a sequence $X_{T}^{\ast}$ computed on the bootstrap data, $X_{T}^{\ast}\overset{p^{\ast}}{\rightarrow}_{p}0$ or $X_{T}^{\ast}=o_{p}^{\ast}(1)$, in probability, denote that $P^{\ast} (|X_{T}^{\ast}|>\epsilon)\rightarrow0$ in probability for any $\epsilon>0$; $X_{T}^{\ast}=O_{p}^{\ast}(1)$, in probability, denotes that there exists a $c>0$ such that $P^{\ast}(|X_{T}^{\ast}|>c)\rightarrow0$ in probability; with $X_{T}^{\ast}\overset{d^{\ast}}{\rightarrow}_{p}X$ (weak convergence in probability) we mean that $E^{\ast}(g(X_{T}^{\ast}))\overset{p}{\rightarrow }E(g(X))$ for all continuous bounded functions $g$, in each case as $T\rightarrow\infty$. Finally, $\mathcal{N}$ denotes a Gaussian random variable and, for $\mu>0$, $\mathcal{E}(\mu)$ denotes an exponential random variable with mean $1/\mu$.

Likelihood-based analysis of the point process

We discuss here likelihood-based estimation for a general class of point process models. For later use when establishing asymptotic validity of the bootstrap, we state explicit sufficient conditions for classic likelihood-based asymptotic theory. Precisely, and as in O1978, we establish consistency and limiting distributions of likelihood-based estimators, as well as the related (likelihood ratio) test statistics.

The model

By assumption, the observed event times are realizations from a univariate point process, i.e. a collection $\{t_{i}\}_{i=1}^{\infty}$, $t_{i}>0$, of stochastic event times with associated waiting times (or durations), $w_{i}:=t_{i}-t_{i-1}$, for $i=1,2,...$ with $t_{0}:=0$; see e.g. DVJ2003 for an introduction to point processes. The point process can be equivalently characterized by the continuous-time counting process

for $t\geq0$, with associated filtration $(\mathcal{F}_{t}),t\geq0$ where $\mathcal{F}_{t}$ is the $\sigma$-field generated by $\{N(s),s\leq t\}$.

In addition, and as used here predominantly, a regular point process is uniquely defined by its conditional intensity process, $\lambda(t)$, $t\geq0$, which captures the instantaneous conditional probability of event occurrences\footnote{Note that the limit on the right-hand side of (ref) is assumed to exist, such that the conditional distribution of the waiting times is continuous.} and is defined as

Observe that, as the point process is assumed to be regular and orderly, $\lambda(t)$ essentially captures the instantaneous conditional probability of observing a single event at each time $t$.

A key example used throughout is the `self-exciting' Hawkes point process, where the conditional intensity is given by

where the $\mu>0$ is the baseline intensity and $\gamma(t)$ is the so-called kernel function, which typically is either exponential,

or following a power law,

where $\alpha,\beta,\delta\geq0$.

Note as the sum in ((ref)) is over all events $t_{i}$ prior to $t$, the Hawkes process has infinite (or long) memory. In contrast, if $\gamma(t)=0$, $\lambda(t)=\mu>0$, then the point process reduces to a homogeneous Poisson process which has i.i.d. exponentially distributed waiting times $w_{i}$ with rate $\mu$, that is, the $w_{i}$'s are i.i.d. $\mathcal{E}(\mu)$ distributed. Likewise, an example of counting process with finite memory (or, $q$ `lags'), sometimes referred to as a `Wold process', is given by

where $\gamma(\cdot)$ is a mapping from $\mathbb{R}_{+}^{q}$ to $\mathbb{R} _{+}$. A specific example is given by \[ \gamma(t-t_{N(t-)},...,t-t_{N(t-)-q+1};\theta)=\sum_{i=1}^{q}\gamma _{i}(t-t_{N(t-)-i+1};\theta), \] with $\gamma_{i}(\cdot)$ being exponential or power law kernel functions as in (ref) and (ref) for $i=1,2,...,q$. Notice that for $q=1$, this is an example of a renewal process, with associated i.i.d. waiting times $w_{i}$ which are not exponentially distributed.

The class of self-exciting point process models is also linked to the ACD model of ER1998, which is based on the following dynamic equation for the waiting times $w_i:=t_i-t_{i-1}$ between events:

where the $\varepsilon_{i}$'s are strictly positive i.i.d. random variables with mean one. The ACD model can be given a point process representation; specifically, the conditional intensity associated to the model, see ER1998, takes the form

where $\lambda_{\varepsilon}(\cdot)=p_{\varepsilon}(\cdot)/S_{\varepsilon }(\cdot)$, with $p_{\varepsilon}$ and $S_{\varepsilon}$ denoting the pdf and the survival function of $\varepsilon_{i}$, respectively. A simple example is the ACD(1) with exponential errors, \[ w_{i}=\psi_{i}\varepsilon_{i},\quad \psi_{i}=\omega+\alpha w_{i-1},\quad \varepsilon_{i}\sim\mathcal{E}(1), \] with intensity given by the piecewise constant function \[ \lambda(t)=\frac{1}{\psi_{N(t-)+1}} =\frac{1}{\omega+\alpha w_{N(t-)}} =\frac{1}{\omega+\alpha(t_{N(t-)}-t_{N(t-)-1})}, \] which is a special case of a Wold process with two `lags'. Should $\varepsilon_{i}$ be a continuous, non-exponentially distributed random variable, it follows from ((ref)) that the intensity of the ACD(1) takes the form \[ \lambda(t) =\lambda_{\varepsilon}\left(\frac{t-t_{N(t-)}}{\omega+\alpha w_{N(t-)}}\right) \frac{1}{\omega+\alpha w_{N(t-)}}. \] Finally, with $\alpha=0$ the ACD reduces to a renewal process with intensity $\lambda(t)=\lambda_{\varepsilon}(\omega^{-1}(t-t_{N(t-)}))\omega^{-1}$.

Likelihood-based estimation

For the statistical analysis we assume that the conditional intensity $\lambda(t)$ in ((ref)) is parameterized by a finite-dimensional vector of unknown parameters $\theta\in\Theta\subseteq\mathbb{R}^{d}$, $d:=\dim\theta$. To emphasize the dependence of the intensity on $\theta$, we write $\lambda(t;\theta)$ and, for the associated counting process, $N(t;\theta)$. For notational convenience, when evaluated at the true value, we write $\lambda(t;\theta_{0})=:\lambda(t)$ and $N(t;\theta_{0})=:N(t)$.

Consider a sample of event times $t_{1},t_{2},\ldots,t_{n_{T}}$ observed in a time interval $[0,T]$, with $n_{T}=N(T)$ the total number of events in the interval. Standard arguments as in DVJ2003 imply that the joint log-likelihood function $\ell_{T}(\theta)$ can be written as

where $\Lambda(\cdot;\theta)$ is the so-called integrated intensity, given by

and we assume $t_{n_{T}}=T$ (such that $T$ coincides with the last event time) in deducing the second equality in (ref). The maximum likelihood estimator (MLE) $\hat{\theta}_{T}$ is defined by

For the Hawkes model, $\lambda(t;\theta)$ is given by (ref) with $\theta=(\mu,\alpha,\beta)^{\prime}$ for the exponential kernel in (ref) and $\theta=(\mu,\alpha,\beta,\eta)^{\prime}$ for the power law kernel in (ref). Then the log-likelihood function $\ell_{T}(\theta)$ in (ref) becomes

Note that for the special case of a homogeneous Poisson process where $\lambda(t)=\mu$, then $\theta=\mu$, and the log-likelihood simplifies to $\ell_{T}(\theta)=n_{T}\log\theta-T\theta$. Hence, in this special case the MLE has the closed form $\hat{\theta}_{T}=n_{T}/T$.

Asymptotic theory

For the asymptotic theory of $\hat{\theta}_{T}$ we assume that the information set $\mathcal{F}_{t}$ is defined as the $\sigma$-field generated by $\{N(s,t],-\infty<s\leq t\}$. Under mild requirements O1978, the analysis presented below extends to the case where $\mathcal{F}_{t}=\{N(s),0\leq s\leq t\}$. Likewise, we assume for simplicity that $t_{n_{T}}=T$.

A key role in the asymptotic analysis here -- as well as for the novel bootstrap asymptotics below -- is played by the Doob-Meyer decomposition of $N(t)$ in (ref), which is given by \[ N(t)=M(t)+A(t). \] Here $M$ is a square integrable continuous-time $\mathcal{F}_{t}$-local martingale and $A(t)$ is the compensator of $N(t)$, which in this case is given by the integrated intensity $\Lambda(t;\theta_{0})=\Lambda(t)$ in (ref); that is, $A(t)=\int_{0}^{t} \lambda(s)ds=\Lambda(t)$. By definition,

is a continuous-time martingale, and we may write $E(dM(t)|\mathcal{F} _{t-})=E(dN(t)-\lambda(t)dt|\mathcal{F}_{t-})=0$ for any $t>0$. Since $\lambda(t)$ is $\mathcal{F}_{t-}$-measurable, it follows that $E(dN(t)|\mathcal{F}_{t-})=\lambda(t)dt$ O1978 which will be used repeatedly throughout for both the standard and the bootstrap asymptotic analyses. Furthermore, we make the following technical assumptions.

\noindentAssumption 1

(a) The parameter space $\Theta\subseteq\mathbb{R}^{d}$ is compact with $\theta_{0}\in\Theta_{0}\subset\Theta$;

(b) For $\theta\in\Theta_{0}$, with $N(\cdot;\theta)$ denoting the counting process $N(\cdot)$ indexed by $\theta$, $N(\cdot ;\theta)$ is an orderly point process with stationary and ergodic increments. Moreover, $E\big(\sup_{m\geq1}m(N(t+1/m;\theta)-N(t;\theta ))^{2}\big)<\infty$;

(c) the intensity process $\lambda(t;\theta)$ satisfies the following conditions almost surely: (i) it is predictable (left-continuous) for all $\theta$, continuous in $\theta$ and strictly positive; (ii) for all $\theta$, $|\lambda(t;\theta)|\leq\xi_{1}(\theta)$ with $E(\xi_{1}(\theta )^{2})<\infty$; (iii) $\lambda(t;\theta_{1})=\lambda(t;\theta_{2})$ if and only if $\theta_{1}=\theta_{2}$.

Notice that Assumption 1(c) implies in particular that $\log\lambda(t;\theta)$ has a finite second order moment.

Consistency of the MLE is given in the next theorem from O1978.

For the analysis of the score and the information, and for establishing the asymptotic normality of the MLE, we make use of the Assumption 2 below, where we use the following notation: for any function $f(t;\theta)$ of $\theta$ (and $t$), $f(t):=f(t;\theta_{0})$, $\partial_{\theta}f(t;\theta)=\partial f(t;\theta)/\partial\theta$ and $\partial_{\theta_{0}}f(t)=\left. \partial f(t;\theta)/\partial\theta\right\vert _{\theta=\theta_{0}}$ (and similarly for higher order and partial derivatives).

\noindentAssumption 2

(a) The intensity process $\lambda(t;\theta)$ satisfies the following conditions almost surely: (i) $\lambda(t;\theta)$ is continuously differentiable with respect to $\theta$ up to order three, for all $t\geq0$; (ii) $E((\partial_{\theta_{i}}\lambda(t;\theta))^{2})<\infty$ and $E((\partial_{\theta_{i},\theta_{j}}^{2}\lambda(t;\theta))^{2})<\infty$ for all $\theta$;

(b) With $h(t;\theta):=\lambda(t;\theta)^{-1}(\partial _{\theta}\lambda(t;\theta))(\partial_{\theta}\lambda(t;\theta))^{\prime}$ and $I(\theta):=E(h(t;\theta))$, it holds that $I(\theta_{0})>0$ and each element of $h(t):=h(t;\theta_{0})$ has finite variance;

(c) With $N_{\epsilon}(\vartheta)$ denoting a neighborhood of $\vartheta$, for all $\vartheta\in\Theta$, \[ \sup_{\theta\in N_{\epsilon}(\vartheta)}|\partial_{\theta_{i},\theta _{j},\theta_{k}}^{3}\lambda(t;\theta)|\leq c_{ijk}(t),\quad\sup_{\theta\in N_{\epsilon}(\vartheta)}|\partial_{\theta_{i},\theta_{j},\theta_{k}}^{3} \log\lambda(t;\theta)|\leq d_{ijk}(t) \] where $c_{ijk}(t),d_{ijk}(t)$ are stationary and ergodic processes with $E(c_{ijk}(t))<\infty$ and $E(\lambda(t)^{2}d_{ijk}^{2}(t))<\infty$.

Note that Assumption 2(c) differs from standard requirements as in O1978 which address uniformity of the Hessian.

Next, let $S_{T}(\theta):=\partial_{\theta}\ell_{T}(\theta)$ and $H_{T} (\theta):=\partial_{\theta}^{2}\ell_{T}(\theta)$ denote the score and the Hessian, respectively. Using (ref) and (ref) it holds that

Applying JR2004 the following theorem holds, where we here also consider the distribution of the likelihood ratio test statistic $LR_T(\theta_{0})$ for a simple null hypothesis $H_{0}:\theta=\theta_{0}$.

The bootstrap

We discuss here two bootstrap schemes. The first bootstrap, which is novel, is denoted as the `fixed intensity bootstrap' (FIB). The FIB as proposed here builds on ideas from the `fixed design bootstrap' in regression (and time series) models W1986,GK2004, as well as the so-called `fixed volatility bootstrap' in conditional volatility modelling CPR2018, in the sense that the bootstrap intensity function is fixed across bootstrap repetitions. The second scheme, which has been applied in e.g. ELL2011 and SNCWWSB2011, is here denoted the `recursive intensity bootstrap' (RIB). As will be discussed later, in practice the FIB is simpler and faster to implement than the RIB, in addition to being valid under milder regularity conditions. Since no theory exists for either of the FIB and the RIB schemes, in Section (ref) we establish validity of the bootstrap for both.

Random time change

A key property we employ in defining our bootstrap algorithms is that using the integrated intensity to transform the original event times $\{t_{i}\}$ to another sequence of event times $\{s_{i}\}$ gives a homogeneous Poisson process with unit intensity. Equivalently, the original and non i.i.d. waiting times $\{w_{i}\}$, $w_{i}=t_{i}-t_{i-1}$, can be transformed into new waiting times $\{v_{i}\}$, $v_{i}=s_{i}-s_{i-1}$, which are i.i.d. $\mathcal{E}\left( 1\right) $, see also DVJ2003.

The time change transformation $t_{i}\mapsto s_{i}$ is given by

where the integrated intensity $\Lambda(t;\theta)$ is defined in ((ref)). Moreover, with

the associated transformed waiting times $v_{i}(\theta)$ are given by

for $i=1,2,...$. By definition, at the true value $\theta_{0}$ the transformed waiting times $v_{i}:=v_{i}(\theta_{0})=\Lambda(t_{i},t_{i-1};\theta_{0})$ are i.i.d. $\mathcal{E}\left( 1\right) $, such that the transformed event times, $s_{i}:=s_{i}(\theta_{0})$, form a homogeneous Poisson process with unit intensity. For the Hawkes process in (ref), $\lambda(t)=\mu+\sum_{t_{j}<t}\gamma(t-t_{j})$, and hence

which, for the exponential kernel in ((ref)), reduces to

For the implementation of the bootstrap, the reverse time transformation $s_{i}\mapsto t_{i}$ is of key interest. Specifically, consider initially a sequence $\{v_{i}\}$ of waiting times in the transformed time scale, generated as i.i.d. and $\mathcal{E}\left( 1\right) $-distributed. Then, under the true model, we can (numerically) invert the mapping (ref) and generate the $i$th waiting time $w_{i}$ (or, equivalently, the $i$th event time) recursively in terms of the $i$th waiting time in transformed time scale $v_{i}$ and the past event times $\{t_{j},j=1,...,i-1\}$. The recursion is initiated by generating the first waiting time $w_{1}$ as $w_{1}=\Lambda ^{-1}(v_{1};\theta_{0})$, where $v_{1}$ is the first ($\mathcal{E}\left( 1\right) $-distributed) waiting time in transformed time scale.

As detailed below, for both the FIB and RIB schemes, we generate i.i.d. random event times in the transformed time scale, which are next transformed to the original time scale using the intensity dynamics estimated from the data. The key difference between the two algorithms is whether the transformation from transformed to original event times is fixed or sequential (and hence random) across bootstrap samples.

Fixed intensity bootstrap

Given a sample of event times $\{t_{i}\}_{i=1}^{n_{T}}$ in $[0,T]$, fix the bootstrap true value parameter, $\theta_{T}^{\ast}$. As is standard, one may for example set $\theta_{T}^{\ast}=\hat{\theta}_{T}$, the unrestricted MLE based on $\{t_{i}\}_{i=1}^{n_{T}}$; for hypothesis testing, one may also set $\theta_{T}^{\ast}=\tilde{\theta}_{T}$, the MLE restricted by the null hypothesis.

For the FIB, where the intensity is kept fixed across replications, denote the intensity process implied by the bootstrap true value as $\hat{\lambda}(t):=\lambda(t;\theta_{T}^{*})$, and the corresponding integrated intensity process as

By definition, $\hat{\lambda}(t)$ and $\hat{\Lambda}(t)$ depend on the original data through the observed event times $\{t_{i}\}_{i=1}^{n_{T}}$ and bootstrap true value $\theta_{T}^{*}$. Therefore, by construction, $\hat{\lambda}(t)$ and $\hat{\Lambda}(t)$ are known and fixed conditionally on the data.

\noindentAlgorithm 1 (FIB)

(i) Generate a (conditionally on the original data) i.i.d. sample $\{v_{i}^{\ast}\}$ of bootstrap transformed waiting times from the $\mathcal{E}\left( 1\right) $ distribution; the bootstrap transformed event times are then given by $\{s_{i}^{\ast}\}$ where $s_{i}^{\ast}=\sum_{j=1} ^{i}v_{j}^{\ast}$.

(ii) Construct the bootstrap event times in the original time scale as \[ t_{i}^{\ast}=\hat{\Lambda}^{-1}(s_{i}^{\ast}), \] for $i=1,\ldots,n_{T}^{\ast}$, where $\hat\Lambda$ is defined in (ref) and the number of bootstrap events $n_{T}^{\ast}$ is \[ n_{T}^{\ast}:=\max\{k:s_{k}^{\ast}\leq\hat{\Lambda}(T)\}=\max\left\{ k:t_{k}^{\ast}\leq T\right\} \text{;} \] the associated bootstrap counting process is $N^{\ast}(t):=\sum_{i\geq 1}\mathbb{I}(t_{i}^{\ast}\leq t)$, for $t\in\lbrack0,T]$.

(iii) Define the bootstrap MLE as $\hat{\theta}_{T}^{*} :=\arg\max_{\theta\in\Theta}\ell_{T}^{*}(\theta)$ with bootstrap log-likelihood

Some remarks are in order.

Recursive intensity bootstrap

The RIB resembles the recursive bootstrap in time series models, see e.g. CR2021 for a review. Thus, and in contrast to the FIB, the RIB conditional intensity, denoted here by $\lambda^{*}(t;\theta)$, is constructed using the functional form of the original intensity $\lambda(t;\theta)$, but in terms of recursively obtained bootstrap event times $t_{i}^{*}$. This entails that, for any $\theta\in\Theta$, $\lambda ^{*}(t;\theta)$ is a random process, even conditionally on the original data, and hence differs from the FIB intensity, which is fixed across bootstrap repetitions. Note also that the recursively obtained bootstrap intensity process $\lambda^{*}(t;\theta)$ inherits the same properties, in terms of e.g. differentiability with respect to $\theta$, of the original intensity process $\lambda(t;\theta)$.

We define $\lambda^{*}(t):=\lambda^{*}(t;\theta_{T}^{*})$ and

The RIB is then defined as follows.

\noindentAlgorithm 2 (RIB)

(i) As in Algorithm 1.

(ii) For $i=1,...,n_{T}^{\ast}$, construct the bootstrap event times $t_{i}^{\ast}$ in the original time scale recursively (see Remark (ref) below) as \[ t_{i}^{\ast}=\Lambda^{\ast-1}(s_{i}^{\ast}), \] for $i=1,\ldots,n_{T}^{\ast}$, where the number of bootstrap events $n_{T}^{\ast}$ is \[ n_{T}^{\ast}:=\max\{k:s_{k}^{\ast}\leq\Lambda^{\ast}(T)\}=\max\left\{ k:t_{k}^{\ast}\leq T\right\} \] and $\Lambda^{\ast}(t)$ is defined in (ref); the associated bootstrap counting process is $N^{\ast}(t):=\sum_{i\geq1}\mathbb{I} (t_{i}^{\ast}\leq t)$, for $t\in\lbrack0,T]$.

(iii) Define the bootstrap MLE as $\hat{\theta}_{T}^{*} :=\arg\max_{\theta\in\Theta}\ell_{T}^{*}(\theta)$ with bootstrap log-likelihood

Validity of bootstrap inference

In this section, we establish bootstrap asymptotic validity for the FIB and RIB bootstrap schemes outlined above. As emphasized the bootstrap true parameter is assumed to be consistent, $\theta_{T}^{*}\rightarrow_{p} \theta_{0}$, which holds, e.g., for the particular choices where $\theta_{T}^{*} =\hat{\theta}_{T}$ (unrestricted bootstrap) or $\theta_{T}^{*}=\tilde{\theta }_{T}$ (restricted bootstrap) under the null.

Throughout, we let $\mathcal{F}_{t}^{\ast}$ denote the $\sigma$-field generated by $\{N^{\ast}(s),$ $0\leq s\leq t\}$ and $\mathcal{F}_{t-}^{\ast}$ be its left limit. Notice that, since the distribution of $N^{\ast}$ depends on $T$, formally we have an array $\mathcal{F}_{T,t}^{\ast}:=\{N_{T}^{\ast }(s),$ $0\leq s\leq t\leq T$, $T\geq0\}$; for simplicity, in the following we suppress the dependence on $T$ and write $N_{T}^{\ast}(t)$ and $\mathcal{F} _{T,t}^{\ast}$ simply as $N^{\ast}(t)$ and $\mathcal{F}_{t}^{\ast}$.

Preliminaries

As for the non-bootstrap asymptotic analysis, define the bootstrap martingale \[ M^{\ast}(t)=N^{\ast}(t)-\Lambda_{N^{\ast}}(t). \] Here $\Lambda_{N^{\ast}}(t)$ is the integrated conditional intensity of either the FIB or the RIB bootstrap process $N^{\ast}(t)$ (see also Remark (ref)(i)) and hence it corresponds to the bootstrap compensator of $N^{\ast}(t)$ conditionally on the data. Consequently, $M^{\ast}(t)$ is a continuous-time $\mathcal{F}_{t}^{\ast}$ local martingale conditionally on the data. Moreover, for any process $\xi^{\ast}(t)$ which (conditionally on the original data) is predictable with respect to $\mathcal{F}_{t}^{\ast}$, the (Stieltjes) stochastic integral process

is also (conditionally on the original data) a continuous-time martingale.

Validity of the FIB

We first consider the FIB. From the bootstrap log-likelihood defined in (ref), we derive the corresponding bootstrap score and Hessian,

where $\xi(t;\theta):=\partial_{\theta}\log\lambda(t;\theta)$, $\zeta (t;\theta):=\partial_{\theta}^{2}\log\lambda(t;\theta)$ and $h(t;\theta)$ is defined in Assumption 2.

Notice that $S_{T}^{\ast}(\theta)$ and $H_{T}^{\ast}(\theta)$ depend on the bootstrap data only through $N^{\ast}(t)$ which, conditionally on the original data, is an inhomogeneous Poisson point process with fixed conditional intensity given by $\hat{\lambda}(t)=\lambda(t;\theta_{T}^{\ast})$. With $M^{\ast}(t):=N^{\ast}(t)-\hat{\Lambda}(t)=N^{\ast}(t)-\Lambda(t;\theta _{T}^{\ast})$, the score and Hessian evaluated at the bootstrap true value $\theta_{T}^{\ast}$ can be rewritten as

where $\hat{\xi}(t)=\xi(t;\theta_{T}^{\ast})$, $\hat{h}(t)=h(t;\theta _{T}^{\ast})$ and $\hat{\zeta}(t)=\zeta(t;\theta_{T}^{\ast})$.

Using the fact that $M^{\ast}$ is a martingale, we prove in the appendix the following lemma, which requires only a mild strengthening of the assumptions in Theorem (ref).

The following theorem shows the first-order validity of the FIB and of the associated likelihood ratio test.

Validity of the RIB

For the RIB, the bootstrap score and Hessian at the bootstrap true value $\theta_{T}^{\ast}$ mimic their counterparts on the original data, see (ref)--(ref). Specifically, with $\lambda^{\ast}(t):=\lambda^{\ast}(t;\theta_{T}^{\ast})$ and $M^{\ast }(t)=N^{\ast}(t)-\int_{0}^{t}\lambda^{\ast}(t)dt$,

where $h^{\ast}(t):=h^{\ast}(t;\theta_{T}^{\ast})$. The next lemma shows that the RIB score and Hessian mimic the large sample properties of the original score and Hessian. It requires an additional assumption, see (ref) below, which is not required for the FIB. In order to introduce it, we emphasize that the quantity $h(t;\theta)$ in Assumption 2 depends on the data generating process, and hence on the true parameter $\theta_{0}$. That is, $h(t;\theta)=h_{\theta_{0}}(t;\theta)$.

The proof is based on the fact that for any fixed $T$ and conditionally on the data, the bootstrap sample can be made stationary.

For bootstrap consistency, we modify Assumption 2(c) as follows.

\noindentAssumption 2(c$^{\ast}$)

Assumption 2(c) holds with $c_{ijk}(t)=c_{ijk}(t;\theta_{0})$ and $d_{ijk}(t)=d_{ijk}(t;\theta_{0})$ replaced by $\sup_{\theta\in\Theta_{0} }(c_{ijk}(t;\theta))$ and $\sup_{\theta\in\Theta_{0}}(d_{ijk}(t;\theta))$, respectively.

The modification is necessary in order to bound the third order derivatives of the RIB likelihood.

Non-parametric FIB and RIB

In the presented parametric bootstrap, bootstrap event times are obtained in transformed time scale by cumulating randomly-generated i.i.d. $\mathcal{E}\left( 1\right) $ waiting times. This was motivated by the fact that waiting times $v_{i}=v_{i}(\theta_{0})$ in (ref) are i.i.d. $\mathcal{E}\left( 1\right) $-distributed for $i=1,...,n_{T}$ and, moreover, with $\theta_{T}^{\ast}=\theta_{0}+o_{p}(1)$,

However, in the case of a misspecified model, it may be the case that the transformed waiting times $\hat{v}_{i}$ are not exponentially distributed (asymptotically). Therefore, we consider here the point process bootstrap equivalent of the well-known residual-based i.i.d. bootstrap in discrete time series models. Specifically, after the point process model is fit to data, the residuals to resample from can be taken as the waiting times in transformed time scale, i.e. $\hat{v}_{i}$, $i=1,...,n_{T}$. Then, the bootstrap waiting times in transformed time can be generated as an i.i.d sample from the sample $\{\hat{v}_{i}\}_{i=1}^{n_{T}}$. This algorithm is denoted here as the `non-parametric bootstrap', and can be implemented for both FIB and RIB bootstraps, see below.

For the bootstrap in conditional mean and variance time series models, the residuals are typically centered and/or scaled prior to the implementation of the bootstrap. Similarly, here the waiting times $\hat{v}_{i}$ need to be properly standardized, such that the bootstrap transformed waiting times $v_{i}^{\ast}$ match (as a minimum) the mean of the $\mathcal{E}\left( 1\right) $ distribution, i.e. $E^{\ast}(v_{i}^{\ast})=1$. This is achieved by sampling from $\hat{v}_{i}^{c}$ given by

where $\bar{v}_{T}:=n_{T}^{-1}\sum_{j=1}^{n_{T}}\hat{v}_{j}$. Note that $\hat{v}_{i}^{c}>0$ for all $i$, and therefore a random draw from $\{\hat {v}_{i}^{c}\}_{i=1}^{n_{T}}$ has, conditionally on the original data, unit expected value, i.e. $E^{\ast}(v_{i}^{\ast})=n_{T}^{-1}\sum_{i=1}^{n_{T}} \hat{v}_{i}^{c}=1$.

With the transformed waiting times $\{\hat{v}_{i}^{c}\}_{i=1}^{n_{T}}$ defined in (ref), the proposed non-parametric bootstrap algorithm is as follows.

\noindentAlgorithm 3 (Non-Parametric Bootstrap)

(i) Generate a sample $\{v_{i}^{*}\}$ of bootstrap transformed waiting times by resampling with replacement from $\{\hat{v} _{i}^{c}\}_{i=1}^{n_{T}}$, such that

where $u_{i}^{\ast}$ is an i.i.d. discrete uniformly distributed sequence on $\{1,\dots,n_{T}\}$. The bootstrap transformed event times are then given by $s_{i}^{\ast}=\sum_{j=1}^{i}v_{j}^{\ast}$.

(ii)-(iii) as in Algorithm 1 or Algorithm 2 depending on whether it is a fixed intensity or recursive intensity bootstrap.

To provide an intuition about validity of this bootstrap and about the importance of rescaling, consider a simple Poisson process model with intensity $\lambda(t)=\theta$, where interest is in inference on $\theta$ using the (unrestricted) bootstrap. Recall that the log-likelihood for the original sample is $\ell_{T}(\theta)=\int\log\theta dN(t)-\int\theta dt=n_{T}\log\theta-T\theta$, with associated bootstrap score $\theta^{-1}\int dN(t)-T=\theta^{-1}n_{T}-T$, which leads to the unique MLE, $\hat{\theta} _{T}=n_{T}/T$. To implement the non-parametric bootstrap, consider the transformed waiting times, see Section (ref), which in this case are given by $\hat{v}_{i}=\hat{\theta}_{T}w_{i}$, with $w_{i} =t_{i}-t_{i-1}$ the original observed waiting times. The non-parametric bootstrap generates the $v_{i}^{\ast}$'s by initially resampling from the rescaled $\hat{v}_{i}^{c}$ defined in (ref); next, the $v_{i} ^{\ast}$'s are transformed back in the original time scale using the inverse mapping $w_{i}^{\ast}=v_{i}^{\ast}/\hat{\theta}_{T}$. This leads to the bootstrap event times $t_{i}^{\ast}:=\sum_{j=1}^{i}w_{j}^{\ast}$ with associated bootstrap counting process $N^{\ast}(t):=\sum_{i\geq1} \mathbb{I}(t_{i}^{\ast}\leq t)$. The bootstrap likelihood and score are then given by $\ell_{T}^{\ast}(\theta)=\int\log\theta dN^{\ast}(t)-\int\theta dt=n_{T}^{\ast}\log\theta-T\theta$ and $S_{T}^{\ast}(\theta)=\theta^{-1}\int dN^{\ast}(t)-T=\theta^{-1}n_{T}^{\ast}-T$, respectively, where as earlier $n_{T}^{\ast}$ denotes the total number of events, $n_{T}^{\ast}=\max \{k:\sum_{1}^{k}w_{i}^{\ast}\leq T\}$.

Consider next the bootstrap score at the true value $\theta_{T}^{\ast} =\hat{\theta}_{T}$,

Because of the rescaling in (ref), $E^{\ast}(v_{i}^{\ast}-1)=0$. This is a key feature for the bootstrap score to mimic the large-sample behavior of the original score. In contrast, without rescaling, the bootstrap mean of $v_{i}^{\ast}-1$ would be of order $O_{p}(n_{T}^{-1/2})=O_{p} (T^{-1/2})$ (the order being sharp), thereby introducing an asymptotically non-negligible (random) bias term in the distribution of the bootstrap score.

In order to analyze the large sample properties of the non-parametric bootstrap score, it is important to observe that a standard (bootstrap version of the) CLT cannot be applied to (ref) because the number of terms in the sum is itself random. That is, $S_{T}^{\ast} (\hat{\theta}_{T})$ is a randomly selected partial sum. Its behavior can however be analyzed by considering the following FCLT for i.i.d. waiting times, which for non-bootstrap sequences is due to B1968 (the extension to bootstrap random variables is straightforward and is omitted for brevity).

By using the fact that $\hat{\theta}_{T}$ is consistent and that the sample variance of the transformed waiting times converges to one, an immediate application of Theorem (ref) yields that \[ T^{-1/2}S_{T}^{\ast}(\hat{\theta}_{T})\overset{d^{\ast}}{\rightarrow} _{p}\mathcal{N}(0,\theta_{0}^{-1}), \] which matches the asymptotic distribution of the original score. For the Hessian, \[ T^{-1}H_{T\ }^{\ast}(\hat{\theta}_{T})=T^{-1}n_{T}^{\ast}/\hat{\theta}_{T} ^{2}\overset{p^{\ast}}{\rightarrow}_{p}1/\theta_{0}\text{,} \] as $T^{-1}n_{T}^{\ast}=\hat{\theta}_{T}^{\ast}=\theta_{0}+o_{p^{\ast}}\left( 1\right) $, in probability, applying again Theorem (ref). By standard arguments this implies that \[ T^{1/2}(\hat{\theta}_{T}^{\ast}-\theta_{T}^{\ast})\overset{d^{\ast} }{\rightarrow}_{p}\mathcal{N}(0,\theta_{0})\text{,} \] and similarly, $LR_{T}^{\ast}(\theta_{T}^{\ast})\overset{d^{\ast}} {\rightarrow}_{p}\chi_{1}^{2}$. The general (non-Poisson) case is more involved due to the fact that conditionally on the data the bootstrap waiting times in transformed time scale have a discrete distribution. Although this feature is not crucial in the Poisson case, the general case involves the analysis of random terms of the form $\int\xi(t)dN^{\ast}(t)$ and an explicit calculation of the compensator of $N^{\ast}(t)$.

We conclude by noticing that, as shown in the next section, the non-parametric bootstrap performs as well as the parametric bootstrap.

Relation with the bootstrap for ACD models

In this section we discuss the relation between our proposed bootstrap algorithms and theory and extant results on the bootstrap for ACD models; see in particular FG2005, GKS2015, PHS2016 and PS2017, PS2021 for the related class of multiplicative error models [MEM].

Consider, initially, the exponential ACD process [EACD] which, by ((ref)) and ((ref)), has intensity \[ \lambda\left( t\right) =\frac{1}{\psi_{N\left( t-\right) +1}} =\lambda\left( t-t_{N\left( t-\right) },t-t_{N\left( t-\right) -1},....\right) , \] and associated integrated intensity

Note that, to simplify notation, we omit here the dependence on $\theta$ parametrizing the intensity function and hence $\psi_{i}$.

It follows that our proposed RIB algorithms are related to recursive bootstraps in the ACD framework. To see this, recall that for the parametric RIB, we first generate the sequence of transformed waiting times $\left\{ v_{i}^{\ast}\right\} $ as i.i.d. $\mathcal{E}\left( 1\right) $, while for the non-parametric RIB we resample from the original (standardized) transformed waiting times $v_{i},i=1,...,n_{T}$ which, using ((ref)), are given by $v_{i}=w_{i}/\psi_{i}$ in the case where $\theta_{T}^{\ast}=\theta_{0}$ without loss of generality. Next, the bootstrap waiting times $w_{i}^{\ast}$ are generating recursively as

with

which is equivalent to a recursive bootstrap for the EACD model (either parametric or non-parametric). Therefore, the theory we develop in this paper can also be used to establish bootstrap validity for EACD models. A crucial difference between the recursive bootstrap for MEM is that in ((ref)) the number of event times $n_{T} ^{\ast}$ is random, $\sum\nolimits_{i=1}^{n_{T}^{\ast}}w_{i}^{\ast}\leq T$, such that the event times fall within the interval $[0,T]$. In contrast, in the recursive MEM case $n_{T}^{\ast}=n_{T}$, which implies that $\sum \nolimits_{i=1}^{n_{T}}w_{i}^{\ast}$ can be much smaller or even larger than $T$.

For the case of ACD with non-exponentially distributed errors $\varepsilon _{i}$, the intensity is ((ref)) with corresponding integrated intensity

where $S_{\varepsilon}$ is one minus the cdf of $\varepsilon_{i}$. In the parametric case, with $v_{i}^{\ast}$ drawn as $\mathcal{E}(1)$, using ((ref)) in the bootstrap world, we recursively obtain

where

with the $U_{i}^{\ast}$'s being i.i.d. uniform in $[0,1]$; hence, apart for the random stopping time, our RIB covers the recursive bootstrap for ACD.

For the non-parametric case, existing bootstrap algorithms for ACD generate the errors $\varepsilon_{i}^{\ast}$ by resampling the residuals $\varepsilon _{i}=w_{i}/\psi_{i}$, while the RIB first generates the bootstrap waiting times in transformed scale $v_{i}^{\ast}$ by resampling the estimated $v_{i}$; these are later used to generated the bootstrap errors $\varepsilon_{i}^{\ast }$ and then $w_{i}^{\ast}$, see ((ref)) and ((ref)) above.

Finally, consider our (either parametric or non-parametric) FIB applied to the ACD. It would be tempting to think that our FIB would correspond to a `fixed conditional expected duration' bootstrap in the ACD world, where the bootstrap waiting times are generated as

where $\hat{\psi}_{i}$ is the $i$-th estimated conditional expected duration on the original data. Although this algorithm, which resembles the fixed volatility bootstrap for ARCH processes proposed in CPR2018 and has not been investigated previously in the literature, seems to be an interesting development, it does not correspond to our FIB algorithm. In particular, the FIB uses the inverse of the estimated (integrated) intensity, say $\hat{\Lambda}^{-1}$, to transform the bootstrap $v_{i}^{\ast}$ into the bootstrap waiting times $w_{i}^{\ast}$ and generates a number of event times which is random in the bootstrap world; in contrast, a bootstrap based on ((ref)) generates a number of events, given by $n_{T}$, which is fixed in the bootstrap world.

Monte Carlo Simulations

In this section we consider the finite sample properties of asymptotic and bootstrap-based confidence intervals and hypothesis tests for the well-known and much used case of a Hawkes process. By considering a detailed simulation study based on the exponential kernel, we analyze how the bootstrap compares to asymptotic inference for different values of key quantities such as the `branching ratio' (defined below) and the decaying rate of the memory of past events. We consider both the RIB and the proposed FIB schemes, parametric as well as non-parametric.

Model and implementation

In the simulations, we consider the Hawkes process with exponential kernel function, $\gamma(x;\alpha,\beta)=\alpha e^{-\beta x}$ and conditional intensity \[ \lambda(t;\theta)=\mu+\sum_{t_{i}<t}\gamma(t-t_{i};\alpha,\beta), \] with $\theta=(\mu,\alpha,\beta)^{\prime}$, $\mu,\alpha,\beta>0$, see also (ref). Here $\mu$ is the baseline intensity; $\alpha$ is the jump size of the intensity when a new event occurs; $\beta$ is the exponential decaying rate, which determines how fast the memory of past events declines to zero. In terms of $\alpha$ and $\beta$, a key quantity is the branching ratio

which describes how quickly the number of events increases\footnote{More precisely, in the Poisson cluster representation of the self-exciting point process HO1974, the branching ratio defines the expected number of direct offsprings spawned by an `immigrant' event.}. Moreover, with $\mu,\alpha,\beta>0$, stationarity of the Hawkes process requires the branching ratio to satisfy $0<a<1$, in which case the mean intensity $m$ is well defined and given by \[ m:=E(\lambda(t))=\frac{\mu}{1-a}. \] Hence, the stationary region is given by $\{\theta=(\mu,\alpha,\beta)^{\prime }\in\mathbb{R}\times\mathbb{R}\times\mathbb{R}:\mu>0,0<\alpha<\beta\}$. A few remarks about the simulation scheme are as follows.

We simulate three stationary Hawkes processes (denoted by Models 1--3) with true parameters $\theta_{0}$ set as follows. For all simulated processes, the mean intensity is set to unity ($m_{0}=1$), while different levels of the branching ratio $a_{0}=\alpha_{0}/\beta_{0}$ are considered; specifically, we set $a_{0}\in\{0.2,0.5,0.8\}$. For each simulation, we consider three parameterizations (see A--C below) to allow different jump sizes and decaying behavior of the intensity. In all cases, we consider samples over $[0,T]$ for $T\in\{50,100,200\}$ with initial burn-in period $[-M,0)$ for $M=500$. The number of valid Monte Carlo replications (see Remark (ref)(iv)) is $10,000$, and the number of bootstrap repetitions is $B=199$.

The parameter configurations are summarized in Table (ref) along with the (Monte Carlo) probabilities that the SC fails. It can be noticed that the probabilities of SC failure are severely high only for Model 1A when $T=50$. This is because the number of events generated for $T=50$ is extremely volatile and the likelihood of observing samples with a small number of events (hence, not informative enough for estimating the model reasonably well) is indeed high. Another reason is that, as is known, it is hard to precisely estimate the parameters when the true parameters $\alpha_{0}$ and $\beta_{0}$ are close to the zero boundary and $T$ is small. The reparameterization by branching ratio helps to resolve some numerical issues in estimation, as discussed in Remark (ref)(iii) but the improvement is not sufficient when the branching ratio itself is also low as in the case of Model 1A. Nevertheless, despite the quite extreme parameter setting of Model 1A, we decided to keep it in our Monte Carlo simulation for completion.

For each parameter configuration and sample size, we report the coverage probabilities (estimated over the Monte Carlo replications) of confidence intervals at the $95\%$ nominal level, using both asymptotic and bootstrap methods. Asymptotic confidence intervals for the individual parameters as well as the (joint) confidence ellipsoid are based on the sample Hessian. We also report the coverage of (asymptotic and bootstrap) confidence intervals for the branching ratio, $a=\alpha/\beta$. For bootstrap confidence intervals we consider the naive percentile interval method.

Finally, we also report the (null) empirical rejection probabilities of likelihood ratio tests for the hypothesis $H_{0}:\theta=\theta_{0}$. For the bootstrap tests, we implement the unrestricted bootstrap (i.e., without the null imposed on the bootstrap sample); results for the restricted bootstrap (i.e., with the null imposed on the bootstrap sample) do not differ substantially.

Results

The coverage probabilities of the asymptotic and bootstrap confidence intervals [CI] for individual parameters are presented in Table (ref). We can see that, in general, the asymptotic CIs suffer from the problem of undercoverage for almost all models and sample spans, and this fact is particularly severe for some of the cases. In contrast, the bootstrap methods, especially the FIB, powerfully correct these distortions.

Below we provide a summary of the problems related to the asymptotic CIs for each individual parameter (branching ratio $a$, baseline intensity $\mu$, intensity jump size $\alpha$ and decay rate $\beta$).

(i) The undercoverage of the asymptotic CI for the branching ratio is severe in finite sample for all Models 1--3. The coverage deteriorates as the true value of branching ratio increases (moving from Model 1 to 3), and as the true values of $\alpha$ and $\beta$ decrease (moving from Model C to A). Accordingly, the performance of the asymptotic CI for the branching ratio is the worst for Model 3A, where the coverage probability is $86.5\%$ for $T=50$. Larger $\alpha_{0}$ and $\beta_{0}$ seem to improve the coverage rate of the branching ratio, and this improvement is the most significant for Model 1 where the branching ratio is low.

(ii) The asymptotic CI for the baseline intensity $\mu$ performs poorly in finite samples when $\mu_{0}$ is low. Note that for Model 3, where $\mu_{0}=0.2$, the empirical coverage probabilities are $90.5\%$, $89.5\%$, and $88.1\%$ for Model 3A, 3B and 3C, respectively, when $T=50$. In contrast, these probabilities are all above $90\%$ for Models 1 and 2.

(iii) The problem of undercoverage deteriorates when $\alpha_{0}$ is larger (moving from Model A to C). There are no significant changes in the coverage of $\alpha$ over different values of the branching ratio. Improvements in the coverage of $\alpha$ seem to come only from increasing the sample span $T$. In general, the coverage is acceptable.

(iv) The undercoverage of $\beta$ is severe for Model 1 with small branching ratio, but the coverage rate improves noticeably as branching ratio increases, and as sample span $T$ increases. In particular, the asymptotic CI coverage for $\beta$ is almost perfect for Models 3A--C even when $T=50$. The performance is independent of the value of $\beta$.

In contrast to the coverage of asymptotic CIs, which show evident finite sample distortions, the empirical coverage probabilities of the bootstrap percentile intervals based on the fixed intensity scheme (for both the parametric and non-parametric methods, labelled `PRFB' and `NPFB' in Table (ref)) are very close to the nominal level, for almost all simulation models and even when the sample span is very short ($T=50$). The only exceptions are for the coverage of branching ratio in Model 3, where the coverage probabilities of the parametric FIB and the non-parametric FIB are slightly below $95\%$, the nominal level. Nevertheless, the CIs of the two recursive intensity bootstraps, although performing generally better than asymptotic CIs for the coverage of parameter $\alpha$ and $\beta$, share similar features of finite sample distortion as asymptotic CIs. For instance, the coverage of $\mu$ deteriorates as $\mu_{0}$ decreases (for both the parametric and non-parametric RIBs); the coverage of $\beta$ is much below the nominal level for Model 1 where branching ratio is low, while it converges to the nominal level as the branching ratio increases (for the RIB); finally, we observe that the coverage of the branching ratio deteriorates when the branching ratio increases.

Unreported simulations of average lengths of the $95\%$ asymptotic and bootstrap confidence intervals for each parameter show that bootstrap confidence intervals are not significantly wider than asymptotic confidence intervals, except for Models 1A and 1B in which the parameter settings are relatively more extreme, or when the sample span is short ($T=50$). The wider bootstrap confidence intervals reveal the higher uncertainty associated to parameter estimation, and is in line with existing literature on the bootstrap.

Table (ref) presents the joint coverage rate of the asymptotic and bootstrap confidence ellipsoids [CE], for both parameterizations $\theta=(\mu,\alpha,\beta)^{\prime}$ and $\tilde{\theta}=(\mu,a,\beta )^{\prime}$. Noticeably, here the benefit of using bootstrap methods to improve the finite sample joint coverage is way more than evident. The performance of the asymptotic CEs is clearly unsatisfactory: For all nine models, the empirical coverage probabilities of the asymptotic CEs are below $89\%$ when $T=50$; despite a gradual improvement of the coverage rates as sample span $T$ increases, the coverage rates when $T=200$ are still below the nominal level for all models (the joint coverage probabilities of Model 1B are even less than $88\%$ when $T=200$). On the contrary, all bootstrap methods produce the joint CEs that cover the true parameters with probabilities very close to the nominal level,\footnote{We do observe that there is some tendency of over-coverage of the bootstrap joint CEs for Model 2A and 3A for relatively short sample spans.} across different models and different sample spans.

Finally, in Table (ref) we report the empirical rejection probabilities of the asymptotic and unrestricted bootstrap likelihood-ratio tests for the null hypothesis $H_{0}:\theta=\theta_{0}$. In general, both the asymptotic and bootstrap tests perform satisfactorily well in terms of size, especially when $T=100$ and $200$. Nevertheless, we do notice that the asymptotic test tends to be oversized for larger values of the branching ratio. This can be seen by inspecting the rejection probabilities of the asymptotic test on $H_{0}$ for Model 3 (which has the largest branching ratio, $a=0.8$) for $T=50,100$. In particular, the asymptotic test is severely oversized for all three sub-models of Model 3, particularly so when $T=50$. In contrast, we do not see much variability of the bootstrap empirical rejection probabilities across different models or sample spans -- they are all very close to the nominal level (slightly conservative in some cases).

Empirical illustrations

To illustrate how the proposed bootstrap schemes work in applications, we consider two empirical examples. The first consists of `extreme occurrences' in US stock market data, as measured by empirical quantiles of the Dow Jones Index, see ELL2011. We use this application to compare the four different bootstrap schemes discussed in the paper. Next, we analyze recent Danish COVID-19 tweets using the non-parametric FIB. We illustrate how bootstrap confidence intervals reveal the presence of a structural break in the parameters, whereas confidence intervals based on the asymptotic Gaussian approximation do not.

Dow Jones Index

As in ELL2011, we consider Dow Jones Index (DJI) daily (log) returns observed over the period January 1, 1994 to December 31, 2010. The event times corresponding to extreme returns are given by the trading days where the corresponding daily return is below the $10\%$ empirical quantile (negative occurrences), resulting in $n_{T}=428$ events during the period of $T=6144$ days considered. Figure (ref) (top panel) shows the event times and the associated counting process.

To analyze the data, we consider a Hawkes model with intensity reparameterized as

where $a$ is the branching ratio and $\gamma$ is the (exponential) kernel; that is, $\gamma(t;\theta)=\beta\exp(-\beta t)$, see also ((ref)) and Section (ref). With parameter vector $\theta=(\mu,a,\beta)^{\prime}$, the MLE $\hat{\theta}$ is obtained by maximizing the log-likelihood in ((ref)) subject to $\mu,\beta>0$, $0<a<1$ and with initial values from ELL2011. Estimation results are reported in Table (ref); the estimated intensity is portrayed in the bottom panel of Figure (ref). The MLE $\hat{\theta}$ is very similar to ELL2011, and we observe in particular that the branching ratio $a$ appears to be well inside the stationary region.

As previously emphasized, if the model is correctly specified, the transformed waiting times should be i.i.d. $\mathcal{E}\left( 1\right) $. Therefore, the model fit can be evaluated by considering the estimated transformed waiting times \[ \hat{v}_{i}=\Lambda(t_{i},t_{i-1};\hat{\theta}),\quad i=1,2,...,n_{T}, \] with $\Lambda$ defined in ((ref)). Figure (ref) contains QQ-plots and Kolmogorov-Smirnov (KS) plots, as well as sample autocorrelograms and related tests. Based on these, we see no clear signs of model misspecification. Precisely, the QQ plot of $\hat{v}_{i}$ against a unit exponential distribution has no significant deviations from the identity line, except a few quantiles in the extreme upper tail, as also confirmed by the KS statistic p-value ($0.147$). Moreover, while the observed waiting times $w_{i}$ are autocorrelated, this is not the case for the transformed waiting times $\hat{v}_{i}$ (and its squares, $\hat{v}_{i}^{2}$).

We next compare the different bootstrap algorithms in terms of confidence intervals for the parameters, and compare these with the asymptotic CIs. With $\{\hat{\theta}_{T,i:b}^{\ast}\}_{b=1}^{B}$ the i.i.d. bootstrap realizations of the $i$-th element of $\hat{\theta}_{T}^{\ast}$, the bootstrap CIs reported are based on the empirical $\alpha/2$ and $1-\alpha/2$ quantiles of the empirical distribution function of the $\hat{\theta}_{T,i:b}^{\ast}$'s. In Table (ref), while we find no noticeable difference between the parametric and non-parametric bootstraps, the bootstrapped CIs based on the FIB are less wide when compared to the asymptotic and RIB CIs (recall also from the Monte Carlo results that in general the bootstrap coverage probabilities are better than those associated to the asymptotic CIs). The observed difference between the FIB and RIB CIs is likely to be caused by the added randomness in the sequential computation of the RIB. Interestingly, the FIB and RIB bootstrap CIs are further away from the non-stationary region ($a\geq1$) than the asymptotic CIs.

COVID-19 Tweets

We consider the arrival times of tweets related to the COVID-19 pandemic, recorded on March 11 (06:00-00:00) 2020, when during a press briefing the Danish Prime Minister at 20:30 announced the first lockdown of Denmark. In total, there are $n_{T}=1822$ events from $1166$ unique individuals, with each event time $\{t_{i}\}_{i=0}^{n_{T}}$ ($t_{0}=0$) measured with a time resolution of $1$ second within the $T=18$ hours considered. In order to analyze the effects of the announcement, we analyze the full sample, as well as the pre-press briefing sample (06:00--20:30), and the post-press briefing sample (20:30--00:00). In Figure (ref), we show the observed counting process $N(t)$ for $t\in\lbrack0,T]$ as well as an initial proxy for the intensity given by the number of events per $15$-minute intervals. It is worth noticing that there is a surge in activity after 20:30, visible both in the counting process and the increased intensity.

As for the DJI data, we consider the Hawkes model with exponential kernel. Based on the diagnostics (see Figure (ref)), the model seems to be well specified in all the three (sub)samples. However, we observe a large difference between the estimates reported for the first subsample and for the second subsample, see Table (ref). In particular, the effect of the response to the announcement is a substantial increase in the intensity. One may also note that the estimated memory parameter $\beta$ for the full period is between the estimates for the pre-announcement and post-announcement periods. Table (ref) also reports asymptotic CIs and FIB CIs. As can clearly be seen, the bootstrap CIs indicates the presence of non-overlapping parameter estimates for the samples before and after the announcement. This possibly reflects different types of dynamics in the two samples, and indicates a structural break around the press briefing. We note that this is not detectable by the standard misspecification tests for the full sample, and is much less pronounced from the reported asymptotic CIs for the three samples (in particular so for the baseline $\mu$).

In addition, we have also considered the power law kernel, where $\gamma(t;\theta)$ in (ref) is replaced by a power law, see (ref). Interestingly, unreported results show that, in terms of model misspecification, one is unable to discriminate between the two models, and moreover that the estimates of the baseline $\mu$ and branching ratio $a$ are virtually indistinguishable from those obtained using the exponential kernel. Finally, estimation based on the power law kernel (unlike the exponential kernel) is highly sensitive to initial values, which may reflect the large correlation of the parameter estimators for power law kernels.

Conclusions

In this paper we have discussed the theoretical foundations and practical implementations of bootstrap inference for self-exciting point process models. Applications of the bootstrap in order to improve upon the poor quality of asymptotic approximations are scarce in the literature. Classic `recursive intensity bootstrap' (RIB) schemes have been proposed in the recent literature, although without proof of their first-order validity. RIB schemes can also be quite involved to implement in practice, as they generally require numerical integration for the recursive computation of the intensity for each bootstrap repetition. To improve, we have introduced a new bootstrap scheme, the `fixed intensity bootstrap' (FIB), where the conditional intensity is kept fixed across bootstrap repetitions. By doing so, conditionally on the original data the bootstrap data generating process follows a simple inhomogeneous point process with known intensity; therefore, it is very simple to implement and to use in practice. For both bootstrap schemes, we have provided a new bootstrap (asymptotic) theory, which allows to assess bootstrap validity for both bootstraps. Monte Carlo evidence supports the idea that the bootstrap is a valid inference method when applied to point process models.

The results in the paper could be extended in several directions. On top of the obvious extension to multivariate point process models, an interesting one is how to deal with marked point process models. Marked (self-exciting) processes are particularly useful in applications, as the intensity function can be made dependent on a set of `marks' associated to past events (for financial returns, the trading volumes; for energy prices, the magnitude of price spikes; for tweets, the number of followers; for earthquakes modelling, the magnitude of the earthquakes). In this context the proposed FIB seems very powerful as re-sampling with a fixed intensity, even as a function of marks, is feasible and easy to implement. As an example, consider briefly an extension of the Hawkes model with exponential kernel in (ref). One may include real-valued marks, or covariates, $y_{t}\in\mathbb{R}^{d}$ in the conditional intensity $\lambda(t;\theta)$ as for example, \[ \lambda\left( t;\theta\right) =\mu\left( y_{t}\right) +\alpha\sum _{t_{i}<t}\beta\left( y_{t_{i}}\right) \gamma\left( t-t_{i}\right) , \] where $\mu,\beta:\mathbb{R}^{d}\mathbb{\rightarrow R}_{+}$, see e.g. CHH2015 for an application to price spikes in electricity markets. Under the assumption of `strongly exogenous' (or, ancillary) marks, similar to exogenous covariates in discrete time Poisson autoregressions ACKR2016 and with $\theta$ the parameters parameterizing the extended Hawkes intensity, estimation and inference based on the FIB utilize the original event times and marks, $\{t_{i},y_{t_{i}}\}_{i=1}^{N(T)}$. Thus, in contrast to the RIB and other existing recursive bootstraps, bootstrap inference based on FIB would not require further assumptions (apart from stationarity) of the covariates.

A further extension is to develop model misspecification-robust bootstrap methods. In particular, throughout the paper we have assumed that the model is correctly specified. This assumption implies that the bootstrap can be implemented parametrically by constructing bootstrap waiting times from an i.i.d. sequence of mean one exponential random variables (the waiting times in transformed time scale), as discussed in Sections (ref) and (ref). However, misspecification of the model (in the simplest case, data are modelled as a Poisson process, but the waiting times form a renewal process) may result in i.i.d., but non-exponential (transformed) waiting times. Although in this case the parametric bootstraps could fail, we believe that the non-parametric bootstrap algorithms discussed in Section (ref) could serve as the basis of novel misspecification-robust bootstrap methods. All these extensions are left for future research.

Acknowledgements

We are grateful to Torben Andersen (Co-Editor) and two anonymous referees for many constructive comments and suggestions on an earlier version of the paper. We have also benefited from discussions and feedback from seminar participants at Saint-Petersburg State University (CEBA talks), Singapore Management University, Macquarie University, as well as participants of the 9th Italian Congress of Econometrics and Empirical Economics (University of Cagliari) and the 2021 Virtual Workshop on Financial Econometrics (Durham University).

This research was supported by the Danish Council for Independent Research (DSF Grant 015-00028B), the Center for Information and Bubble Studies, University of Copenhagen, the Italian Ministry of University and Research (PRIN 2017 Grant 2017TA7TYC) and the University of Sydney (Faculty Research Future Fix 2020 Grant). Part of this paper was written while Giuseppe Cavaliere was visiting the School of Economics of the University of Sydney; financial support and hospitality are gratefully acknowledged. Finally, the authors acknowledge the technical assistance provided by the Sydney Informatics Hub of the University of Sydney for the high-performance computing and cloud services.