Decompounding discrete distributions: A non-parametric Bayesian approach
Shota Gugushvili Ester Mariucci Frank van der Meulen


1.25

Decompounding discrete distributions]Decompounding discrete distributions: A non-parametric Bayesian approach

S. Gugushvili]Shota Gugushvili^1
^1Biometris

	Wageningen University & Research
shota@yesdatasolutions.com

E. Mariucci]Ester Mariucci^2
^2Institut für Mathematik

Potsdam Universität
mariucci@ovgu.de

F. H. van der Meulen]Frank van der Meulen^3
^3Delft Institute of Applied Mathematics

	Delft University of Technology
f.h.vandermeulen@tudelft.nl

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'







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übel. On a theoretical side, we study the posterior from the frequentist point of view and prove that as the sample size n→∞, it contracts around the `true', data-generating parameters at  rate 1/√(n), up to a log n factor. 





[
    [
    
=====





§ INTRODUCTION





 §.§ Problem formulation


Let N=(N_t : t≥ 0) be a Poisson process with a constant intensity λ>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≥ 0) is 

    X_t=∑_j=1^N_t Y_j,

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 <cit.>.

Assume that the process X is observed at discrete times 0<t_1<t_2<…<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,…,X_t_n, our goal is to estimate the jump size distribution P and the intensity λ.  We specifically restrict our attention to the case where P is a discrete distribution, P(ℕ)=1, and we will write p=(p_k)_k∈ℕ 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^-λ, which does not allow consistent estimation of the intensity λ. For a discussion of applications of this CPP model in risk theory, see <cit.>.

Define the increments Z_i = X_t_i - X_t_i-1, i=1,…,n. Then 𝒵_n = (Z_i:i = 1,…,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(ℕ_0)=1. As 𝒵_n carries as much information as (X_t_i:i=1,…,n) does, we can base our estimation procedure directly on the increments 𝒵_n. Since summing up the jumps Y_j's amounts to compounding their distributions, the inverse problem of recovering P and λ from Z_i can be referred to as decompounding; see <cit.>.

There are two natural ways to parametrise the CPP model: either in terms of the pair (λ,p), or in terms of the Lévy measure ν=(ν_k)_k∈ℕ of the process X, see <cit.>. A relationship between the two is λ=∑_k=1^∞ν_k and p=ν/λ. 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.



 §.§ Approach and results


In this paper, we take a non-parametric Bayesian approach to estimation of the Lévy measure ν of X. See <cit.> and <cit.> 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 <cit.>). The posterior contraction rate, up to a practically insignificant log n factor, turns out to be 1/√(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.



 §.§ 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 <cit.>, <cit.>, <cit.> and <cit.>. 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 <cit.> and <cit.>, 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, <cit.> 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. <cit.> fix this problem by truncation and renormalisation. This works, but looks artificial. As noted in <cit.>, in practice the latter approach breaks down if no zero values are observed among Z_i's. <cit.> 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 <cit.> and <cit.> do not allow obvious generalisations to non-equidistant observation times {t_i}. In <cit.>, 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 <cit.>  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 <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.> and <cit.>. We refer to <cit.>, <cit.>, <cit.> and <cit.> for more extensive literature surveys.




 §.§ Outline

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



 §.§ Notation

For two sequences {a_n} and {b_n} of positive real numbers, the notation a_n≲ b_n  (or b_n≳ a_n) means that there exists a constant C>0 that is independent of n and such that a_n≤ C b_n. We write a_n≍ b_n if both a_n≲ b_n and a_n≳ b_n hold.
We denote a prior (possibly depending on the sample size n) by Π_n. The corresponding posterior measure is denoted by  Π_n(·|𝒵_n). The Gamma distribution with shape parameter a and rate parameter b (a,b >0) is denoted by Gamma(a,b). Its density is given by
x ↦b^a/Γ(a) x^a-1e^-bx,    x>0,
where Γ is the Gamma function. The inverse Gamma distribution with shape parameter a and scale parameter b is denoted by (a,b). Its density is 
x ↦b^a/Γ(a) x^-a-1e^-b / x,    x>0. We use the notation Exp(a) for an exponential distribution with mean 1/a.
Finally, given a metric d on a set 𝒬 and ϵ > 0, the covering number N(ϵ,𝒬, d) is defined as the minimal number of balls of radius ϵ needed to cover 𝒬.


§ 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.



 §.§ Prior


We define a prior Π on ν through a hierarchical specification 

    {ν_k}_k=1^∞| a, m, β_k    Gamma(a,1/β_k)	·1_{1≤ k≤ m}, 
    β_1,…, β_m |γ   (c,γ),
    γ   ∼Exp(1).
	
 Note that the (fixed) hyperparameters m∈ℕ, a, c >0 are denoted by  Latin letters. 

The hyperparameter m incorporates our  a priori opinion on the  support of the Lévy measure ν, 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,…,n Z_i; this latter is the maximal value suggested by the data 𝒵_n at hand. Nevertheless, we may simultaneously expect that the `true', data-generating ν charges full mass only to a proper, perhaps even a small subset of the set {1,…,m}. In other words, ν 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 ν_0, or Z_i is predominantly formed of many small jumps, which in turn indicates that the intensity λ of the Poisson arrival process N may be large. To achieve accurate estimation results, a prior should take a possible sparsity of ν into account. This is precisely the reason of our hierarchical definition of the prior Π: a small β_k encourages a priori the shrinkage of the components ν_k of ν towards zero.





 §.§ Data augmentation


Assume temporarily t_i=i, i=1,…,n, and write q=(q_k)_k∈ℕ_0 for q_k=q({k}). Then Z_i have the distribution

    q = e^-λ∑_j=0^∞λ^j/j!p^∗ j,

with * denoting convolution.  The compounding mapping (λ,p)↦ q can be expressed explicitly via the Panjer recursion (see <cit.>):

    q_0 = e^-λ,    q_k = λ/k∑_j=1^k j p_j q_k-j,    k∈ℕ.

This recursion can be inverted to give the inverse mapping q↦(λ,p) via

    λ = - log q_0,    p_k = - q_k/q_0log q_0 - 1/kq_0∑_j=1^k-1jp_jq_k-j,    k∈ℕ.


In view of (<ref>), the likelihood in the CPP model is explicit. Nevertheless, an attempt to directly use (<ref>) 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 (<ref>) and the Panjer recursion, we will employ  data augmentation, see <cit.>. 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∈[0,T]) of CPP, the likelihood is tractable and is proportional to

    e^-T∑_k=1^mν_k∏_k=1^m ν_k^μ_k,

see <cit.>, p. 498. Here

    μ_k = #{Y_j=k},

i.e. the total number of jumps of size k.

Then the  prior Π from Subsection <ref> leads to conjugate posterior computations. In fact, the full conditionals are

    ν_k |{μ_k}, {β_k}Gamma(a+μ_k,1/β_k+T),    k=1,…,m,
    β_k |{ν_k}, γ(a + c,γ + ν_k),    k=1,…,m,
    γ|{β_k}∼Gamma(c m+1, 1+∑_k=1^m β_k^-1).

Therefore, the Gibbs sampler for posterior inference on ν 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.,  <cit.> and Section 24.5 in <cit.> 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

    Z_i = ∑_j=1^m j μ_ij,

where (μ_ij:i=1,…,n, j=1,…,m) are independent, and μ_ij∼Poisson(Δ_iν_j) for ν_j = λ p_j and Δ_i = t_i - t_i-1; see Corollary 11.3.4 in <cit.>. Furthermore, for μ_k as in (<ref>) we trivially have μ_k=∑_i=1^n μ_ik. Data augmentation iterates the following two steps:

	
  * Draw (μ_ij) conditional on the data 𝒵_n and the parameter ν.
	
  * Draw ν conditional on (μ_ij).

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

Thus, keeping ν fixed, for each i we want to compute the conditional distribution
(μ_ij:j=1,…,m) | Z_i, and furthermore, we want to be able to simulate from this distribution. In turn, this will immediately allow us to simulate μ_k conditional on the data 𝒵_n. Now, with (·) referring to probability under the parameter ν, it holds that

    (μ_i1=k_1,…,μ_im=k_m | Z_i = z_i) 
    
     = 1/(Z_i = z_i) e^-Δ_i ∑_j=1^m ν_j∏_j=1^m(Δ_iν_j)^k_j/k_j!{∑_j=1^m j k_j = z_i}.

Knowledge of the normalising constant (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 <cit.>, 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 <cit.>. We start by observing that for a fixed i, the support of ( ·|  Z_i = z_i) is precisely the set 𝒮_i of non-negative solutions (k_1,…,k_m) of the Diophantine equation
∑_j=1^m j k_j = z_i.
The R package nilde (see <cit.>) implements an algorithm from <cit.> that finds all such solutions for given integers m and z_i. By Markovianity of the process X, we can simulate the vectors (μ_i1,…,μ_im)  independently for each i=1,…,n. If z_i=0 or 1, there is only one solution to the Diophantine equation: the trivial solution (0,…,0) in the first case, and the solution (1,0,…,0) in the second case; for such z_i, no simulation is required, as (μ_i1,…,μ_im) is known explicitly. We thus only need to consider each i∈ℐ = {i:z_i ≠ 0  or  1} in turn, and design a Metropolis-Hastings move on the set of the corresponding solutions 𝒮_i. Fix once and for all an ordering of elements in 𝒮_i (this could be, e.g., lexicographic, or reverse lexicographic); we use the notation |𝒮_i| for the cardinality of 𝒮_i. Let μ = (μ_i1,…,μ_im) be the current state of the chain, corresponding to the ℓth element s_ℓ of 𝒮_i. A proposal μ^∘ = (μ_i1^∘,…,μ_im^∘) is obtained as follows:





	
  * If ℓ = 1, draw μ^∘ uniformly at random among the elements {s_2,s_|𝒮_i|} of 𝒮_i.
	
  * If ℓ = |𝒮_i|, draw  μ^∘ uniformly at random among the elements {s_1,s_ |𝒮_i| -1 } of 𝒮_i.
	
  * If ℓ≠ 1 or |𝒮_i|, draw  μ^∘ uniformly at random among the elements {s_ℓ-1,s_ℓ+1} of 𝒮_i.


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

	
	
  * Draw μ^∘ = (μ_i1^∘,…,μ_im^∘) uniformly at random from 𝒮_i.

The two proposals lead to reversible moves, and one may also alternate them with probabilities π and 1-π, e.g. π=0.8. The logarithm of the acceptance probability of a move from (μ_i1,…,μ_im) to (μ_i1^∘,…,μ_im^∘) is computed as

    log A = ∑_k=1^m (μ_ik^∘ - μ_ik) log(Δ_i  ν_k) + ∑_k=1^m {log(μ_ik!) - log(μ_ik^∘!) }.

The move is accepted if log U ≤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 (μ_i1^∘,…,μ_im^∘). Otherwise the chain stays in (μ_i1,…,μ_im). 












§ 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 <cit.>. Two real data examples are given in Section <ref>. 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 β_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 <cit.>. The code and datasets for replication of our examples are available on GitHub[See <https://github.com/fmeulen/Bdd>] and Zenodo, see <cit.>.



 §.§ Uniform base distribution


This simulation example follows with some extensions that in <cit.>. Let λ_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, Δ_i=1 for 1≤ i≤ n;
  
  (b) n=500, Δ_i=1 for 1≤ i≤ n (the data under (a) are augmented with 400 extra observations);
  
  (c) n=500, Δ_i=Unif(0,2) for 1≤ i≤ n.  



We set m= min(15,Z_(n)), where Z_(n) = max_1≤ i ≤ 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 ν_k. The results for the first 10 coefficients are shown in Figure  <ref>. For comparison, the estimator from <cit.> is also included in the figure. 

For setting (b), traceplots of every 50th iteration for a couple of coefficients ν_k are shown in Figure <ref>. 








We measure the error of an estimate {ν̂_k} by Err(ν, ν̂) = ∑_k=1^∞ |ν̂_k - ν_k|. The errors are reported in Table <ref>. In all  settings, for these particular realisations of the simulated data, the Bayesian procedure outperformed the truncated estimator from <cit.>. For setting (c), the latter produces a poor result, as was to be expected, given that it is derived under the assumption Δ_i=1 for all i. An advantage of the Bayesian procedure is the included measure of uncertainty, namely the credible intervals for ν_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 <cit.> and <cit.>, 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., <cit.> and <cit.>, Chapter 23). Here instead we opt for a numerical illustration, which is reported in Appendix <ref>.












 §.§ Geometric base distribution



The setup of this synthetic data example likewise follows that in <cit.>. Assume q is a geometric distribution with parameter α, i.e. q_k=(1-α)^k α for 0 < α <1, k∈_0. Then λ=-logα, and

    p_k = - (1-α)^k/klogα,    k∈ℕ.

Hence, ν_k=(1-α)^k/k.

We consider two simulation setups:


  (a) n=500, Δ_i=1 for 1≤ i ≤ n and  α=1/3;

  (b) n=500, Δ_i=1 for 1≤ i ≤ n and α=1/6.

We set m= min(15,Z_(n)) and ran the sampler according to the settings of Section <ref>. 
The results for both scenarios (a) and (b) are reported in Figure <ref>. In Table <ref> 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 α = 1/3 is smaller than that for α = 1/6. This appears intuitive, as a smaller value of α corresponds to a larger value of λ. 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.















 §.§ 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:

  
  * The setting from Section <ref> with n=250. We took m=min(15, Z_(n)).   
  
  * The setting from Section <ref> with α=1/3. We took m=min(20, Z_(n)).    

In both cases we assumed Δ_i=1 for all 1≤ i ≤ 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 <ref> 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 ν_k and the estimate within a particular simulation run  by ν̂_k, the error is defined by  Err(ν, ν̂) = ∑_k=1^∞ |ν̂_k - ν_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 <ref> 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. 





 §.§ 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 ∑_j=1^m j k_j = z for all unique values  z in the observation set. For moderate values of z, say z≤ 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 μ_ij do not take  into account the current values of ν_1,…, ν_m; secondly, the size of the state space that needs to be explored  increases exponentially with m. 



§ REAL DATA EXAMPLES



 §.§ 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 <cit.>. Each observation is an integer from 0 to 4, giving the number of deaths for a given year and a given cavalry corps, with overall counts reported in Table <ref>. The data are extracted from the table on p. 25 in <cit.>. 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 <ref>. Clearly, point estimates of both methods are in agreement and lend support to the Poisson model for this dataset. 






 §.§ Plant data


Our second real example is the one used in <cit.>. Consider the data in Table <ref>, taken from <cit.>. 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 0 to 12. The second row of Table <ref> 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., <cit.> 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 <ref>. There are some small differences between the posterior mean and the Buchmann-Grübel estimate, but overall they are very similar. 






§ FREQUENTIST ASYMPTOTICS



In this section we assume that the observation times {t_i} are equidistant: t_i = i, i=1,…,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 ν_0; see <cit.> for a thorough treatment of Bayesian asymptotics from the frequentist point of view. From now on the subscript 0 in various quantities will refer to the data-generating distribution.

Our strategy consists in proving that the posterior contraction rate for ν_0, given the sample 𝒵_n=(Z_1,…, Z_n), can be derived from the posterior contraction rate for q_0 given 𝒵_n, which is mathematically easier since Z_1,…, 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 λ of the Poisson process with a prior distribution.
  
  * Independently, endow the vector (p_1,…, p_m) with a Dirichlet distribution with parameter (α_1,…,α_m).
  
  * Set a priori p_k=0 for all k > m. 

This is a somewhat simplified version of the prior we used in Section <ref>, 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 ξ_1,…,ξ_m are independent Gamma distributed random variables, ξ_i ∼Gamma(α_i,1), then for η_i = ξ_i / ∑_j=1^m ξ_j, the vector (η_1,…,η_m) follows the Dirichlet distribution with parameter (α_1,…,α_m); furthermore, we have that ∑_j=1^m ξ_j ∼Gamma(∑_j=1^m α_j,1), and η_i are independent of ∑_j=1^m ξ_j.

In our asymptotic setting, we will make m=m_n dependent on n and let m_n→∞ at a suitable rate as n→∞.

Recall that we write Q=(q_k)_k∈ℕ_0 for q_k=Q({k}). Let Q denote the collection of all probability measures supported on ℕ.

	Suppose there exists  α, such that 0<α≤α_i ≤ 1 for all 1≤ i ≤ m_n. Suppose λ∼Gamma(a,b) with a∈ (0,1] and that ν_0 has a compact support. 
Then, for any γ >1,

    Π_n( ν - ν_0 _1 ≥log^γ n/√(n) | 𝒵_n ) → 0
 in Q_0^n-probability, as n→∞.



	Note that since the support of ν_0 is not assumed to be known, our CPP model is still naturally non-parametric. The assumption of the compact support of ν_0 does not cover the simulation example of Section <ref>. 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 <ref>. 






 §.§ Basic posterior inequality via the stability estimate

A key step of the proof of Theorem <ref> is the stability estimate in Equation (<ref>) below, which bounds the total variation distance between the Lévy measures ν, ν^' in terms of the total variation distance between the corresponding compound distributions q,q^'.


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 <cit.>, that will lead to the stability estimate we are after.

First we introduce some notation, and also recall a few useful facts summarised in <cit.>. The space of absolutely summable sequences is defined as
ℓ_1 { a∈^ℕ_0: ∑_j=0^∞ |a_j| < ∞},
with a norm given by
a_1 = ∑_j=0^∞ |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∈ℓ_1, we have the inequality

    a∗ b_1 ≤a_1 b_1,

where ∗ denotes convolution of a and b. We define a mapping a↦exp(a) from ℓ_1 into ℓ_1 via

    exp(a)=∑_j=0^∞a^∗ j/j!.


The exponential has the following two useful properties:

    exp(a+b)=exp(a)∗exp(b),    a,b∈ℓ_1,

and

    exp(a)=exp(b) ⟹ a=b,    a,b∈ℓ_1.

We define a sequence δ_0=(δ_0,k)_k∈ℕ_0, such that δ_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 (<ref>) in terms of convolutions of infinite sequences: q=exp(λ(p-δ_0)). Its convolution inverse, i.e. q^∗(-1) such that q^∗(-1)∗ q = δ_0, is given by r=q^∗(-1)=exp(-λ(p-δ_0)). Note that r∈ℓ_1. We have the following recursive expressions

    r_0 = 1/q_0,    r_k = - 1/q_0∑_j=1^k q_j r_k-j,    k∈ℕ.





Let q,q^' correspond to two pairs (λ,p) and (λ^',p^'), respectively (and r correspond to q, i.e. the pair (λ,p)). 
Then, in accordance with the notation introduced above and provided that q^'-q_1<r_1^-1, it holds that

    ν^' - ν_1 = λ^'p^' - λ p _1 ≤r_1  q^' - q _1/1 - r_1  q^' - q _1 .



The result is a direct consequence of Lemma 3 in <cit.>, which states that 

    (λ^'-λ)δ_0 + λ p - λ^'p^' = ∑_j=1^∞1/j(r∗(q-q^'))^∗ j

whenever q^'-q_1<r_1^-1.
Taking the ·_1-norm on both sides and some elementary bounding via (<ref>) imply that

    |λ^'-λ|+ λ^'p^' - λ p _1 ≤r_1  q^' - q _1/1 - r_1  q^' - q _1 ,

and thus Equation (<ref>) follows.

We will use Equation (<ref>) to establish the key inequality for the posterior measure Π(·|𝒵_n). We recall once again that the subscript 0 refers to `true', data-generating quantities.


For any prior Π on ν, for any ε∈ (0,1] and for any n≥ 1, the following posterior inequality holds:

    Π( ν - ν_0 _1 ≥ε | 𝒵_n )≤ 2 Π(  q - q_0 _1 ≥ε/2r_0_1 | 𝒵_n ).




Write {ν: ν - ν_0 _1 ≥ε} as a union of the sets

    {ν: ν - ν_0 _1 ≥ε}∩{ν: r_0_1  q - q_0 _1 < 1/2 }

and

    {ν: ν - ν_0 _1 ≥ε}∩{ν: r_0_1  q - q_0 _1 ≥  1/2 }.

Thanks to Lemma <ref>, the set

    {ν: ν - ν_0 _1 ≥ε}∩{ν: r_0_1  q - q_0 _1 < 1/2 }

is a subset of {ν:  q - q_0 _1 ≥ε/(2r_0_1) }. The proof is concluded by observing that {ν: ν - ν_0 _1 ≥ε}∩{ν: r_0_1  q - q_0 _1 ≥  1/2 } is a subset of {ν:  q - q_0 _1 ≥ε/(2r_0_1) }, too, since ε≤ 1. 

In general, stability estimates like the one in Equation (<ref>) 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., <cit.>.




 §.§ Proof of Theorem <ref>

The usefulness of Proposition <ref> above lies in the fact that the posterior contraction rate in the inverse problem of estimating the Lévy measure ν_0 from indirect observations 𝒵_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 ν_0 and is therefore fixed in the proofs). The general machinery developed in <cit.> can be applied to handle the latter, and also several inequalities obtained in <cit.> are useful in that respect.
In particular, we make use of the following inequality for the Hellinger distance, 

    h(q_, p, q_^', p^') ≤√(λ) h(p,p^') +|√(λ) - √(λ^')|,

Cf. Lemma 1 in <cit.>. 

To ease our notation, in the sequel we will often write  q and q^' instead of q_λ, p and q_λ^', p^', respectively.

Denote

    KL(q_0,q) = Q_0(logq_0/q),    V(q_0,q) = Q_0(logq_0/q)^2.

Another two inequalities we will use are the following: let λ, λ_0 ∈ [λ, λ]. Then there exists a positive constant C, such that 

    KL(q_0, q) ≤C(KL(p_0,p) +|λ_0- λ|^2), 
    
     V(q_0, q) ≤C(V(p_0,p)+ KL(p_0,p) +|λ_0- λ|^2);



cf. equations (14) and (15) in Lemma 1 in <cit.>.

These three inequalities can be obtained by adjustment of the arguments used in <cit.>. However, we opted to give their direct proofs in Lemma <ref> from Appendix  <ref> under slightly weaker conditions than required in <cit.>.

Our proof of Theorem <ref> proceeds via verification of the conditions for posterior contraction in Theorem 2.1 in <cit.>.  In our setting, the latter result reads  as follows.

Assume Z_n =(Z_1,…, Z_n), where Z_1,…, Z_n are independent and identically distributed with distribution q_0. Let h denote the Hellinger metric on Q, a collection of all measures with support in ℕ. 
Suppose that for a sequence {_n} with _n → 0 and n_n^2→∞, a constant C>0 and sets Q_n⊂Q, we have

    log N(_n, Q_n, h) ≤ n_n^2,


    Π_n(Q∖Q_n) ≤exp(-n_n^2(C+4)),


    Π_n( q KL(q_0,q) ≤_n^2,  V(q_0,q) ≤_n^2) ≥exp(-C n_n^2).

Then, for sufficiently large M>0, we have that Π_n (Q h(q,q_0)≥ M_n |Z_n ) → 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 {Λ_n}, {Λ_n}, and 
define the sieves

    Q_n ={ q_,p∈ [Λ_n, Λ_n], supp p ⊆{1,…,m_n}}.





  §.§.§ Entropy

We start with bounding the entropy of the sieve Q_n for h-balls of radius _n. 

	
Assume that as n→∞,

    m_n→∞,   ϵ_n → 0,   Λ_n → 0,   Λ_n →∞.

Then 

    log N(ϵ_n, Q_n, h) ≲ m_n {log (m_n) + log (Λ_n) +  log( 1/_n) } + log(1/Λ_n).






For , ' ≥Λ_n,

    |√()-√(')| = |-'|/√()+√(')≤1/2√(Λ_n) | -'|.

Furthermore, from Section 3.3 in <cit.>,

    h(p, p')   ≤√( p - p^'_1 )≤√( m_n  p-p' _∞).

Combining the preceding two displays and Equation (<ref>), we get 

    h(q_, p, q_', p') ≤√(Λ_n m_n p-p'_∞) + 1/2√(Λ_n) | -'|.

Hence, if

    p-p'_∞≤_n^2/4 Λ_n m_n,     |-'| ≤√(Λ_n)_n,

then  the Hellinger distance between q_, p and q_', p' is bounded by _n. 
To cover [Λ_n, Λ_n], we need at most ⌊Λ_n/2_n √(Λ_n)⌋ +1 intervals of length 2 √(Λ_n)_n.
To cover discrete distributions with support in {1,…,m_n}, we need at most 

    (  ⌊2Λ_n m_n /_n^2⌋ +1  )^m_n

L_∞-balls of radius _n^2 / (4Λ_n m_n). Under assumption (<ref>), 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 Q_n is of order

    ( Λ_n m_n /_n^2)^m_n×Λ_n/_n √(Λ_n) .

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




  §.§.§ Remaining prior mass

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

 For λ∼Gamma(a,b) with 0<a ≤ 1,
	
    Π_n(Q∖Q_n) ≲Λ_n^a - 1 e^- b Λ_n +  Λ_n .




We have (with a slight abuse of notation)
	
    Π_n(Q∖Q_n) = Π_n( [Λ_n, ∞)) +Π_n( [0,Λ_n)).

Now,

    Π_n( λ≥Λ_n ) =  b^a /Γ(a) ∫_Λ_n ^∞λ^a - 1 e^-b λλ≲Λ_n^a - 1 e^-b Λ_n.

Furthermore,

    Π_n( [0,Λ_n)) = b^a/Γ(a)  ∫_0^Λ_n λ^a - 1 e^- bλλ≲Λ_n^a.

The proof is concluded.





  §.§.§ 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 (0,∞) × (0,1) → (0,∞) by 
    g(,c) = C ^2/2[log(e/c)]^2,
 where C is the constant appearing in the statement of Lemma <ref> below.

Assume that

  
  * there exists  α, such that 0<α≤α_i ≤ 1 for all 1≤ i ≤ m_n;
  
  * strictly positive sequences p_n → 0, _n→ 0 and m_n→∞ satisfy the inequalities m_n  g(_n,p_n)<1 and p_n <g(_n,p_n)^2.

Define 

    B_n() = {q ∈Q_n   KL(q_0,q) ≤^2,  V(q_0,q) ≤^2}.

Then 

    Π_n(B_n(_n))  ≳Π_n(  |λ_0-λ| ≤_n)
    ×Γ(∑_i=1^m_nα_i)  exp(-m_nlog(1/(g(_n,p_n)^2-p_n)) - m_n log(1/α)).

Here _n = _n/√(3C), with a constant C > 0 not depending on n. 


Define 

    B_n() = {(λ, p)λ∈ [Λ_n, Λ_n], min_1≤ i ≤ m_n p_i ≥p_n, supp p ⊆{1,…,m_n}, 
      KL(p_0,p) ≤^2,  V(p_0,p) ≤^2,  |λ_0-λ| ≤} .

For all n large enough and  small, we have {λ: |λ_0-λ| ≤}⊆ [Λ_n, Λ_n]. Then by inequalities in Lemma <ref>,  B_n()⊂ B_n(√(3C)), with a constant C that can be taken the same for all large enough n; see the arguments in Section 4.2 in <cit.>. Hence, using the a priori independence of p and λ,

    Π_n(B_n(_n)) ≥Π_n(B_n(_n))  = Π_n(   |λ_0-λ| ≤_n) × U_n,
 where 

    U_n=  Π_n( {p   KL(p_0, p) ≤_n^2,    V(p_0,p) ≤_n^2,  min_1≤ i ≤ m_n p_i ≥p_n}).
 
Furthermore, by Lemma <ref> from Appendix <ref>, we have 

    U_n ≥Π_n({p ∑_i=1^m_n |p_0i-p_i| ≤ 2 g(_n, p_n),  min_1≤ i ≤ m_n p_i ≥p_n}).

The statement of the lemma now follows upon  applying Lemma <ref> from Appendix <ref> with η=p_n and =g(_n, p_n). 






  §.§.§ Using bounds in Theorem <ref>

We take

    m_n ≍log n,   _n ≍log^γ n/√(n),  p_n ≍1/n^2,
    Λ_n ≍log^2γ n,   Λ_n ≍exp(-const·log^2γn)

with appropriately selected proportionality constants, and verify the conditions in Theorem <ref>.

Firstly, condition (<ref>) is trivially satisfied. Therefore, we can invoke Lemma <ref> and conclude that the entropy is upper bounded by a multiple of log^2γ n, since γ > 1. Now log^2γ n ≲ n_n^2, and this verifies the entropy condition in Theorem <ref>.

Be Lemma <ref>, 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 <ref> we need that the term

    Π_n(  |λ_0-λ| ≤_n)  exp(-m_nlog(1/(g(_n,p_n)^2-p_n)) - m_n log(1/α))

is lower bounded by
 exp( - C n_n^2) for some large enough C>0. Now, Π_n(   |λ_0-λ| ≤_n)  ≍_n. Take  the logarithm on both sides of the above display and note that by our conditions 
 
    log( Π_n(  |λ_0-λ| ≤_n) ) ≳log( _n) ≳ -n_n^2.

 Likewise,
 
    m_nlog(1/(g(_n,p_n)^2-p_n)) + m_n log(1/α) ≲ n_n^2,

so that the prior mass condition holds.
 
Thus we have verified all the conditions of Theorem <ref>. The resulting posterior contraction rate is _n ≍log^γ n/√(n).




§ OUTLOOK



In this paper we introduced a non-parametric Bayesian approach to estimation of the Lévy measure ν of a discretely observed CPP, when the support of ν is a subset of ℕ. We constructed an algorithm for sampling from the posterior distribution of ν, and showed that in practice our procedure performs well and measures up to a benchmark frequentist plug-in approach from <cit.>. 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→∞, the posterior concentrates around the `true', data-generating distribution. The corresponding posterior contraction rate is the (nearly) optimal rate log^γn/√(n) for an arbitrary γ > 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 ℤ. The corresponding model has garnered substantial interest in financial applications, see <cit.>. 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.


italic-apa-good





§ TECHNICAL RESULTS







	Let q (resp. q^') be the law at time 1 of a compound Poisson process with intensity λ (resp. λ^') and jump distribution p (resp. p^'). Suppose that p and p^' are distributions concentrated on . Then,
	
    KL(q,q')   ≤λ KL(p,p')+λ'-λ+λlogλ/λ',
    
    V(q,q')   ≤ 2λ (V(p,p')+2KL(p,p'))+2KL(p,p')^2λ^2 
         +2λlog(λ/λ')(2(λ'-λ)+(λ-1)log(λ/λ'))+2(λ'-λ)^2,
    
    		h(q,q')   ≤√(λ)h(p,p')+√(1-e^-1/2(√(λ)-√(λ'))^2)≤√(λ)h(p,p')+|√(λ)-√(λ')|.

	In particular, if λ,λ^'∈[Λ, Λ] with 0<Λ≤Λ<∞, then there exists a positive constant C, that depends on Λ, Λ, such that
	
    KL(q,q^')   ≤C( KL(p,p^')+|λ - λ^'|),
    
    	V(q,q^')   ≤C(V(p,p')+KL(p,p')+(λ - λ^')^2),
    
    	h(q,q^')   ≤C h(p,p^') + |√(λ) - √(λ^')|.




	If KL(p,p^') and V(p,p^') 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≥ 1,
	
    KL(p^*n,p'^*n)   ≤ n KL(p,p'),
    
            V(p^*n,p'^*n)   ≤ n V(p,p') +4nKL(p,p')+n(n-1)KL(p,p')^2,
    
            h^2(p^*n,p'^*n)   ≤ n h^2(p,p').

	The assertions are trivial for n=1. Assuming that the first one holds for n-1 with n≥ 2, we will now show that it holds for n as well. Using the notation p^*n(i) for the ith element of p^*n and similarly in the case of p^'*n, we have
	
    KL(p^*n,p'^*n)   =∑_i∈ p^*n(i)log(p^*n(i)/p'^*n(i))
       =∑_i∈∑_k∈ p^*(n-1)(k) p(i-k)log(∑_k∈ p^*(n-1)(k) p(i-k)/∑_k∈ p'^*(n-1)(k) p'(i-k))
       ≤∑_i∈∑_k∈p^*(n-1)(k) p(i-k)log(p^*(n-1)(k) p(i-k)/p'^*(n-1)(k) p'(i-k)) 
       =∑_i∈∑_k∈ p^*(n-1)(k) p(i-k)(log(p^*(n-1)(k)/p'^*(n-1)(k))+log(p(i-k)/p'(i-k)))
       =KL(p^*(n-1),p'^*(n-1))+KL(p,p'),

	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
	
    ∑_k∈ p^*(n-1)(k)=1,   ∑_i∈p(i-k)=1,   ∀ k∈.

	By induction, we deduce that KL(p^*n,p'^*n)≤ n KL(p,p'). 
	
	The proof of the inequality for V is similar: we assume the inequality is true for n-1 with n≥ 2, and will show it holds for n as well. Write V(p^*n,p'^*n)=R_1+R_2 for
	
    R_1   =∑_i∈ p^*n(i)log^2(p^*n(i)/p'^*n(i)) _{p^*n(i)/p'^*n(i)≥ 1 },
    
    	R_2   =∑_i∈ p^*n(i)log^2(p^*n(i)/p'^*n(i)) _{p^*n(i)/p'^*n(i)<1 }.

	Observe that the function x ↦ (xlog^2x)_{ x≥ 1 } is convex. By Jensen's inequality we have for positive a_k,b_k that
	
    ( ∑_k a_k ) log^2( ∑_k a_k /∑_k b_k ) _{∑_k a_k ≥∑_k b_k }≤∑_k a_k log^2( a_k/b_k) _{ a_k/b_k ≥ 1 }.

	Using this inequality and
	
    p^*n(i) = ∑_k∈ℕ p^*(n-1)(k) p(i-k),    p'^*n(i) = ∑_k∈ℕ p'^*(n-1)(k) p'(i-k),

	we get that
	
    R_1
    	
    	   ≤∑_i∈∑_k∈ p(i-k)p^*(n-1)(k)log^2(p(i-k)p^*(n-1)(k)/p'(i-k)p'^*(n-1)(k))
       =∑_i∈∑_k∈ p(i-k)p^*(n-1)(k)(log(p(i-k)/p'(i-k))+log(p^*(n-1)(k)/p'^*(n-1)(k)))^2
       =V(p^*(n-1),p'^*(n-1))+V(p,p')+2KL(p,p')KL(p^*(n-1),p'^*(n-1))
       ≤ V(p^*(n-1),p'^*(n-1))+V(p,p') +2(n-1)KL(p,p')^2
       ≤ nV(p,p^') + 4(n-1)KL(p,p^') + n(n-1)KL(p,p^')^2,

	where in the last inequality we used the induction hypothesis.
	
	
	
	
	Now recall an elementary inequality
	
    e^-xx^2 ≤ 4 ( e^-x/2 - 1 )^2

	valid for x ≥ 0; see p. 12 in <cit.>. Applying this inequality to
	
    x = - log(p^*n(i)/p'^*n(i))

	such that p^*n(i)/p'^*n(i)<1, we get
	
    p^*n(i)/ p'^*n(i) log^2(p^*n(i)/p'^*n(i)) ≤ 4 ( √(p^*n(i)/p'^*n(i)) - 1 )^2.

	By multiplying both sides of the above inequality with p'^*n(i), summing the result through i and recalling the definition of the Hellinger distance, we get that
	
    R_2    ≤ 4 ∑_i∈ℕ( √( p^*n(i) ) - √( p'^*n(i) ))^2 = 4 h^2( p^*n,p'^*n) 
       ≤ 4KL(p^*n,p'^*n) ≤ 4 n KL(p,p').

	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
	
    h^2(p^*n,p'^*n)=∑_i∈ p^*n(i)g(p'^*n(i)/p^*n(i))

	for a convex function g(x)=(1-√(x))^2 _[0,∞)(x). Then, by the same reasoning as above, we have
	
    h^2(p^*n,p'^*n)   ≤∑_i∈∑_k∈ p^*(n-1)(k)p(i-k)g(p'^*(n-1)(k)p'(i-k)/p^*(n-1)(k)p(i-k))
       =∑_i∈∑_k∈((√(p^*(n-1)(k))-√(p'^*(n-1)(k)))√(p(i-k))
       =∑_i,k∈((+√(p'^*(n-1)(k))(√(p(i-k))-√(p'(i-k))))^2
       =h^2(p^*(n-1),p'^*(n-1))+h^2(p,p')
          +2 ∑_i∈∑_k∈(√(p^*(n-1)(k))-√(p'^*(n-1)(k)))
         +2∑_i,k∈(×(√(p(i-k))-√(p'(i-k)))√(p(i-k)p^*(n-1)(k)).

	Now note that the last summand satisfies
	
    2(∑_k∈√(p^'*(n-1)(k)p^*(n-1)(k))-1)(1-∑_k∈√(p(k)p'(k)))
     =-1/2h^2( p^*(n-1) , p^'*(n-1)(k) ) h^2(p,p^') ≤ 0.

	We therefore conclude that
	 
    h^2(p^*n,p'^*n)≤ h^2(p^*(n-1),p'^*(n-1))+h^2(p,p'),

	 which leads to the desired inequality, by an induction argument. 
	
	Step 2:  Now we prove the inequalities
	
    KL(q,q')   ≤∑_n=0^∞(N=n)KL(p^*n,p'^*n) +KL(N,N'),
    
    	V(q,q')   ≤2∑_n=0^∞(N=n)(V(p^*n,p'^*n)+2KL(p^*n,p'^*n)) +2V(N,N'),
    
    	h(q,q')   ≤√(∑_n=0^∞(N=n)h^2(p^*n,p'^*n))+h(N,N').

	Here N and N' are Poisson random variables with means λ and λ', respectively, and with a slight abuse of notation, KL(N,N'), V(N,N') and h(N,N') are the KL and V divergences and the Hellinger distance between the corresponding laws.
	
	Note that
	
    q(i) = ∑_n=0^∞ p^*n(i)P(N=n),    q^'(i) = ∑_n=0^∞ p^'*n(i)P(N'=n).

	Using this and the log-sum inequality,
	
    KL(q,q')   =∑_i∈ q(i)log(q(i)/q'(i))
    	
    	≤∑_i∈∑_n∈ p^*n(i) P (N=n)log( p^*n(i) P(N=n)/ p'^*n(i) P(N'=n))
       =∑_i∈∑_n∈  p^*n(i)P(N=n)(log( p^*n(i)/ p'^*n(i))+log(P(N=n)/P(N'=n)))
    
    	   =∑_n=0^∞ P(N=n)KL( p^*n, p'^*n) +KL(N,N').

	
	
	
	
	
	For the divergence V, write V(q,q')=B_1+B_2 for
	
    B_1   =∑_i∈ q(i) log^2(q(i)/q'(i))_{ q(i)/ q'(i)≥ 1}
       ≤∑_i∈∑_n∈ p^*n(i)(N=n)log^2(p^*n(i)(N=n)/p'^*n(i)(N'=n))
       ≤ 2 ∑_i∈∑_n∈ p^*n(i)(N=n)(log^2(p^*n(i)/p'^*n(i))+log^2((N=n)/(N'=n)))
       =2∑_n=0^∞ V(p^*n,p'^*n)(N=n)+2V(N,N'),
    
            B_2   =∑_i∈ q(i) log^2(q(i)/q'(i))_{ q(i)/ q'(i)<1}.

	To control B_2, we use the same arguments as in the proof of inequalities (12) and (15) in <cit.>, getting 

    B_2≤ 4KL(q,q')≤ 4 ∑_n=0^∞ KL(p^*n,p'^*n)(N=n).


This gives the required inequality for the V divergence.

	Finally, we prove the inequality for the Hellinger distance. Denoting the law of ∑_j=1^N Y_j by q̃, it holds by the triangle inequality that
	h(q,q')≤ h(q,q̃)+h(q̃,q'). Since g(x)=(1-√(x))^2_[0,∞)(x) is a convex function,
	
    h^2(q,q̃)   ≤∑_i∈∑_n∈ p^*n(i)(N=n) g(p'^*n(i)/p^*n(i))=∑_n∈(N=n)h^2( p^*n, p'^*n).

	It remains to prove that h(q̃,q')≤ h(N,N'). This again follows by convexity of g, since
	
    h^2(q̃,q')   =∑_i∈∑_n∈ p'^*n(i)(N=n)g(∑_n=0^∞ p'^*n(i)(N'=n)/∑_n=0^∞ p'^*n(i)(N=n))
       ≤∑_i∈∑_n∈ p'^*n(i)(N=n)g((N'=n)/(N=n))=h^2(N,N').

	
	Step 3: From Steps 1 and 2 we derive that
	
    KL(q,q')   ≤ KL(p,p') [N]+KL(N,N'),
    
    	V(q,q')   ≤ 2 [N](V(p,p')+2KL(p,p'))+2V(N,N')
          +2([N^2]-[N])KL(p,p')^2,
    
    	h(q,q')   ≤ h(p,p')√([N])+h(N,N').

	
	Now the proof of the lemma follows from these three inequalities upon noticing that KL(p,p')^2≤ V(p,p'), and recalling that
	
    [N]   =λ, 
    [N^2]   =λ+λ^2,
    
    	KL(N,N')     =λ'-λ+λlogλ/λ'≲ (λ-λ')^2,
    
    	h^2(N,N')   =1-e^-1/2(√(λ)-√(λ'))^2≤ |√(λ)-√(λ')|,
    
    	V(N,N')   =λlog(λ/λ')(2(λ'-λ)+(λ-1)log(λ/λ'))+(λ'-λ)^2≲ (λ-λ')^2.




	Let >0. 
	Suppose p=(p_1,…, p_m) and  p^'=(p^'_1,…, p^'_m) are points in the m-dimensional unit simplex, and let min_1≤ i ≤ m p_i ≥ c for some c ∈ (0,1). Then there exists a universal constant C>0, such that the inequality
	
    ∑_i=1^m | p^'_i-p_i| ≤ C ^2/[log (e/c)]^2

	implies that KL(p^', p) ≤^2 and V(p^',p) ≤^2 hold. 


	Lemma 8 in <cit.> assures that there exists a  constant C (not depending on either p or p^'), such that 
	
    KL(p^', p) ≤C h^2(p^',p) [ 1 + log( p^'/p_∞) ]

	and
	
    V(p^', p) ≤C h^2(p^',p)  [ 1 + log( p^'/p_∞) ]^2.

	From Section 3.3 in <cit.> we have h^2(p^',p)≤∑_i=1^m |p^'_i-p_i|. 
	
	
	
	
	Furthermore, since 0<c<1 and min_1≤ i ≤ m p_i ≥ c,
	
    1 ≤ 1 + log( p^'/p_∞)  ≤ 1  +log (1/c)=log (e/c).

	Therefore,
	
    max(KL(p^',p), V(p^',p)) ≤C[log(e/c)]^2 ∑_i=1^m |p^'_i-p_i| ,

	from which the assertion of the lemma follows trivially.





	Let m≥ 2 be an integer.	Suppose (p_1,…,p_m) ∼Dir(α_1,…, α_m). Let p_0=(p_01,…, p_0m) be an arbitrary point in the m-dimensional unit simplex. Assume there exists  α, such that 0<α≤α_i ≤ 1 for all 1≤ i ≤ m. Let >0, and let η be such that η<^2. Then if m<1,
	
    Π_n(∑_i=1^m |p_i - p_0i| ≤ 2 , min_1≤ i ≤ m p_i≥η)  
               ≥Γ(∑_i=1^m α_i)  exp(-mlog(1/(^2-η)) - m log(1/α)).



	By arguments analogous to those in the proofs of Lemma 6.1 in <cit.> and Lemma 10 in <cit.>, we obtain that the left-hand side of (<ref>) can be lower bounded by 
	
    Γ(∑_i=1^m α_i)/∏_i=1^m Γ(α_i)∏_i=1^m-1∫_max(p_0i-^2, η^2)^min(p_0i+^2,1) x_i^α_i-1 x_i.

	The length of the integration interval in each of the integrals in the above product is lower bounded by ^2-η. Using that α≤α_i≤ 1, we deduce that the preceding display is lower bounded by 
	
    Γ(∑_i=1^m α_i) α^m exp(-(m-1)log(1/^2-η)).

	This entails the result. 






§ 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 <ref>. 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 ν. Finally, for each k, the α/2-th and (1-α/2)-th sample quantiles were obtained to yield 1-α bootstrap confidence intervals for ν_k. The results with α=0.05 are displayed in Figure <ref>. 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.


