Decompounding discrete distributions: A non-parametric Bayesian approach

TL;DR

This paper introduces a non-parametric Bayesian method for estimating the jump distribution and intensity of a discretely observed compound Poisson process, with theoretical guarantees and practical applications.

Contribution

It develops a novel Bayesian approach with MCMC sampling for decompounding discrete distributions, including theoretical contraction results.

Findings

The method performs well on simulated data.

It outperforms the frequentist plug-in estimator in experiments.

Posterior contraction rate is established as 1/√n up to a log factor.

Abstract

Suppose that a compound Poisson process is observed discretely in time and assume that its jump distribution is supported on the set of natural numbers. In this paper we propose a non-parametric Bayesian approach to estimate the intensity of the underlying Poisson process and the distribution of the jumps. We provide a MCMC scheme for obtaining samples from the posterior. We apply our method on both simulated and real data examples, and compare its performance with the frequentist plug-in estimator proposed by Buchmann and Gr\"ubel. On a theoretical side, we study the posterior from the frequentist point of view and prove that as the sample size $n\rightarrow\infty$, it contracts around the `true', data-generating parameters at rate $1/\sqrt{n}$, up to a $\log n$ factor.

Figures

Tables

Table 1. Table 1. Results for scenarios (a)–(c) from Section  3.1 .

Simulation setting	(a)	(b)	(c)
Buchmann-Grübel estimator	$1.40$	$0.32$	$1.44$
Posterior mean	$0.15$	$0.07$	$0.12$

Table 2. Table 2. Results for scenarios (a)–(b) from Section  3.2 .

Simulation setting	(a)	(b)
Buchmann-Grübel estimator	$0.28$	$0.60$
Posterior mean	$0.52$	$1.05$

Table 3. Table 3. Data on the number of soldiers in the Prussian cavalry killed by horse kicks. See the table on p. 25 in von Bortkewitsch ( 1898 ) .

Deaths	$0$	$1$	$2$	$3$	$4$
Counts	$109$	$65$	$22$	$3$	$1$

Table 4. Table 4. Plant population data from Evans ( 1953 ) .

Plants	$0$	$1$	$2$	$3$	$4$	$5$	$6$	$7$	$8$	$9$	$10$	$11$	$12$
Counts	$274$	$71$	$58$	$36$	$20$	$12$	$10$	$7$	$6$	$3$	$0$	$2$	$1$

Equations

$$X_t=\sum _{j=1}^{N_t}{Y}_j,$$

$$\{{\nu }_k{\}}_{k=1}^{\infty }\mid a,m,{\beta }_k$$

$${\beta }_1,\dots ,{\beta }_m\mid \gamma$$

$$\gamma$$

$$q=e^{-\lambda }\sum _{j=0}^{\infty }\frac{{\lambda }^j}{j!}{p}^{\ast j},$$

$$q_0=e^{-\lambda },q_k=\frac{\lambda }{k}\sum _{j=1}^{k}j{p}_j{q}_{k-j},k\in \mathbb{N}.$$

$$\lambda =-\log q_0,p_k=-\frac{q_k}{q_0\log q_0}-\frac{1}{k{q}_0}\sum _{j=1}^{k-1}j{p}_j{q}_{k-j},k\in \mathbb{N}.$$

$$e^{-T{\sum }_{k=1}^m{\nu }_k}\prod _{k=1}^m{\nu }_k^{{\mu }_k},$$

$${\mu }_k=\#\{Y_j=k\},$$

$${\nu }_k\mid \{{\mu }_k\},\{{\beta }_k\}\stackrel{\text{i.i.d.}}{\sim }\mathrm{Gamma}(a+{\mu }_k,1/{\beta }_k+T),k=1,\dots ,m,$$

$${\beta }_k\mid \{{\nu }_k\},\gamma \stackrel{\text{i.i.d.}}{\sim }\mathrm{IG}(a+c,\gamma +{\nu }_k),k=1,\dots ,m,$$

$$\gamma \mid \{{\beta }_k\}\sim \mathrm{Gamma}\left(cm+1,1+\sum _{k=1}^m{\beta }_k^{-1}\right).$$

$$Z_i=\sum _{j=1}^{m}j{\mu }_{ij},$$

$$\operatorname{Pr}(\mu_{i1}=k_{1},\ldots,\mu_{im}=k_{m}\mid Z_{i}=z_{i})\\ =\frac{1}{\operatorname{Pr}(Z_{i}=z_{i})}e^{-\Delta_{i}\sum_{j=1}^{m}\nu_{j}}\prod_{j=1}^{m}\frac{(\Delta_{i}\nu_{j})^{k_{j}}}{k_{j}!}\mathbf{1}\Big{\{}\sum_{j=1}^{m}jk_{j}=z_{i}\Big{\}}.$$

$$\log A=\sum _{k=1}^m({\mu }_{ik}^{\circ }-{\mu }_{ik})\log ({\Delta }_i{\nu }_k)+\sum _{k=1}^m\left\{\log ({\mu }_{ik}!)-\log ({\mu }_{ik}^{\circ }!)\right\}.$$

$$p_k=-\frac{(1-\alpha {)}^k}{k\log \alpha },k\in \mathbb{N}.$$

$${\Pi }_n\left(\Vert \nu -{\nu }_0{\Vert }_1\ge \frac{{\log }^{\gamma }n}{\sqrt{n}}|{\mathcal{Z}}_n\right)\to 0$$

$$\Vert a\ast b{\Vert }_1\le \Vert a{\Vert }_1\Vert b{\Vert }_1,$$

$$\exp (a)=\sum _{j=0}^{\infty }\frac{a^{\ast j}}{j!}.$$

$$\exp (a+b)=\exp (a)\ast \exp (b),a,b\in {\ell }_1,$$

$$\exp (a)=\exp (b)⟹a=b,a,b\in {\ell }_1.$$

$$r_0=\frac{1}{q_0},r_k=-\frac{1}{q_0}\sum _{j=1}^k{q}_j{r}_{k-j},k\in \mathbb{N}.$$

$$\Vert {\nu }^{\prime }-\nu {\Vert }_1=\Vert {\lambda }^{\prime }{p}^{\prime }-\lambda p{\Vert }_1\le \frac{\Vert r{\Vert }_1\Vert q^{\prime }-q{\Vert }_1}{1-\Vert r{\Vert }_1\Vert q^{\prime }-q{\Vert }_1}.$$

$$({\lambda }^{\prime }-\lambda ){\delta }_0+\lambda p-{\lambda }^{\prime }{p}^{\prime }=\sum _{j=1}^{\infty }\frac{1}{j}(r\ast (q-q^{\prime }){)}^{\ast j}$$

$$|{\lambda }^{\prime }-\lambda |+\Vert {\lambda }^{\prime }{p}^{\prime }-\lambda p{\Vert }_1\le \frac{\Vert r{\Vert }_1\Vert q^{\prime }-q{\Vert }_1}{1-\Vert r{\Vert }_1\Vert q^{\prime }-q{\Vert }_1},$$

$$\Pi \left(\Vert \nu -{\nu }_0{\Vert }_1\ge \epsilon |{\mathcal{Z}}_n\right)\le 2\Pi \left(\Vert q-q_0{\Vert }_1\ge \frac{\epsilon }{2\Vert r_0{\Vert }_1}|{\mathcal{Z}}_n\right).$$

$$\left\{\nu :\Vert \nu -{\nu }_0{\Vert }_1\ge \epsilon \right\}\cap \left\{\nu :\Vert r_0{\Vert }_1\Vert q-q_0{\Vert }_1<1/2\right\}$$

$$\left\{\nu :\Vert \nu -{\nu }_0{\Vert }_1\ge \epsilon \right\}\cap \left\{\nu :\Vert r_0{\Vert }_1\Vert q-q_0{\Vert }_1\ge 1/2\right\}.$$

$$\left\{\nu :\Vert \nu -{\nu }_0{\Vert }_1\ge \epsilon \right\}\cap \left\{\nu :\Vert r_0{\Vert }_1\Vert q-q_0{\Vert }_1<1/2\right\}$$

$$h(q_{\lambda ,p},q_{{\lambda }^{\prime },p^{\prime }})\le \sqrt{\lambda }h(p,p^{\prime })+|\sqrt{\lambda }-\sqrt{{\lambda }^{\prime }}|,$$

$$KL(q_0,q)=Q_0\left(\log \frac{q_0}{q}\right),V(q_0,q)=Q_0{\left(\log \frac{q_0}{q}\right)}^2.$$

$$\begin{array}{r}KL(q_0,q)\le \overline{C}\left(KL(p_0,p)+|{\lambda }_0-\lambda {|}^2\right),\\ V(q_0,q)\le \overline{C}\left(V(p_0,p)+KL(p_0,p)+|{\lambda }_0-\lambda {|}^2\right);\end{array}$$

Full text

The research leading to the results in this paper has received funding from the European Research Council under ERC Grant Agreement 320637, from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – 314838170, GRK 2297 MathCoRe, and from the Deutsche Forschungsgemeinschaft (DFG) through the grant CRC 1294 ‘Data Assimilation’

1. Introduction

1.1. Problem formulation

Let $N=(N_{t}:t\geq 0)$ be a Poisson process with a constant intensity $\lambda>0$, and let $Y_{i}$ be a sequence of independent random variables, each with distribution $P$, that are also independent of $N$. By definition, a compound Poisson process (abbreviated CPP) $X=(X_{t}:t\geq 0)$ is

[TABLE]

where here and below the sum over an empty index set is understood to be equal to zero. In particular, $X_{0}=0$. CPP constitutes a classical model in, e.g., risk theory, see Embrechts et al. (1997).

Assume that the process $X$ is observed at discrete times $0<t_{1}<t_{2}<\ldots<t_{n}=T$, where the instants $t_{i}$ are not necessarily equidistant on $[0,T]$. Based on the observations $X_{t_{1}},X_{t_{2}},\ldots,X_{t_{n}}$, our goal is to estimate the jump size distribution $P$ and the intensity $\lambda$. We specifically restrict our attention to the case where $P$ is a discrete distribution, $P(\mathbb{N})=1$, and we will write $p=(p_{k})_{k\in\mathbb{N}}$ for the probability mass function corresponding to $P$, where $p_{k}=P(\{k\})$. A similar notation will be used for any other discrete law. The distribution $P$ is called the base distribution. Abusing terminology, we will at times identify it with the corresponding probability mass function $p.$ An assumption that $P$ has no atom at zero is made for identifiability: otherwise this atom gets confounded with $e^{-\lambda}$, which does not allow consistent estimation of the intensity $\lambda$. For a discussion of applications of this CPP model in risk theory, see Zhang et al. (2014).

Define the increments $Z_{i}=X_{t_{i}}-X_{t_{i-1}}$, $i=1,\ldots,n$. Then $\mathcal{Z}_{n}=(Z_{i}:i=1,\ldots,n)$ is a sequence of independent random variables. When $\{t_{i}\}$ are equidistant on $[0,T]$, the random variables $Z_{i}$ have in fact a common distribution $Q$ satisfying $Q(\mathbb{N}_{0})=1$. As $\mathcal{Z}_{n}$ carries as much information as $(X_{t_{i}}:i=1,\dots,n)$ does, we can base our estimation procedure directly on the increments $\mathcal{Z}_{n}$. Since summing up the jumps $Y_{j}$’s amounts to compounding their distributions, the inverse problem of recovering $P$ and $\lambda$ from $Z_{i}$ can be referred to as decompounding; see Buchmann & Grübel (2003).

There are two natural ways to parametrise the CPP model: either in terms of the pair $(\lambda,p)$, or in terms of the Lévy measure $\nu=(\nu_{k})_{k\in\mathbb{N}}$ of the process $X$, see Sato (2013). A relationship between the two is $\lambda=\sum_{k=1}^{\infty}\nu_{k}$ and $p=\nu/\lambda$. Inferential conclusions in one parametrisation can be easily translated into inferential conclusions into another parametrisation. However, for our specific statistical approach the Lévy measure parametrisation turns out to be more advantageous from the computational point of view.

1.2. Approach and results

In this paper, we take a non-parametric Bayesian approach to estimation of the Lévy measure $\nu$ of $X.$ See Ghosal & van der Vaart (2017) and Müller et al. (2015) for modern expositions of Bayesian non-parametrics. A case for non-parametric Bayesian methods has already been made elsewhere in the literature, and will not be repeated here. On the practical side, we implement our procedure via the Gibbs sampler and data augmentation, and show that it performs well under various simulation setups. On the theoretical side, we establish its consistency and derive the corresponding posterior contraction rate, which can be thought of as an analogue of a convergence rate of a frequentist estimator (see Ghosal & van der Vaart (2017)). The posterior contraction rate, up to a practically insignificant $\log n$ factor, turns out to be $1/\sqrt{n}$, which is an optimal rate for non-parametric estimation of cumulative distribution functions. Our contribution thus nicely bridges practical and theoretical aspects of Bayesian non-parametrics.

1.3. Related literature

To provide a better motivation for our model and approach, in this subsection we briefly survey the existing literature. A Bayesian approach to non-parametric inference for Lévy processes is a very recent and emerging topic, with references limited at the moment to Belomestny et al. (2019), Gugushvili et al. (2015), Gugushvili et al. (2018) and Nickl & Söhl (2017). These deal exclusively with the case when the Lévy measure is absolutely continuous with respect to the Lebesgue density. At least from the computational point of view, these works are of no help in our present context.

Related frequentist papers for CPP models with discrete base distributions are Buchmann & Grübel (2003) and Buchmann & Grübel (2004), which, after earlier contributions dating from the previous century, in fact revived interest in non-parametric techniques for Lévy processes. To estimate the base distribution $p$, Buchmann & Grübel (2003) employ a frequentist plug-in approach relying on the Panjer recursion (i.e., an empirical cumulative distribution estimate of $q$ is plugged into the Panjer recursion equations to yield an estimate of $p$; see below on the Panjer recursion). The drawback is that the parameter estimates are not guaranteed to be non-negative. Buchmann & Grübel (2004) fix this problem by truncation and renormalisation. This works, but looks artificial. As noted in Buchmann & Grübel (2004), in practice the latter approach breaks down if no zero values are observed among $Z_{i}$’s. Buchmann & Grübel (2004) establish weak convergence of their modified estimator, but on the downside its asymptotic distribution is unwieldy to give confidence statements on $p$. Most importantly, the plug-in approaches in Buchmann & Grübel (2003) and Buchmann & Grübel (2004) do not allow obvious generalisations to non-equidistant observation times $\{t_{i}\}.$ In Lindo et al. (2018), another frequentist estimator of the jump measure is introduced, that is obtained via the steepest descent technique as a solution to an optimisation problem over the cone of positive measures. The emphasis in Lindo et al. (2018) is on numerical aspects; again, no obvious generalisation to the case of non-equidistant $\{t_{i}\}$ is available.

Finally, some important, predominantly theoretical references on inference for Lévy processes are Comte & Genon-Catalot (2011), Duval & Hoffmann (2011), Kappus (2014), Neumann & Reiß (2009), Nickl & Reiß (2012), van Es et al. (2007) and Trabs (2015). We refer to Belomestny et al. (2015), Coca (2018), Coca (2018) and Duval & Mariucci (2017) for more extensive literature surveys.

1.4. Outline

The rest of the paper is organised as follows: in Section 2 we introduce our approach and describe an algorithm for drawing from the posterior distribution. In Sections 3 and 4 we study its performance on synthetic and real examples. Section 5 is devoted to the examination of asymptotic frequentist properties of our procedure. An outlook on our results is given in Section 6. Finally, in Appendix A technical lemmas used in the proofs of Section 5 are collected, whereas Appendix B contains some additional simulation results.

1.5. Notation

For two sequences $\{a_{n}\}$ and $\{b_{n}\}$ of positive real numbers, the notation $a_{n}\lesssim b_{n}$ (or $b_{n}\gtrsim a_{n}$) means that there exists a constant $C>0$ that is independent of $n$ and such that $a_{n}\leq Cb_{n}.$ We write $a_{n}\asymp b_{n}$ if both $a_{n}\lesssim b_{n}$ and $a_{n}\gtrsim b_{n}$ hold. We denote a prior (possibly depending on the sample size $n$) by $\Pi_{n}$. The corresponding posterior measure is denoted by $\Pi_{n}(\cdot\mid\mathcal{Z}_{n}).$ The Gamma distribution with shape parameter $a$ and rate parameter $b$ ($a,b>0$) is denoted by $\operatorname{Gamma}(a,b)$. Its density is given by $x\mapsto\frac{b^{a}}{\Gamma(a)}x^{a-1}e^{-bx},\quad x>0$, where $\Gamma$ is the Gamma function. The inverse Gamma distribution with shape parameter $a$ and scale parameter $b$ is denoted by $\operatorname{IG}(a,b)$. Its density is $x\mapsto\frac{b^{a}}{\Gamma(a)}x^{-a-1}e^{-b/x},\quad x>0$. We use the notation $\operatorname{Exp}(a)$ for an exponential distribution with mean $1/a$. Finally, given a metric $d$ on a set $\mathcal{Q}$ and $\epsilon>0$, the covering number $N(\epsilon,\mathcal{Q},d)$ is defined as the minimal number of balls of radius $\epsilon$ needed to cover $\mathcal{Q}$.

2. Algorithm for drawing from the posterior

A Bayesian statistical approach relies on the combination of the likelihood and the prior on the parameter of interest through Bayes’ formula. We start with specifying the prior. As far as the likelihood is concerned, although explicit, it is intractable from the computational point of view for non-parametric inference in CPP models. We will circumvent the latter problem by means of data augmentation, as detailed below.

2.1. Prior

We define a prior $\Pi$ on $\nu$ through a hierarchical specification

[TABLE]

Note that the (fixed) hyperparameters $m\in\mathbb{N}$, $a,c>0$ are denoted by Latin letters.

The hyperparameter $m$ incorporates our a priori opinion on the support of the Lévy measure $\nu$, or equivalently, the base measure $p$. In applications, the support of $p$ may be unknown, which necessitates the use of a large $m$, e.g. $m=\max_{i=1,\ldots,n}Z_{i}$; this latter is the maximal value suggested by the data $\mathcal{Z}_{n}$ at hand. Nevertheless, we may simultaneously expect that the ‘true’, data-generating $\nu$ charges full mass only to a proper, perhaps even a small subset of the set $\{1,\ldots,m\}$. In other words, $\nu$ may form a sparse sequence, with many components equal to zero. In fact, there are at least two plausible explanations for an occurrence of a large increment $Z_{i}$ in the data: either a few large jumps $Y_{j}$’s occurred, which points towards a large right endpoint of the support of $\nu_{0}$, or $Z_{i}$ is predominantly formed of many small jumps, which in turn indicates that the intensity $\lambda$ of the Poisson arrival process $N$ may be large. To achieve accurate estimation results, a prior should take a possible sparsity of $\nu$ into account. This is precisely the reason of our hierarchical definition of the prior $\Pi$: a small $\beta_{k}$ encourages a priori the shrinkage of the components $\nu_{k}$ of $\nu$ towards zero.

2.2. Data augmentation

Assume temporarily $t_{i}=i,$ $i=1,\ldots,n$, and write $q=(q_{k})_{k\in\mathbb{N}_{0}}$ for $q_{k}=q(\{k\})$. Then $Z_{i}$ have the distribution

[TABLE]

with $*$ denoting convolution. The compounding mapping $(\lambda,p)\mapsto q$ can be expressed explicitly via the Panjer recursion (see Panjer (1981)):

[TABLE]

This recursion can be inverted to give the inverse mapping $q\mapsto(\lambda,p)$ via

[TABLE]

In view of (2), the likelihood in the CPP model is explicit. Nevertheless, an attempt to directly use (2) or the Panjer recursions in posterior computations results in a numerically intractable procedure. Equally important is the fact that a Panjer recursion based approach would not apply to non-equidistant observation times $\{t_{i}\}$. Therefore, instead of (2) and the Panjer recursion, we will employ data augmentation, see Tanner & Wong (1987). We switch back to the case when $\{t_{i}\}$ are not necessarily uniformly spaced. The details of our procedure are as follows: when the process $X$ is observed continuously over the time interval $[0,T]$, so that our observations are a full sample path $X^{(T)}=(X_{t}:t\in[0,T])$ of CPP, the likelihood is tractable and is proportional to

[TABLE]

see Shreve (2004), p. 498. Here

[TABLE]

i.e. the total number of jumps of size $k$. Then the prior $\Pi$ from Subsection 2.1 leads to conjugate posterior computations. In fact, the full conditionals are

[TABLE]

Therefore, the Gibbs sampler for posterior inference on $\nu$ can be implemented. The Gibbs sampler cycles through the above conditionals a large number of times, generating approximate (dependent) samples from the posterior. See, e.g., Gelfand & Smith (1990) and Section 24.5 in Wasserman (2004) on the Gibbs sampler and its use in Bayesian statistics.

As we do not observe the process $X$ continuously, we will combine the above with the data augmentation device. First note that we have

[TABLE]

where $(\mu_{ij}:i=1,\ldots,n,j=1,\ldots,m)$ are independent, and $\mu_{ij}\sim\operatorname{Poisson}(\Delta_{i}\nu_{j})$ for $\nu_{j}=\lambda p_{j}$ and $\Delta_{i}=t_{i}-t_{i-1}$; see Corollary 11.3.4 in Shreve (2004). Furthermore, for $\mu_{k}$ as in (3) we trivially have $\mu_{k}=\sum_{i=1}^{n}\mu_{ik}$. Data augmentation iterates the following two steps:

1. (i)

Draw $(\mu_{ij})$ conditional on the data $\mathcal{Z}_{n}$ and the parameter $\nu$. 2. (ii)

Draw $\nu$ conditional on $(\mu_{ij})$.

Once the algorithm has been run long enough, this gives approximate (dependent) samples from the posterior of $\nu$. We already know how to deal with step (ii); now we need to handle step (i).

Thus, keeping $\nu$ fixed, for each $i$ we want to compute the conditional distribution $(\mu_{ij}:j=1,\ldots,m)\mid Z_{i}$, and furthermore, we want to be able to simulate from this distribution. In turn, this will immediately allow us to simulate $\mu_{k}$ conditional on the data $\mathcal{Z}_{n}$. Now, with $\operatorname{Pr}(\cdot)$ referring to probability under the parameter $\nu$, it holds that

[TABLE]

Knowledge of the normalising constant $\operatorname{Pr}(Z_{i}=z_{i})$ will not be needed in our approach.

In general, simulation from a discrete multivariate distribution is non-trivial; some general options are discussed in Devroye (1986), Chapter XI, Section 1.5, but are unlikely to work easily for a large $m$. We will take an alternative route and use the Metropolis-Hastings algorithm, see, e.g., Section 24.4 in Wasserman (2004). We start by observing that for a fixed $i$, the support of $\operatorname{Pr}(\cdot\mid Z_{i}=z_{i})$ is precisely the set $\mathcal{S}_{i}$ of non-negative solutions $(k_{1},\ldots,k_{m})$ of the Diophantine equation $\sum_{j=1}^{m}jk_{j}=z_{i}.$ The R package nilde (see Pya Arnqvist et al. (2018)) implements an algorithm from Voinov & Nikulin (1997) that finds all such solutions for given integers $m$ and $z_{i}$. By Markovianity of the process $X$, we can simulate the vectors $(\mu_{i1},\ldots,\mu_{im})$ independently for each $i=1,\ldots,n$. If $z_{i}=0$ or $1$, there is only one solution to the Diophantine equation: the trivial solution $(0,\ldots,0)$ in the first case, and the solution $(1,0,\ldots,0)$ in the second case; for such $z_{i}$, no simulation is required, as $(\mu_{i1},\ldots,\mu_{im})$ is known explicitly. We thus only need to consider each $i\in\mathcal{I}=\{i:z_{i}\neq 0\textrm{ or }1\}$ in turn, and design a Metropolis-Hastings move on the set of the corresponding solutions $\mathcal{S}_{i}$. Fix once and for all an ordering of elements in $\mathcal{S}_{i}$ (this could be, e.g., lexicographic, or reverse lexicographic); we use the notation $|\mathcal{S}_{i}|$ for the cardinality of $\mathcal{S}_{i}$. Let $\mu=(\mu_{i1},\ldots,\mu_{im})$ be the current state of the chain, corresponding to the $\ell$th element $s_{\ell}$ of $\mathcal{S}_{i}$. A proposal $\mu^{\circ}=(\mu_{i1}^{\circ},\ldots,\mu_{im}^{\circ})$ is obtained as follows:

1. (i)

If $\ell=1$, draw $\mu^{\circ}$ uniformly at random among the elements $\{s_{2},s_{|\mathcal{S}_{i}|}\}$ of $\mathcal{S}_{i}$. 2. (ii)

If $\ell=|\mathcal{S}_{i}|$, draw $\mu^{\circ}$ uniformly at random among the elements $\{s_{1},s_{|\mathcal{S}_{i}|-1}\}$ of $\mathcal{S}_{i}$. 3. (iii)

If $\ell\neq 1$ or $|\mathcal{S}_{i}|$, draw $\mu^{\circ}$ uniformly at random among the elements $\{s_{\ell-1},s_{\ell+1}\}$ of $\mathcal{S}_{i}$.

Occasionally, one may want to propose another type of a move too.

1. (iv)

Draw $\mu^{\circ}=(\mu_{i1}^{\circ},\ldots,\mu_{im}^{\circ})$ uniformly at random from $\mathcal{S}_{i}$.

The two proposals lead to reversible moves, and one may also alternate them with probabilities $\pi$ and $1-\pi$, e.g. $\pi=0.8$. The logarithm of the acceptance probability of a move from $(\mu_{i1},\ldots,\mu_{im})$ to $(\mu_{i1}^{\circ},\ldots,\mu_{im}^{\circ})$ is computed as

[TABLE]

The move is accepted if $\log U\leq\log A$ for $U$ an independently generated uniform random variate on $[0,1]$, and in that case the current state of the chain is reset to $(\mu_{i1}^{\circ},\ldots,\mu_{im}^{\circ})$. Otherwise the chain stays in $(\mu_{i1},\ldots,\mu_{im})$.

3. Simulation examples

In this section, we test performance of our approach in a range of representative simulation examples. For benchmarking, we use the frequentist plug-in estimator from Buchmann & Grübel (2004). Two real data examples are given in Section 4. Unless otherwise stated, we took $c=2$ and $a=0.01$ as hyperparameters in our prior specification. As can be seen from the update formulae for the Gibbs sampler, as long as $a$ is not taken too large, its precise value is not very influential on the posterior, given a reasonable sample size. The value $c=2$ ensures that the update step for $\beta_{k}$ has finite variance. At each step of updating the imputed data for increment size $z$ we have chosen with probability $0.2$ to propose uniformly from all solutions to the Diophantine equation (for that particular value of $z$).

We implemented our procedure in Julia, see Bezanson et al. (2017). The code and datasets for replication of our examples are available on GitHub111See https://github.com/fmeulen/Bdd and Zenodo, see Gugushvili et al. (2019).

3.1. Uniform base distribution

This simulation example follows with some extensions that in Buchmann & Grübel (2004). Let $\lambda_{0}=2$, and let $P_{0}$ be the discrete uniform distribution on $\{1,4,6\}$. We simulated data according to the following settings:

• (a)

$n=100$, $\Delta_{i}=1$ for $1\leq i\leq n$;

• (b)

$n=500$, $\Delta_{i}=1$ for $1\leq i\leq n$ (the data under (a) are augmented with $400$ extra observations);

• (c)

$n=500$, $\Delta_{i}=\operatorname{Unif}(0,2)$ for $1\leq i\leq n$.

We set $m=\min(15,Z_{(n)})$, where $Z_{(n)}=\max_{1\leq i\leq n}Z_{i}.$ In all cases this led to $m=15$, as the value of $Z_{(n)}$ was equal to $30$, $35$ and $40$ for the simulated data under settings (a), (b) and (c), respectively. The Gibbs sampler was run for $500{,}000$ iterations, of which the first $250{,}000$ were discarded as burn-in. From the remaining samples, the posterior mean and $2.5\%$ and $97.5\%$ percentiles were computed for each coefficient $\nu_{k}$. The results for the first $10$ coefficients are shown in Figure 1. For comparison, the estimator from Buchmann & Grübel (2004) is also included in the figure.

For setting (b), traceplots of every $50$th iteration for a couple of coefficients $\nu_{k}$ are shown in Figure 2.

We measure the error of an estimate $\{\hat{\nu}_{k}\}$ by $\operatorname{Err}(\nu,\hat{\nu})=\sum_{k=1}^{\infty}|\hat{\nu}_{k}-\nu_{k}|$. The errors are reported in Table 1. In all settings, for these particular realisations of the simulated data, the Bayesian procedure outperformed the truncated estimator from Buchmann & Grübel (2004). For setting (c), the latter produces a poor result, as was to be expected, given that it is derived under the assumption $\Delta_{i}=1$ for all $i$. An advantage of the Bayesian procedure is the included measure of uncertainty, namely the credible intervals for $\nu_{k}$. On the other hand, for the Buchmann-Grübel estimator it is hardly possible to derive confidence intervals via an asymptotic method, since the limiting distribution of the estimator is fairly complicated. Although not considered in the original publications Buchmann & Grübel (2003) and Buchmann & Grübel (2004), a natural alternative is the bootstrap. A detailed examination of the performance of the latter and its comparison to that of the Bayesian method lies beyond the scope of the present paper. Indeed, any thorough study would require, on one hand, the asymptotic justification of bootstrap confidence intervals, and on another hand establishing frequentist coverage properties of our Bayesian procedure. In that respect, good performance of neither method is automatically warranted (see, e.g., van der Pas et al. (2017) and van der Vaart (1998), Chapter 23). Here instead we opt for a numerical illustration, which is reported in Appendix B.

3.2. Geometric base distribution

The setup of this synthetic data example likewise follows that in Buchmann & Grübel (2004). Assume $q$ is a geometric distribution with parameter $\alpha$, i.e. $q_{k}=(1-\alpha)^{k}\alpha$ for $0<\alpha<1$, $k\in\mathbb{N}_{0}$. Then $\lambda=-\log\alpha$, and

[TABLE]

Hence, $\nu_{k}=(1-\alpha)^{k}/k$.

We consider two simulation setups:

• (a)

$n=500$, $\Delta_{i}=1$ for $1\leq i\leq n$ and $\alpha=1/3$;

• (b)

$n=500$, $\Delta_{i}=1$ for $1\leq i\leq n$ and $\alpha=1/6$.

We set $m=\min(15,Z_{(n)})$ and ran the sampler according to the settings of Section 3.1. The results for both scenarios (a) and (b) are reported in Figure 3. In Table 2 we also report estimation errors in one simulation run. For this example and these generated data, the Bayesian procedure gives less precise point estimates than the Buchmann-Grübel method. Note that estimation error for $\alpha=1/3$ is smaller than that for $\alpha=1/6.$ This appears intuitive, as a smaller value of $\alpha$ corresponds to a larger value of $\lambda.$ The latter implies that on average each $Z_{i}$ is a superposition of a larger number of jumps, which renders the decompounding problem more difficult. However, this argument is hard to formalise.

3.3. Monte Carlo study

For a more thorough comparison of the Buchmann-Grübel estimator and our Bayesian method, we performed a small Monte Carlo experiment. We considered two settings:

1. (i)

The setting from Section 3.1 with $n=250$. We took $m=\min(15,Z_{(n)})$. 2. (ii)

The setting from Section 3.2 with $\alpha=1/3$. We took $m=\min(20,Z_{(n)})$.

In both cases we assumed $\Delta_{i}=1$ for all $1\leq i\leq n$. The number of Monte Carlo repetitions was taken equal to $50$. We took $400{,}000$ MCMC iterations and discarded the first half of these as burn-in. In Figure 4 we give a graphical display of the results by means of boxplots of the errors. Here, as before, if the true values are denoted by $\nu_{k}$ and the estimate within a particular simulation run by $\hat{\nu}_{k}$, the error is defined by $\operatorname{Err}(\nu,\hat{\nu})=\sum_{k=1}^{\infty}|\hat{\nu}_{k}-\nu_{k}|$ (we truncated the infinite summation to $50$). The results agree with our earlier findings, in that there is no clear “winner” in the comparison. Note that for the setting (ii) we considered both $c=2$ and $c=0.01$ in the prior specification. Both values give similar performance of the Bayesian method. This provides insight into sensitivity of our results with respect to the prior specification. A minor difference between the middle and righmost panel of Figure 4 may be attributed to Monte Carlo error: the $50$ simulated datasets on which these panels are based are not the same. Note that the prior promotes sparsity, and in that respect it is not surprising it does better when the true data-generating Lévy measure is sparse.

3.4. Computing time

In terms of computational effort, the time it takes to evaluate the Buchmann-Grübel estimator is negligible compared to our algorithm for sampling from the posterior. This is not surprising, as that frequentist estimator relies on a plug-in approach, whereas in our case an approximation to the posterior is obtained by MCMC simulation. However, if proper uncertainty quantification is desired, then the Bayesian method is advantageous in the sense that it does not solely produce a point estimate.

Note that the proposed MCMC scheme requires determination of the solutions to the Diophantine equation $\sum_{j=1}^{m}jk_{j}=z$ for all unique values $z$ in the observation set. For moderate values of $z$, say $z\leq 30,$ this is rather quick, but for large values of $z$ the computing time increases exponentially, as does the amount of the allocated memory. The computing time of each Metropolis-Hastings step is then very small, but we potentially need a very large number of iterations to reach stationarity. The latter is caused firstly by the fact that at a particular iteration, our proposals for $\mu_{ij}$ do not take into account the current values of $\nu_{1},\ldots,\nu_{m}$; secondly, the size of the state space that needs to be explored increases exponentially with $m$.

4. Real data examples

4.1. Horse kick data

To further illustrate our procedure, we will use the von Bortkewitsch data on the number of soldiers in the Prussian cavalry killed by horse kicks (available by year and by cavalry corps); this example was also employed in Buchmann & Grübel (2003). Each observation is an integer from [math] to $4$, giving the number of deaths for a given year and a given cavalry corps, with overall counts reported in Table 3. The data are extracted from the table on p. 25 in von Bortkewitsch (1898). Note that von Bortkewitsch corrects for the fact that the Guards and I, VI and XI cavalry corpses of the Prussian army had a different organisation from other units, and justifiably omits the corresponding counts from consideration.

It has been demonstrated by von Bortkewitsch that the Poisson distribution fits the horse kick data remarkably well. Assuming instead that observations follow a compound Poisson distribution is a stretch of imagination, as that would correspond to a horse running amok and killing possibly more than one soldier in one go. Nevertheless, this example provides a good sanity check for our statistical procedure.

The estimation results are graphically depicted in Figure 5. Clearly, point estimates of both methods are in agreement and lend support to the Poisson model for this dataset.

4.2. Plant data

Our second real example is the one used in Buchmann & Grübel (2004). Consider the data in Table 4, taken from Evans (1953). The data were collected as follows: the area was divided into plots of equal size and in each plot the number of plants was counted; the number of plants in each plot ranges from [math] to $12$. The second row of Table 4 gives the counts of plots containing a given number of plants; thus, there were $274$ plots that contained no plant, $71$ that contained $1$ plant, etc. It is customary in the ecological literature to model such count data as i.i.d. realizations from a compound Poisson distribution. Thus, e.g., Neyman (1939) advocated the use of a Poisson base distribution in this context; another option here is a geometric base distribution. Given existence of several distinct modelling possibilities, performing an exploratory non-parametric analysis appears to be a sensible strategy.

The estimation results are graphically depicted in Figure 6. There are some small differences between the posterior mean and the Buchmann-Grübel estimate, but overall they are very similar.

5. Frequentist asymptotics

In this section we assume that the observation times $\{t_{i}\}$ are equidistant: $t_{i}=i,i=1,\ldots,n$. To evaluate our Bayesian method from a theoretical point of view, we will verify that it is consistent, and we will establish the rate at which the posterior contracts around the ‘true’, data-generating Lévy measure $\nu_{0}$; see Ghosal & van der Vaart (2017) for a thorough treatment of Bayesian asymptotics from the frequentist point of view. From now on the subscript [math] in various quantities will refer to the data-generating distribution.

Our strategy consists in proving that the posterior contraction rate for $\nu_{0}$, given the sample $\mathcal{Z}_{n}=(Z_{1},\ldots,Z_{n})$, can be derived from the posterior contraction rate for $q_{0}$ given $\mathcal{Z}_{n}$, which is mathematically easier since $Z_{1},\ldots,Z_{n}$ is a sequence of independent and identically distributed random variables with distribution $q_{0}$. We therefore effectively avoid dealing directly with the inverse nature of the problem of estimating $p_{0}$.

The prior we consider in this section is defined as follows:

• •

Endow the rate $\lambda$ of the Poisson process with a prior distribution.

• •

Independently, endow the vector $(p_{1},\ldots,p_{m})$ with a Dirichlet distribution with parameter $(\alpha_{1},\ldots,\alpha_{m})$.

• •

Set a priori $p_{k}=0$ for all $k>m$.

This is a somewhat simplified version of the prior we used in Section 2, which allows us to concentrate on essential features of the problem, without need to clutter the analysis with extra and unenlightening technicalities. Also remember the well-known relationship between the Gamma and Dirichlet distributions: if $\xi_{1},\ldots,\xi_{m}$ are independent Gamma distributed random variables, $\xi_{i}\sim\operatorname{Gamma}(\alpha_{i},1)$, then for $\eta_{i}=\xi_{i}/\sum_{j=1}^{m}\xi_{j}$, the vector $(\eta_{1},\ldots,\eta_{m})$ follows the Dirichlet distribution with parameter $(\alpha_{1},\ldots,\alpha_{m})$; furthermore, we have that $\sum_{j=1}^{m}\xi_{j}\sim\operatorname{Gamma}\left(\sum_{j=1}^{m}\alpha_{j},1\right)$, and $\eta_{i}$ are independent of $\sum_{j=1}^{m}\xi_{j}$.

In our asymptotic setting, we will make $m=m_{n}$ dependent on $n$ and let $m_{n}\rightarrow\infty$ at a suitable rate as $n\rightarrow\infty$.

Recall that we write $Q=(q_{k})_{k\in\mathbb{N}_{0}}$ for $q_{k}=Q(\{k\})$. Let ${\mathcal{Q}}$ denote the collection of all probability measures supported on $\mathbb{N}$.

Theorem 1**.**

Suppose there exists $\underline{\alpha}$, such that $0<\underline{\alpha}\leq\alpha_{i}\leq 1$ for all $1\leq i\leq m_{n}$. Suppose $\lambda\sim\operatorname{Gamma}(a,b)$ with $a\in(0,1]$ and that $\nu_{0}$ has a compact support. Then, for any $\gamma>1$,

[TABLE]

in $Q_{0}^{n}$-probability, as $n\to\infty$.

Remark 1*.*

Note that since the support of $\nu_{0}$ is not assumed to be known, our CPP model is still naturally non-parametric. The assumption of the compact support of $\nu_{0}$ does not cover the simulation example of Section 3.2. However, its relaxation appears to pose very difficult technical challenges and is not attempted in this work.

The remainder of this section is devoted to the proof of Theorem 1.

5.1. Basic posterior inequality via the stability estimate

A key step of the proof of Theorem 1 is the stability estimate in Equation (5) below, which bounds the total variation distance between the Lévy measures $\nu,\nu^{\prime}$ in terms of the total variation distance between the corresponding compound distributions $q,q^{\prime}$.

In principle, it is conceivable that the Panjer recursion should allow one to bound probability distances between $P$-probabilities via distances between $Q$-probabilities; we call such a bound a stability estimate. Nevertheless, explicit as the equations of the Panjer recursion are, they are still somewhat unwieldy for that purpose. Hence we will use another inversion formula from Buchmann & Grübel (2003), that will lead to the stability estimate we are after.

First we introduce some notation, and also recall a few useful facts summarised in Buchmann & Grübel (2003). The space of absolutely summable sequences is defined as $\ell_{1}\coloneqq\left\{a\in\mathbb{R}^{\mathbb{N}_{0}}:\sum_{j=0}^{\infty}|a_{j}|<\infty\right\}$, with a norm given by $\|a\|_{1}=\sum_{j=0}^{\infty}|a_{j}|$. For probability vectors $a,b$, the norm $\|a-b\|_{1}$ is (twice) the total variation distance between $a$ and $b$. For any $a,b\in\ell_{1}$, we have the inequality

[TABLE]

where $\ast$ denotes convolution of $a$ and $b.$ We define a mapping $a\mapsto\exp(a)$ from $\ell_{1}$ into $\ell_{1}$ via

[TABLE]

The exponential has the following two useful properties:

[TABLE]

and

[TABLE]

We define a sequence $\delta_{0}=(\delta_{0,k})_{k\in\mathbb{N}_{0}}$, such that $\delta_{0,0}=1$ and its all other entries are equal to zero. Then, using the above properties of the exponential, we can write concisely the compounding mapping in (2) in terms of convolutions of infinite sequences: $q=\exp(\lambda(p-\delta_{0}))$. Its convolution inverse, i.e. $q^{\ast(-1)}$ such that $q^{\ast(-1)}\ast q=\delta_{0}$, is given by $r=q^{\ast(-1)}=\exp(-\lambda(p-\delta_{0}))$. Note that $r\in\ell_{1}$. We have the following recursive expressions

[TABLE]

Lemma 1**.**

Let $q,q^{\prime}$ correspond to two pairs $(\lambda,p)$ and $(\lambda^{\prime},p^{\prime})$, respectively (and $r$ correspond to $q$, i.e. the pair $(\lambda,p)$). Then, in accordance with the notation introduced above and provided that $\|q^{\prime}-q\|_{1}<\|r\|_{1}^{-1}$, it holds that

[TABLE]

Proof.

The result is a direct consequence of Lemma 3 in Buchmann & Grübel (2003), which states that

[TABLE]

whenever $\|q^{\prime}-q\|_{1}<\|r\|_{1}^{-1}$. Taking the $\|\cdot\|_{1}$-norm on both sides and some elementary bounding via (4) imply that

[TABLE]

and thus Equation (5) follows. ∎

We will use Equation (5) to establish the key inequality for the posterior measure $\Pi(\cdot\mid\mathcal{Z}_{n})$. We recall once again that the subscript [math] refers to ‘true’, data-generating quantities.

Proposition 1**.**

For any prior $\Pi$ on $\nu$, for any $\varepsilon\in(0,1]$ and for any $n\geq 1$, the following posterior inequality holds:

[TABLE]

Proof.

Write $\{\nu:\|\nu-\nu_{0}\|_{1}\geq\varepsilon\}$ as a union of the sets

[TABLE]

and

[TABLE]

Thanks to Lemma 1, the set

[TABLE]

is a subset of $\{\nu:\|q-q_{0}\|_{1}\geq{\varepsilon}/(2\|r_{0}\|_{1})\}$. The proof is concluded by observing that $\left\{\nu:\|\nu-\nu_{0}\|_{1}\geq\varepsilon\right\}\cap\left\{\nu:\|r_{0}\|_{1}\|q-q_{0}\|_{1}\geq 1/2\right\}$ is a subset of $\{\nu:\|q-q_{0}\|_{1}\geq{\varepsilon}/(2\|r_{0}\|_{1})\}$, too, since $\varepsilon\leq 1$. ∎

In general, stability estimates like the one in Equation (5) are unknown in the literature on Lévy processes. Consequently, studying Bayesian asymptotics for Lévy models, even in the CPP case, necessitates the use of very intricate arguments under restrictive assumptions; see, e.g., Nickl & Söhl (2017).

5.2. Proof of Theorem 1

The usefulness of Proposition 1 above lies in the fact that the posterior contraction rate in the inverse problem of estimating the Lévy measure $\nu_{0}$ from indirect observations $\mathcal{Z}_{n}$ can be now deduced from the posterior contraction rate in the direct problem of estimating the compound distribution $q_{0}$, which is easier (observe that $r_{0}$ is determined by $\nu_{0}$ and is therefore fixed in the proofs). The general machinery developed in Ghosal et al. (2000) can be applied to handle the latter, and also several inequalities obtained in Gugushvili et al. (2015) are useful in that respect. In particular, we make use of the following inequality for the Hellinger distance,

[TABLE]

Cf. Lemma 1 in Gugushvili et al. (2015). To ease our notation, in the sequel we will often write $q$ and $q^{\prime}$ instead of $q_{\lambda,p}$ and $q_{\lambda^{\prime},p^{\prime}},$ respectively.

Denote

[TABLE]

Another two inequalities we will use are the following: let $\lambda,\lambda_{0}\in[\underline{\lambda},\overline{\lambda}].$ Then there exists a positive constant $\overline{C}$, such that

[TABLE]

cf. equations (14) and (15) in Lemma 1 in Gugushvili et al. (2015).

These three inequalities can be obtained by adjustment of the arguments used in Gugushvili et al. (2015). However, we opted to give their direct proofs in Lemma 5 from Appendix  A under slightly weaker conditions than required in Gugushvili et al. (2015).

Our proof of Theorem 1 proceeds via verification of the conditions for posterior contraction in Theorem 2.1 in Ghosal et al. (2000). In our setting, the latter result reads as follows.

Theorem 2**.**

Assume ${\mathcal{Z}}_{n}=(Z_{1},\ldots,Z_{n}),$ where $Z_{1},\ldots,Z_{n}$ are independent and identically distributed with distribution $q_{0}$. Let $h$ denote the Hellinger metric on ${\mathcal{Q}}$, a collection of all measures with support in $\mathbb{N}$. Suppose that for a sequence $\{\epsilon_{n}\}$ with $\epsilon_{n}\to 0$ and $n\epsilon_{n}^{2}\to\infty$, a constant $C>0$ and sets ${\mathcal{Q}}_{n}\subset{\mathcal{Q}}$, we have

[TABLE]

[TABLE]

[TABLE]

Then, for sufficiently large $M>0$, we have that $\Pi_{n}\left(Q\colon h(q,q_{0})\geq M\epsilon_{n}\mid{\mathcal{Z}}_{n}\right)\to 0$ in $Q_{0}^{n}$-probability.

We will now verify the three conditions of this theorem, which we refer to as the entropy condition, the remaining mass condition, and the prior mass condition, respectively. To that end, fix strictly positive sequences $\{\underline{\Lambda}_{n}\}$, $\{\overline{\Lambda}_{n}\}$, and define the sieves

[TABLE]

5.2.1. Entropy

We start with bounding the entropy of the sieve ${\mathcal{Q}}_{n}$ for $h$-balls of radius $\epsilon_{n}$.

Lemma 2**.**

Assume that as $n\to\infty,$

[TABLE]

Then

[TABLE]

Proof.

For $\lambda,\lambda^{\prime}\geq\underline{\Lambda}_{n}$,

[TABLE]

Furthermore, from Section 3.3 in Pollard (2002),

[TABLE]

Combining the preceding two displays and Equation (6), we get

[TABLE]

Hence, if

[TABLE]

then the Hellinger distance between $q_{\lambda,p}$ and $q_{\lambda^{\prime},p^{\prime}}$ is bounded by $\epsilon_{n}$. To cover $[\underline{\Lambda}_{n},\overline{\Lambda}_{n}]$, we need at most $\lfloor\frac{\overline{\Lambda}_{n}}{2\epsilon_{n}\sqrt{\underline{\Lambda}_{n}}}\rfloor+1$ intervals of length $2\sqrt{\underline{\Lambda}_{n}}\epsilon_{n}$. To cover discrete distributions with support in $\{1,\ldots,m_{n}\}$, we need at most

[TABLE]

$L_{\infty}$-balls of radius $\epsilon_{n}^{2}/(4\overline{\Lambda}_{n}m_{n})$. Under assumption (8), the summand $1$ in the above display is asymptotically negligible and can be omitted. In that case, the number of $h$-balls that we need to cover ${\mathcal{Q}}_{n}$ is of order

[TABLE]

Taking the logarithm and next a straightforward rearrangement of the terms gives the statement of the lemma. ∎

5.2.2. Remaining prior mass

Now we will derive an inequality for the remaining prior mass.

Lemma 3**.**

For $\lambda\sim\operatorname{Gamma}(a,b)$ with $0<a\leq 1$,

[TABLE]

Proof.

We have (with a slight abuse of notation)

[TABLE]

Now,

[TABLE]

Furthermore,

[TABLE]

The proof is concluded. ∎

5.2.3. Prior mass

Finally, we lower bound the prior mass in a small Kullback-Leibler neighbourhood of the data-generating compound distribution $q_{0}$. Define the function $g\colon(0,\infty)\times(0,1)\to(0,\infty)$ by

[TABLE]

where $C$ is the constant appearing in the statement of Lemma 6 below.

Lemma 4**.**

Assume that

1. (i)

there exists $\underline{\alpha}$, such that $0<\underline{\alpha}\leq\alpha_{i}\leq 1$ for all $1\leq i\leq m_{n}$; 2. (ii)

strictly positive sequences $\underline{p}_{n}\rightarrow 0$, $\epsilon_{n}\rightarrow 0$ and $m_{n}\rightarrow\infty$ satisfy the inequalities $m_{n}g(\epsilon_{n},\underline{p}_{n})<1$ and $\underline{p}_{n}<g(\epsilon_{n},\underline{p}_{n})^{2}.$

Define

[TABLE]

Then

[TABLE]

Here $\widetilde{\epsilon}_{n}=\epsilon_{n}/\sqrt{3\overline{C}}$, with a constant $\overline{C}>0$ not depending on $n$.

Proof.

Define

[TABLE]

For all $n$ large enough and $\epsilon$ small, we have $\{\lambda:|\lambda_{0}-\lambda|\leq\epsilon\}\subseteq[\underline{\Lambda}_{n},\overline{\Lambda}_{n}].$ Then by inequalities in Lemma 5, $\widetilde{B}_{n}(\epsilon)\subset B_{n}\left(\sqrt{3\overline{C}}\epsilon\right)$, with a constant $\overline{C}$ that can be taken the same for all large enough $n$; see the arguments in Section 4.2 in Gugushvili et al. (2015). Hence, using the a priori independence of $p$ and $\lambda$,

[TABLE]

where

[TABLE]

Furthermore, by Lemma 6 from Appendix A, we have

[TABLE]

The statement of the lemma now follows upon applying Lemma 7 from Appendix A with $\eta=\underline{p}_{n}$ and $\epsilon=g(\widetilde{\epsilon}_{n},\underline{p}_{n})$. ∎

5.2.4. Using bounds in Theorem 2

We take

[TABLE]

with appropriately selected proportionality constants, and verify the conditions in Theorem 2.

Firstly, condition (8) is trivially satisfied. Therefore, we can invoke Lemma 2 and conclude that the entropy is upper bounded by a multiple of $\log^{2\gamma}n$, since $\gamma>1.$ Now $\log^{2\gamma}n\lesssim n\epsilon_{n}^{2}$, and this verifies the entropy condition in Theorem 2.

Be Lemma 3, for a suitable choice of the constant $C$ the remaining prior mass condition is likewise satisfied.

Finally, for the prior mass condition in a small Kullback-Leibler neighbourhood to hold, by Lemma 4 we need that the term

[TABLE]

is lower bounded by $\exp(-Cn\epsilon_{n}^{2})$ for some large enough $C>0$. Now, $\Pi_{n}\left(|\lambda_{0}-\lambda|\leq\widetilde{\epsilon}_{n}\right)\asymp\widetilde{\epsilon}_{n}.$ Take the logarithm on both sides of the above display and note that by our conditions

[TABLE]

Likewise,

[TABLE]

so that the prior mass condition holds.

Thus we have verified all the conditions of Theorem 2. The resulting posterior contraction rate is $\epsilon_{n}\asymp{\log^{\gamma}n}/{\sqrt{n}}.$

6. Outlook

In this paper we introduced a non-parametric Bayesian approach to estimation of the Lévy measure $\nu$ of a discretely observed CPP, when the support of $\nu$ is a subset of $\mathbb{N}.$ We constructed an algorithm for sampling from the posterior distribution of $\nu$, and showed that in practice our procedure performs well and measures up to a benchmark frequentist plug-in approach from Buchmann & Grübel (2004). Although computationally more demanding and slower than the latter, our method has an added benefit of providing uncertainty quantification in parameter estimates through the spread of the posterior distribution. On the theoretical side we show that our procedure is consistent, in that asymptotically, as the sample size $n\rightarrow\infty,$ the posterior concentrates around the ‘true’, data-generating distribution. The corresponding posterior contraction rate is the (nearly) optimal rate $\log^{\gamma}n/\sqrt{n}$ for an arbitrary $\gamma>1$, if we are to ignore a practically insignificant $\log n$ factor.

Among several generalisations of our results, the one that looks the most promising is extension of our methodology to CPP processes with jump size distributions supported on the set of integers $\mathbb{Z}.$ The corresponding model has garnered substantial interest in financial applications, see Barndorff-Nielsen et al. (2012). We leave this extension as a topic of possible future research.

Acknowledgements

The authors would like to thank the Associate Editor and the referee for their detailed and constructive comments on the paper.

Appendix A Technical results

Lemma 5**.**

Let $q$ (resp. $q^{\prime}$) be the law at time $1$ of a compound Poisson process with intensity $\lambda$ (resp. $\lambda^{\prime}$) and jump distribution $p$ (resp. $p^{\prime}$). Suppose that $p$ and $p^{\prime}$ are distributions concentrated on $\mathbb{N}$. Then,

[TABLE]

In particular, if $\lambda,\lambda^{\prime}\in[\underline{\Lambda},\overline{\Lambda}]$ with $0<\underline{\Lambda}\leq\overline{\Lambda}<\infty$, then there exists a positive constant $\overline{C}$, that depends on $\underline{\Lambda}$, $\overline{\Lambda}$, such that

[TABLE]

Proof.

If $KL(p,p^{\prime})$ and $V(p,p^{\prime})$ are infinite, then the above inequalities are trivially satisfied. Therefore, we can assume these two divergences are finite. With this in mind, the proof of the lemma is divided into three steps.

Step 1: We begin by proving that for any $n\geq 1,$

[TABLE]

The assertions are trivial for $n=1$. Assuming that the first one holds for $n-1$ with $n\geq 2$, we will now show that it holds for $n$ as well. Using the notation $p^{*n}(i)$ for the $i$th element of $p^{*n}$ and similarly in the case of $p^{\prime*n}$, we have

[TABLE]

where the inequality follows from the log-sum inequality, and the last equality is obtained by means of Fubini’s theorem combined with the facts that

[TABLE]

By induction, we deduce that KL\big{(}p^{*n},p^{\prime*n}\big{)}\leq nKL(p,p^{\prime}).

The proof of the inequality for $V$ is similar: we assume the inequality is true for $n-1$ with $n\geq 2$, and will show it holds for $n$ as well. Write $V\left(p^{*n},p^{\prime*n}\right)=R_{1}+R_{2}$ for

[TABLE]

Observe that the function $x\mapsto(x\log^{2}x)\mathbf{1}_{\{x\geq 1\}}$ is convex. By Jensen’s inequality we have for positive $a_{k},b_{k}$ that

[TABLE]

Using this inequality and

[TABLE]

we get that

[TABLE]

where in the last inequality we used the induction hypothesis. Now recall an elementary inequality

[TABLE]

valid for $x\geq 0$; see p. 12 in Gugushvili et al. (2015). Applying this inequality to

[TABLE]

such that ${p^{*n}(i)}/{p^{\prime*n}(i)}<1,$ we get

[TABLE]

By multiplying both sides of the above inequality with $p^{\prime*n}(i),$ summing the result through $i$ and recalling the definition of the Hellinger distance, we get that

[TABLE]

To conclude the proof of the inequality for $V$, we combine the bounds derived for $R_{1}$ and $R_{2}.$

As far as the inequality for the Hellinger distance is concerned, we observe that

[TABLE]

for a convex function $g(x)=(1-\sqrt{x})^{2}\mathbf{1}_{[0,\infty)}(x).$ Then, by the same reasoning as above, we have

[TABLE]

Now note that the last summand satisfies

[TABLE]

We therefore conclude that

[TABLE]

which leads to the desired inequality, by an induction argument.

**Step 2: ** Now we prove the inequalities

[TABLE]

Here $N$ and $N^{\prime}$ are Poisson random variables with means $\lambda$ and $\lambda^{\prime}$, respectively, and with a slight abuse of notation, $KL(N,N^{\prime})$, $V(N,N^{\prime})$ and $h(N,N^{\prime})$ are the $KL$ and $V$ divergences and the Hellinger distance between the corresponding laws.

Note that

[TABLE]

Using this and the log-sum inequality,

[TABLE]

For the divergence $V,$ write $V(q,q^{\prime})=B_{1}+B_{2}$ for

[TABLE]

To control $B_{2}$, we use the same arguments as in the proof of inequalities (12) and (15) in Gugushvili et al. (2015), getting

[TABLE]

This gives the required inequality for the $V$ divergence.

Finally, we prove the inequality for the Hellinger distance. Denoting the law of $\sum_{j=1}^{N}Y_{j}$ by $\tilde{q}$, it holds by the triangle inequality that $h(q,q^{\prime})\leq h(q,\tilde{q})+h(\tilde{q},q^{\prime}).$ Since $g(x)=(1-\sqrt{x})^{2}\mathbf{1}_{[0,\infty)}(x)$ is a convex function,

[TABLE]

It remains to prove that $h(\tilde{q},q^{\prime})\leq h(N,N^{\prime}).$ This again follows by convexity of $g$, since

[TABLE]

Step 3: From Steps 1 and 2 we derive that

[TABLE]

Now the proof of the lemma follows from these three inequalities upon noticing that $KL(p,p^{\prime})^{2}\leq V(p,p^{\prime})$, and recalling that

[TABLE]

∎

Lemma 6**.**

Let $\epsilon>0$. Suppose $p=(p_{1},\ldots,p_{m})$ and $p^{\prime}=(p^{\prime}_{1},\ldots,p^{\prime}_{m})$ are points in the $m$-dimensional unit simplex, and let $\min_{1\leq i\leq m}p_{i}\geq c$ for some $c\in(0,1)$. Then there exists a universal constant $C>0$, such that the inequality

[TABLE]

implies that $KL(p^{\prime},p)\leq\epsilon^{2}$ and $V(p^{\prime},p)\leq\epsilon^{2}$ hold.

Proof.

Lemma 8 in Ghosal & van der Vaart (2007) assures that there exists a constant $\overline{C}$ (not depending on either $p$ or $p^{\prime}$), such that

[TABLE]

and

[TABLE]

From Section 3.3 in Pollard (2002) we have $h^{2}(p^{\prime},p)\leq\sum_{i=1}^{m}|p^{\prime}_{i}-p_{i}|$. Furthermore, since $0<c<1$ and $\min_{1\leq i\leq m}p_{i}\geq c$,

[TABLE]

Therefore,

[TABLE]

from which the assertion of the lemma follows trivially. ∎

Lemma 7**.**

Let $m\geq 2$ be an integer. Suppose $(p_{1},\ldots,p_{m})\sim\operatorname{Dir}(\alpha_{1},\ldots,\alpha_{m})$. Let $p_{0}=(p_{01},\ldots,p_{0m})$ be an arbitrary point in the $m$-dimensional unit simplex. Assume there exists $\underline{\alpha}$, such that $0<\underline{\alpha}\leq\alpha_{i}\leq 1$ for all $1\leq i\leq m$. Let $\epsilon>0$, and let $\eta$ be such that $\eta<\epsilon^{2}$. Then if $m\epsilon<1$,

[TABLE]

Proof.

By arguments analogous to those in the proofs of Lemma 6.1 in Ghosal et al. (2000) and Lemma 10 in Ghosal & van der Vaart (2007), we obtain that the left-hand side of (10) can be lower bounded by

[TABLE]

The length of the integration interval in each of the integrals in the above product is lower bounded by $\epsilon^{2}-\eta$. Using that $\underline{\alpha}\leq\alpha_{i}\leq 1$, we deduce that the preceding display is lower bounded by

[TABLE]

This entails the result. ∎

Appendix B Bootstrap confidence intervals

Here we report a small comparison between the Bayesian credible intervals and the bootstrap confidence intervals for the Buchmann-Grübel estimator. The setup and the simulated dataset that we used are the same as in Subsection 3.1. The bootstrap confidence intervals were computed as follows: $B=9999$ bootstrap samples were generated from the compound Poisson model under the Buchmann-Grübel estimates computed from the observed data. These bootstrap samples were then fed back to the Buchmann-Grübel procedure to yield $B$ bootstrap estimates of the Lévy measure $\nu$. Finally, for each $k$, the $\alpha/2$-th and $(1-\alpha/2)$-th sample quantiles were obtained to yield $1-\alpha$ bootstrap confidence intervals for $\nu_{k}.$ The results with $\alpha=0.05$ are displayed in Figure 7. One observes that while both methods result in a good coverage for this specific dataset, the bootstrap appears to be noticeably more conservative than the Bayesian approach, as evidenced by the width of the intervals.