Séminaire Lotharingien de Combinatoire 67 (2015), Article B67a

SYMBOLIC CALCULUS IN MATHEMATICAL STATISTICS: A REVIEW

ELVIRA DI NARDO

Abstract. In the last ten years, the employment of symbolic methods has substantially extended both the theory and the applications of statistics and probability. This survey reviews the development of a symbolic technique arising from classical umbral calculus, as introduced by Rota and Taylor in 1994. The usefulness of this symbolic technique is twofold. The first is to show how new algebraic identities drive in discovering insights among topics apparently very far from each other and related to probability and statistics. One of the main tools is a formal generalization of the convolution of identical probability distributions, which allows us to employ compound Poisson random variables in various topics that are only somewhat interrelated. Having got a different and deeper viewpoint, the second goal is to show how to set up algorithmic processes performing efficiently algebraic calculations. In particular, the challenge of finding these symbolic procedures should lead to a new method, and it poses new problems involving both computational and conceptual issues. Evidence of efficiency in applying this symbolic method will be shown within statistical inference, parameter estimation, Lévy processes, and, more generally, problems involving multivariate functions. The symbolic representation of Sheffer polynomial sequences allows us to carry out a unifying theory of classical, Boolean and free cumulants. Recent connections within random matrices have extended the applications of the symbolic method.

Contents

1. Introduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2. The symbolic method of moments . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.1. Moments of sampling distributions . . . . . . . . . . . . . . . . . . . 9 2.2. Sheppard’s corrections. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

3. Compound Poisson processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.1. Time-space harmonic polynomials . . . . . . . . . . . . . . . . . . . . 23

4. Cumulants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 4.1. Randomized compound Poisson r.v.’s. . . . . . . . . . . . . . . . . 30 4.2. Lévy processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

5. Sheffer polynomial sequences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 5.1. Solving linear recurrences . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 5.2. Lagrange inversion formula . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

5.3. Generalized Bell polynomials . . . . . . . . . . . . . . . . . . . . . . . . . 40 6. Fast symbolic computation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

6.1. k-statistics and polykays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 6.2. Computing cumulants through generalized Abel

polynomials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

7. Multivariate Faà di Bruno formula. . . . . . . . . . . . . . . . . . . . . . . . . 53 7.1. Multivariate cumulants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 7.2. Multivariate Lévy processes . . . . . . . . . . . . . . . . . . . . . . . . . . 58 7.3. Multivariate Hermite polynomials.. . . . . . . . . . . . . . . . . . . . 59 7.4. Multivariate Bernoulli polynomials.. . . . . . . . . . . . . . . . . . . 61

8. Appendix 1: Tables of computation times. . . . . . . . . . . . . . . . . . 62 9. Appendix 2: Families of time-space harmonic polynomials . 64 10. Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 References. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

1. Introduction

Most practical problems do not have exact solutions, but numerical computing helps in finding solutions to desired precision. Numerical methods give data that can be analyzed to reveal trends or approximations. In contrast to numerical methods, symbolic methods treat objects that are either formal expressions, or are algebraic in nature. In many cases, numerical methods will not give sufficient information about the nature of the problem; this is where symbolic methods fit in. When used properly, they can give us more insight into the problem we are trying to solve. Another advantage is that, while numerical methods may simply fail to compute correct results, symbolic methods yield closed or explicit formulas. Although the term “symbolic method” is also used in different contexts (see for example [42]), here we refer to a set of manipulation techniques aimed at performing algebraic calculations, preferably through an algorithmic approach, in order to find efficient mechanical processes that may be passed to a computer. This is usually called “symbolic computation”. Despite the fact that its development took place later than that of numerical or graphical algorithms, symbolic computation has had a comparable impact, especially in theoretical and applied statistical research. Topics range from asymptotic expansions to Edgeworth series, from likelihood functions to saddlepoint approximations.

Since symbolic computation is exact, its accuracy is not a fundamental issue, as it is for numerical methods. Replacing pencil and paper by a computer is not the only aim of symbolic methods. Computational cost is a question naturally arising when problems of complex nature are addressed. The advantage of symbolic methods is that, by approaching a problem from a different mathematical point of view, efficiency is sometimes a by-product.

The symbolic method we are going to review arises from the socalled “classical umbral calculus”, a creation of Gian-Carlo Rota, who spent most of his scientific life looking for an elegant and as general as possible formalization of this calculus.

It is certainly due to Rota’s efforts that the heuristic method, invented and largely employed by Rev. John Blissard between 1861 and 1868 (see references in [33]), now has a solid foundation and widespread application within the mathematical community. Despite the impressive amount of scientific literature following the introduction of this calculus, Rota, in 1994, had the courage to turn his umbral calculus upside down, outlining what he considered to be a new and correct syntax for this matter. The classical umbral calculus we use consists essentially in a moment calculus, since its basic device is to represent a unital sequence of numbers by symbols α, called umbra, i.e., to associate the sequence 1, a_{1}, a_{2}, … to the sequence 1, α, α^{2}, … of powers of α via an operator 𝔼 that looks like the expectation of random variables (r.v.’s) [35, 75]. The i-th element of the sequence is then called the ith moment of α. This symbolic approach shares with free probability [58, 90] the use of moments as a tool to characterize r.v.’s. An analog of freeness^{1} is defined in a manner resembling independence of tensors. In contrast to free probability, where different linear operators may be applied to the same non-commutative r.v., in the symbolic method of moments just one operator is employed, but the same sequence of moments may correspond to more than one symbol. For a pair of commutative r.v.’s, freeness is equivalent to at least one of them having vanishing variance. So freeness is a pure non-commutative device, and the symbolic approach we are going to review can be considered as the natural analog of a non-commutative field.

“As sometimes happens in the practice of the mathematical investigation, the subject we deal with here does not develop the original idea from which our research started in the spring of 1997, but this paper is closely related to it. In that period, Gian-Carlo Rota was visiting professor at the University of Basilicata and, during a conversation just before he left, he shared with us his keen interest in a research project: to develop a combinatorial random variable theory. The delicate question arising on the underlying foundational side and the short time available led us to continue the discussion via email, weaving it with different activities over a period of several months. The following year, Gian-Carlo Rota held his last course in Cortona; we seized the opportunity to spend some time with him. We resumed the thread of our conversations and presented him with the doubts that gradually took hold of us. As usually, his contribution disclosed new horizons that have led us to write these pages.” This quotation is a part of the introduction of the paper “Umbral nature of the Poisson random variables” written by the author together with Domenico Senato and published in a volume dedicated to Gian-Carlo Rota after his death [33]. This paper was the first of a substantial list of articles devoted to the applications of classical umbral calculus in probability and statistics. In these papers, more definitions have been added, and the method has been restyled according to the needs faced when tackling some particular problems. The more these new tools have been developed, the closer the method has come to the method of moments employed in the study of random matrices, which the author believes to be a natural milestone of this theory and at the same time a good starting point for further developments and applications. This is why the version here reviewed is called symbolic method of moments in a talk [18] on new challenges and future developments of mathematical statistics given by the author in Dublin in 2011. Two papers have to be mentioned in support of this personal belief: a paper published in the Annals of Statistics devoted to the computation of spectral statistics for random matrices [26] and a second one, published in the Journal of Multivariate Analysis, where the symbolic method of moments allows us to compute very general joint moments of non-central Wishart distributions [19]. The author chose to not include these topics in this review, since the overall setting is not yet completely settled, and a serious rethinking of the matter in terms of cumulants^{2} instead of moments seems necessary.

1In free probability, free r.v.’s correspond to classical independent r.v.’s.

The first goal of this paper is to give a general overview of how the method has been developed and where it has been used with success. While sticking to the original idea, the current version of the method offers a much richer set of tools. Therefore the introduction of these tools is done gradually along the paper according to the outline of applications. The second goal is to prove its effectiveness in some applications that over the years have contributed to refining the theory. An example particularly meaningful is the computation of k-statistics (Section 6.1), a challenging problem since their first introduction in 1929 as unbiased estimators of cumulants [41]. Many authors have proposed different techniques aiming at performing this computation in a reasonable time. By using the symbolic method here reviewed, an efficient computation of k-statistics can be performed by means of a suitable generalization of randomized compound Poisson r.v.’s.

The paper is organized as follows. Notation is introduced in Section 2 using the semantics of r.v.’s. In this way, the intention is to permit the reader with no prior knowledge of symbolic calculations to become comfortable with them. Section 2.1 shows how to perform these computations in dealing with moments of sampling distributions, reducing the underlying combinatorics of symmetric functions to a few relations covering a great variety of calculations. Section 2.2 is devoted to a special statistical device, called Sheppard’s correction, useful in correcting estimations of moments when sample data are grouped into classes. The resulting closed-form formula gives an efficient algorithm to perform computations. In the same section, its multivariate version is sketched, showing what we believe is one of the most interesting feature of this symbolic method: its multivariate generalization. More details are given in Section 7, where a symbolic version of the multivariate Faà di Bruno formula is also introduced. The umbral calculus is particularly suited not only in dealing with number sequences, but also with polynomial sequences. Time-space harmonic polynomials are an example in which the symbolic method allows us to simplify the proofs of many properties, as well as their computation. Introduced in Section 3, these polynomials generate stochastic processes that result in a special family of martingales employed in mathematical finance. Time-space harmonic polynomials rely on the symbolic representation of Lévy processes, as given in Section 4.2, obtained by using a generalization of compound Poisson processes, as given in Section 4.1. A Lévy process is better described by its cumulant sequences than by its moments, as happens for many other r.v.’s, the Gaussian or the Poisson r.v.’s. Section 4 is devoted to introducing umbrae whose moments are cumulants of a given sequence of numbers. Besides classical cumulants, quite recently [90] new families of cumulants have been introduced, the Boolean and the free cumulants by using the algebra of multiplicative functions on suitable lattices. This link is exploited in Section 6.2, where the computational efficiency of the symbolic method of moments is highlighted. In particular, by means of the symbolic representation of Abel polynomials, all families of cumulants can be parametrized by the same formula. This parametrization allows us to write an algorithm giving moments in terms of classical, Boolean and free cumulants, and vice-versa. Other topics on Sheffer sequences are treated in Section 5, including Lagrange inversion and Riordan arrays, and solutions of difference equations. Section 6.1 shows how to compute unbiased estimators of cumulants by using exponential polynomials, which are special Sheffer sequences. The resulting algorithm is faster than others given in the literature, and it has been generalized to the multivariate case. In the last section, the multivariate umbral calculus is introduced, and various applications are given involving special families of multivariate polynomials.

2For the definition of cumulants of a number sequence, see Section 4 and in particular Equation (4.1).

This text is the written version of lectures given at the 67-th Seminaire Lotharingien de Combinatoire in September 2011; joint session with XVII Incontro Italiano di Combinatoria Algebrica. Due to the length and to avoid making dull reading, the formalism is kept to the minimum, and references are given step by step for those wishing to go into details and to understand the technical proofs of formal matters.

2. The symbolic method of moments

The method we are going to introduce in this section consists of shifting powers a^{i} to indexes a_{i} by using a linear functional. The idea of dealing with subscripts as they were powers was extensively used by the mathematical community since the nineteenth century, although with no rigorous foundation. The method, sometimes called Blissard’s symbolic method, is often attributed to J. Sylvester. He was the first who introduced the term umbral from latin to denote the shadowy techniques used to prove certain polynomial equations. Roman and Rota have equipped these “tricks” with the right formalism in 1978, by using the theory of linear operators [71]. Despite the impressive amount of applications following this rigorous foundation (see references in [33]), in 1994 Rota proposed a new version of the method in a paper entitled “The classical umbral calculus” [75]. The reason of this rethinking is well summarized in the following quotation [75]: “Although the notation of Hopf algebra satisfied the most ardent advocate of spic-and-span rigor, the translation of classical umbral calculus into the newly found rigorous language made the method altogether unwieldy and unmanageable. Not only was the eerie feeling of witchcraft lost in the translation, but, after such a translation, the use of calculus to simplify computation and sharpen our intuition was lost by the wayside.” The term classical is strictly related to which sequence is chosen as the unity of the calculus. Rota has always referred to the sequence {1}, and the same is done in the version we are going to introduce, that from now on will be called symbolic method of moments.

Let ℝ be the real field^{3}. The symbolic method of moments consists of a set 𝒜 = {α, γ, …} of elements called umbrae, and a linear functional 𝔼 : ℝ[𝒜] → ℝ, called evaluation, such that 𝔼[1] = 1 and

(2.1)  𝔼[α^{i}γ^{j} ⋯] = 𝔼[α^{i}]E[γ^{j}] ⋯ (uncorrelation property)

for non-negative integers i, j, . . ..

Definition 2.1. [33] The sequence {1, a_{1}, a_{2}, …} ∈ ℝ is said to be umbrally represented by an umbra α if 𝔼[α^{i}] = a_{i} for all positive integers i.

In the following, in place of {1, a_{1}, a_{2}, …} we will use the notation {a_{i}}_{i≥0} assuming a_{0} = 1. Special umbrae are given in Table 1.

3The umbral calculus given in [75] considers a commutative integral domain whose quotient field is of characteristic zero. For the applications given in the following, the real field ℝ is sufficient.

4The i-th Bell number is the number of partitions of a finite non-empty set with i elements or the i-th coefficient times i! in the Taylor series expansion of the function exp(e^{t} − 1).

Table 1. Special umbrae

Umbrae  Moments for positive integers i
Augmentation umbra 𝔼[ε^{i}] = 0
Unity umbra  𝔼[u^{i}] = 1
Singleton umbra  𝔼[χ^{i}] = δ_{i,1}, the Kronecker delta
Bell umbra  𝔼[β^{i}] = 𝔅_{i}, the i-th Bell number^{4}
Bernoulli umbra  𝔼[ι^{i}] = B_{i}, the i-th Bernoulli number^{5}
Euler umbra  𝔼[ξ^{i}] = 𝔈_{i}, the i-th Euler number^{6}

Associating a measure with a sequence of numbers is familiar in probability theory, when the i-th term of a sequence can be considered as the i-th moment of a r.v., under suitable hypotheses (the so-called Hamburger moment problem [94]). As Rota underlines in Problem 1: the algebra of probability [73], all of probability theory could be done in terms of r.v.’s alone by taking an ordered commutative algebra over the reals, and endowing it with a positive linear functional. Following this analogy, a_{i} is called the i-th moment of the umbra α.

Definition 2.2. A r.v. X is said to be represented by an umbra α if its sequence of moments {a_{i}}_{i≥0} is umbrally represented by α.

In order to avoid misunderstandings, the expectation of a r.v. X will be denoted by E, and its i-th moment by E[X^{i}].

Example 2.3. If P (X = 0) = 1, the r.v. X is represented by the augmentation umbra ε. If P (X = 1) = 1, the r.v. X is represented by the unity umbra u. The Poisson r.v. Po(1) of parameter 1 is represented by the Bell umbra β. More examples of umbrae representing classical r.v.’s can be found in [28].

Not all r.v.’s can be represented by umbrae. For example, the Cauchy r.v. does not admit moments. Moreover, not all umbrae represent r.v.’s. For example, the singleton umbra in Table 1 does not represent a r.v. since its variance is negative 𝔼[χ^{2}] − 𝔼[χ]^{2} = −1, even though this umbra will play a fundamental role in the whole theory. If {a_{i}}_{i≥0}, {g_{i}}_{i≥0} and {b_{i}}_{i≥0} are sequences umbrally represented by the umbrae α, γ and ζ, respectively, then the sequence {h_{i}}_{i≥0} with

h_{i} = ∑_{k_{1}+k_{2}+···+k_{m}=i} \binom{i}{k_{1}, k_{2}, …, k_{m}}a_{k_{1}} g_{k_{2}} ⋯ b_{k_{m}}

5Many characterizations of Bernoulli numbers can be found in the literature, and each one may be used to define these numbers. Here we refer to the sequence of numbers satisfying B_{0} = 1 and ∑_{k=0}^{i} \binom{i+1}{k}B_{k} = 0 for i = 1, 2, ….

6The Euler numbers are the coefficients of the formal power series 2e^{z}/[e^{2z} + 1].

is umbrally represented by \underbrace{α + γ + ⋯ + ζ}_{m}, that is,

(2.2)

𝔼[(\underbrace{α + γ + ⋯ + ζ}_{m})^{i}] = ∑_{k_{1}+k_{2}+···+k_{m}=i} \binom{i}{k_{1}, k_{2}, …, k_{m}}a_{k_{1}} g_{k_{2}} ⋯ b_{k_{m}}.

The right-hand side of Equation (2.2) represents the i-th moment of a summation of independent r.v.’s.

The second fundamental device of the symbolic method of moments is to represent the same sequence of moments by different umbrae.

Definition 2.4. Two umbrae α and γ are said to be similar if and only if 𝔼[α^{i}] = 𝔼[γ^{i}] for all positive integers i. In symbols α ≡ γ.

If we replace the set {α, γ, …, ζ} by a set of m distinct and similar umbrae {α, α′, …, α′′}, Equation (2.2) gives

(2.3)  𝔼[(\underbrace{α + α′ + ⋯ + α′′}_{m})^{i}] = ∑_{λ⊢i}(m)_{ν_{λ}}d_{λ}a_{λ},

where the summation is over all partitions^{7} λ of the integer i, a_{λ} = a_{1}^{r_{1}} a_{2}^{r_{2}} ⋯ and

d_{λ} = \frac{i!}{(1!)^{r_{1}}(2!)^{r_{2}} ⋯ r_{1}! r_{2}! ⋯}

is the number of i-set partitions with block sizes given by the parts of λ. The new symbol m.α denotes the summation α + α′ + ⋯ + α′′ and is called dot-product of m and α. This new symbol is referred to as auxiliary umbra, in order to underline that it is not in 𝒜 . So Equation (2.3) may be rewritten as

(2.4)  𝔼[(m.α)^{i}] = ∑_{λ⊢i}(m)_{ν_{λ}}d_{λ}a_{λ}.

The introduction of new symbols is not only a way of lightening notation, but also to represent sequences of moments obtained by suitable generalizations of the right-hand side of Equation (2.4). This will become clearer as we move along further. According to Definition 2.2, the auxiliary umbra m.α represents a summation of m independent and identically distributed (i.i.d.) r.v.’s. Indeed, identically distributed r.v.’s share the same sequence of moments, if they have.

If we denote by 𝒳 the alphabet of auxiliary umbrae, the evaluation 𝔼 may be defined using the alphabet 𝒜 ∪ 𝒳 , and so we can deal with auxiliary umbrae as they were elements of the base alphabet [75]. In the following, we still refer to the alphabet 𝒜 as the overall alphabet of umbrae, assuming that auxiliary umbrae are included.

7Recall that a partition of an integer i is a sequence λ = (λ_{1}, λ_{2}, …, λ_{t}), where λ_{j} are weakly decreasing integers and ∑_{j=1}^{t} λ_{j} = i. The integers λ_{j} are called parts of λ. The length of λ is the number of its parts and will be indicated by ν_{λ}. A different notation is λ = (1^{r_{1}}, 2^{r_{2}}, …), where r_{j} is the number of parts of λ equal to j and r_{1} + r_{2} + ⋯ = ν_{λ}. We use the classical notation λ ⊢ i to denote that “λ is a partition of i”.

Summations of i.i.d.r.v.’s appear in many sample statistics. So the first application we give in statistics concerns sampling distributions.

2.1. Moments of sampling distributions. It is recognized that an appropriate choice of language and notation can simplify and clarify many statistical calculations [54]. In managing algebraic expressions such as the variance of sample means, or, more generally, moments of sampling distributions, the main difficulty is the manual computation. Symbolic computation has removed many of such difficulties. The need of efficient algorithms has increased the attention on how to speed up the computational time. Many of these algorithms rely on the algebra of symmetric functions [1]: U-statistics are an example which arise naturally when looking for minimum-variance unbiased estimators. A U-statistic of a random sample has the form

(2.5)  U = \frac{1}{(n)_{k}} ∑ Φ(X_{j_{1}}, X_{j_{2}}, …, X_{j_{k}}),

where X_{1}, X_{2}, …, X_{n} are n independent r.v.’s, and the sum ranges over the set of all permutations (j_{1}, j_{2}, …, j_{k}) of k integers chosen in [n] := {1, 2, …, n}. If X_{1}, X_{2}, …, X_{n} are identically distributed r.v.’s with common cumulative distribution function F (x), U is an unbiased estimator of the population parameter

Θ(F) = ∫ ∫ ⋯ ∫ Φ(x_{1}, x_{2}, …, x_{k}) dF (x_{1}) dF (x_{2}) ⋯ dF (x_{k}).

In this case, the function Φ may be assumed to be a symmetric function of its arguments. Often, in applications, Φ is a polynomial in the X_{j}’s, so that U-statistics are symmetric polynomials. By virtue of the fundamental theorem of symmetric polynomials, U-statistics can be expressed in terms of elementary symmetric polynomials. The symbolic method proposed here provides a way to find this expression [22].

The starting point is to replace ℝ by ℝ[x_{1}, x_{2}, …, x_{n}], where x_{1}, x_{2}, …, x_{n} are indeterminates. This device allows us to deal with multivariable umbral polynomials p ∈ ℝ[x_{1}, x_{2}, …, x_{n}]. The uncorrelation property (2.1) is then updated to

𝔼[x_{i}^{k_{1}}x_{j}^{k_{2}} ⋯ α^{k_{3}}γ^{k_{4}} ⋯] = x_{i}^{k_{1}}x_{j}^{k_{2}} ⋯ 𝔼[α^{k_{3}}]𝔼[γ^{k_{4}}] ⋯

for any set of distinct umbrae in 𝒜 , for i, j ∈ [n], with i ≠ j, and for all positive integers k_{1}, k_{2}, k_{3}, k_{4}.

For umbral polynomials p, q ∈ ℝ[x_{1}, x_{2}, …, x_{n}][𝒜], we may relax the similarity equivalence and introduce another umbral equivalence:

p ≃ q if and only if 𝔼[p] = 𝔼[q].

A polynomial sequence {p_{i}}_{i≥0} ∈ ℝ[x_{1}, x_{2}, …, x_{n}], with p_{0} = 1 and degree(p_{i}) = i, is represented by an umbra ψ, called polynomial umbra, if 𝔼[ψ^{i}] = p_{i} for all positive integers i. The classical bases of the algebra of symmetric polynomials can be represented by polynomial umbrae [21]. For example, the i-th elementary symmetric polynomial e_{i}(x_{1}, x_{2}, …, x_{n}) satisfies

(2.6) e_{i}(x_{1}, x_{2}, …, x_{n}) ≃ \frac{(χ_{1}x_{1} + χ_{2}x_{2} + ⋯ + χ_{n}x_{n})^{i}}{i!} , for i ≤ n,

with χ_{1}, χ_{2}, …, χ_{n} uncorrelated singleton umbrae. If the indeterminates x_{1}, x_{2}, …, x_{n} are replaced by n uncorrelated umbrae α_{1}, α_{2}, …, α_{n} similar to an umbra α, then Equivalence (2.6) becomes

(2.7)  e_{i}(α_{1}, α_{2}, …, α_{n}) ≃ \frac{[n.(χα)]^{i}}{i!} , for i ≤ n.

The umbral polynomial e_{i}(α_{1}, α_{2}, …, α_{n}) is called the i-th umbral elementary symmetric polynomial. The umbrae {α_{1}, α_{2}, …, α_{n}} in (2.7) may also be replaced by some powers such as {α_{1}^{j}, α_{2}^{j}, …, α_{n}^{j}}. Let us consider the monomial symmetric polynomial

𝔪_{λ}(x_{1}, x_{2}, …, x_{n}) = ∑ x_{1}^{λ_{1}}x_{2}^{λ_{2}} ⋯ x_{n}^{λ_{n}}

with λ ⊢ i ≤ n, where the notation is symbolic and must be interpreted in the way that all different images of the monomial x_{1}^{λ_{1}}x_{2}^{λ_{2}} ⋯ x_{n}^{λ_{n}} under permutations of the variables have to be summed. If the indeterminates x_{1}, x_{2}, …, x_{n} are replaced by n uncorrelated umbrae α_{1}, α_{2}, …, α_{n} similar to an umbra α, then

(2.8)  𝔪_{λ}(α_{1}, α_{2}, …, α_{n}) ≃ \frac{[n.(χα)]^{r_{1}}}{r_{1}!} \frac{[n.(χα^{2})]^{r_{2}}}{r_{2}!} ⋯,

which is a product of umbral elementary polynomials in (2.7). In many statistical calculations, the so-called augmented monomial symmetric polynomials [94] are involved:

𝔪̃_{λ}(x_{1}, x_{2}, …, x_{n})

= ∑_{j_{1}≠...≠j_{r_{1}}≠j_{r_{1}+1}≠...≠j_{r_{1}+r_{2}}≠···}x_{j_{1}} ⋯ x_{j_{r_{1}}} x_{j_{r_{1}+1}}^{2} ⋯ x_{j_{r_{1}+r_{2}}}^{2} · · · .

These polynomials are obtained from 𝔪_{λ}(x_{1}, x_{2}, …, x_{n}) as follows:

(2.9)  𝔪_{λ}(x_{1}, x_{2}, …, x_{n}) = \frac{𝔪̃_{λ}(x_{1}, x_{2}, …, x_{n})}{r_{1}!r_{2}! ⋯}.

Then, from Equivalence (2.8) and Equation (2.9), we have

(2.10)  𝔪̃_{λ}(α_{1}, α_{2}, …, α_{n}) ≃ [n.(χα)]^{r_{1}}[n.(χα^{2})]^{r_{2}} ⋯,

which is again a product of umbral elementary polynomials (2.7).

Example 2.5. If λ = (1^{2}, 3) and n ≥ 5 then

[n.(χα)]^{2}[n.(χα^{3})] = (χ_{1}α_{1} + ⋯ + χ_{n}α_{n})^{2}(χ_{1}α_{1}^{3} + ⋯ + χ_{n}α_{n}^{3})

= ( ∑_{1≤j_{1}≠j_{2}≤n} χ_{j_{1}}α_{j_{1}}χ_{j_{2}}α_{j_{2}}) ( ∑_{1≤j_{3}≤n} χ_{j_{3}}α_{j_{3}}^{3}) ≃ ∑_{1≤j_{1}≠j_{2}≠j_{3}≤n} α_{j_{1}}α_{j_{2}}α_{j_{3}}^{3}.

This last equivalence follows by observing that 𝔼[χ_{j_{1}}χ_{j_{2}}χ_{j_{3}}] vanishes whenever there is at least one pair of equal indices among {j_{1}, j_{2}, j_{3}}.

More details on symmetric polynomials and polynomial umbrae are given in [21]. Due to Equivalence (2.10), the fundamental expectation result, at the core of unbiased estimation of moments, may be restated as follows.

Theorem 2.6. [21] If λ = (1^{r_{1}}, 2^{r_{2}}, …) ⊢ i ≤ n, then

𝔼 {[n.(χα)]^{r_{1}}[n.(χα^{2})]^{r_{2}} ⋯ } = (n)_{ν_{λ}}a_{λ}.

From Equation (2.4) and Theorem 2.6, the next corollary follows.

Corollary 2.7. If i ≤ n, then

(2.11)  (m.α)^{i} ≃ ∑_{λ⊢i} \frac{(m)_{ν_{λ}}}{(n)_{ν_{λ}}} d_{λ} [n.(χα)]^{r_{1}}[n.(χα^{2})]^{r_{2}} · · · .

From Equivalence (2.11), if m = n, then

(2.12)  (n.α)^{i} ≃ ∑_{λ⊢i} d_{λ}[n.(χα)]^{r_{1}}[n.(χα^{2})]^{r_{2}} · · · .

Theorem 2.6 discloses a more general result: products of moments are moments of products of umbral elementary symmetric polynomials (2.7).

By analogy with Equation (2.5), the umbral symmetric polynomial on the right-hand side of Equivalence (2.11),

\frac{1}{(n)_{ν_{λ}}} [n.(χα)]^{r_{1}}[n.(χα^{2})]^{r_{2}} ⋯,

is called the U-statistic of uncorrelated and similar umbrae α_{1}, α_{2},…, α_{n}. Then Theorem 2.6 says that U-statistics can be expressed as products of umbral elementary symmetric polynomials, and Corollary 2.7 says that any statistic involving a summation of i.i.d.r.v.’s can be expressed as a linear combination of U-statistics.

A more general problem consists in computing population moments of sample statistics. These are symmetric functions of not independent r.v.’s. Within the symbolic method of moments, the dependence between r.v.’s corresponds to working with umbral monomials with non-disjoint supports. Indeed, the support of an umbral polynomial p is the set of all umbrae occurring in p. If μ and ν are umbral monomials with non-disjoint supports, then 𝔼[μ^{i} ν^{j}] ≠ 𝔼[μ^{i}]𝔼[ν^{j}] for all positive integers i and j. In the following, unless otherwise specified, we work with umbral monomials with non-disjoint supports and give an example of how to use these monomials to perform computations with not independent r.v.’s.

Example 2.8. [22] Let us consider a bivariate random sample of r.v.’s (X_{1}, Y_{1}), …, (X_{n}, Y_{n}), that is, (X_{1}, …, X_{n}) is a random sample with parent distribution X not necessarily independent from the parent distribution Y of the random sample (Y_{1}, …, Y_{n}). In order to compute the expectation of

(2.13)  (∑_{i=1}^{n} X_{i})^{2} (∑_{i=1}^{n} Y_{i}) or (∑_{i≠j}^{n} X_{i}^{2}X_{j}) (∑_{i=1}^{n} X_{i}^{2}Y_{i})^{2}

by the symbolic method of moments, the main tool is to represent the sample statistics in (2.13) as a suitable linear combination of umbral monomials μ_{1} and μ_{2} representing the populations X and Y, respectively. For the former sample statistic (2.13), we have

E [(∑_{i=1}^{n} X_{i})^{2} (∑_{i=1}^{n} Y_{i})] = 𝔼[(n.μ_{1})^{2} (n.μ_{2})].

For the latter, let us observe that E[(∑_{i≠j}^{n} X_{i}^{2}X_{j})] = 𝔼[(n.χμ_{1}^{2})(n.χμ_{1})] and E[(∑_{i=1}^{n} X_{i}^{2}Y_{i})] = 𝔼[n.(χμ_{1}^{2} μ_{2})], so that the overall expectation of the product is given by

(2.14)  𝔼[n.(χ_{1}μ_{1}^{2}) n.(χ_{1}μ_{1}) n.(χ_{2}μ_{1}^{2}μ_{2}) n.(χ_{3}μ_{1}^{2}μ_{2})].

Both evaluations can be performed by using a suitable generalization of Theorem 2.6 involving multisets of umbral monomials.

Let us recall that a multiset of umbral monomials {μ_{1}, μ_{2}, …, μ_{k}} is

(2.15)  M = \{\underbrace{μ_{1}, …, μ_{1}}_{f(μ_{1})}, \underbrace{μ_{2}, …, μ_{2}}_{f(μ_{2})}, …, \underbrace{μ_{k}, …, μ_{k}}_{f(μ_{k})}\}.

The length of the multiset M in (2.15) is |M| = f(μ_{1}) + f(μ_{2}) + ⋯ + f(μ_{k}). The support of the multiset M is M̄ = {μ_{1}, …, μ_{k}}. A subdivision of a multiset M is a multiset S of l non-empty submultisets M_{i} = (M̄_{i}, h_{i}) of M with l ≤ |M| that satisfy [21]

a) ⋃_{i=1}^{l} M̄_{i} = M̄;
b) ∑_{i=1}^{l} h_{i}(μ) = f(μ) for every μ ∈ M̄.

Roughly speaking, a subdivision is a multiset which plays the same role for M as a partition for a set. A subdivision of the multiset M in (2.15) is denoted by

S = \{\underbrace{M_{1} …, M_{1}}_{g(M_{1})}, \underbrace{M_{2}, …, M_{2}}_{g(M_{2})}, …, \underbrace{M_{l}, …, M_{l}}_{g(M_{l})}\}.

Example 2.9. For the multiset M = {μ_{1}, μ_{1}, μ_{2}}, a subdivision is S = {{μ_{1}}, {μ_{1}, μ_{2}}} or S = {{μ_{1}}, {μ_{1}}, {μ_{2}}}. In the latter case, we

write S = \{\underbrace{\{μ_{1}\}}_{2}, \underbrace{\{μ_{2}\}}_{1}\}.

A MAPLE algorithm to find all multiset subdivisions is described in [23]. This procedure is faster than the iterated full partitions of Andrews and Stafford [1], given that it takes into account the multiplicity of all elements of M. The higher this multiplicity is, the more the procedure is efficient. The number of set partitions of [|M|] corresponding to the same subdivision S is

(2.16) d_{S} = \frac{1}{g(M_{1})! g(M_{2})! ⋯ g(M_{l})!} ∏_{i=1}^{k} \frac{f(μ_{i})!}{h_{1}(μ_{i})! h_{2}(μ_{i})! ⋯ h_{l}(μ_{i})!},

where h_{j}(μ_{i}) is the multiplicity of μ_{i} in the submultiset M_{j}.

Example 2.10. For M = {μ_{1}, μ_{1}, μ_{2}} let us consider the subdivision S = {{μ_{1}}, {μ_{1}, μ_{2}}}. If we label the elements of M = {μ_{1}, μ_{1}, μ_{2}} by the elements of {1, 2, 3}, the partitions of {1, 2, 3} corresponding to S = {{μ_{1}}, {μ_{1}, μ_{2}}} are 1|23 and 2|13, so d_{S} = 2. The same result arises from computing (2.16), since for M_{1} = {μ_{1}} and M_{2} = {μ_{1}, μ_{2}} we have g(M_{1}) = 1, g(M_{2}) = 1 and h_{1}(μ_{1}) = 1, h_{2}(μ_{1}) = h_{2}(μ_{2}) = 1. Moreover f(μ_{1}) = 2 and f(μ_{2}) = 1.

If we set (n.μ)_{M} = ∏_{i}(n.μ_{i})^{f(μ_{i})}, then, from Equivalence (2.12), we obtain

(2.17)  (n.μ)_{M} ≃ ∑_{S} d_{S} [n.(χμ)]_{S},

where the summation is over all subdivisions S of the multiset M, the integer d_{S} is given in (2.16), and

(2.18) [n.(χμ)]_{S} = ∏_{M_{i}∈S̄}[n.(χμ_{M_{i}})]^{g(M_{i})} with μ_{M} = ∏_{μ∈M̄} μ^{f(μ)}.

Equivalence (2.17) allows us to generalize Theorem 2.6 to umbral monomials with non-disjoint supports.

Theorem 2.11. [21] We have

𝔼 ([n.(χμ)]_{S}) = (n)_{|S|} ∏_{j=1}^{l} m_{h_{j}(μ_{1}),h_{j}(μ_{2}),...,h_{j}(μ_{k})}^{g(M_{j})},

where

m_{h_{j}(μ_{1}),h_{j}(μ_{2}),...,h_{j}(μ_{k})} = 𝔼 [μ_{1}^{h_{j}(μ_{1})}μ_{2}^{h_{j}(μ_{2})} ⋯ μ_{k}^{h_{j}(μ_{k})}] = 𝔼[μ_{M_{j}}].

Example 2.8 (continued). In order to compute 𝔼[(n.μ_{1})^{2}(n.μ_{2})], we use Equivalence (2.17). For M = {μ_{1}, μ_{1}, μ_{2}} we have

(n.μ_{1})^{2}(n.μ_{2}) ≃ [n.(χμ_{1}^{2}μ_{2})] + 2 [n.(χμ_{1})] [n.(χμ_{1}μ_{2})]

+ [n.(χμ_{1}^{2})][n.(χμ_{2})] + [n.(χμ_{1})]^{2}[n.(χμ_{2})],

and, from Theorem 2.11 with S = {{μ_{1}, μ_{1}}, {μ_{2}}}, we get

𝔼[(n.μ_{1})^{2}(n.μ_{2})] = n m_{2,1}+2(n)_{2}m_{1,0}m_{1,1}+(n)_{2}m_{2,0}m_{0,1}+(n)_{3}m_{1,0}^{2}m_{0,1},

with m_{i,j} = E[X^{i} Y^{j}] = 𝔼[μ_{1}^{i} μ_{2}^{j}]. In order to compute (2.14), Equivalence (2.17) and Theorem 2.11 have to be applied to

M = {χ_{1}μ_{1}^{2}, χ_{1}μ_{1}, χ_{2}μ_{1}^{2}μ_{2}, χ_{3}μ_{1}^{2}μ_{2}}.

Since M is a set, the subdivisions S are just set partitions. In difference to the strategy proposed by [105] to solve the same problem, the singleton umbrae in M allow us to speed up the computation. For example, consider the partition π_{1} = {{χ_{1}μ_{1}^{2}, χ_{1}μ_{1}}, {χ_{2}μ_{1}^{2}μ_{2}, χ_{3}μ_{1}^{2}μ_{2}}}. We have

[n.(χμ)]_{π_{1}} = [n.(χ χ_{1}^{2} μ_{1}^{3})] [n.(χ χ_{2}χ_{3} μ_{1}^{4}μ_{2}^{2})].

Since χ , χ_{2} , χ_{3} are uncorrelated umbrae, we have 𝔼[χ χ_{2}χ_{3} μ_{1}^{4}μ_{2}^{2}] = 𝔼[μ_{1}^{4}μ_{2}^{2}], but 𝔼[χ χ_{1}^{2} μ_{1}^{3}] = 0 since 𝔼[χ_{1}^{2}] = 0. Therefore we have

𝔼[n.(χμ)]_{π_{1}}] = 0.

This is the trick to speed up the computation: whenever, in one — or more — blocks of the subdivision S, there are at least two umbral monomials involving correlated singleton umbrae, the auxiliary umbra [n.(χμ)]_{S} has evaluation equal to zero. In the iteration procedure employed to build the subdivision S, this allows us to discard these blocks of the subdivision in the subsequent steps, see [22] for more details.

Further applications of Theorems 2.6 and 2.11 can be found in [22]. See Table 3 of Appendix 1 for a comparison of computation times with procedures existing in the literature.

2.2. Sheppard’s corrections. In the real world, continuous variables are observed and recorded in finite precision through a rounding or coarsening operation, i.e., a grouping rule. A compromise between the desire to know and the cost of knowing is then a necessary consequence. The literature on grouped data spans different research areas, see for example [46]. Due to the relevance of the method of moments as a tool to learn about a distribution, great attention has been paid in the literature to the computation of moments when data are grouped into classes. The moments computed by means of the resulting grouped frequency distribution are a first approximation to the estimated moments of the parent distribution, but they suffer from the error committed during the process of grouping. The correction for grouping is a sum of two terms, the first depending on the length of the grouping interval, the second being a periodic function of the position. For a continuous r.v., this very old problem was first discussed by Thiele [97], who studied the second term without taking into consideration the first, and then by Sheppard [86] who studied the first term, without taking into consideration the second. Compelling reasons have been given in the literature for neglecting the second term and thus for using Sheppard’s corrections, that are nowadays still employed. Quite recently the appropriateness of Sheppard’s corrections was re-examined in connection with some applications, see, for example, [15] and [102].

Let us consider a r.v. X with probability density function f, representing the continuous parent distribution underlying data. Denote by a_{i} its i-th moment (also called raw moment to distinguish it from the estimated one),

(2.19)  a_{i} = ∫_{−∞}^{∞} x^{i}f(x) dx.

Let ã_{i} be the i-th moment of the grouped distribution,

(2.20)  ã_{i} = \frac{1}{h} ∫_{−∞}^{∞} x^{i} ∫_{− \frac{1}{2} h}^{\frac{1}{2} h} f(x + z) dz dx,

with h the length of the grouping intervals. Then the i-th raw moment a_{i} can be reconstructed via Sheppard’s corrections as

(2.21)  a_{i} = ∑_{j=0}^{i} \binom{i}{j} (2^{1−j} − 1) B_{j} h^{j} ã_{i−j},

where B_{j} is the j-th Bernoulli number.

The derivation of Sheppard’s corrections was a popular topic in the first half of the last century, see [44] for a historical account. This is because Sheppard deduced Equation (2.21) by using the Euler–Maclaurin summation formula and by assuming that the density function had high order contact with the x-axis at both ends. So there was a considerable controversy on the set of sufficient conditions to be required in order to use formula (2.21). If the rounding lattice is assumed to be random, these sufficient conditions can be removed.

Grouping includes also censoring or splitting data into categories during collection or publication, and so it does not only involve continuous variables. The derivation of corrections for raw moments of a discrete parent distribution followed a different path. They were first given in the Editorial of Vol. 1, no. 1, of Annals of Mathematical Statistics (page 111). The method used to develop a general formula was extremely laborious. Some years later, Craig [14] considerably reduced and simplified the derivation of these corrections without requiring any other condition and stating these corrections on the average. Except for the papers of Craig [14] and Baten [7], no attention has been paid to multivariate generalizations of Sheppard’s corrections, probably due to the complexity of the resulting formulas.

The symbolic method has allowed us to find such corrections both for continuous and discrete parent distributions in closed-form formulas, which have been naturally extended to multivariate parent distributions [17]. The result is an algorithm to implement all these formulas in any symbolic package [20]. A similar employment of the symbolic method can be performed within wavelet theory. Indeed, the reconstruction of the full moment information a_{i} in (2.21) through the grouped moments ã_{i} can be seen as some kind of multilevel analysis, see [79, 85].

If the sequence {ã_{i}}_{i≥0} in (2.20) is umbrally represented by the umbra α̃ and the sequence {a_{i}}_{i≥0} in (2.19) is umbrally represented by the umbra α, Sheppard’s corrections (2.21) can be obtained by expanding the i-th moment of

(2.22)  α ≡ α̃ + h (ι + \frac{1}{2}) ,

where ι is the Bernoulli umbra. Equivalence (2.22) still holds, if the moments refer to a parent distribution over a finite interval. When a_{i} is the i-th moment of a discrete population, and ã_{i} is the i-th moment calculated from the grouped frequencies, Equivalence (2.22) needs to be modified to

(2.23)  α ≡ α̃ + h (ι + \frac{1}{2}) + \frac{h}{m} (−1.ι − \frac{1}{2}) ,

where we assume that m consecutive values of the discrete population are grouped in a frequency class of width h. In Equivalence (2.23), −1.ι denotes one more auxiliary umbra, the inverse umbra of the Bernoulli umbra.^{8} Sheppard’s corrections (2.22) can be recovered from (2.23) letting m → ∞ in a suitable way.

As done in Section 2.1, the generalization of Equivalences (2.22) and (2.23) to the multivariate setting follows by introducing suitable umbral monomials.

Let (μ_{1}, μ_{2}, …, μ_{k}) be a k-tuple of umbral monomials with nondisjoint supports.

Definition 2.12. The sequence {m_{t_{1}, t_{2},..., t_{k}}} is umbrally represented by the k-tuple (μ_{1}, μ_{2}, …, μ_{k}) if E[μ_{1}^{t_{1}}μ_{2}^{t_{2}} ⋯ μ_{k}^{t_{k}}] = m_{t_{1}, t_{2},..., t_{k}} for all positive integers t_{1}, t_{2}, …, t_{k}.

8Given two umbrae α and γ, they are said to be inverse to each other if α+γ ≡ ε.

We denote the inverse of the umbra α by −1.α. An umbra and its inverse are uncorrelated. The inverse of an umbra is not unique, but any two umbrae inverse to any given umbra are similar.

Following the analogy with random vectors, m_{t_{1}, t_{2},..., t_{k}} is called the multivariate moment of the k-tuple (μ_{1}, μ_{2}, …, μ_{k}) of order t = (t_{1}, t_{2}, …, t_{k}) .

According to the latter notation in (2.18), if M is a multiset with support M̄ = {μ_{1}, μ_{2}, …, μ_{k}} and multiplicities f(μ_{j}) = t_{j} for j = 1, 2, …, k, then we may write E[μ_{M}] = m_{t_{1}, t_{2},..., t_{k}}. Moreover, assume there exists a multivariate random vector X = (X_{1}, X_{2}, …, X_{k}) with joint density f_{X}(x) over ℝ^{k} such that

(2.24)  m_{t_{1}, t_{2},..., t_{k}} = ∫_{ℝ^{k}} x^{t}f_{X}(x) dx

is its multivariate moment of order t. Note that, as in the univariate case, without loss of generality we may consider joint densities with range^{9} of bounded rectangle type. The moments calculated from the grouped frequencies are

(2.25) m̃_{t_{1}, t_{2},..., t_{k}} = \frac{1}{h_{1}h_{2} ⋯ h_{k}} ∫_{A_{k}} ∫_{ℝ^{k}} ∏_{j=1}^{k}(x_{j} + z_{j})^{t_{j}}f_{X}(x) dx dz,

where

A_{k} = {z = (z_{1}, z_{2}, …, z_{k}) ∈ ℝ^{k} :

z_{j} ∈ (−\frac{1}{2}h_{j}, \frac{1}{2}h_{j}) for all j ∈ \{1, 2, …, k\}\} ,

and h_{j} ∈ ℝ − {0} are the window widths for any component. Then we have the following theorem.

Theorem 2.13. If the sequence {m_{t_{1}, t_{2},..., t_{k}}} in (2.24) is umbrally represented by the k-tuple (μ_{1}, μ_{2}, …, μ_{k}) and the sequence {m̃_{t_{1}, t_{2},..., t_{k}}} in (2.25) is umbrally represented by the k-tuple (μ̃_{1}, μ̃_{2}, …, μ̃_{k}), then

(2.26)  μ_{M} ≡ [μ̃ + h (ι + \frac{1}{2})]_{M} ,

where M is the multiset given in (2.15),

[μ̃ + h (ι + \frac{1}{2})]_{M} = ∏_{j=1}^{k} [μ̃_{j} + h_{j} (ι_{j} + \frac{1}{2})]^{t_{j}} ,

and {ι_{j}} are uncorrelated Bernoulli umbrae.

Equivalence (2.26) is implemented in [20] by means of the following steps. In order to recover expressions of raw multivariate moments in terms of grouped moments, we need to multiply summations like

9The range of a joint (probability) density (function) f_{X}(x) is the subset of ℝ^{k} where f_{X}(x) is different from 0.

∑_{s_{j}=1}^{i} \binom{i}{s_{j}}μ̃_{j}^{s_{j}}h_{j}^{i−s_{j}} (2^{1−i+s_{j}} − 1) B_{i−s_{j}},

corresponding to the i-th power of μ̃_{j} + h_{j} (ι_{j} + \frac{1}{2}) , and then to replace occurrences of products like μ̃_{1}^{s_{1}}μ̃_{2}^{s_{2}} ⋯ μ̃_{k}^{s_{k}} by m̃_{s_{1}, s_{2},..., s_{k}}. If the multivariate parent distribution is discrete, Equivalence (2.26) has to be updated to

μ_{M} ≡ [μ̃ + h (ι + \frac{1}{2}) + \frac{h}{m} (−1.ι − \frac{1}{2})]_{M}

where the symbol on the right-hand side denotes the product

∏_{j=1}^{k} [μ̃_{j} + h_{j} (ι_{j} + \frac{1}{2}) + \frac{h_{j}}{m_{j}} (−1.ι_{j} − \frac{1}{2})]^{t_{j}} ,

and m_{j} denotes the number of consecutive values grouped in a frequency class of width h_{j}.

3. Compound Poisson processes

In 1973, Rota [76] claimed that compound Poisson processes are related to polynomial sequences {p_{i}(x)}_{i≥0} of binomial type. Such a sequence of polynomials is defined by p_{i}(x) being of degree i for all positive integers i, p_{0}(x) = 1, and

(3.1)  p_{i}(x + y) = ∑_{j=0}^{i} \binom{i}{j}p_{j}(x)p_{i−j}(y),  i = 0, 1, ….

This connection has ramifications into several other areas of analysis and probability: Lagrange expansions, renewal theory, exponential family distributions and infinitely divisible processes. Due to its relevance, different authors have tried to formalize this connection. Stam [91] studied polynomial sequences of binomial type in terms of an integer-valued compound Poisson process. He was especially interested in the asymptotic behavior of the probability generating function p_{i}(x)/p_{i}(1) for i → ∞. Partial results are available on the radius of convergence under suitable conditions, involving the power series with coefficients {p_{i}(x)}_{i≥0}. The resulting theory relies on a complicated system of notations. Constantine and Savit [12] have derived a generalization of Dobińsky’s formula by means of compound Poisson processes. Pitman [66] has investigated some probabilistic aspects of Dobińsky’s formula.

The theory of Bell umbrae gives a natural way to relate compound Poisson processes to polynomial sequences of binomial type. Historically, this connection was the first application in probability of classical umbral calculus as proposed in 1994. More details are given in [33].

So let us come back to the Bell umbra β, introduced in Section 2. Instead of using Bell numbers, we may characterize this umbra by using its factorial moments. The factorial moments of an umbra are

(3.2)  a_{(i)} = 𝔼[(α)_{i}] = \begin{cases}1, & if i = 0, \\ E[α(α − 1) ⋯ (α − i + 1)], & if i > 0. \end{cases}

Theorem 3.1. [33] The Bell umbra β is the unique umbra (up to similarity) with all factorial moments equal to 1.

In Equation (2.4), set 𝔼[(m.α)^{i}] = q_{i}(m) and observe that q_{i}(m) is a polynomial of degree i in m. Suppose that we replace m by an indeterminate^{10} x. The symbol having sequence {q_{i}(x)}_{i≥0} as moments is denoted by x.α and called the dot-product of x and α, that is,

(3.3)  q_{i}(x) = 𝔼[(x.α)^{i}] = ∑_{λ⊢i}(x)_{ν_{λ}}d_{λ}a_{λ}

for all positive integers i. A Poisson r.v. Po(x) is represented by x.β, which is called the Bell polynomial umbra.

By using the exponential (partial) Bell polynomials B_{i,k} (cf. [70]), a different way to write q_{i}(x) in (3.3) is

(3.4)  q_{i}(x) = 𝔼 [(x.α)^{i}] = ∑_{k=1}^{i}(x)_{k}B_{i,k}(a_{1}, a_{2}, …, a_{i−k+1}).

Equation (3.4) with α replaced by the Bell umbra β returns the socalled exponential polynomials Φ_{i}(x),

(3.5)  𝔼 [(x.β)^{i}] = Φ_{i}(x) = ∑_{k=1}^{i}(x)_{k}S(i, k),

where S(i, k) are Stirling numbers of the second kind. According to [9], these polynomials were introduced by S. Ramanujan in his unpublished notebooks. Later, exponential polynomials were studied by Bell [8] and Touchard [100]. Rota, Kahaner and Odlyzko [76] have stated their basic properties via umbral operators. If, in (3.5), we replace x by a non-negative integer n, the resulting umbra n.β is the sum of n similar and uncorrelated Bell umbrae, likewise in probability theory where a Poisson r.v. Po(n) is the sum of n i.i.d. Poisson r.v.’s Po(1). More generally, the closure under convolution of Poisson distributions F_{t}, that is, F_{s} ⋆ F_{t} = F_{s+t}, is simply reformulated as^{11}

(t + s).β ≡ t.β + s.β

10We replace ℝ by ℝ[x, y]. The indeterminate y is added because of (3.1).

11According to property iv) of Corollary 1 in [33], the umbra (t + s).β is similar to t.β + s.β′, with β and β′ uncorrelated Bell umbrae. If s ≠ t then s.β and t.β are two distinct symbols, and so uncorrelated. Then we can use s.β instead of s.β′.

with t, s ∈ ℝ. If, in (3.5), we replace x by a generic umbra α, the auxiliary umbra α.β satisfies

(α.β)^{i} ≃ Φ_{i}(α) = ∑_{k=1}^{i}(α)_{k}S(i, k).

The umbra α.β represents a random sum of independent Poisson r.v.’s Po(1) indexed by an integer r.v. Y represented by α, that is, a randomized Poisson r.v. Po(Y) with parameter Y.

If we swap the umbrae α and β, the resulting umbra is the α-partition umbra β.α representing a compound Poisson r.v. with parameter 1. Let us recall that a compound Poisson r.v. with parameter 1 is a random sum S_{N} = X_{1}+X_{2}+⋯+X_{N}, where N ∼ Po(1) and {X_{j}} are i.i.d.r.v.’s. In the α-partition umbra β.α, β represents N and α represents one of {X_{j}}. Moreover, since the Poisson r.v. Po(x) is umbrally represented by the Bell polynomial umbra x.β, a compound Poisson r.v. with parameter x is represented by the polynomial α-partition umbra x.β.α. This encodes the connection between compound Poisson processes and polynomial sequences of binomial type, which is what we have claimed at the beginning of this section. The name “partition umbra” has a probabilistic background. Indeed, the parameter of a Poisson r.v. is usually denoted by x = λt, with t representing time. When the (time) interval in which t ranges is partitioned into non-overlapping ones, their contributions are stochastically independent and add up to S_{N}. This circumstance is umbrally expressed by the relation

(3.6)  (x + y).β.α ≡ x.β.α + y.β.α,

paralleling the property (3.1) for the binomial sequences p_{i}(x) = 𝔼[(x.β.α)^{i}]. If {a_{i}}_{i≥0} is umbrally represented by the umbra α, moments of x.β.α are

(3.7) 𝔼[(x.β.α)^{i}] = ∑_{k=1}^{i} x^{k}B_{i,k}(a_{1}, a_{2}, …, a_{i−k+1}) = ∑_{λ⊢i} d_{λ} x^{ν_{λ}}a_{λ},

where B_{i,k} are the partial Bell exponential polynomials. For x = 1, we infer from (3.7) that

(3.8) 𝔼[(β.α)^{i}] = ∑_{k=1}^{i} B_{i,k}(a_{1}, a_{2}, …, a_{i−k+1}) = 𝒴_{i}(a_{1}, a_{2}, …, a_{i}),

where 𝒴_{i} = 𝒴_{i}(a_{1}, a_{2}, …, a_{i}) is the i-th complete Bell exponential polynomial [70]. Then the α-partition umbra represents the binomial sequence {𝒴_{i}}.

The following example proves that partition umbrae are employed in representing r.v.’s with distribution different from the Poisson distribution.

Example 3.2. A Gaussian r.v. N (m, s) is represented by the umbra m + β.(sη), where η is an umbra satisfying E[η^{i}] = δ_{i,2}, for all positive integers i. Indeed, the i-th moment of m + β.(sη) is

𝔼\{[m + β.(sη)]^{i}\} = ∑_{k=0}^{⌊i/2⌋} (\frac{s^{2}}{2} )^{k} \frac{(i)_{2k}}{k!} m^{i−2k},

which gives the i-th element of a well-known sequence of orthogonal polynomials, the Hermite polynomials H_{i}^{(ν)}(x) [70], with m = x and s^{2} = −ν.

Theorem 3.3. [33] The α-partition umbra satisfies the equivalence

(3.9)  (β.α)^{i} ≃ α (β.α + α)^{i−1}, i = 1, 2, …,

and conversely every umbra satisfying Equivalence (3.9) is an α-part-ition umbra.

One more property of the α-partition umbra involves disjoint sums. The disjoint sum α ˙+ γ of α and γ is an auxiliary umbra representing the sequence {a_{i} +g_{i}}_{i≥0} with {a_{i}}_{i≥0} and {g_{i}}_{i≥0} umbrally represented by the umbra α and γ, respectively.

Theorem 3.4. [34] β.(α ˙+ γ) ≡ β.α + β.γ.

If in the α-partition umbra β.α, we replace the umbra β by a different umbra γ, we get a new auxiliary umbra, the dot-product γ.α. More formally, let us consider the polynomial q_{i}(x) in (3.3), and suppose to replace x by an umbra γ with moments {g_{i}}_{i≥0}. The polynomial

(3.10)  q_{i}(γ) = ∑_{λ⊢i}(γ)_{ν_{λ}}d_{λ}a_{λ}

is an umbral polynomial with support {γ}. The auxiliary umbra γ.α is an auxiliary umbra with moments

(3.11)  𝔼[(γ.α)^{i}] = 𝔼[q_{i}(γ)],

for all positive integers i. A very special dot-product umbra is χ.α, with χ the singleton umbra. This umbra is the keystone in dealing with sequences of cumulants. More details will be given in Section 4. Indeed we will show that any umbra is the partition umbra of χ.α, that is, the symbolic method of moments is like a calculus of measures on Poisson algebras. In the following, we give an example of how to use the umbra χ.α to disclose a not well-known connection between partition polynomials and Kailath–Segall polynomials.

Example 3.5 (The Kailath–Segall formula). Let {X_{t}}_{t≥0} be a centered Lévy process^{12} with moments of all orders. The variations of the process {X_{t}^{(i)}}_{t≥0} are defined as X_{t}^{(1)} = X_{t} and

12A centered Lévy process is a zero mean stochastic process with independent and stationary increments [80], see also Section 3.1.

X_{t}^{(2)} = [X, X]_{t} = σ^{2}t + ∑_{0<s≤t}(ΔX_{s})^{2}, X_{t}^{(i)} = ∑_{0<s≤t}(ΔX_{s})^{i}, for i ≥ 3,

where ΔX_{s} = X_{s} − X_{s-}. The iterated stochastic integrals

P_{t}^{(0)} = 1, P_{t}^{(1)} = X_{t}, P_{t}^{(i)} = ∫_{0}^{t} P_{s-}^{(i−1)}dX_{s}, for i ≥ 2,

are related to the variations {X_{t}^{(i)}}_{t≥0} by the Kailath–Segall for-

mula [47]

(3.12) P_{t}^{(i)} = \frac{1}{i} (P_{t}^{(i−1)}X_{t}^{(1)} − P_{t}^{(i−2)}X_{t}^{(2)} + ⋯ + (−1)^{i+1}P_{t}^{(0)}X_{t}^{(i)}) .

It follows by induction that P_{t}^{(i)} is a polynomial in X_{t}^{(1)}, X_{t}^{(2)}, …, X_{t}^{(i)}, called the i-th Kailath–Segall polynomial. Theorem 3.6. [30] If {Υ_{t}}_{t≥0} and {ψ_{t}}_{t≥0} are two families of umbrae with 𝔼[Υ_{t}^{i}] = i!E [P_{t}^{(i)}] and 𝔼[ψ_{t}^{i}] = E [X_{t}^{(i)}] , for all positive integers i, then

(3.13)  Υ_{t} ≡ β.[(χ.χ)ψ_{t}] and (χ.χ)ψ_{t} ≡ χ.Υ_{t}.

Equivalences (3.13) give the umbral version of the Kailath–Segall formula (the former) and its inversion (the latter). The proof essentially relies on Theorem 3.3, since Equation (3.12) is equivalent to 𝔼[Υ_{t}^{i}] = 𝔼 [ψ_{t}(Υ_{t} + ψ_{t})^{i−1}] . By the former Equivalence (3.13), the Kailath–Segall polynomials result in the complete Bell exponential polynomials (3.8) in the variables

{(−1)^{i−1}(i − 1)!E [X_{t}^{(i)}]} .

The latter similarity in Equivalence (3.13) is a generalization of

χ.(χ_{1}x_{1} + ⋯ + χ_{n}x_{n}) ≡ (χ.χ)σ,

which gives the elementary symmetric polynomials in terms of power sum symmetric polynomials, umbrally represented by the umbra σ, see [21]. That is, if we replace the jumps {ΔX_{s}} in X_{t}^{(i)} by suitable indeterminates {x_{s}}, then the Kailath–Segall polynomials reduce to the polynomials given in [4]. These polynomials deserve further investigations in connection with stochastic integrals [62, 78], for which we believe that an approach via the symbolic method of moments is fruitful.

Theorem 3.6 gives the symbolic representation of Lévy process variations. The next subsection is devoted to the symbolic representation of Lévy processes. This representation extends the claimed connection between binomial sequences and compound Poisson processes to the more general class of infinite divisible distributions.

3.1. Time-space harmonic polynomials. Lévy processes are stochastic processes playing a fundamental role in various fields such as physics, engineering and actuarial science [6]. They are particularly popular in mathematical finance as models for jumps or spikes of asset prices in financial markets.

The martingale property is one of the fundamental mathematical properties which underlies many finance models, when no knowledge of the past price process can help to predict the future asset price. In particular, a market is said to be efficient if the price process is modeled by a martingale (for more details see [82]).

Lévy processes are not martingales if they have an increasing or decreasing drift. One tool to compensate this fault is to use Wald’s martingale, which has various applications [53]. If {X_{t}}_{t≥0} = {X_{t}} is a Lévy process,^{13} then the Wald exponential martingale is defined as

(3.14)  M_{t} = \frac{exp\{zX_{t}\}}{E[exp\{zX_{t}\}]}.

Definition (3.14) requires the existence of the moment generating function E[exp{zX_{t}}] in a suitable neighborhood of the origin. A way to overcome this gap is to look for polynomials P (x, t), with the property that, when the indeterminate x is replaced by a Lévy process {X_{t}}, the resulting stochastic process is a martingale. These polynomials have been introduced first by Neveu [56] for random walks and then by Sengupta [43] for Lévy processes. They are called time-space harmonic polynomials. Definition 3.7. A family of polynomials {P (x, t)}_{t≥0} is called timespace harmonic with respect to a stochastic process {X_{t}}_{t≥0} if

E[P (X_{t}, t) | 𝔉_{s}] = P (X_{s}, s),

for 0 ≤ s ≤ t, where 𝔉_{s} = σ (X_{τ} : 0 ≤ τ ≤ s) is the natural filtration^{14} associated with {X_{t}}_{t≥0}.

Wald’s martingale (3.14) has been recently used in order to characterize time-space harmonic polynomials [84], but without reaching a closed expression.

The aim of this subsection is to show how to recover a very general class of time-space harmonic polynomials by using a symbolic representation of Lévy processes. In particular, our attention focuses on special properties of Lévy processes, namely [80]

a) X_{0} = 0 a.s. (almost surely);

b) for all n ≥ 1 and for all 0 ≤ t_{1} ≤ t_{2} ≤ … ≤ t_{n} < ∞, the r.v.’s X_{t_{2}} − X_{t_{1}}, X_{t_{3}} − X_{t_{2}}, …, X_{t_{n}} − X_{t_{n−1}} are independent; c) for 0 ≤ s ≤ t, X_{t+s} − X_{s} and X_{t} are identically distributed

13When no misunderstanding is possible, we simply use {X_{t}} to denote a stochastic process.

14A natural filtration 𝔉_{t} is the σ-algebra generated by the pre-images X_{s}^{−1}(B)

for Borel subsets B of ℝ and times s with 0 ≤ s ≤ t.

(stationarity). Due to these properties, at time t the information of the process {X_{t}} may be split into i.i.d. subprocesses:

(3.15)  X_{t} =^{d} \underbrace{ΔX_{t/n} + ΔX_{t/n} + ⋯ + ΔX_{t/n}}_{n} .

Property (3.15) characterizes the so-called infinitely divisible r.v.’s, which are in one-to-one correspondence with Lévy processes. This is also why Lévy himself refers to these processes as a sub-class of additive processes. Starting from (2.4), a symbolic version of a Lévy process can be obtained by means of the same steps as employed in the construction of q_{i}(x) at the beginning of Section 3. If the positive integer m is replaced by t ∈ ℝ in (2.4), then the auxiliary umbra t.α satisfies

(3.16)  𝔼[(t.α)^{i}] = q_{i}(t) = ∑_{λ⊢i}(t)_{ν_{λ}} d_{λ} a_{λ},

for all positive integers i. The umbra t.α is the dot-product of t and α. Since m.α is a summation of m uncorrelated and similar umbrae, t.α parallels the decomposition (3.15). Its properties are the same as the ones of Lévy processes when they admit finite moments. Here there are some examples.

Binomial property. The i-th moment q_{i}(t) = E(X_{t}^{i}) of a Lévy process is a polynomial function of t which satisfies the property (3.1) of binomial sequences. In terms of umbrae, this property is equivalent to

(3.17)  (t + s).α ≡ t.α + s.α.

We cannot miss the analogy between Equivalences (3.6) and (3.17). Indeed, an umbra representing a Lévy process is similar to an umbra representing a compound Poisson process. This is because for any umbra α there exists an umbra γ satisfying α ≡ β.γ, so that t.α ≡ t.β.γ. The proof of this statement is given in Section 4.

Additivity property. If {W_{t}} and {Z_{t}} are two independent Lévy processes, then the process {X_{t}} with X_{t} = W_{t} + Z_{t} is a Lévy process. In terms of umbrae, the additivity property corresponds to

t.(α + γ) ≡ t.α + t.γ, with α, γ ∈ 𝒜 .

Homogeneity property. If {X_{t}} is a Lévy process and c ∈ ℝ, then {c X_{t}} is a Lévy process. This homogeneity property is equivalent to

t.(cα) ≡ c(t.α), with α ∈ 𝒜 .

Nested processes. If {X_{t}} is a Lévy process, then {(X_{t})_{s}} is a Lévy process, that is, a nested Lévy process is again a Lévy process. The analog for t.α is t.(s.α) ≡ s.(t.α) ≡ (st).α, with t, s ∈ ℝ.

According to Definition 3.7, in order to verify that an umbral polynomial is a time-space harmonic polynomial, a suitable notion of conditional evaluation has been introduced in the polynomial ring ℝ[x][𝒜], see [30]. Denote by X the set X = {α}.

Definition 3.8. The linear operator 𝔼( · α) : ℝ[x][𝒜] ⟶ ℝ[X] satisfying

i) 𝔼(1 α) = 1; ii) 𝔼(x^{m}α^{n}γ^{i}ς^{j} ⋯ α) = x^{m}α^{n}𝔼(γ^{i})𝔼(ς^{j}) ⋯ for uncorrelated um-

brae α, γ, ς, … and for positive integers m, n, i, j, …

is called conditional evaluation with respect to α.

In other words, the conditional evaluation with respect to α deals with the umbra α as an indeterminate. As it happens for r.v.’s, the conditional evaluation is an element of ℝ[x][𝒜], and the overall evaluation of 𝔼(p α) gives 𝔼(p), with p ∈ ℝ[x][𝒜].

The conditional evaluation needs to be carefully handled when dot products are involved. Let us observe that the conditional evaluation with respect to n.α satisfies

𝔼[(n + 1).α | n.α] = 𝔼(n.α + α′ | n.α) = n.α + 𝔼(α′),

with α′ an umbra similar to α. By similar arguments, for all positive integers n and m, since 𝔼{[(n + m).α]^{i} | n.α} = 𝔼{[n.α + m.α′]^{i} | n.α}, we have

(3.18)  𝔼\{[(n + m).α]^{i} | n.α\} = ∑_{j=0}^{i} \binom{i}{j}(n.α)^{j}𝔼[(m.α′)^{i−j}].

Equation (3.18) suggests how to define the conditional evaluation with respect to the auxiliary umbra s.α.

Definition 3.9. The conditional evaluation of t.α with respect to the auxiliary umbra s.α is

𝔼[(t.α)^{i} | s.α] = ∑_{j=0}^{i} \binom{i}{j}(s.α)^{j}𝔼\{[(t − s).α′]^{i−j}\}.

Theorem 3.10. [17] For all non-negative integers i, the family of polynomials

(3.19)  Q_{i}(x, t) = 𝔼[(−t.α + x)^{i}] ∈ ℝ[x]

is time-space harmonic with respect to {t.α}_{t≥0}, that is,

(3.20)  𝔼 [Q_{i}(t.α, t) s.α] = Q_{i}(s.α, s)  for s, t ≥ 0.

The umbra −t.α in (3.19) is a special auxiliary umbra with the remarkable property

(3.21)  −t.α + t.α′ ≡ ε,

where α′ ≡ α. Due to Equivalence (3.21), the umbra −t.α is the inverse^{15} umbra of t.α. Its moments can be obtained via Equation (3.16) by replacing t by −t. The family {Q_{i}(x, t)}_{i≥0} is a basis for the space of time-space harmonic polynomials [29] so that many families of known time-space harmonic polynomials can be recovered as a linear combination of {Q_{i}(x, t)}_{i≥0}. For example, Hermite polynomials are time-space harmonic with respect to Wiener processes, and Poisson–Charlier polynomials are time-space harmonic with respect to Poisson processes. More examples are given in Appendix 2.

Remark 3.11. Theorem 3.10 may be proved when the Lévy process admits moments up to a finite order m. In this case, the polynomials Q_{i}(x, t) in Equation (3.19) are defined up to i ≤ m. Then formal power series are replaced by polynomials of degree m, and operations like summation and multiplication among formal power series can be performed in terms of polynomials [5].

Example 3.12 (Random walks). When the parameter t is replaced by a positive integer n, Theorem 3.10 still holds. The sequence {n.α}_{n≥0} represents a random walk {X_{n}}_{n≥0} with X_{0} = 0, X_{n} = M_{1} + M_{2} + ⋯ + M_{n}, and {M_{n}}_{n≥0} a sequence of i.i.d.r.v.’s. Note that the sequence {M_{n}}_{n≥0} satisfies M_{0} = X_{0} = 0 and M_{n} = X_{n} − X_{n−1}, for all positive integers n.

To recover the connection between Wald’s martingale (3.14) and the polynomials Q_{i}(x, t), we have to introduce the generating function of an umbra α, that is, the formal power series [34]

f(α, z) = 1 + ∑_{i≥1} a_{i} \frac{z^{i}}{i!} ∈ ℝ[x][[z]],

whose coefficients are the moments of the umbra α. Formal power series allow us to work with generating functions [106] which do not have a positive radius of convergence or have indeterminate coefficients [92]. Similar umbrae α ≡ γ have the same generating function f(α, z) = f(γ, z). This property allows us to recover the generating functions of the auxiliary umbrae. Table 2 gives some examples.

In particular, the generating function of the time-space harmonic polynomial umbra −t.α + x is

(3.22)  f(−t.α + x, z) = \frac{exp\{xz\}}{f(α, z)^{t}} = 1 + ∑_{i≥1} Q_{i}(x, t)\frac{z^{i}}{i!} .

15Since −t.α and t.α are two distinct symbols, they can be considered uncorrelated. When no confusion arises, we will use the symbol t.α instead of t.α′.

Table 2. Special generating functions

Umbrae  Generating functions
Unity umbra u  f(u, z) = 1
Singleton umbra χ f(χ, z) = 1 + z
Bell umbra β  f(β, z) = exp[e^{z} − 1]
Summation α + γ f(α + γ, z) = f(α, z)f(γ, z)
Dot product n.α f(n.α, z) = f(α, z)^{n}
Dot product t.α f(t.α, z) = f(α, z)^{t}
Inverse −t.α  f(−t.α, z) = 1/f(α, z)^{t}

By replacing x by t.α in (3.22), we can compare Equation (3.22) with the following series expansion for Wald’s exponential martingale (3.14)

\frac{exp\{zX_{t}\}}{E[exp\{zX_{t}\}]} = 1 + ∑_{i≥1} R_{i}(X_{t}, t)\frac{z^{i}}{i!} .

The algebraic structure of formal power series is isomorphic to sequences endowed with the convolution product, each series corresponding to the sequence of its coefficients [92]. Thus, equality of two formal power series is interpreted as equality of their corresponding coefficients, that is, E[R_{i}(X_{t}, t)] = 𝔼[Q_{i}(t.α, t)].

Since 𝔼[Q_{i}(t.α, t)] = 0, we have E[R_{i}(X_{t}, t)] = 0 for i ≥ 1 and

1 + ∑_{i≥1} E[R_{i}(X_{t}, t)]\frac{z^{i}}{i!} = 1

in agreement with Wald’s identity [56]. Therefore, the sequence {E[R_{i}(X_{t}, t)]}_{t≥0} is umbrally represented by the augmentation umbra ε, with f(ε, z) = 1. But this is exactly what happens if in −t.α + x we replace x by t.α, due to (3.21).

Time-space harmonic polynomials with respect to {t.α}_{t≥0} are characterized by the following theorem.

Theorem 3.13. [17] A polynomial P (x, t) = ∑_{j=0}^{k} p_{j}(t) x^{j}, of degree k for all t ≥ 0, is time-space harmonic with respect to a Lévy process {X_{t}}, represented by t.α, if and only if

p_{j}(t) = ∑_{i=j}^{k} \binom{i}{j} p_{i}(0) 𝔼[(−t.α)^{i−j}], for j = 0, …, k.

An open problem is the problem of extending the symbolic method of moments to matrix-valued stochastic processes by using Lévy processes and time-space harmonic polynomials. Matrix-value time-space harmonic polynomials include the matrix-valued counterparts of classical polynomials, as the Hermite, Laguerre or Jacobi polynomials, which have gained increasing interest in multivariate statistics. Different methods have been proposed in the literature, see [51] and references therein, by using eigenvalues of random matrices or hypergeometric functions of matrix argument. It would be interesting to check if dealing with hypergeometric functions by means of the symbolic method leads to more feasible expressions as happens for the multivariate generating functions, see Section 7.

4. Cumulants

A probability measure ℙ is infinitely divisible if for any positive integer n there is a probability measure ℙ_{n} such that ℙ = (ℙ_{n})^{n∗}, where (ℙ_{n})^{n∗} denotes the n-fold convolution of ℙ_{n}. This property parallels Equation (3.15) for Lévy processes {X_{t}} and can be generalized. Indeed, if we denote by ℙ_{t} the probability measure of X_{t}, then ℙ_{t} = (ℙ_{1})^{t∗} with ℙ_{1} the probability measure of X_{1}. Since the convolution of probability measures corresponds to a summation of r.v.’s, the r.v. X_{t} can be expressed as a summation of t times the r.v. X_{1}. The symbolic method of moments allows us to generalize this decomposition by replacing t ∈ ℝ by the umbra γ. The resulting auxiliary umbra is the dot-product of α and γ in (3.10), that represents a summation of γ times the umbra α. A special dot-product is the umbra χ.α, where χ is the singleton umbra, since its moments are cumulants of the sequence {a_{i}}_{i≥0} umbrally represented by α.

Among the number sequences related to r.v.’s, cumulants play a central role characterizing many r.v.’s occurring in classical stochastic processes. For example, a Poisson r.v. is the unique probability distribution for which all cumulants are equal to the parameter. A Gaussian r.v. is the unique probability distribution for which all cumulants vanish beyond the second.

(Formal) cumulants {c_{j}} of a sequence {a_{i}} are usually defined by the identity of formal power series [73]

(4.1)  1 + ∑_{i≥1} a_{i} \frac{t^{i}}{i!} = exp (∑_{j≥1} c_{j} \frac{t^{j}}{j!}),

which uniquely determines the non-linear functionals {c_{j}} as polynomials in {a_{i}}. In this definition, any questions concerning convergence of involved series may be disregarded [92]. Moreover, many difficulties connected to the so-called problem of cumulants smooth out. Here, with the problem of cumulants, we refer to the characterization of sequences that are cumulants of some probability distribution. The simplest example is that the second cumulant of a probability distribution must always be non-negative, and is zero if and only if all of the higher cumulants are zero. Cumulants are not subject to such constraints when they are analyzed from a symbolic point of view.

Example 4.1 (The variance). The variance Var(X) = E(X^{2}) − E(X)^{2} of a r.v. is the cumulant of second order, an index very useful in statistics. It possesses remarkable properties that may be extended to cumulants of any order:

i) it is invariant under translation: Var(X + a) = Var(X) for any constant a;

ii) Var(X + Y) = Var(X) + Var(Y), if X and Y are independent

r.v.’s; iii) Var(X) is a polynomial in the moments of the r.v. X.

A first approach to the theory of cumulants via the classical umbral calculus was given by Rota and Shen in [74]. By using the symbolic method of moments, the formulas commonly used to express moments of a r.v. in terms of cumulants, and vice-versa, are encoded in a dotproduct involving the singleton umbra [34].

Definition 4.2. The α-cumulant umbra κ_{α} is the auxiliary umbra satisfying κ_{α} ≡ χ.α.

The three main algebraic properties of cumulants can be stated as follows:

i) additivity property:

χ.(α + γ) ≡ χ.α ˙+ χ.γ,

that is, if {c_{i}(X + Y)} is the sequence of cumulants of the summation of two independent r.v.’s X and Y, then c_{i}(X+Y) = c_{i}(X) + c_{i}(Y); ii) homogeneity property:

χ.(aα) ≡ a(χ.α), for all a ∈ ℝ;

that is, c_{i}(aX) = a^{i} c_{i}(X) if {c_{i}(aX)} denotes the sequence of cumulants of aX.

iii) semi-invariance property:

χ.(α + a u) ≡ χ.α ˙+ aχ,

that is, for all a ∈ ℝ, we have c_{1}(X + a) = c_{1}(X) + a and c_{i}(X + a) = c_{i}(X) for all positive integers i greater than 1.

The next proposition follows from Equivalence (3.9) and paves the way to a symbolic handling of Abel sequences, as will become clearer in Section 6.

Proposition 4.3. [32] If κ_{α} is the α-cumulant umbra, then, for all positive integers i,

(4.2)  α^{i} ≃ κ_{α}(κ_{α} + α)^{i−1} and κi_{α} ≃ α(α − 1.α)^{i−1}.

The former Equivalence (4.2) has been assumed by Rota and Shen as definition of the cumulant umbra [74]. In terms of moments, this equivalence reads

a_{i} = ∑_{j=0}^{i−1} \binom{i}{j}a_{j}c_{i−j},

largely used in statistic framework [65, 87].

Since 𝔼[(χ)_{i}] = (−1)^{i−1}(i − 1)! for all positive integers i [34], the formula expressing the cumulants {c_{j}} in terms of the moments {a_{i}} is recovered from Equations (3.10) and (3.11), that is,

(4.3)  c_{j} = ∑_{λ⊢j} d_{λ}(−1)^{ν_{λ}−1}(ν_{λ} − 1)! a_{λ}.

Theorem 4.4 (Inversion theorem). [34] If κ_{α} is the α-cumulant umbra, then α ≡ β.κ_{α}.

In particular, from (3.8) we have

(4.4)  a_{i} = 𝒴_{i}(c_{1}, c_{2}, …, c_{i}) and a_{i} = ∑_{λ⊢i} d_{λ}c_{λ}.

The inversion theorem in Theorem 4.4 justifies Definition 4.2 since in terms of generating functions we have f(α, z) = exp[f(κ_{α}, z) − 1], corresponding to Equation (4.1). Moreover, since, from Inversion Theorem 4.4, any umbra α could be seen as the partition umbra of its cumulant umbra κ_{α}, it is possible to prove a more general result: every polynomial sequence of binomial type is completely determined by the sequence of its (formal) cumulants. On the other hand, we have already seen in the previous section that the polynomial umbra x.β.κ_{α} admits a sequence of moments of binomial type, see Equivalence (3.6). Due to Inversion Theorem 4.4, any polynomial sequence of binomial type is umbrally represented by a suitable polynomial umbra x.α. This statement gives the connection between Lévy processes and compound Poisson processes. We will get back to this connection in Example 4.5.

If we swap the umbra α and the umbra χ in χ.α, the resulting auxiliary umbra has moments equal to the factorial moments {a_{(i)}}_{i≥0} of α given in (3.2), that is, E[(α.χ)^{i}] = a_{(i)}. The umbra α.χ is called the α-factorial umbra. Factorial moments provide very concise expressions for moments of some discrete distributions, such as the binomial distribution [34], for example.

4.1. Randomized compound Poisson r.v.’s. The way to characterize the umbra representing a r.v. X is to determine the sequence of its moments {a_{i}}. When this sequence exists, this can be done by comparing its moment generating function E[exp(zX)] with the generating function of the umbra. Recall that, when E[exp(zX)] = f(z), for z in a suitable neighborhood of the origin, then f(z) admits an exponential expansion in terms of moments, which are completely determined by the corresponding distribution function (and vice-versa). In this case, the moment generating function encodes all the information of X, and the notion of identically distributed r.v.’s corresponds to the similarity among umbrae. In the symbolic method of moments, instead, the convergence of f(α, z) is not relevant [92]. This means that we can define the umbra whose moments are the same as the moments of a log-normal r.v.,^{16} even if this r.v. does not admit a moment generating function. In Section 3, generating functions of umbrae have been used to recover Wald’s martingale. But generating functions have more advantages. One of which is the encoding of symbolic operations among umbrae in terms of suitable operations among generating functions. For example, Equation (4.1) is encoded in the equivalence α ≡ β.κ_{α} of Inversion Theorem 4.4, and, more generally, the exponential of a generating function corresponds to the generating function of the polynomial α-partition umbra

(4.5)  f(x.β.α, z) = exp [x (f(α, z) − 1)].

Example 4.5 (Lévy processes). If {X_{t}} is a Lévy process, then its moment generating function is E[e^{zX_{t}}] = [ϕ(z)]^{t}, where ϕ(z) is the moment generating function of X_{1} [80]. This property parallels the infinite divisibility property of ℙ_{t} given at the beginning of this section. In particular,

E[e^{zX_{t}}] = exp[t log ϕ(z)] = exp[t k(z)],

where k(z) is the cumulant generating function of X_{1}, satisfying k(0) = 0. Comparing exp[t k(z)] with exp[t(f(κ_{α}, z) − 1)] given in Equation (4.5), with x replaced by t, we infer that

i) a Lévy process is umbrally represented by t.β.κ_{α} representing in turn a compound Poisson process of parameter t;

ii) the umbra κ_{α} represents the coefficients of k(z), which are the

cumulants of X_{1}; iii) the compound Poisson process of parameter t involves i.i.d.r.v.’s

{X_{i}}_{i≥0} whose moments are the cumulants of X_{1}.

The composition of generating functions corresponds to an iterated dot-product. Indeed, moments of the auxiliary umbra (γ.β).α can be obtained from (3.7), by replacing x by γ:

(4.6)  𝔼[Φ_{i}(γ.β)] = ∑_{λ⊢i} d_{λ} g_{ν_{λ}}a_{λ}.

The auxiliary umbra (γ.β).α is called the composition umbra of α and γ, since its i-th moment corresponds to the i-th coefficient of the composition f(γ, f(α, z) − 1).

16If X is a standard Gaussian r.v. N (0, 1), then the log-normal r.v. is Y = exp(X).

Theorem 4.6. [33] If α, γ ∈ 𝒜, then (γ.β).α ≡ γ.(β.α).

Due to Theorem 4.6, we may write the composition umbra as γ.β.α. The composition umbra (γ.β).α represents a compound randomized Poisson r.v., i.e., a random sum indexed by a randomized Poisson r.v. represented by the umbra γ.β. The composition umbra γ.(β.α) generalizes the random sum of i.i.d. compound Poisson r.v.’s with parameter 1, indexed by an integer r.v. X, i.e., a randomized compound Poisson r.v. with random parameter X represented by the umbra γ.

Paralleling the composition of generating functions, the compositional inverse of α is the auxiliary umbra α^{<−1>} having generating function f(α^{<−1>}, z) satisfying

f[α^{<−1>}, f(α, z) − 1] = f[α, f(α^{<−1>}, z) − 1] = 1 + z

or α.β.α^{<−1>} ≡ α^{<−1>}.β.α ≡ χ, since f(χ, z) = 1 + z.

Theorem 4.7. [27] An umbra α has a compositional inverse α^{<−1>} if and only if 𝔼[α] = a_{1} ≠ 0.

A very special umbra is the compositional inverse u^{<−1>} of the unity umbra u which satisfies

(4.7)  u.β.u^{<−1>} ≡ u^{<−1>}.β.u ≡ χ.

The following result is very useful in performing computations with dot-products.

Proposition 4.8. [34] χ ≡ u^{<−1>}.β ≡ β.u^{<−1>} and β.χ ≡ u ≡ χ.β.

The compositional inverse of an umbra is closely related to the Lagrange inversion formula, but we will use Sheffer polynomial sequences in Section 5 to compute its moments.

We end this section by introducing the connection between composition umbrae and convolutions of multiplicative functions. Multiplicative functions have been used to define free and Boolean cumulants by means of the lattice of non-crossing and interval partitions, respectively [58, 90]. This is an issue that we will return to again in Section 6.

Example 4.9. [32] Let Π_{n} be the set of all partitions of [n] equipped with the usual refinement order ≤, and denote the minimum and the maximum by 0_{n} and 1_{n}, respectively. Denote by |π| the number of blocks of π ∈ Π_{n}. If σ ∈ Π_{n} and σ ≤ π, then there exists a unique sequence of non-negative integers (k_{1}, k_{2}, …, k_{n}) with k_{1} + 2 k_{2} + ⋯ + n k_{n} = |σ| and k_{1} + k_{2} + ⋯ + k_{n} = |π| such that

[σ, π] = Π_{1}^{k_{1}} × Π_{1}^{k_{1}} × ⋯ × Π_{n}^{k_{n}},

with [σ, π] = {τ ∈ Π_{n} : σ ≤ τ ≤ π}. The sequence (k_{1}, k_{2}, …, k_{n}) is called the type of the interval [σ, π], and k_{i} is the number of blocks of π that are the union of i blocks of σ. A function 𝔣 : Π_{n} × Π_{n} → ℂ is said to be multiplicative if 𝔣(σ, π) = f_{1}^{k_{1}}f_{2}^{k_{2}} ⋯ f_{n}^{k_{n}} with σ ≤ π and f_{n} := 𝔣(0_{n}, 1_{n}). The Möbius function μ is multiplicative with μ_{n} = (−1)^{n−1}(n−1)!. Since the sequence {(−1)^{n−1}(n−1)!} is represented by the compositional inverse of the unity umbra u, the umbral counterpart of the Möbius function μ is u^{<−1>}. The zeta function ζ is multiplicative with ζ_{n} = 1. The umbral counterpart of the zeta function ζ is u. As the umbra u^{<−1>} is the compositional inverse of u, the Möbius function corresponds to the inverse of the zeta function. Indeed, the convolution ⋆ between two multiplicative functions 𝔣 and 𝔤 is defined by

(𝔣 ⋆ 𝔤)(σ, π) := ∑_{σ≤τ≤π} 𝔣(σ, τ) 𝔤(τ, π).

The function 𝔥 = 𝔣 ⋆ 𝔤 is also multiplicative with h_{n} = (𝔣 ⋆ 𝔤)(0_{n}, 1_{n}) and

(4.8)  h_{n} = ∑_{π∈Π_{n}} f_{π} g_{ℓ(π)},

where f_{π} := f_{1}^{k_{1}}f_{2}^{k_{2}} ⋯ f_{n}^{k_{n}} and k_{i} is the number of blocks of π of cardinality i. The convolution of two multiplicative functions (4.8) corresponds to a composition umbra since Equation (4.6) can be indexed by set partitions [21], that is,

𝔼[(γ.β.α)^{n}] = ∑_{π∈Π_{n}} g_{π} a_{|π|}.

The identity with respect to the convolution ⋆ is the Delta function Δ with Δ_{n} = δ_{1,n}. Its umbral counterpart is the singleton umbra χ. The Möbius function μ and the zeta function ζ are inverses of each other with respect to ⋆, that is, μ ⋆ ζ = ζ ⋆ μ = Δ. So, this last equality parallels Equivalences (4.7).

4.2. Lévy processes. In Section 3.1, we have stressed the central role played by Lévy processes in modeling financial markets. They are preferable to Brownian motion since deviation from normality can be explained by the jump or spike component. The symbolic representation of Lévy processes {X_{t}} by means of the family of auxiliary umbrae {t.α}_{t≥0} takes into account the strong interplay between these processes and infinitely divisible distributions, but does not reveal their jumpdiffusion structure. The same considerations hold when α is replaced by β.κ_{α}, which discloses the role played by cumulants of X_{1} in the distribution of X_{t}, but does not add more information on paths of a Lévy process. Instead, Lévy processes are also called jump-diffusion processes since they result from a summation of a Brownian motion and a compensated compound Poisson process,^{17} a property that characterizes unexpected changes in the paths, if any. The following theorem gives the jump-diffusion decomposition of a Lévy process.

17A compensated Poisson process is {N_{t} − λt} with {N_{t}} a Poisson process Po(λt).

Theorem 4.10. [80] If {X_{t}} is a Lévy process, then

X_{t} = m t + s B_{t} + ∑_{k=1}^{N_{t}} J_{k} − t a λ,

where m ∈ ℝ, s > 0, {B_{t}} is a standard Brownian motion, {N_{t}} is a Poisson process Po(λt), and {J_{k}}_{k≥1} are i.i.d.r.v.’s with E(J) = a < ∞.

All sources of randomness are mutually independent. In Theorem 4.10, the diffusion component of a Lévy process is represented by the Wiener process m t+s B_{t}, while the jump component is represented by the compensated compound Poisson process (∑_{k=1}^{N_{t}} J_{k} − taλ) . The symbolic method reduces the jump-diffusion decomposition to a compound Poisson process. This symbolic representation is obtained by taking into consideration the so-called Lévy–Khintchine formula, which splits the increment X_{1} in a summation according to the jump-diffusion property.

Theorem 4.11 (Lévy–Khintchine formula). [82] If {X_{t}} is a Lévy process, then E[e^{zX_{t}}] = [ϕ(z)]^{t}, with

(4.9) ϕ(z) = exp \{z m + \frac{1}{2} s^{2} z^{2} + ∫_{ℝ} (e^{zx} − 1 − zx1_{\{|x|≤1\}}) d ν(x)\} .

The triple (m, s^{2}, ν) is called Lévy triple, and ν is the Lévy measure. Since in Example 4.5 we have seen that f(t.β.κ_{α}, z) = [ϕ(z)]^{t}, with f(κ_{α}, z) = ϕ(z), what we need is to characterize the umbra κ_{α} according to (4.9). First, let us observe that the Lévy measure ν contains much information on X_{t}. If ν(ℝ) < ∞, almost all paths of X_{t} have a finite number of jumps on every compact interval. Under suitable normalization, this corresponds to the condition ∫_{ℝ} d ν(x) = 1, in agreement with 𝔼[1] = 1. If ν admits all moments, then Equation (4.9) may be rewritten as

(4.10) ϕ(z) = exp (c_{0}z + \frac{1}{2}s^{2}z^{2}) exp \{∫_{ℝ} (e^{zx} − 1 − zx) d ν(x)\}

where c_{0} = m + ∫_{{|x|≥1}} x d ν(x). Equation (4.10) allows us to split the umbra κ_{α} in a disjoint sum, as the following theorem shows.

Theorem 4.12. A Lévy process {X_{t}} is umbrally represented by the family of auxiliary umbrae

{t.β.[c_{0}χ ˙+ sη ˙+ γ]}_{t≥0},

where γ is the umbra associated to the Lévy measure, that is, f(γ, z) = 1 + ∫_{ℝ} (e^{zx} − 1 − zx) d ν(x), and η is an umbra with f(η, z) = 1 + z^{2}/2.

The jump-diffusion decomposition in Theorem 4.10 parallels the symbolic representation in Theorem 4.12, since from Theorem 3.4 we have

t.β.[c_{0}χ ˙+ sη ˙+ γ] ≡ t.β.[c_{0}χ ˙+ sη] + t.β.γ,

and the umbra (see Example 3.2)

β.[c_{0}χ ˙+ sη] ≡ c_{0} + β.(sη)

represents a Gaussian r.v. with mean c_{0} (called drift) and variance s^{2}. The compensated Poisson process is instead umbrally represented by the auxiliary umbra t.β.γ. Indeed, the umbra γ can be further decomposed into a pure point part, corresponding to the moment generating function ∫_{ℝ} e^{zx} d ν(x), and a singular part, corresponding to the moment generating function ∫_{ℝ}(1 + z x) d ν(x). The singular part of the Lévy measure is encoded in the umbra χ, which indeed does not have a probabilistic counterpart.

The singleton umbra also plays a special role with respect to the martingale property of a Lévy process. A centered Lévy process is a process with E[X_{t}] = 0 for all t ≥ 0, and, in this case, almost all paths of X_{t} have finite variation.

Theorem 4.13. [3] A Lévy process is a martingale if and only if its drift is c_{0} = 0.

A centered Lévy process is umbrally represented by the auxiliary umbra {t.β.(sη ˙+ γ)}_{t≥0}, where the explicit contribution of the singleton umbra is not visible anymore.

5. Sheffer polynomial sequences

Many characterizations of Sheffer polynomial sequences are given in the literature. We refer to the one involving generating functions.

Definition 5.1 (Sheffer sequences). A sequence {s_{i}(x)} of polynomials is called a Sheffer sequence if and only if its generating function is given by

1 + ∑_{i≥1} s_{i}(x)\frac{t^{i}}{i!} = h(t) exp\{x g(t)\},

with h(t) = h_{0} +h_{1}t+h_{2}t^{2} +⋯, g(t) = g_{1}t+g_{2}t^{2} +⋯, and h_{0}, g_{1} ≠ 0.

Applications of Sheffer sequences range from analysis to statistics, from combinatorics to physics. A good account of fields where these polynomials are employed is given in [16]. Their importance stems from the fact that many polynomial sequences, such as Laguerre polynomials, Meixner polynomials of first and second kind, Bernoulli polynomials, Poisson–Charlier polynomials, and Stirling polynomials are Sheffer polynomial sequences. Moreover, the time-space harmonic polynomials introduced in Section 3 are special Sheffer sequences as well.

Sheffer sequences can be considered as the core of the umbral calculus as introduced by Roman and Rota in [71]. There, one finds several chapters devoted to identities, properties and recurrences involving Sheffer sequences. Following these properties, Taylor has characterized their coefficients by using the classical umbral calculus [96]. In order to handle Sheffer polynomials by the symbolic method of moments, a special auxiliary umbra has been introduced in [27]. One of the main advantages of this symbolic representation is that properties of Sheffer sequences are reduced to a few equivalences.

In order to keep the length of the paper within reasonable bounds, here we limit ourselves to presenting some applications of the symbolic method involving Sheffer polynomial sequences. Whenever necessary, we recall the umbral equivalences we need. For a more complete treatment of this subject, the reader should consult [27].

The first application we deal with concerns the finding of solutions to linear recurrences. In many special combinatorial problems, the hardest part in their solution may be the discovery of an effective recursion. Once a recursion has been established, Sheffer polynomials are often a simple and general tool for finding solutions in closed form. The main contributions in this respect are due to Niederhausen [59] (see references therein), by means of finite operator calculus [76]. Using the symbolic method of moments, computations are simplified [36]. We give several examples below.

A second application involves the Lagrange inversion formula. This formula gives the coefficients of the compositional inverse of a formal power series. Indeed, in Definition 5.1 the formal power series g(t) has a compositional inverse, since g_{1} ≠ 0. There are many different proofs of this formula, many of them assuming that the involved power series are convergent and represent analytic functions on a disk. Using the symbolic representation of Sheffer sequences, the proof reduces to a simple computation [27].

A third and last application concerns Riordan arrays. A(n exponential) Riordan array is a pair (g(t), f(t)) of (exponential) formal power series, where g(t) is an invertible series and f(0) = 0. The pair defines an infinite lower triangular array according to the rule

d_{i,k} = i-th coefficient of g(t)\frac{[f(t)]^{k}}{k!} , for 0 ≤ k ≤ i < ∞.

Riordan arrays are used in a constructive way to find the generating function of many combinatorial sums, to characterize families of orthogonal polynomials, and to compute determinants of Hankel or Toeplitz matrices. Using the symbolic method of moments, we give a characterization of Riordan arrays. This characterization is very general since the same holds for the coefficients of Sheffer polynomials and connection constants [27]. If {s_{i}(x)} and {r_{i}(x)} are Sheffer sequences, the constants b_{i,k} in the expression

(5.1)  r_{i}(x) = ∑_{k=0}^{i}b_{i,k}s_{k}(x)

are known as connection constants. Several changes of bases among polynomial sequences may be efficiently computed by using connection constants. Connection constants, Riordan arrays, and coefficients of Sheffer polynomials share the same expression in terms of generalized Bell polynomials [27]. These polynomials are a generalization of complete Bell exponential polynomials (3.8) and also play a fundamental role in the multivariate version of the symbolic method.

The starting point of the symbolic representation of Sheffer sequences is the following definition.

Definition 5.2. A polynomial umbra ς_{x} is called a Sheffer umbra for (α, γ) if and only if ς_{x} ≡ α + x.γ^{∗}, where γ^{∗} ≡ β.γ^{<−1>} is called the adjoint umbra of γ.

Of course, to be well posed, in Definition 5.2 the umbra γ has to admit the compositional inverse γ^{<−1>}, that is, E[γ] = g_{1} ≠ 0. From Definition 5.2, the generating function of a Sheffer umbra is

f(ς_{x}, z) = f(α, z) exp{x [f^{<−1>}(γ, z) − 1]},

in agreement with Definition 5.1.

5.1. Solving linear recurrences. In order to show how to solve linear recurrences by using the symbolic method, we present two examples, the former involving a uniform r.v., the latter Dyck paths [27].

Example 5.3. Suppose we are asked to solve the difference equation

(5.2)  q_{i}(x + 1) = q_{i}(x) + q_{i−1}(x)

under the condition ∫_{0}^{1} q_{i}(x) dx = 1 for all positive integers i. This initial condition is equivalent to 𝔼[q_{i}(U)] = 1, where U is a uniform r.v. on the interval (0, 1). Set q_{i}(x) = s_{i}(x)/i!. Then, from (5.2), we have

(5.3)  s_{i}(x + 1) = s_{i}(x) + i s_{i−1}(x).

A polynomial umbra ς_{x} is a Sheffer umbra if and only if there exists an umbra γ, equipped with compositional inverse, which satisfies ς_{γ+x.u} ≡ χ + ς_{x} [27], that is, if and only if s_{i}(γ + x.u) = s_{i}(x) + i s_{i−1}(x) for all positive integers i. By comparing this last equality with (5.3), we have to choose as umbra γ the unity umbra u. Since u^{∗} ≡ χ, we have q_{i}(x) ≃ (x.χ + α)^{i}/i! for α depending on the initial condition. By recalling that ∫_{0}^{1} p(x) dx = 𝔼[p(−1.ι)] for any p(x) ∈ ℝ[x] [17], we infer

∫_{0}^{1} q_{i}(x) dx = 1 if and only if 𝔼[q_{i}(−1.ι)] = 1

if and only if 𝔼[s_{i}(−1.ι)] = i!.

Let the sequence {i!} be represented by the umbra ū, called the Boolean unity umbra.^{18} Then, from 𝔼[q_{i}(−1.ι)] = 1, we have

−1.ι.χ + α ≡ ū if and only if α ≡ ū + ι.χ.

Solutions of the recursion (5.2) are moments of the Sheffer umbra (ι + β + x.u).χ divided by i! □

Example 5.4. A ballot path takes up-steps (u) and right-steps (r), starting at the origin and staying weakly above the diagonal. For example, ururuur is a ballot path (starting at the origin and) terminating at the point (3, 4) . Let D (i, m) be the number of ballot paths terminating at (i, m) which do not contain the pattern (substring) urru.

It is not difficult to see that D(i, m) satisfies the recurrence relation

D(i, m) − D(i − 1, m) = D(i, m − 1) − D(i − 2, m − 1) + D(i − 3, m − 1)

under the initial conditions D(i, i) = D(i − 1, i) and D (0, 0) = 1. Set D(i, m) = s_{i}(m)/i!. The previous recurrence relation may be written as

(5.4) s_{i}(m)−i s_{i−1}(m) = s_{i}(m−1)−(i)_{2} s_{i−2}(m−1)+(i)_{3} s_{i−3}(m−1),

with initial condition s_{i}(i) = i s_{i−1}(i). Replace m by x in Equation (5.4). The umbral counterpart of the Sheffer identity^{19} [27] is ς_{x+y} ≡ ς_{x}+y.γ^{∗}, which yields ς_{x} ≡ ς_{x−1}+γ^{∗} for y = −1. Since s_{i}(x)−i s_{i−1}(x) ≃ (ς_{x}−χ)^{i}, we may replace ς_{x} by ς_{x−1} + γ^{∗} so that

(5.5)  (ς_{x} − χ)^{i} ≃ (ς_{x−1} + γ^{∗} − χ)^{i}.

Since

(ς_{x−1} + γ^{∗} − χ)^{i} ≃ ς_{x−1}^{i} + i ∑_{j=0}^{i−1} \binom{i − 1}{j}ς_{x−1}^{i−j−1} [\frac{(γ^{∗})^{j+1}}{j + 1} − (γ^{∗})^{j}] ,

from Equivalence (5.5) and Equation (5.4), we infer

f(γ^{∗}, t) = 1 + t + ∑_{k≥3} t^{k} = (1 − t^{2} + t^{3})/(1 − t).

18Its moments are called Boolean cumulants, introduced in Section 6. For Boolean cumulants, the umbra ū plays the same role as the unity umbra u for classical cumulants.

19Sheffer identity: A sequence {s_{i}(x)} of Sheffer polynomials satisfies s_{i}(x+y) = ∑_{k=0}^{i} \binom{i}{k}s_{k}(x)s_{i−k}(y).

The initial condition s_{i}(i) = i s_{i−1}(i) is equivalent to looking for an umbra α satisfying

(5.6)  (α + i .γ^{∗})^{i} ≃ i (α + i .γ^{∗})^{i−1}.

Assume α ≡ ū + ζ. In this case the initial condition (5.6) gives (ζ + i .γ^{∗})^{i} ≃ ε^{i}, and then ζ^{i} ≃ −i [γ^{<−1>}]^{i} for all positive integers i. The solution of the recurrence relation (5.4) is

\frac{s_{i}(x)}{i!} ≃ \frac{(ū + ζ + x.γ^{∗})^{i}}{i!} . □

Two special classes of Sheffer polynomials are often employed in applications: polynomial sequences of binomial type and Appell sequences. In Section 3, we have already seen that polynomial sequences of binomial type are umbrally represented by the polynomial umbra x.β.γ, which is a Sheffer umbra for (ε, γ).

An Appell sequence is a sequence of polynomials {p_{i}(x)} satisfying the identity

\frac{d}{dx} p_{i}(x) = i p_{i−1}(x), for i = 1, 2, ….

This sequence of polynomials is represented by a Sheffer umbra for (α, χ). The resulting umbra α+x.u is called Appell umbra. Examples of Appell polynomials are the Bernoulli polynomials q_{i}(x) = E[(ι + x.u)^{i}] introduced in Section 2, the Euler polynomials

q_{i}(x) = E [x.u + (\frac{1}{2}[ξ − u])^{i}] ,

with ξ the Euler umbra, and the time-space harmonic polynomials Q_{i}(x, t) = E[(x + t.α)^{i}] introduced in Section 3.

5.2. Lagrange inversion formula. The Lagrange inversion formula gives the coefficients of the compositional inverse of a formal power series and thus of the compositional inverse of an umbra γ. Following Rota, Shen and Taylor [77], the proof we propose relies on the circumstance that any sequence of binomial type can be represented by Abel polynomials {x (x − i a)^{i−1}}, a ∈ ℝ. We use umbral Abel polynomials {x (x − i.γ)^{i−1}}. An exhaustive treatment of umbral Abel polynomials can be found in [63]. In contrast to [77], the employment of the singleton umbra reduces the proof of the Lagrange inversion formula to a simple computation. We need a preliminary definition.

Definition 5.5. If {g_{i}} is a sequence umbrally represented by an umbra γ, then its derivative umbra γ_{D} has moments g_{D,i} = i g_{i−1} for all positive integers i.

The derivative umbra has moments obtained from a symbolic derivation with respect to γ, that is, i g_{i−1} = 𝔼[i γ^{i−1}] = 𝔼[∂_{γ}γ^{i}].

Theorem 5.6 (Abel representation of binomial sequences).

[27] If γ is an umbra equipped with a compositional inverse, then

(5.7)  (x.γ∗_{D})^{i} ≃ x (x − i.γ)^{i−1}, for i = 1, 2, ….

Theorem 5.6 includes the well-known Transfer Formula [71].

Theorem 5.7 (Lagrange inversion formula). [27]

(γ_{D}^{<−1>})^{i} ≃ (−i.γ)^{i−1}, for i = 1, 2, ….

Proof. In (5.7), replace x by the singleton umbra χ and observe that

(χ.β.γ_{D}^{<−1>})^{i} ≃ (γ_{D}^{<−1>})^{i} ≃ χ (χ − i .γ)^{i−1},

due to χ.β ≡ u (see Proposition 4.8). The result follows by applying the binomial expansion to χ (χ−i.γ)^{i−1} and by recalling that χ^{k+1} ≃ 0 for k = 1, 2, ..., i − 1. □

As −i.γ ≡ i.(−1.γ), moments of γ_{D}^{<−1>} can be computed by using the inverse of γ. The umbra γD has first moment equal to 1. It is possible to generalize Theorem 5.7 to umbrae having first moment different from 1, see [27].

5.3. Generalized Bell polynomials. By Theorem 5.6, the umbral analog of the well-known Abel identity^{20} is:

(5.8)  (x + y)^{i} ≃ ∑_{k=0}^{i} \binom{i}{k}(x + k.γ)^{i−k}y(y − k.γ)^{k−1}.

The polynomials

𝔹_{i,k}^{(γ)}(x) = \binom{i}{k}(x + k.γ)^{i−k},

given in Equation (5.8), are called generalized Bell umbral polynomials since E[𝔹_{i,k}^{(γ)}(0)] = B_{i,k} (g_{D,1}; . . . . . .; g_{D,i−k+1}) , where B_{i,k} are the partial Bell exponential polynomials given in Equation (3.7), and {g_{D,i}} are the moments of the umbra γ_{D}. Paralleling the complete exponential Bell polynomials, the generalized complete Bell umbral polynomials are

𝕐_{i}^{(γ)}(x) = (x + β.γ_{D})^{i} ≃ ∑_{k=0}^{i}\binom{i}{k}(x + k.γ)^{i−k} = ∑_{k=0}^{i}𝔹_{i,k}^{(γ)}(x).

Generalized Bell polynomials are useful to give a unifying representation of the coefficients of Sheffer polynomials, of Riordan arrays, and of connection constants, as we show in the following. Assume 𝔼[γ] = 1. Generalizations to the case 𝔼[γ] ≠ 1 are given in [27]. We need a preliminary definition.

Definition 5.8. If {g_{i}} is a sequence umbrally represented by an umbra γ, then its primitive umbra γ_{P} has moments g_{P,i} = g_{i+1}/(i + 1) for all positive integers i.

20The Abel identity is: (x + y)^{i} = ∑_{k=0}^{i} \binom{i}{k}(x + k a)^{i−k}y(y − k a)^{k−1}, a ∈ ℝ.

The primitive umbra has moments obtained from a symbolic integration with respect to γ, that is,

\frac{g_{i+1}}{i + 1} = 𝔼 (\frac{γ^{i+1}}{i + 1}) = 𝔼 (∫ γ^{i}dγ) .

For 0 ≤ k ≤ i, let c_{i,k} denote the k-th coefficient of the i-th Sheffer polynomial s_{i}(x).

Theorem 5.9 (Coefficients of Sheffer polynomials). [27] If {s_{i}(x)} is a polynomial sequence umbrally represented by a Sheffer umbra for the pair (α, γ), then c_{i,k} = 𝔼 [𝔹_{i,k}^{(γ_{P})}(α)] for all positive integers i and 0 ≤ k ≤ i.

Theorem 5.10 (Riordan array). [27] The elements of an exponential Riordan array for the pair (f(α, t), f^{<−1>}(γ, t)−1) are umbrally represented by d_{i,k} = 𝔼 [𝔹_{i,k}^{(γ_{P})}(α)] for all positive integers i and 0 ≤ k ≤ i.

So, the theory of Riordan arrays and that of Sheffer sequences are the two sides of the same coin.

Theorem 5.11 (Connection constants). [27] If {s_{i}(x)} is a polynomial sequence umbrally represented by a Sheffer umbra for (α, γ) and {r_{i}(x)} is a polynomial sequence umbrally represented by a Sheffer umbra for (ξ, ζ), then b_{i,k} = 𝔼 [𝔹_{i,k}^{(ϕ_{P})}(ϱ)] with {b_{i,k}} given in (5.1) and ϱ ≡ (−1.α.β.γ + ξ.β.ζ).ζ^{∗}, ϕ ≡ ζ.γ^{∗}.

6. Fast symbolic computation

There are many packages devoted to numerical/graphical statistical toolsets but not doing algebraic/symbolic computations. Quite often, the packages filling this gap are not open source. One of the softwares largely used in the statistical community is R.^{21} R is a much stronger numerical programming environment, and the procedures including symbolic software are not yet specifically oriented for statistical calculations. The availability of a widely spread open source symbolic platform will be of great interest, especially with interface capabilities to external programs. The conceptual aspects related to symbolic methods involve more strictly algebraic and combinatorial issues, and the topics introduced in the previous sections are an example.

What we regard as symbolic computation is now evolving towards a universal algebraic language which aims at combining syntactic elegance and computational efficiency. The classical umbral calculus is no doubt an old topic. However, the version here reviewed has not only the virtue of improving the syntactic elegance of the umbral theory but also to offer a different and deeper viewpoint, for example on the role played by cumulants in working with compound Poisson processes.

21Available at (http://www.r-project.org/)

In this section, we show how this different and deeper viewpoint allows us to reach computational efficiency as by-product. These examples should push statisticians to climb the steep learning curve of a new method and mathematicians to believe what Terence Tao claimed in his blog: “A different choice of foundations can lead to a different way of thinking about the subject, and thus to ask a different set of questions and to discover a different set of proofs and solutions. Thus it is often of value to understand multiple foundational perspectives at once, to get a truly stereoscopic view of the subject.”

Two examples are proposed here. The first concerns efficient algorithms aiming at the computation of experimental measurements of cumulants, known in the literature as polykays. The sampling behavior of polykays is much simpler than sample moments [54], but their employment was not so widespread in the statistical community due to the past computational complexity in recovering their expression. Only recently, the symbolic method of moments has contributed in speeding up the procedures for their computations [23]. The second example shows how to recover classical, Boolean, and free cumulants from moments and vice-versa, by using only one family of umbral polynomials. Quite surprisingly this algorithm is an application of the symbolic representation of Sheffer sequences making any of the existing procedures obsolete. This part has a consistent theoretical background devoted to the symbolic representation of classical, Boolean and free cumulants, their generalization to Abel-type cumulants, and some open problems which are still under investigation.

6.1. k-statistics and polykays. In dealing with computations, the approach used to compute U-statistics and estimators of multivariate moments is not always efficient. So, further manipulations are required. As an example, we show how to generate k-statistics, which are Ustatistics for cumulants, by using a suitable compound Poisson r.v. The i-th k-statistic k_{i} is the unique symmetric unbiased estimator of the cumulant c_{i} of a given statistical distribution [94], that is, E[k_{i}] = c_{i}.

A frequently asked question is: why is it so relevant to get efficiency in these calculations? Usually higher order objects require enormous amounts of data in order to get a good accuracy. Nevertheless, there are different areas, such as astronomy [69], astrophysics [40], biophysics [55], and neuroscience [93], in which there is need to compute high order k-statistics in order to detect a Gaussian population. Undoubtedly, a challenging target is to have efficient procedures to deal with the involved huge amount of algebraic and symbolic computations.

The theory of k-statistics has a long history beginning with Fisher [41], who rediscovered the theory of half-invariants of Thiele [98]. Fisher introduced k-statistics (single and multivariate) as new symmetric functions of a random sample, with the aim of estimating cumulants without using moment estimators. Dressel [38] developed a theory of more general functions, revisited later by Tukey [101], who called them polykays. The whole subject is described by Stuart and Ord [94] in great detail. In the 1980s, tensor notation was used by Speed [88] and extended to polykays and single k-statistics. This extension uncovers the coefficients defining polykays to be values of the Möbius function on the lattice of set partitions. As a consequence, Speed used the set-theoretic approach to symmetric functions introduced by Doubilet [37]. In the same period, McCullagh [54] simplified the tensor notation of Kaplan [48] by introducing the notion of generalized cumulants. Symbolic operators for expectation and the derivation of unbiased estimators for multiple sums were introduced by Andrews et al. [1].

A classical way to recover k-statistics is to use Theorem 2.6 in (4.3). The resulting U-statistics are expressed in terms of augmented symmetric polynomials. Then, to get k-statistics, we have to express augmented symmetric polynomials in terms of power sum symmetric polynomials.

Theorem 6.1. [21] For i ≤ n, we have^{22}

(6.1)  (χ.α)^{i} ≃ ∑_{λ⊢i} \frac{(χ.χ)^{ν_{λ}}}{(n)_{ν_{λ}}} d_{λ} ∑_{π∈Π_{ν_{λ}}}(χ.χ)^{.π}(n.α)_{S_{π}},

where S_{π} is the subdivision of the multiset P_{λ} = {α^{(r_{1})}, α^{2}^{(r_{2})}, …} corresponding to the partition π ∈ Π_{ν_{λ}}.

A significant computational cost is required to find set partitions or multiset subdivisions, although these procedures are already optimized [22]. By using the symbolic method, an improvement in the performance is achievable. The main idea is to recover k-statistics as cumulants of compound Poisson r.v.’s replacing the augmented symmetric polynomials by exponential polynomials (3.5). The implementation is then sped up since there is no need to compute set partitions.

To this aim, we introduce the multiplicative inverse of an umbra.

Definition 6.2. The multiplicative inverse of an umbra α is the umbra γ satisfying αγ ≡ u.

In dealing with a saturated umbral calculus, the multiplicative inverse of an umbra is not unique, but any two multiplicatively inverse umbrae of the umbra α are similar. From Definition 6.2, we have a_{n}g_{n} = 1 for all non-negative integers n, i.e., g_{n} = 1/a_{n}, where {a_{n}} and {g_{n}} are umbrally represented by the umbrae α and γ, respectively. In the following, the multiplicative inverse of an umbra α will be denoted by the symbol 1/α. The following proposition shows how to express cumulants of an umbra by suitable polynomials linked to moments of compound Poisson r.v.’s.

22For the notation, see Subsection 2.1.

Theorem 6.3. [21] If c_{i}(y) = E[(n.χ.y.β.α)^{i}], then

(χ.α)^{i} ≃ c_{i} (\frac{χ.χ}{n.χ}) , for i = 1, 2, ….

Thus, moments of the α-cumulant umbra can be computed by evaluating the umbral polynomials c_{i}(y), after having replaced the indeterminate y by the umbra χ.χ/n.χ. The polynomials c_{i}(y) are moments of a summation of n uncorrelated cumulant umbrae of y.β.α, the latter in its turn being a polynomial α-partition umbra representing a compound Poisson r.v. of parameter y. The last step is to express the polynomials c_{i}(y) in terms of power sum symmetric polynomials n.α^{r} for r ≤ n.

Theorem 6.4. [21] If

p_{n}(y) = ∑_{k=1}^{n}(−1)^{k−1}(k − 1)! S(n, k) y^{k},

where S(n, k) are the Stirling numbers of the second kind, then

(6.2)  (χ.α)^{i} ≃ ∑_{λ⊢i} d_{λ} p_{λ} (\frac{χ.χ}{n.χ}) (n.α)^{r_{1}}(n.α^{2})^{r_{2}} ⋯,

with p_{λ}(y) = [p_{1}(y)]^{r_{1}}[p_{2}(y)]^{r_{2}} · · · .

Comparing Equivalence (6.2) with Equivalence (6.1), the reduction of complexity is evident since no set partitions are involved.

Theorem 6.4 has been generalized to polykays. Polykays are symmetric statistics k_{r,..., t} satisfying

E[k_{r,..., t}] = c_{r} ⋯ c_{t},

where c_{r}, …, c_{t} are cumulants. For simplicity, in the following we just deal with two subindices k_{r,t}, the generalization to more than two being straightforward. As a product of uncorrelated cumulants, (χ.α)^{r}(χ′.α′)^{t} is the umbral counterpart of the polykay k_{r,t}, with χ, χ′ uncorrelated singleton umbrae, and α, α′ uncorrelated and similar umbrae. For polykays, the polynomial

p_{r,t}(y) = ∑_{(λ⊢r, η⊢t)} y^{ν_{λ}+ν_{η}} (n)_{ν_{λ}+ν_{η}} d_{λ+η} a_{λ+η}

plays the same role as the polynomial c_{i}(y) in Theorem 6.3. The following theorem is the analog of Theorem 6.4 for polykays.

Theorem 6.5. [21] If q_{r,t} is the umbral polynomial obtained via p_{r,t}(y) by replacing y^{ν_{λ}+ν_{η}} by

\frac{(χ.χ)^{ν_{λ}}(χ′.χ′)^{ν_{η}}}{(n.χ)^{ν_{λ}+ν_{η}}} \frac{d_{λ}d_{η}}{d_{λ+η}} ,

then (χ.α)^{r}(χ′.α′)^{t} ≃ q_{r,t}.

Now, let us consider a polynomial umbra ρ_{y,k} whose first k moments are all equal to y, that is, ρ_{y,k}^{i} ≃ (χ.y.β)^{i} for i = 0, 1, 2, …, k, and all the remaining moments are zero. If k = max{r, t}, then

[n.(ρ_{y,k} α)]^{r+t} ≃ ∑_{λ⊢(r+t)} d_{ξ} (χ.ρ_{y,k})_{ξ} (n.α)^{s_{1}} (n.α^{2})^{s_{2}} ⋯,

which in turn gives the polynomials p_{r,t}(y) in terms of power sum symmetric polynomials, as [n.(ρ_{y,k} α)]^{r+t} ≃ p_{r,t}(y). This device together with Theorem 6.5 speeds up the computation of polykays.

The resulting algorithm can be fruitfully employed also in the estimation of products of moments. Indeed, unbiased estimators of products of moments can be computed by using polykays, as we show in the following. Let us consider a set {α_{1}, α_{2}, …, α_{n}} of n uncorrelated umbrae similar to α. Set α^{.π} = α_{i_{1}}^{|B_{1}|}α_{i_{2}}^{|B_{2}|} ⋯ α_{i_{k}}^{|B_{k}|}, where π = {B_{1}, B_{2}, …, B_{k}} is a partition of [n] and i_{1}, i_{2}, …, i_{k} are distinct integers chosen in [n]. We have E[α^{.π}] = a_{λ}, where λ is the type of the partition π. The formula giving products of moments in terms of polykays is [21]

(6.3)  α^{.π} ≃ ∑_{\begin{matrix}τ∈Π_{n} \\ τ≤π \end{matrix}}(χ.α)^{.τ}.

By using the Möbius inversion

(6.4)  (χ.α)^{.π} ≃ ∑_{\begin{matrix}τ∈Π_{n} \\ τ≤π \end{matrix}}μ(τ, π) α^{.π},

we get polykays in terms of products of moments. Equivalences (6.3) and (6.4) have also a different meaning. Taking the evaluation of both sides, formulas connecting classical cumulants and moments are obtained in terms of multiplicative functions. These formulas are

a_{π} = ∑_{τ≤π}c_{τ} for all π, if and only if c_{π} = ∑_{τ≤π}μ(τ, π)a_{τ} for all π,

where a_{π} = a_{|B_{1}|}a_{|B_{2}|} ⋯ a_{|B_{k}|} (the same for c_{π}), and

μ(τ, π) = (−1)^{n−r}(2!)^{r_{3}}(3!)^{r_{4}} ⋯ ((n − 1)!)^{r_{n}}

is the Möbius function for the classical partition lattice, with n = |τ|, r = |π|, and λ(τ, π) = (1^{r_{1}}, 2^{r_{2}}, …, n^{r_{n}}) the class of [τ, π].

Since all these formulas can be generalized to correlated umbral monomials (see Definition 2.12), the next section is devoted to the multivariate umbral calculus and its applications in probability and statistics. Since any umbra is a composition umbra, the key to generalize this symbolic calculus to k-tuples of umbral monomials is the multivariate Faà di Bruno formula, which is related to the composition of multivariate formal power series. For the fast algorithms computing multivariate k-statistics and multivariate polykays by using the symbolic method, the reader is referred to [23]. In Appendix 1, Table 5 and 6, we show the advantage of these fast algorithms compared to existing procedures.

6.2. Computing cumulants through generalized Abel polynomials.

6.2.1. Classical, Boolean, and free cumulants. In Section 4, cumulants have been introduced as a sequence of numbers providing an alternative to moments of a distribution. Lévy processes are an example because their symbolic representation involves cumulants of X_{1}. Since cumulants linearize the convolution of probability measures in classical, Boolean, and free probability, they have one more advantage: cumulants are a tool to recognize independent r.v.’s. The linearity of classical convolutions corresponds to tensor independent r.v.’s [39]. The linearity of Boolean convolutions corresponds to Boolean independent r.v.’s for which the factorization rule holds, but without the commutative property of the resulting moments. The Boolean convolution [90] was constructed starting from the notion of partial cumulants, which were extensively used in the context of stochastic differential equations. The linearity of free convolutions corresponds to free r.v.’s [104]. The latter belong to the non-commutative probability theory introduced by Voiculescu [103] with a view to tackle problems in the theory of operator algebras and random matrix theory. The combinatorics underlying classical, Boolean, and free cumulants is closely related to the algebra of multiplicative functions (see Example 4.9). Roughly speaking, one can work with classical, Boolean, and free cumulants using the lattice of set, interval, and non-crossing partitions,^{23} respectively.

Since the linearity of cumulants gives a quick access to test whether a given probability measure is a convolution, the availability of a fast procedure for computing cumulants — classical, Boolean, and free — from the sequence of moments turns out to be useful. The symbolic method of moments replaces the combinatorics of partition lattices with the algebra of umbral Abel polynomials given in Theorem 5.6. Indeed, all these families of cumulants share the same parametrization in terms of umbral Abel polynomials. The result is an algorithm providing just one formula for computing all types of cumulants in terms of moments and vice-versa.

23The combinatorics underlying non-crossing partitions has been first studied by Kreweras [50] and Poupard [68]. Within free probability, non-crossing partitions are extensively used by Speicher [90].

The aim of this section is to introduce this parametrization, to describe the resulting algorithm and some problems that this parametrization lives open. Moreover we define an Abel-type cumulant and explain why the connection among classical, Boolean, and free cumulants cannot be encoded in a straight way in the formal power series language, as pointed out also in [58].

Let us consider again the umbral Abel polynomials {x(x − i.α)^{i−1}}.

Definition 6.6. The umbral polynomial 𝔞_{i}^{(m)}(x, α) = x (x + m.α)^{i−1}, with m ∈ ℤ, is called the generalized umbral Abel polynomial.

We have already used generalized umbral Abel polynomials in Section 4 for classical cumulants. Indeed, from the latter (4.2), moments of the α-cumulant umbra are umbrally equivalent to 𝔞_{i}^{(−1)}(α, α)

(6.5)  κi_{α} ≃ α(α − 1.α)^{i−1} ≃ 𝔞_{i}^{(−1)}(α, α).

The former (4.2) inverts Equivalence (6.5) and can be restated as

α^{i} ≃ κ_{α}(κ_{α} + α)^{i−1} ≃ 𝔞_{i}^{(1)}(κ_{α}, α).

The connection between these two types of generalized umbral Abel polynomials holds in greater generality as the following proposition shows.

Theorem 6.7 (First Abel inversion theorem). [64]

(6.6) α^{i} ≃ 𝔞_{i}^{(m)}(γ, γ) ≃ γ (γ + m.γ)^{i−1}

if and only if γ^{i} ≃ 𝔞_{i}^{(−m)}(α, γ) ≃ α (α − m.γ)^{i−1}.

A natural question arising here is the following: does the sequence {𝔞_{i}^{(−m)}(α, α)}, resulting from (6.5) for positive integers m, have the three main algebraic properties of cumulants?

Let us consider the umbra ᾱ ≡ α ū whose moments are 𝔼[ᾱ^{i}] = i! a_{i}, and the sequence {𝔞_{i}^{(−2)}(ᾱ, ᾱ)}. By the First Abel Inversion Theorem, given in Theorem 6.7, we have

(6.7) η̄_{α}^{i} ≃ 𝔞_{i}^{(−2)}(ᾱ, ᾱ) ≃ ᾱ (ᾱ − 2.ᾱ)^{i−1}

if and only if ᾱ^{i} ≃ 𝔞_{i}^{(2)}(η̄_{α}, ᾱ) ≃ η̄_{α} (η̄_{α} + 2.ᾱ)^{i−1}.

The umbra η̄_{α} is called the ᾱ-Boolean cumulant umbra. The three main algebraic properties of cumulants can be stated as follows:

i) (additivity) η̄_{ξ} ≡ η̄_{α} +^{.} η̄_{γ} if and only if −1.ξ̄ ≡ −1.ᾱ +^{.} −1.γ̄; ii) (homogeneity) η̄_{cα} ≡ cη̄_{α}, c ∈ ℝ; iii) (semi-invariance) η̄_{ξ} ≡ η̄_{α} +^{.} c χ if and only if −1.ξ̄ ≡ −1.ᾱ +^{.}

−1.(cū).

Its moments are {𝔼[𝔞_{i}^{(−2)}(ᾱ, ᾱ)]} whose expansion in terms of moments is

𝔼[ᾱ(ᾱ − 2.ᾱ)^{i−1}] = ∑_{λ⊢i} d_{λ}(−1)^{ν_{λ}−1}ν_{λ}! ā_{λ}, with ā_{λ} = r_{1}! a_{1}^{r_{1}} r_{2}! a_{2}^{r_{2}} ⋯

by Definition 6.6.

Theorem 6.8. [32] If η̄_{α} is the ᾱ-Boolean cumulant umbra, then

(6.8)  η̄_{α} ≡ ū^{<−1>}.β.ᾱ and ᾱ ≡ ū.β.η̄_{α}.

Comparing Equivalences (6.8) with the equivalences characterizing classical cumulants,

(6.9)  κ_{α} ≡ χ.α ≡ u^{<−1>}.β.α and α ≡ β.κ_{α} ≡ u.β.κ_{α},

the analogy is evident. The same analogy has a counterpart within formal power series, as the following remark shows.

Remark 6.9. Let us observe that the ordinary formal power series

f(z) = 1 + ∑_{i≥1} f_{i}z^{i}

is the generating function of the umbra ᾱ, that is, f(ᾱ, z) = f(z). If we consider the umbra γ satisfying f(γ̄, z) = g(z) = 1 + ∑_{i≥1} g_{i}z^{i}, then the auxiliary umbra ᾱ.β.γ̄ has generating function f[g(t) − 1]. This umbra is the counterpart of the convolution of multiplicative functions on the lattice of interval partitions. Indeed, let us consider the lattice of interval partitions (I_{i}, ≤). A partition π of Π_{i} is called an interval partition if each block B_{j} of π is an interval [a_{j}, b_{j}] of [i], that is, B_{j} = [a_{j}, b_{j}] = {x ∈ [j] | a_{j} ≤ x ≤ b_{j}}, where a_{j}, b_{j} ∈ [i]. The type of each interval [σ, τ] in I_{i} is the same as the one given in Example 4.9. A convolution ⋄ is defined on the algebra of multiplicative functions by

(6.10)  (𝔤 ⋄ 𝔣)(σ, π) := ∑_{{σ≤τ≤π \atop τ∈I_{i}}} 𝔤(σ, τ)𝔣(τ, π),

with 𝔣 and 𝔤 multiplicative functions. Then 𝔥 = 𝔤 ⋄ 𝔣 is a multiplicative function and

h_{i} = ∑_{τ∈I_{i}} f_{ℓ(τ)}g_{τ}.

The convolution (6.10) of two multiplicative functions corresponds to a composition umbra since E[(ᾱ.β.γ̄)^{i}] = ∑_{π∈I_{i}} a_{|π|} g_{π}.

Remark 6.10. There is an interesting parallelism between the umbra having all classical cumulants equal to 1 and the umbra having all Boolean cumulants equal to 1. The Bell umbra β is the unique umbra (up to similarity) having {1}_{i≥1} as sequence of classical cumulants. Its i-th moment is the i-th Bell number, that is, the number of partitions of an i-set. The umbra α satisfying ᾱ ≡ (2ū)_{D} is the unique umbra (up to similarity) having {1}_{i≥1} as sequence of Boolean cumulants. The i-th moment of this umbra is 2^{i−1}, that is, the number of interval partitions |I_{i}|.

Free cumulants are closely related to the Lagrange inversion formula. They were introduced by Voiculescu [104] as semi-invariants via the Rtransform. Indeed, let us consider a non-commutative r.v. X, i.e., an element of a unital non-commutative algebra 𝔸. Suppose ϕ : 𝔸 → ℂ is a unital linear functional. The i-th moment of X is the complex number m_{i} = ϕ(X^{i}), while its generating function is the formal power series M(z) = 1 + ∑_{i≥1} m_{i}z^{i}. The non-crossing (or free) cumulants of X are the coefficients r_{i} of the ordinary power series R(z) = 1 + ∑_{i≥1} r_{i}z^{i} which satisfies the identity

(6.11)  M(z) = R[zM(z)].

If we set M(z) = f(ᾱ, z) and R(z) = f(𝔎̄_{α}, z), then f(ᾱ_{D}, z) = 1 + zM(z), and Equation (6.11) gives ᾱ ≡ 𝔎̄_{α}.β.ᾱ_{D}, which justifies the following definition.

Definition 6.11. The umbra 𝔎̄_{α} satisfying (−1.𝔎̄_{α})_{D} ≡ ᾱ_{D}^{<−1>} is called the free cumulant umbra of α.

Denote the auxiliary umbra (ᾱ_{D}^{<−1>})_{P} by 𝔏ᾱ. We call this umbra Lagrange involution of α. Then, from Definition 6.11, we have

𝔎̄_{α} ≡ −1.𝔏ᾱ,

which is an alternative definition of the free cumulant umbra. The three main algebraic properties of cumulants can be stated as follows:

i) (additivity) 𝔎̄_{ξ} ≡ 𝔎̄_{α} +^{.} 𝔎̄_{γ} if and only if − 1.𝔏ξ̄ ≡ −1.𝔏ᾱ +^{.} −1.𝔏γ̄;

ii) (homogeneity) 𝔎̄_{cα} ≡ c𝔎̄_{α}, c ∈ ℝ; iii) (semi-invariance) 𝔎̄_{ξ} ≡ 𝔎̄_{α} +^{.} cχ if and only if − 1.𝔏ξ̄ ≡

−1.𝔏ᾱ +^{.} −1.𝔏_{cū}.

From Definition 6.11 we have

𝔎̄_{α} ≡ ᾱ.β.ᾱ_{D}^{<−1>} and ᾱ ≡ 𝔎̄_{α}.β.(−1.𝔎̄_{α})_{D}^{<−1>}.

If we compare these two equivalences with those corresponding to the α-cumulant umbra (6.9) and the α-Boolean cumulant umbra (6.8), respectively, the difference is quite obvious. Similarly, while the characterization of classical cumulants and Boolean cumulants is closely related to the composition of generating functions, this is not true for free cumulants. To be more precise, there does not exist a generating function U(t) such that, if R(t) is the free cumulant generating function of M(t), then M(t) = U[R(t)−1], see also [32]. The same difference arises for convolutions of multiplicative functions on the lattice of non-crossing partitions, as the following remark shows.

Remark 6.12. A non-crossing partition π = {B_{1}, B_{2}, …, B_{s}} of the set [i] is a partition such that, if 1 ≤ h < l < k < m ≤ i, with h, k ∈ B_{j}, and l, m ∈ B_{j′}, then j = j′. Let N C_{i} denote the set of all non-crossing partitions of [i]. Since N C_{i} ⊆ Π_{i}, we consider the induced subposet (N C_{i}, ≤) of (Π_{i}, ≤). If π ∈ N C_{i}, then the decomposition of the intervals [0_{i}, π] = {τ ∈ N C_{i} | τ ≤ π} is the same as the one given in Example 4.9, while the intervals [π, 1_{i}] = {τ ∈ N C_{i} | π ≤ τ} have a quite different decomposition, which does not depend on ℓ(π) (see [89]). The convolution ∗ defined on the multiplicative functions is given by

(𝔤 ∗ 𝔣)(σ, π) := ∑_{{σ≤τ≤π \atop τ∈N C_{i}}} 𝔤(σ, τ) 𝔣(τ, π).

Hence, if 𝔥 = 𝔤 ∗ 𝔣, then h_{i} = ∑_{τ∈N C_{i}} g_{τ} 𝔣(τ, 1_{i}). Let us observe that the generating function of f(𝔏ᾱ, z) is the Fourier transform (F𝔣)(z) of multiplicative functions [57]. Indeed, 𝔏ω̄ ≡ 𝔏γ̄ + 𝔏ᾱ corresponds to [F(𝔤 ∗ 𝔣)](z) = (F𝔤)(z)(F𝔣)(z) for all unital multiplicative functions 𝔣 and 𝔤 on the lattice of non-crossing partitions.

Remark 6.13. The umbra having all free cumulants equal to 1 has a special meaning. Indeed, the umbra ϑ̄ ≡ ū.(−1.ū)_{D}^{∗} is the unique umbra satisfying 𝔎̄ϑ̄ ≡ ū. Since ϑ^{i} ≃ C_{i}, where C_{i} denotes the i-th Catalan number, the umbra ϑ is called Catalan umbra. The parallelism with the Bell umbra relies on the fact that C_{i} = |N C_{i}|, so the moments of ϑ give the numbers of non-crossing partitions [90].

Example 6.14 (Wigner semicircle distribution). In free probability, the Wigner semicircle r.v.^{24} is analogous to the Gaussian r.v. in classical probability theory. Indeed, its free cumulants of degree higher than 2 are zero as for classical cumulants of a Gaussian r.v. By using Equivalence (4.6), it is straightforward to prove that the Wigner semicircle r.v. is represented by the umbra ϑ̄.β.η̄, where ϑ is the Catalan umbra and η is the umbra introduced in Example 3.2.

Theorem 6.15 (Abel parametrization). [32] If 𝔎̄ᾱ is the free cumulant umbra of α, then

(6.12)  ᾱ^{i} ≃ 𝔎̄ᾱ(𝔎̄ᾱ + i.𝔎̄ᾱ)^{i−1} and 𝔎̄i_{ᾱ} ≃ ᾱ(ᾱ − i.ᾱ)^{i−1}.

In terms of generalized Abel polynomials, Equivalences (6.12) satisfy

ᾱ^{i} ≃ 𝔞_{i}^{(i)}(𝔎̄ᾱ, 𝔎̄ᾱ) and 𝔎̄i_{ᾱ} ≃ 𝔞_{i}^{(−i)}(ᾱ, ᾱ).

Example 6.16. Complete Bell polynomials a_{i} = 𝒴_{i}(c_{1}, c_{2}, …, c_{i}) given in (4.4) return moments in terms of classical cumulants. For free cumulants, the same role is played by the volume polynomials introduced by Pitman and Stanley [67]. Let us recall that the i-th volume polynomial V_{i}(x) is the following homogeneous polynomial of degree i:

24The Wigner semicircle r.v. has probability density function given by a semicircle of radius R centered at (0, 0) and suitably normalized, that is, f(x) = 2\sqrt{R^{2} − x^{2}} /(πR^{2}) for x ∈ [−R, R] and f(x) = 0 otherwise.

(6.13)  V_{i}(x) = \frac{1}{i!} ∑_{p∈ park(i)}x_{p},

where x_{p} = x_{p_{1}} ⋯ x_{p_{i}} and park(i) is the set of all parking functions of length i.^{25} The volume polynomial (6.13) in the uncorrelated umbrae α_{1}, …, α_{i} similar to α satisfy i! V_{i}(α_{1}, …, α_{i}) ≃ α(α + i.α)^{i−1}, see [63]. In particular, if {r_{i}} denotes the sequence of free cumulants, then a_{i} = i! V_{i}(r_{1}, …, r_{i}), in parallel with (4.4).

6.2.2. Abel-type cumulants. A more general class of cumulants can be defined by using the generalized Abel polynomials. To the best of our knowledge, a previous attempt to give a unifying approach to cumulants was given in [2], but the Boolean case does not seem to fit in. Set 𝔞_{i}^{(−m)}(α) = 𝔞_{i}^{(−m)}(α, α).

Definition 6.17. An Abel-type cumulant c_{i,m}(α) of α is

c_{i,m}(α) ≃ 𝔞_{i}^{(−m)}(α) ≃ α(α − m.α)^{i−1}

if i and m are positive integers.

Let us consider the infinite matrix

(6.14)  C_{α} = ( \begin{matrix}𝔞_{1}^{(−1)}(α) & 𝔞_{1}^{(−2)}(α) & … \\ 𝔞_{2}^{(−1)}(α) & 𝔞_{2}^{(−2)}(α) & … \\ ⋮ & ⋮ \\ 𝔞_{n}^{(−1)}(α) & 𝔞_{n}^{(−2)}(α) & … \\ ⋮ & ⋮ \end{matrix} ) .

Classical cumulants can be recovered from the first column, and Boolean cumulants from the second column. The sequence corresponding to the main diagonal of the matrix (6.14) picks out free cumulants.

Theorem 6.18 (Homogeneity property).

𝔞_{i}^{(−m)}(cα) ≃ c^{i}𝔞_{i}^{(−m)}(α), c ∈ ℝ.

To prove the semi-invariance property, a suitable normalization of 𝔞_{i}^{(−m)}(α) is necessary. For example, for Boolean and free cumulants, the normalization coefficient is equal to {i!}, umbrally represented by the Boolean unity ū. For classical cumulants, no normalization is needed.

Each sequence of cumulants in (6.14) — chosen along columns or along diagonals — linearizes a suitable convolution of umbrae (i.e., of moments). More precisely, it is possible to define the convolution

25A parking function of length i is a sequence p = (p_{1}, …, p_{i}) of i positive integers, whose non-decreasing arrangement p^{↑} = (p_{j_{1}}, …, p_{j_{i}}) satisfies p_{j_{k}} ≤ k.

α +^{(m)} γ induced by an integer m as

(6.15)  𝔞_{i}^{(−m)}(α +^{(m)} γ) ≃ 𝔞_{i}^{(−m)}(α) + 𝔞_{i}^{(−m)}(γ).

These convolutions are commutative. It is possible to prove the existence of the convolution (6.15) induced by any integer m for any pair α, γ of umbrae [32]. It will be interesting to investigate if further sequences of cumulants known in the literature can be detected in the matrix (6.14), as for example the binary cumulants [95, 99], the tree-cumulants [108] and, more generally, the L-cumulants [109]; or if there are more sequences of cumulants that can be extracted from the matrix (6.14).

6.2.3. An algorithm to compute cumulants from moments and viceversa. A procedure to compute free cumulants from moments has been given in [11] by using nested derivatives of (6.11). The symbolic parametrization of cumulants via Abel polynomials suggests a completely different approach to perform such a computation [28]. The main advantage is that the resulting algorithm does not only allow us to compute free cumulants, but also classical and Boolean cumulants.

The algorithm relies on an efficient expansion of the umbral polynomial γ(γ + ρ.γ)^{i−1} for all positive integers i, with ρ, γ ∈ 𝒜 , that is,

(6.16)  𝔼 [γ(γ + ρ.γ)^{i−1}] = ∑_{λ⊢i} 𝔼[(ρ)_{ν_{λ}−1}] d_{λ} g_{λ},

where {g_{i}} is umbrally represented by the umbra γ. To compute the right-hand side of (6.16), we need the factorial moments of ρ. Recall that, if we just know the moments {𝔼[ρ^{i}]}, the factorial moments can be computed by using the well-known change of basis

𝔼[(ρ)_{i}] = ∑_{k=1}^{i} s(i, k)𝔼[ρ^{k}],

where {s(i, k)} are the Stirling numbers of the first kind.

Expansion (6.16) can be fruitfully used to compute cumulants and moments as follows:

i) for classical cumulants in terms of moments: we use κi_{α} ≃ α(α − 1.α)^{i−1} in (4.2); then in Equation (6.16) choose ρ = −1.u, with 𝔼[(−1.u)_{i}] = (−1)_{i} = (−1)^{i}i!, and γ = α;

ii) for moments in terms of classical cumulants: we use α^{i} ≃ κ_{α}(κ_{α} + α)^{i−1} in (4.2); then in Equation (6.16) choose ρ = β, with 𝔼[(β)_{i}] = 1, and γ = κ_{α};

iii) for Boolean cumulants in terms of moments: we use η̄_{α}^{i} ≃ ᾱ (ᾱ − 2.ᾱ)^{i−1} in (6.7); then in Equation (6.16) choose ρ = −2.u, with 𝔼[(−2.u)_{i}] = (−1)^{i}(i + 1)!, and γ = ᾱ;

iv) for moments in terms of Boolean cumulants: we use ᾱ^{i} ≃ η̄_{α} (η̄_{α} + 2.ᾱ)^{i−1} in (6.7); then in Equation (6.16) choose ρ = 2.ū.β, with 𝔼[(2.ū.β)_{i}] = (i + 1)!, and γ = η̄_{α};

v) for free cumulants in terms of moments: we use 𝔎̄i_{ᾱ} ≃ ᾱ(ᾱ − i.ᾱ)^{i−1} in (6.12); then in Equation (6.16) choose ρ = −i.u, with 𝔼[(−i.u)_{j}] = (−i)_{j}, and γ = ᾱ;

vi) for moments in terms of free cumulants: we use ᾱ^{i} ≃ 𝔎̄ᾱ(𝔎̄ᾱ + i.𝔎̄ᾱ)^{i−1} in (6.12); then in Equation (6.16) choose ρ = i.u, with 𝔼[(i.u)_{j}] = (i)_{j}, and γ = 𝔎̄ᾱ.

The symbolic algorithm allows us to compute also Abel-type cumulants by choosing ρ = −m.u, with 𝔼[(−m.u)_{i}] = (−m)_{i}, and γ = α in (6.16).

7. Multivariate Faà di Bruno formula

The multivariate Faà di Bruno formula is a chain rule to compute higher derivatives of multivariable functions. This formula has been recently addressed by the following two approaches: combinatorial methods are used by Constantine and Savits [13], and (only in the bivariate case) by Noschese and Ricci [60]; analytical methods involving Taylor series are proposed by Leipnik and Pearce [52]. We refer to this last paper for a detailed list of references on this subject and for a detailed account of its applications. For statistical applications, a good account is the paper of Savits [81]. A comprehensive survey of univariate and multivariate series approximation in statistics is given in [49].

Computing the multivariate Faà di Bruno formula by means of a symbolic software can be done by recursively applying the chain rule. Despite its conceptual simplicity, applications of the chain rule become impractical already for small values, because the number of additive terms becomes enormous and their computation unwieldy. Moreover, the output is often untidy and further manipulations are required to simplify the result. Consequently, a “compressed” version of the multivariate Faà di Bruno formula is more attractive, as the number of variables or the order of derivatives increase. By using the symbolic method, the “compressed” version of the multivariate Faà di Bruno formula is nothing else but a suitable generalization of the well-known multinomial theorem:

(7.1)  (x_{1} + ⋯ + x_{k})^{i} = ∑_{j_{1}+···+j_{k}=i} \binom{i}{j_{1}, …, j_{k}}x_{1}^{j_{1}} ⋯ x_{k}^{j_{k}},

where the indeterminates are replaced by auxiliary umbrae. Suitable choices of these auxiliary umbrae give rise to an efficient computation of the following compositions: univariate with multivariate, multivariate with univariate, multivariate with the same multivariate, multivariate with different multivariates in an arbitrary number of components [25].

Various attempts to construct a multivariate generalization of the umbral calculus can be found in the literature, see, for example, [61] and [10].

Here, the starting point is the notion of multi-index. If, in the univariate case, the main device is to replace a_{n} by α^{n} via the linear functional 𝔼, in the multivariate case, the main device is to replace sequences like {g_{i_{1},...,i_{k}}}, where i = (i_{1}, …, i_{k}) ∈ ℕ_{0}^{k} is a multi-index, by a product of powers μ_{1}^{i_{1}} ⋯ μ_{n}^{i_{k}}, where {μ_{1}, …, μ_{k}} are umbral monomials in ℝ[𝒜] with not necessarily disjoint support. The i-th power of μ = (μ_{1}, …, μ_{k}) is defined by μ^{i} = μ_{1}^{i_{1}} ⋯ μ_{k}^{i_{k}}.

Products of powers of umbral monomials have already been employed in Theorem 2.11 and in Definition 2.12, with umbrae indexed by multisets. Indeed there is a connection between the combinatorics of multisets and the Faà di Bruno formula, as highlighted in [45]. But, in facing the problem to set up an algorithm to efficiently compute the multivariate Faà di Bruno formula, multi-index notations have simplified formulae and procedures.

Definition 7.1. A sequence {g_{i}}_{i∈ℕ_{0}^{k}} ∈ ℝ, with g_{i} = g_{i_{1},...,i_{k}} and g_{0} = 1, is umbrally represented by the k-tuple μ if

𝔼[μ^{i}] = g_{i}, for all i ∈ ℕ_{0}^{k}.

The elements {g_{i}}_{i∈ℕ_{0}^{k}} are called multivariate moments of μ.

Definition 7.2. Two k-tuples ν and μ of umbral monomials are said to be similar, if they represent the same sequence of multivariate moments. In symbols: ν ≡ μ.

Definition 7.3. Two k-tuples ν and μ of umbral monomials are said to be uncorrelated if, for all i, j ∈ ℕ_{0}^{k}, we have E[ν^{i}μ^{j}] = E[ν^{i}]E[μ^{j}].

If the sequence {g_{i}}_{i∈ℕ_{0}^{k}} is umbrally represented by the k-tuple μ, then the multivariate formal power series

f(μ, t) = 1 + ∑_{i≥1} g_{i} \frac{t^{i}}{i!}

is the multivariate moment generating function of μ.

Definition 7.4. A random vector X is said to be represented by a ktuple μ of umbral monomials, if its sequence of multivariate moments {g_{i}} is umbrally represented by μ.

Thus, a summation of m uncorrelated k-tuples similar to μ represents a summation of m i.i.d. random vectors and can be denoted by m.μ. In order to compute its multivariate moments, we need the notion of composition of a multi-index and the notion of partition of a multi-index.

Definition 7.5 (Composition of a multi-index). A composition λ of a multi-index i. In symbols λ ⊧ i, is a matrix λ = (λ_{ij}) of nonnegative integers, with no zero columns, satisfying λ_{r1} + λ_{r2} + ⋯ = i_{r} for r = 1, 2, …, k.

The number of columns of λ is called the length of λ and is denoted by l(λ).

Definition 7.6 (Partition of a multi-index). A partition of a multi-index i is a composition λ, whose columns are in lexicographic order. In symbols λ ⊢ i.

As for integer partitions, the notation λ = (λ_{1}^{r_{1}}, λ_{2}^{r_{2}}, …) means that in the matrix λ there are r_{1} columns equal to λ_{1}, r_{2} columns equal to λ_{2}, and so on, with λ_{1} < λ_{2} < ⋯ in lexicographic order. The integer r_{i} is the multiplicity of λ_{i}. We set 𝔪(λ) = (r_{1}, r_{2}, …).

Proposition 7.7. [23] If i ∈ ℕ_{0}^{k}, then

𝔼[(m.μ)^{i}] = ∑_{λ⊢i} \frac{i!}{𝔪(λ)! λ!}(m)_{l(λ)}E[μ^{λ^{1}}]^{r_{1}}E[μ^{λ^{2}}]^{r_{2}} ⋯,

where the sum is over all partitions λ of the multi-index i.

Multivariate moments of the auxiliary umbra (α.β.μ)^{i} can be computed from Proposition 7.7:

(7.2)  𝔼[(α.β.μ)^{i}] = ∑_{λ⊢i} \frac{i!}{m(λ)! λ!} a_{l(λ)}E[μ^{λ^{1}}]^{r_{1}}E[μ^{λ^{2}}]^{r_{2}} · · · .

To obtain Equation (7.2), the arguments are the same as the ones employed to get Equation (4.6). The umbra α.β.μ represents a multivariate randomized compound Poisson r.v., that is, a random sum S_{N} = X_{1} + ⋯ + X_{N} of i.i.d. random vectors {X_{i}} indexed by a Poisson r.v. N. The right-hand side of Equivalence (7.2) gives moments of S_{N} if we replace the i-th multivariate moment of μ by the i-th multivariate moment of {X_{i}} and the moments of α by those of N. This result is the same as Theorem 4.1 of [13], obtained with a different proof and different methods.

A way to compute partitions of a multi-index i is by employing subdivisions of the multiset M, when i is its vector of multiplicities.

Example 7.8. The multiset M = {μ_{1}, μ_{1}, μ_{2}} corresponds to the multi-index (2, 1), and the subdivisions

{{μ_{1}, μ_{1}, μ_{2}}}, {{μ_{1}, μ_{1}}, {μ_{2}}}, {{μ_{1}, μ_{2}}, {μ_{1}}}, {{μ_{1}}, {μ_{1}}, {μ_{2}}} correspond to the multi-index partitions

\binom{2}{1}, \binom{2, 0}{0, 1}, \binom{1, 1}{1, 0}, \binom{1, 1, 0}{0, 0, 1}.

By using the auxiliary umbra α.β.μ, we can summarize the different types of composition as follows:

a) (univariate composed with multivariate) The composition of a univariate formal power series f(α, z) with a multivariate formal power series f(μ, z) is defined by

f(α.β.μ, z) = f[α, f(μ, z) − 1].

Equivalence (7.2) gives its i-th coefficient.

b) (multivariate composed with univariate) The composition of a multivariate formal power series f(μ, z) with a univariate formal power series f(γ, z) is defined by

f[(μ_{1} + ⋯ + μ_{k}).β.γ, z] = f[μ, (f(γ, z) − 1, …, f(γ, z) − 1)].

The i-th coefficient of f[(μ_{1} + ⋯ + μ_{k}).β.γ, z] can be computed by replacing γ by (μ_{1} + ⋯ + μ_{k}) in (4.6) and α by γ.

c) (multivariate composed with multivariate) The composition of a multivariate formal power series f(μ, z) with a multivariate formal power series f(ν, z) is defined by

f[(μ_{1} + ⋯ + μ_{k}).β.ν, z] = f[μ, (f(ν, z) − 1, …, f(ν, z) − 1)].

The i-th coefficient of f[(μ_{1} + ⋯ + μ_{k}).β.ν, z] can be computed by evaluating Equivalence (7.2) after having replaced μ by ν and the umbra α by (μ_{1} + ⋯ + μ_{k}).

d) (multivariate composed with multivariate) The composition of a multivariate formal power series f(μ, z) with a multivariate formal power series f(ν, z_{(m)}), with a different number of indeterminates, is defined by

f[(μ_{1} + ⋯ + μ_{k}).β.ν, z_{(m)}] = f[μ, (f(ν, z_{(m)}) − 1, …, f(ν, z_{(m)}) − 1)].

e) (multivariate composed with different multivariates) The composition of a multivariate formal power series f(μ, z) with different multivariate formal power series f(ν_{i}, z_{(m)}), with a different number of indeterminates, is defined by

(7.3) f[μ_{1}.β.ν_{1} + ⋯ + μ_{k}.β.ν_{k}, z_{(m)}] = f[μ, (f(ν_{1}, z_{(m)}) − 1, …, f(ν_{k}, z_{(m)}) − 1)].

The i-th coefficient of the compositions d) and e) can be obtained by using generalized Bell polynomials.

Definition 7.9. The generalized Bell polynomial of order i is

B_{i}^{(ν_{1},...,ν_{k})}(x_{1}, …, x_{k}) = (x_{1}.β.ν_{1} + ⋯ + x_{k}.β.ν_{k})^{i}.

The formal power series whose coefficients are generalized Bell polynomials is exp {∑_{i=1}^{k} x_{i}[f(ν_{i}, z_{(m)}) − 1]} . These polynomials give the multivariate Faà di Bruno formula for all cases listed before:

a) (univariate composed with multivariate) B_{i}^{ν}(α) = (α.β.μ)^{i};

b) (multivariate composed with univariate)

B_{i}^{γ}(μ_{1} + ⋯ + μ_{k}) = [(μ_{1} + ⋯ + μ_{k}).β.γ]^{i};

c);d) (multivariate composed with multivariate)

B_{i}^{ν}(μ_{1} + ⋯ + μ_{k}) = [(μ_{1} + ⋯ + μ_{k}).β.ν]^{i};

e) (multivariate composed with different multivariates)

B_{i}^{(ν_{1},...,ν_{k})}(μ_{1}, …, μ_{k}) ≃ (μ_{1}.β.ν_{1} + ⋯ + μ_{k}.β.ν_{k})^{i}.

A MAPLE algorithm for computing the evaluation of B_{i}^{(ν_{1},...,ν_{k})}(μ_{1}, …, μ_{k}) is given in [25]. The main idea is to apply the multivariate version of (7.1) to get

(7.4) 𝔼 [B_{i}^{(ν_{1},...,ν_{k})}(μ_{1}, …, μ_{k})]

= ∑ ; (i_{1},...,i_{k}):∑_{j=1}^{k} i_{j}=i \binom{k}{i_{1}, …, i_{k}}𝔼 [(μ_{1}.β.ν_{1})^{i^{1}} ⋯ (μ_{k}.β.ν_{k})^{i^{k}}] ,

to expand powers (μ_{i}.β.ν_{i})^{k_{i}}, to multiply the resulting expansions by (μ_{1}.β.ν_{1})^{i^{1}} ⋯ (μ_{d}.β.ν_{d})^{i^{d}}, and then to apply evaluation lowering powers to indexes.

7.1. Multivariate cumulants. Given the sequence {g_{i}}_{i∈ℕ_{0}^{n}}, its multivariate cumulants are the coefficients of the formal power series in the left-hand side of the following equation

1 + ∑_{i≥1} g_{i} \frac{t^{i}}{i!} = exp (∑_{j≥1} c_{j} \frac{t^{j}}{i!} ) .

Proposition 7.10. The auxiliary umbra χ.μ represents the multivariate cumulants of μ.

Proof. The result follows from a), by choosing as umbra α the umbra χ.χ, and observing that χ.μ ≡ (χ.χ).β.μ and

f(χ.χ, t) = 1 + log(t + 1).  □

The multivariate analog of the well-known semi-invariance property of cumulants is χ.(μ+ν) ≡ χ.μ ˙+ χ.ν, where μ and ν are uncorrelated n-tuples of umbral monomials, and the disjoint sum μ ˙+ ν is defined by the multivariate moments 𝔼[(μ ˙+ ν)^{i}] = 𝔼[μ^{i}] + 𝔼[ν^{i}]. Thus, multivariate cumulants linearize convolutions of uncorrelated n-tuples of umbral monomials. By using Equation (7.2) with α replaced by χ.χ, that is, by expanding the i multivariate moments of (χ.χ).β.μ ≡ χ.μ, we get multivariate cumulants in terms of multivariate moments. Still using Equation (7.2) with α replaced by u and μ replaced by χ.μ, that is, expanding the i-th multivariate moment of u.β.(χ.μ) ≡ μ, we get multivariate moments in terms of multivariate cumulants. The same equations are given in [54] using tensor notations. Umbrae indexed by multisets are employed in [21].

7.2. Multivariate Lévy processes. Let X = {X_{t}} be a Lévy pro-

cess on ℝ^{k}, with X_{t} = (X_{1}^{(t)}, …, X_{k}^{(t)}) , where {X_{j}^{(t)}}_{j=1}^{k} are univari-

ate r.v.’s all defined on the same probability space for all t ≥ 0.

In Section 3, we have shown that every Lévy process X = {X_{t}} on ℝ is umbrally represented by the family of auxiliary umbrae {t.α}_{t≥0} defined by 𝔼 [(t.α)^{i}] = E [X_{t}^{i}] , for all non-negative integers i. This result can be generalized to the multivariate case as follows.

Definition 7.11. A Lévy process X = {X_{t}}_{t≥0} on ℝ^{k} is umbrally represented by {t.μ}_{t≥0}, if μ is a k-tuple of umbral monomials representing the k-tuple of random vectors X_{1} = (X_{1}^{(1)}, …, X_{k}^{(1)}) with the property that

E[μ^{i}] = E [{X_{1}^{(1)}}^{i_{1}} ⋯ {X_{k}^{(1)}}^{i_{k}}] ,

for all i = (i_{1}, …, i_{k}) ∈ N_{0}^{k}.

The following theorem gives the multivariate analog of the Lévy– Khintchine formula (4.9).

Theorem 7.12 (Multivariate Lévy–Khintchine formula). [80] If {X_{t}} is a multivariate Lévy process, then E [e^{X^{t}z′}] = [ϕ_{X_{1}}(z)]^{t} with^{26}

(7.5) ϕ_{X_{1}}(z) = exp \{[\frac{1}{2}zΣz′ + mz′

+ ∫_{ℝ^{k}} (e^{zx′} − 1 − z x′1_{{|x|≤1}}(x)) ν(dx)]} ,

where m ∈ ℝ^{k} and Σ is a symmetric, positive-definite k × k matrix.

The triple (m, Σ, ν) is called Lévy triple, Σ is called covariance matrix, and ν is the Lévy measure on ℝ^{k}. As in the univariate case, if we set m_{2} = ∫_{ℝ^{k}} x 1_{{|x|>1}}(x) ν(dx) and m_{1} = m + m_{2} in Equation (7.5), then

ϕ_{X_{1}}(z) = exp \{[\frac{1}{2}zΣz′ + m_{1}z′ + ∫_{ℝ^{k}}(e^{zx′} − 1 − zx′)ν(dx)]\} . (7.6)

Due to the above equation, a multivariate Lévy process is represented by a family of auxiliary umbrae, as given in the following theorem.

26Here, the notation z′ denotes the transpose of the row vector z.

Theorem 7.13. A multivariate Lévy process {X_{t}}_{t≥0} is represented by the family of auxiliary umbrae

(7.7)  {t.β. (ηC′ ˙+ χ.m_{1} ˙+ γ)}_{t≥0} ,

where η = (η_{1}, …, η_{k}) is a k-tuple of uncorrelated umbrae satisfying f(η, z) = 1 + zz′/2, the matrix C is the square root of the covariance matrix Σ in Equation (7.6), and γ is a k-tuple of umbral monomials satisfying

f(γ, z) = 1 + ∫_{ℝ^{k}} (e^{zx′} − 1 − zx′) ν(dx).

As observed previously, every auxiliary umbra, such as t.β.ν, is the symbolic version of a multivariate compound Poisson process of parameter t. Therefore, for t fixed, the symbolic representation (7.7) of a multivariate Lévy process is a multivariate compound Poisson r.v. of parameter t, and the k-tuple (ηC′ ˙+ χ.m_{1} ˙+ γ) umbrally represents the multivariate cumulants of X_{1}.

The auxiliary umbra χ.m_{1} does not have a probabilistic counterpart. If m is not equal to the zero vector, this parallels the well-known difficulty to interpret the Lévy measure as a probability measure in Subsection 4.2.

7.3. Multivariate Hermite polynomials. In Equation (7.6), assume the Lévy measure ν to be zero. Then the moment generating function of {X_{t}}_{t≥0} is E [e^{X^{t}z′}] = [ϕ_{X_{1}}(z)]^{t} with

(7.8)  ϕ_{X_{1}}(z) = exp \{(\frac{1}{2}zΣz′ + m_{1}z′)\} .

Theorem 7.14 (Wiener process). [80] The stochastic process {X_{t}}_{t≥0} having moment generating function (7.8) is a multivariate Wiener process, that is, X_{t} = m_{1} t + CB_{t}, where C is a k × k matrix whose determinant is non-zero satisfying Σ = CC′, and {B_{t}}_{t≥0} is the multivariate standard Brownian motion in ℝ^{k}.

When m_{1} = 0, the Wiener process reduces to a multivariate nonstandard Brownian motion. From Theorem 7.14, we see that the Wiener process X_{t} is a special multivariate Lévy process, and the following result is a corollary of Theorem 7.13.

Corollary 7.15. The Wiener process {X_{t}}_{t≥0} is represented by the family of auxiliary umbrae

{t.β.(ηC′ ˙+ χ.m_{1})}_{t≥0}.

There is a special family of multivariate polynomials which is closely related to the multivariate non-standard Brownian motion: the (multivariate) Hermite polynomials. The i-th Hermite polynomial H_{i}(x, Σ) is defined by

H_{i}(x, Σ) = (−1)^{|i|} \frac{D_{x}^{(i)}ϕ(x; 0, Σ)}{ϕ(x; 0, Σ)},

where ϕ(x; 0, Σ) denotes the multivariate Gaussian density with mean 0 and covariance matrix Σ of full rank k. The polynomials H_{i}(xΣ^{−1}, Σ^{−1}) are orthogonal with respect to ϕ(x; 0, Σ), where Σ^{−1} is the inverse matrix of Σ, and is denoted by H̃_{i}(x, Σ).

Theorem 7.16. [107] If Z ∼ N(0, Σ) and Y ∼ N(0, Σ^{−1}) are multivariate Gaussian vectors with mean 0 and covariance matrix Σ and Σ^{−1}, respectively, then, for all i ∈ ℕ_{0}^{k}, we have

(7.9) H_{i}(x, Σ) = E[(xΣ^{−1} + 𝔦Y )^{i}], and H̃_{i}(x, Σ) = E[(x + 𝔦Z)^{i}],

where 𝔦 is the imaginary unit.

In the literature, Equations (7.9) are called moment representation of the Hermite polynomials. Quite recently, many authors have obtained moment representations for various families of polynomials, especially in the univariate case, see references in [31]. The multivariate symbolic method allows us to widen the field of applicability of this moment representation. For the multivariate Hermite polynomials, we have

(7.10) H_{i}(x, Σ) = 𝔼 [(−1.β.ν + xΣ^{−1})^{i}] ,

and H̃_{i}(x, Σ) = 𝔼 [(−1.β.μ + x)^{i}] ,

with ν a k-tuple of umbral monomials satisfying f(ν, t) = 1 + \frac{1}{2}tΣ^{−1}t′ and μ a k-tuple of umbral monomials satisfying f(μ, t) = 1+ \frac{1}{2}tΣt′. In Section 5, we have shown that Appell polynomials are umbrally represented by the polynomial umbra x.u + α. It is well known that univariate Hermite polynomials are Appell polynomials. Equations (7.10) then show that also multivariate Hermite polynomials are of Appell type. Let us stress that the umbra −1.β enables us to obtain a simple expression for multivariate Hermite polynomials, without the employment of the imaginary unit as done in Equations (7.9). Moreover, since the moments of Z ∼ N(0, Σ) are umbrally represented by the umbra β.μ, the k-tuple μ ≡ χ.(β.μ) in H̃_{i}(x, Σ) umbrally represents the multivariate cumulants of Z.

Proposition 7.17. 𝔼 {H̃_{i}[t.β(ηC′), Σ]} = 0 for all t ≥ 0.

Proof. The result follows by observing that the k-tuple μ in (7.10) satisfies μ ≡ ηC′ and

H̃_{i}[t.β(ηC′), Σ] = 𝔼 {[−t.β.μ + t.β(ηC′)]^{i}} .  □

For the family of polynomials {H̃_{i}^{(t)}(x, Σ)}_{t≥0}, Proposition 7.17 parallels the main property of time-space harmonic polynomials introduced in Section 3.1. Since, in the literature, to the best of our knowledge, a theory of multivariate time-space harmonic polynomials has not yet been given, we believe that the setting introduced here could be a fruitful way to build this theory with respect to multivariate Lévy processes.

More properties on multivariate Hermite polynomials can be found in [25]. We end this subsection by recalling that multivariate Hermite polynomials are special generalized Bell polynomials.

Proposition 7.18. If ν_{x} is a k-tuple of umbral polynomials satisfying f(ν_{x}, t) = 1 + f(ν, t + x) − f(ν, t), with f(ν, t) = 1 + \frac{1}{2}tΣ^{−1}t′, then

H_{i}(x, Σ) = (−1)^{|i|}𝔼 [(−1.β.ν_{x})^{i}] .

Finally we remark that the efficient computation of the multivariate Hermite polynomials as generalized Bell polynomials can help us in constructing an efficient algorithm for a multivariate Edgeworth approximation of multivariate density functions [49].

7.4. Multivariate Bernoulli polynomials. By the symbolic method, it is possible to introduce multivariate Bernoulli polynomials as powers of polynomials whose coefficients involve multivariate Lévy processes. According to [31], this means to give a moment representation for these polynomials.

Definition 7.19. The multivariate Bernoulli umbra ι is the k-tuple (ι, …, ι), with ι the Bernoulli umbra.

Definition 7.20. If v ∈ ℕ_{0}^{k}, the multivariate Bernoulli polynomial of order v is B_{v}^{(t)}(x) = 𝔼[(x + t.ι)^{v}], where ι is the multivariate Bernoulli umbra and t ≥ 0.

As corollary, the multivariate Bernoulli polynomials are

B_{v}^{(t)}(x) = ∑_{j≤v}(\begin{matrix}v \\ j \end{matrix}) x^{j} 𝔼[(t.ι)^{v−j}].

The symbolic representation of the coefficients of B_{v}^{(t)}(x) by 𝔼[(t.ι)^{v}] simplifies the calculus and is computationally efficient, see [31].

Proposition 7.21. [31] B_{v}^{(t)}(−t.ι) = B_{v}^{(t)}[t.(−1.ι)] = 0.

Let us observe that, when in B_{v}^{(t)}(x) we replace x by the auxiliary umbra −t.μ, this means to replace x^{j} by 𝔼[(−t.μ)^{j}]. The auxiliary umbra −t.ι is the symbolic version of a multivariate Lévy process. Indeed, −1.ι is the umbral counterpart of a k-tuple identically distributed with (U, …, U), where U is a uniform r.v. on the interval (0, 1), and corresponds to the random vector X_{1} = (X_{1}^{(1)}, …, X_{k}^{(1)}) given in Definition 7.11. By generalizing the definition of conditional evaluation given in Section 3 to the multivariate case, the multivariate Bernoulli polynomials should become time-space harmonic with respect to the multivariate Lévy processes {−t.ι}_{t≥0}. Indeed, Proposition 7.21 shows that these polynomials share one of the main properties of time-space harmonic polynomials: when the vector of indeterminates is replaced by the corresponding Lévy process, their overall mean is zero. Similar considerations apply for multivariate Euler polynomials, see [31] for details and further connections between these two families of polynomials.

8. Appendix 1: Tables of computation times

In Table 3, computation times are given for computing the mean of statistics such as the second expression in (2.13). The first column refers to the Mathematica routine proposed in [105]. The second column refers to a Maple routine relying on the symbolic method of moments.

Table 3. Comparison of computation times.

[1^{i}2^{j}3^{k} ⋯]  SIP  MAPLE
[5^{3} 9 10][1 2 3 4 5]  0.7 0.1
[5^{3} 8 9 10][1 2 3 4 5]  5.6 0.4
[6 7 8 9 10][1 2 3 4 5]  2.2 0.1
[6 7 8 9 10][1 2][3 4 5]  3.1 0.4
[6 7][8 9 10][1 2][3 4 5]  4.7 1.3
[5 6 7 8 9 10][1 2 3 4 5]  16.7 0.3
[5 6 7 8 9 10][1 2 3 4 5 6] 348.7 1.5
[6 7 8 9 10][6 7][3 4 5][1 2] 125.6 16.4

Table 4 lists computation times (in seconds) observed by using the polynomial expansion in (6.16) and Bryc’s procedure [11], both implemented in MAPLE, release 7, for the conversion from free cumulants to moments.^{27}

Table 4. Comparison of computation times needed to compute free cumulants in terms of moments. Tasks performed on Intel (R) Pentium (R), CPU 3.00 GHz, 512 MB RAM.

i MAPLE (umbral) MAPLE (Bryc)
15  0.015  0.016
18  0.031  0.062
21  0.078  0.141
24  0.172  0.266
27  0.375  0.703

27The output is in the same form as the one given by Bryc’s procedure.

Table 5 shows computation times for three different software packages. The first one has been implemented in Mathematica and refers to the procedures in [1], here called AS algorithms, see also

http://fisher.utstat.toronto.edu/david/SCSI/chap.3.nb. The

second one is the package MathStatica from [72]. The third package, Fast algorithms, has been implemented in Maple 10.x by using the results of Section 6. The source code is given in [24]. Let us remark that, for all procedures, the results are in the same output form and have been performed by the authors on the same platform. To the best of our knowledge, there is no R implementation for k-statistics and polykays.

Table 5. Computation times in sec. for k-statistics and polykays. Missing computation times mean “greater than 20 hours”.

k_{t,..., l} AS Algorithms MathStatica  Fast algorithms
k_{5}  0.06  0.01  0.01
k_{7}  0.31  0.02  0.01
k_{9}  1.44  0.04  0.01
k_{11}  8.36  0.14  0.01
k_{14}  396.39  0.64  0.02
k_{16}  57982.40  2.03  0.08
k_{18}  -  6.90  0.16
k_{20}  -  25.15  0.33
k_{22}  -  81.70  0.80
k_{24}  -  359.40  1.62
k_{26}  -  1581.05  2.51
k_{28}  -  6505.45  4.83
k_{3,2}  0.06  0.02  0.01
k_{4,4}  0.67  0.06  0.02
k_{5,3}  0.69  0.08  0.02
k_{7,5}  34.23  0.79  0.11
k_{7,7}  435.67  2.52  0.26
k_{9,9}  -  27.41  2.26
k_{10,8}  -  30.24  2.98
k_{4,4,4}  34.17  0.64  0.08

Table 6 is the same as Table 5, but for multivariate k-statistics and multivariate polykays. Missing computation times in the second column are due to the fact that no procedures are devoted to multivariate polykays in MathStatica (release 1).

Table 6. Computation times in sec. for multivariate k-statistics and multivariate polykays. In the first column, missing computation times mean “greater than 20 hours”.

k_{t_{1}... t_{r}; l_{1}...l_{m}} AS Algorithms MathStatica  Fast-algorithms
k_{3 2}  0.20  0.02  0.01
k_{4 4}  18.20  0.14  0.02
k_{5 5}  269.19  0.56  0.08
k_{6 5}  1023.33  1.02  0.16
k_{6 6}  -  2.26  0.33
k_{7 6}  -  4.06  0.59
k_{7 7}  -  8.66  1.23
k_{8 6}  -  7.81  1.16
k_{8 7}  -  15.89  2.59
k_{8 8}  -  30.88  5.53
k_{3 3 3}  1211.05  0.92  0.44
k_{4 3 3}  -  2.09  0.34
k_{4 4 3}  -  4.98  1.02
k_{4 4 4}  -  13.97  2.78
k_{1 1; 1 1}  0.05  -  0.01
k_{2 1; 1 1}  0.20  -  0.01
k_{2 2; 1 1}  1.30  -  0.03
k_{2 2; 2 1}  6.50  -  0.06
k_{2 2; 2 2}  34.31  -  0.11
k_{2 1; 2 1; 2 1}  81.13  -  0.17
k_{2 2; 1 1; 1 1}  30.34  -  0.14
k_{2 2; 2 1; 1 1}  127.50  -  0.22
k_{2 2; 2 1; 2 1}  406.78  -  0.47
k_{2 2; 2 2; 1 1}  467.88  -  0.55
k_{2 2; 2 2; 2 1}  1402.55  -  1.14
k_{2 2; 2 2; 2 2}  3787.41  -  2.96

9. Appendix 2: Families of time-space harmonic polynomials

Hermite polynomials

The generalized Hermite polynomials {H_{k}^{(s^{2}t)}(x)} [71] are time-space harmonic polynomials with respect to the Wiener process with zero drift and variance s^{2}, whose umbral counterpart is the family of auxiliary umbrae {t.β.(sη)}_{t≥0}. Indeed,

H_{k}^{(s^{2}t)}(x) = 𝔼[(x − t.β.(sη))^{k}], k = 1, 2, …,

with s ∈ ℝ^{+} and η the umbra introduced in Example 3.2.

Poisson–Charlier polynomials

The Poisson–Charlier polynomials {C̃_{k}(x, λt)} [71] are time-space harmonic with respect to the Poisson process {N_{t}}_{t≥0} with intensity parameter λ > 0, whose umbral counterpart is the family of auxiliary umbrae {(tλ).β}_{t≥0}. Indeed,

C̃_{k}(x, λt) = ∑_{j=1}^{k} s(k, j) 𝔼[(x − (tλ).β)^{j}], k = 1, 2, …,

where {s(k, j)} are the Stirling numbers of the first kind.

Lévy–Sheffer systems

The Lévy–Sheffer system [83] is a sequence of polynomials {V_{k}(x, t)}_{t≥0} time-space harmonic with respect to the Lévy process umbrally represented by {−t.α.β.γ^{<−1>}}_{t≥0}, since

V_{k}(x, t) = ∑_{i=0}^{k} 𝔼[(x − t.α.β.γ^{<−1>})^{i}] B_{k,i}(g_{1}, …, g_{k−i+1}), k = 1, 2, …,

where B_{k,i} are partial Bell exponential polynomials and {g_{i}} is the sequence umbrally represented by γ.

Laguerre polynomials

The Laguerre polynomials {L_{k}^{t−k}(x)}_{t≥0} [71] are time-space harmonic with respect to the gamma process {G_{t}(1, 1)}_{t≥0}, with scale parameter 1 and shape parameter 1, whose umbral counterpart is the family of auxiliary umbrae {t.ū}_{t≥0}. Indeed,

(−1)^{k} k! L_{k}^{t−k}(x) = 𝔼[(x − t.ū)^{k}], k = 1, 2, …,

where ū is the Boolean unity.

Actuarial polynomials

The actuarial polynomials {g_{k}(x, λt)}_{t≥0} [71] are time-space harmonic with respect to the gamma process {G_{t}(λ, 1)}_{t≥0}, with scale parameter λ > 0 and shape parameter 1, whose umbral counterpart is the family of auxiliary umbrae {(λt).ū}_{t≥0}. Indeed,

g_{k}(x, λt) = ∑_{i=1}^{k} 𝔼[(x − (λt).ū)^{i}] B_{k,i}(m_{1}, …, m_{k−i−1}), k = 1, 2, …,

where B_{k,i} are partial Bell exponential polynomials and m_{i} = 𝔼[(χ.(−χ))^{i}], for i = 1, …, k.

Meixner polynomials

The Meixner polynomials of the first kind {M_{k}(x, t, p)}_{t≥0} [71] are time-space harmonic with respect to the Pascal process {P a(p, t)}_{t≥0}, whose umbral counterpart is the family of auxiliary umbrae {t.ū.d.β}_{t≥0}. Indeed,

(−1)^{k} (t)_{k} M_{k}(x, t, p) = ∑_{i=1}^{k} 𝔼[(x − t.ū.d.β)^{i}] B_{k,i}(m_{1}, …, m_{k−i−1}),

where B_{k,i} are partial Bell exponential polynomials and for i = 1, …, k

m_{i} = 𝔼 [(χ. (−1.χ + \frac{χ}{p}))^{i}] .

Bernoulli polynomials

The Bernoulli polynomials {B_{k}(x, n)} [71] are time-space harmonic with respect to the random walk {n.(−1.ι)}_{n≥0}. Indeed,

B_{k}(x, n) = 𝔼 [(x − n. (−1.ι))^{k}] , k = 1, 2, ….

The n-th increment M_{n} in Example 3.12 is represented by the umbra −1.ι corresponding to a uniform r.v. on the interval [0, 1].

Euler polynomials

The Euler polynomials {E_{k}(x, n)} [71] are time-space harmonic with respect to the random walk \{n. [\frac{1}{2} (−1.ξ + u)]\}_{n≥0}. Indeed,

E_{k}(x, n) = 𝔼 [(x − n. [\frac{1}{2} (−1.ξ + u)])^{k}] , k = 1, 2, ….

The n-th increment M_{n} in Example 3.12 is represented by the umbra \frac{1}{2} (−1.ξ + u) corresponding to a Bernoulli r.v. with parameter 1/2.

Krawtchouk polynomials

The Krawtchouk polynomials {K_{k}(x, p, n)} [71] are time-space harmonic with respect to the random walk {n.(−1.υ)}_{n≥0}. Indeed,

\frac{n!}{(n − k)!} K_{k}(x, p, n) = ∑_{i=1}^{k} 𝔼[(x + n.υ)^{i}] B_{k,i}(m_{1}, …, m_{k−i+1}),

where B_{k,i} are partial Bell exponential polynomials,

m_{i} = 𝔼 [(χ. (−1.χ − \frac{χ}{d}))^{i}]

for i = 1, …, k, and p/q = d, p + q = 1. The n-th increment M_{n} in Example 3.12 is represented by the umbra −1.υ corresponding to a Bernoulli r.v. with parameter p.

Pseudo–Narumi polynomials

The pseudo–Narumi polynomials {N_{k}(x, an)}, with a ∈ ℕ [71], are time-space harmonic with respect to the random walk {(an).(−1.ι)}_{n≥0}, with a.(−1.ι) the umbral counterpart of a sum of a ∈ ℕ r.v.’s with uniform distribution on the interval [0, 1]. Indeed,

k! N_{k}(x, an) = ∑_{i=1}^{k} 𝔼[(x − (an).(−1.ι))^{i}]B_{k,i}(m_{1}, …, m_{k−i+1}),

where B_{k,i} are partial Bell exponential polynomials and m_{i} = 𝔼[(u^{<−1>})^{i}], for i = 1, …, k.

10. Acknowledgments

I would like to thank Domenico Senato and Henry Crapo for helpful discussions and suggestions during the writing of this survey. I would like to also thank Pasquale Petrullo and Rosaria Simone for carefully reading the paper, leading to an improvement of its technical quality. A special thank to Christian Krattenthaler whose remarks and criticisms have stimulated the author in looking for a simplified exposure devoted to readers coming in new to the matter.

References

[1] David F. Andrews and Jamie E. Stafford, Symbolic Computation for Statistical Inference, Oxford Statistical Science Series, vol. 21, Oxford University Press, 2000.

[2] Michael Anshelevich, Partition-dependent stochastic measures and qdeformed cumulants, Doc. Math. 6 (2001), 343–384.

[3] David Applebaum, Lévy Processes and Stochastic Calculus, Cambridge University Press, 2004.

[4] Florin Avram and Murad S. Taqqu, Symmetric polynomials of random variables attracted to an infinitely divisible law, Probab. Theory Rel. Fields 27 (1986), 491–500.

[5] Florin Avram and Murad S. Taqqu, Noncentral limit theorems and Appell polynomials, Ann. Probab. 15 (1987), 767–775.

[6] Ole E. Barndorff-Nielsen, Thomas Mikosch and Sidney I. Resnick, Lévy Processes: Theory and Applications, Birkhäuser, 2001.

[7] William D. Baten, Correction for the moments of a frequency distribution in two variables, Ann. Math. Statist. 2 (1931), 309–319.

[8] Eric T. Bell, Exponential polynomials, Ann. of Math. 35 (1934), 258–277.

[9] Bruce C. Berndt, Ramanujan’s Notebooks, Part 1 , Springer-Verlag, New York, 1985.

[10] James W. Brown, On multivariable Sheffer sequences, J. Math. Anal. Appl. 69 (1979), 398–410.

[11] Włodzimierz Bryc, Computing moments of free additive convolution of measures, Appl. Math. Comput. 194 (2007), 561–567.

[12] Gregory M. Constantine and Thomas H. Savits, A stochastic representation of partition identities, SIAM J. Discrete Math. 7 (1994), 194–202.

[13] Gregory M. Constantine and Thomas H. Savits, A multivariate Faà di Bruno formula with applications, Trans. Amer. Math. Soc. 348 (1996), 503–520.

[14] Cecil C. Craig, Sheppard’s corrections for a discrete variable, Ann. Math. Statist. 7 (1936), 55–61.

[15] Arthur P. Dempster and Donald B. Rubin, Rounding error in regression: the appropriateness of Sheppard’s corrections, J. R. Statist. Soc. Ser. B: Statist. Methodol. 45 (1983), 51–59.

[16] Alessandro Di Bucchianico and Daniel S. Loeb, A selected Survey of Umbral Calculus, Electron. J. Combin. 2 Dynamic Survey 3, updated 2000, 28 pp.

[17] Elvira Di Nardo, A new approach to Sheppard’s corrections, Math. Methods Statist. 12 (2010), 151–162.

[18] Elvira Di Nardo, Moment symbolic calculus in statistics, 58-th World Statistics Congress of the International Statistical Institute, Dublin, 2011,

http://www.youtube.com/watch?v=bf7k1NGLQSY&feature=youtu.be.

[19] Elvira Di Nardo, On a symbolic representation of non-central Wishart random matrices with applications, J. Mult. Anal. 125 (2014), 121–135.

[20] Elvira Di Nardo and Giuseppe Guarino, New algorithms for computing Sheppard’s corrections, Worksheet Maple Software, 2010, available at

http://www.maplesoft.com/applications/view.aspx?SID=59515.

[21] Elvira Di Nardo, Giuseppe Guarino and Domenico Senato, A unifying framework for k-statistics, polykays and their multivariate generalizations, Bernoulli 14 (2008), 440–468.

[22] Elvira Di Nardo, Giuseppe Guarino and Domenico Senato, Symbolic computation of moments of sampling distributions, Comput. Statist. Data Anal. 52 (2008), 4909–4922.

[23] Elvira Di Nardo, Giuseppe Guarino and Domenico Senato, A new method for fast computing unbiased estimators of cumulants, Statist. Comput. 19 (2009), 155–165.

[24] Elvira Di Nardo, Giuseppe Guarino and Domenico Senato, Fast algorithms for k-statistics, polykays and their multivariate generalizations, Worksheet Maple Software, 2009, available at http://www.maplesoft.comappli

cations/view.aspx?SID=33041.

[25] Elvira Di Nardo, Giuseppe Guarino and Domenico Senato, A new algorithm for computing the multivariate Faà di Bruno’s formula, Appl. Math. Comp. 217 (2011), 6286–6295.

[26] Elvira Di Nardo, Peter McCullagh and Domenico Senato, Natural statistics for spectral samples, Ann. Stat. 41, (2013), 982–1004.

[27] Elvira Di Nardo, Heinrich Niederhausen and Domenico Senato, A symbolic handling of Sheffer sequences, Ann. Mat. Pura Appl. 190 (2011), 489–506.

[28] Elvira Di Nardo and Imma Oliva, On the computation of classical, boolean and free cumulants, Appl. Math. Comput. 208 (2009), 347–354.

[29] Elvira Di Nardo and Imma Oliva, On some applications of a symbolic representation of non-centered Lévy processes, Comm. Statist. Theory Methods 42 (2013), 3974–3988.

[30] Elvira Di Nardo and Imma Oliva, A new family of time-space harmonic polynomials with respect to Lévy processes, Ann. Mat. Pura Appl. (4) 192 (2013), 917–929.

[31] Elvira Di Nardo and Imma Oliva, Multivariate Bernoulli and Euler polynomials via Lévy processes, Appl. Math. Lett. 25 (2012), 1179–1184.

[32] Elvira Di Nardo, Pasquale Petrullo and Domenico Senato, Cumulants and convolutions via Abel polynomials, Europ. J. Combin. 31 (2010), 1792–1804.

[33] Elvira Di Nardo and Domenico Senato, Umbral nature of the Poisson random variables, in: Algebraic Combinatorics and Computer Science: a tribute to Gian-Carlo Rota, Henry Crapo and Domenico Senato (eds.), Springer–Verlag, Milano, 2001, pp. 245–266.

[34] Elvira Di Nardo and Domenico Senato, An umbral setting for cumulants and factorial moments, Europ. J. Combin. 27 (2006), 394–413.

[35] Elvira Di Nardo and Domenico Senato, The eleventh and twelveth problems of Rota’s Fubini Lectures: from cumulants to free probability theory, in: From Combinatorics to Philosophy: the Legacy of G.-C. Rota, E. Damiani, O. D’Antona, V. Marra, F. Palombi (eds.), Springer–Verlag, 2009, pp. 91–130.

[36] Elvira Di Nardo and Domenico Senato, Symbolic solutions of some linear recurrences, J. Statist. Plann. Inference 142(2), 423–429.

[37] Peter Doubilet, On the foundations of combinatorial theory VII: Symmetric functions through the theory of distribution and occupancy, Stud. Appl. Math. 51 (1972), 377–396.

[38] Paul L. Dressel, Statistical seminvariants and their estimates with particular emphasis on their relation to algebraic invariants, Ann. Math. Statist. 11 (1940), 33–57.

[39] William Feller, An Introduction to Probability Theory and Its Applications, vol. II, John Wiley & Sons, New York, 1966.

[40] Pedro G. Ferreira, Joao Magueijo and Joseph Silk, Cumulants as non-Gaussian qualifiers, Phys. Rev. D 56 (1997), 4592–4603.

[41] Ronald A. Fisher, Moments and product moments of sampling distributions, Proc. Lond. Math. Soc. 30 (1929), 199–238.

[42] Philippe Flajolet and Robert Sedgewick, Analytic Combinatorics, Cambridge Univ. Press, 2009.

[43] Anindya Goswami and Arindam Sengupta, Time-space polynomial martingales generated by a discrete-time martingale, J. Theoret. Probab. 8 (1995), 417–431.

[44] Anders Hald, On the history of the correction for grouping 1873−1922, Scand. J. Statist. 28 (2001), 417–428.

[45] Michael Hardy, Combinatorics of partial derivatives, Electron. J. Combin. 13 (2006), Article #R1, 13 pp.

[46] Daniel F. Heitjan, Inference from grouped continuous data: a review, Statist. Sci. 4 (1989), 164–183.

[47] Thomas Kailath and Adrian Segall, Orthogonal functionals of independent increment processes, IEEE Trans. Inform. Theory 22 (1976), 287–298.

[48] Edward L. Kaplan, Tensor notation and the sampling cumulants of kstatistics, Biometrika 39 (1952), 319–323.

[49] John E. Kolassa, Series Approximation Methods in Statistics, Lecture Notes in Statistics, vol. 88, Springer–Verlag, 1994.

[50] Germain Kreweras, Sur les partitions non croisées d’un cycle, Discrete Math. 1 (1972), 333–350.

[51] Stephan Lawi, Hermite and Laguerre polynomials and matrix valued stochastic processes, Elect. Comm. Probab. 13 (2008), 67–84.

[52] Roy B. Leipnik and Charles E. M. Pearce, The multivariate Faà di Bruno formula and multivariate Taylor expansions with explicit integral remainder term, Austral. New Zealand Industr. Appl. Math. J. 48 (2007), 327–341.

[53] Yuh-Dauh Lyuu, Financial Engineering and Computation: Principles, Mathematics, Algorithms, Cambridge University Press, Cambridge, 1999.

[54] Peter McCullagh, Tensor Methods in Statistics, Monographs on Statistics and Applied Probability, Chapman & Hall, London, 1987.

[55] Joachim D. Müller, Cumulant analysis in fluorescence fluctuation spectroscopy, Biophys J. 86 (2004), 3981–3992.

[56] Jacques Neveu, Discrete-Parameter Martingales, North-Holland/American Elsevier, 1975.

[57] Alexandru Nica and Roland Speicher, A “Fourier transform” for multiplicative functions on non-crossing partitions, J. Algebraic Combin. 6 (1997), 141– 160.

[58] Alexander Nica and Roland Speicher, Lectures on the Combinatorics of Free Probability, London Mathematical Society, Lecture Note Series, vol. 335, Cambridge University Press, 2006.

[59] Heinrich Niederhausen, Finite Operator Calculus, project book, 2010,

available at http://math.fau.edu/niederhausen/HTML/Research/Umbral
Calculus/bookS2010.pdf.

[60] Silvia Noschese and Paolo E. Ricci, Differentiation of multivariable composite functions and Bell polynomials, J. Comput. Anal. Appl. 5 (2003), 333–340.

[61] Christopher Parrish, Multivariate umbral calculus, Linear Multilinear Algebra 6 (1978), 93–109.

[62] Giovanni Peccati and Murad S. Taqqu, Wiener Chaos: Moments, Cumulants and Diagrams: A Survey with Computer Implementation, Bocconi and Springer, Milano, 2011.

[63] Pasquale Petrullo, A symbolic treatment of Abel polynomials, in: From Combinatorics to Philosophy: The Legacy of G.-C. Rota, E. Damiani, O. D’Antona, V. Marra, F. Palombi (eds.), Springer–Verlag, Milano, 2009, 210–221.

[64] Pasquale Petrullo, A Symbolic Approach to Cumulants and Their Connections with The Abel Identity and The Kerov Polynomials., Ph.D. thesis, Università degli Studi della Basilicata, Potenza, 2011.

[65] Giovanni Pistone and Henry P. Wynn, Cumulant varieties, J. Symbolic Comput. 41 (2006), 210–221.

[66] Jim Pitman, Some probabilistic aspects of set partitions, Amer. Math. Monthly 104 (1997), 201–209.

[67] Jim Pitman and Richard P. Stanley, A polytope related to empirical distributions, plane trees, parking functions and the associahedron, Discrete Comput. Geom. 27, (2002), 603–634.

[68] Yves Poupard, Étude et dénombrement parallèles des partitions non-croisées d’un cycle et des découpages d’un polygone convexe, Discrete Math. 2 (1972), 279–288.

[69] Sudhakar Prasad and Nicolas Menicucci, Fisher information with respect to cumulants, I.E.E.E. Trans. Inform. Theory 50 (2004), 638–642.

[70] John Riordan, An Introduction to Combinatorial Analysis, John Wiley and Sons, Inc., New York, 1958.

[71] Steven M. Roman and Gian-Carlo Rota, The umbral calculus, Adv. Math. 27 (1978), 95–188.

[72] Colin Rose and Murray D. Smith, Mathematical Statistics with Mathematica, Springer Text in Statistics, Springer–Verlag, 2002.

[73] Gian-Carlo Rota, Twelve problems in probability no one likes to bring up, in: Algebraic Combinatorics and Computer Science: a tribute to Gian-Carlo Rota, Henry Crapo and Domenico Senato (eds.), Springer–Verlag, Milano, 2001, pp. 57–93.

[74] Gian-Carlo Rota and Jianhong Shen, On the combinatorics of cumulants, J. Combin. Theory Ser. A 91 (2000), 283–304.

[75] Gian-Carlo Rota and Brian D. Taylor, The classical umbral calculus, SIAM J. Math. Anal. 25 (1994), 694–711.

[76] Gian-Carlo Rota, David Kahaner and Andrew Odlyzko, On the foundations of combinatorial theory. VIII. Finite operator calculus, J. Math. Anal. Appl. 42 (1973), 684–760.

[77] Gian-Carlo Rota, Jianhong Shen and Brian D. Taylor, All polynomials of binomial type are represented by Abel polynomials, Ann. Sci. Norm. Super. Pisa Cl. Sci. 25 (1998), 731–738.

[78] Gian-Carlo Rota and Timothy C. Wallstrom, Stochastic integrals: a combinatorial approach, Ann. Probab. 25 (1997), 1257–1283.

[79] Sandra Saliani and Domenico Senato, Compactly supported wavelets through the classical umbral calculus, J. Fourier Anal. Appl. 12 (2006), 27–36.

[80] Ken-Iti Sato, Lévy Processes and Infinitely Divisible Distributions, Cambridge Studies in Advanced Mathematics, vol. 68, Cambridge University Press, Cambridge, 1999.

[81] Thomas H. Savits, Some statistical applications of Faà di Bruno, J. Multivariate Anal. 97 (2006), 2131–2140.

[82] Wim Schoutens, Lévy Processes in Finance: Pricing Financial Derivatives, Wiley Series in Probability and Statistics, 2003.

[83] Wim Schoutens, Stochastic Processes and Orthogonal Polynomials, Lecture Notes in Statistics, vol. 146, Springer–Verlag, 2000.

[84] Arindam Sengupta, Markov processes, time-space harmonic functions and polynomials, Statist. Probab. Lett. 78 (2008), 3277–3280.

[85] Jianhong Shen, Combinatorics for wavelets: the umbral refinement equation, Stud. Appl. Math. 103 (1999), 121–147.

[86] William F. Sheppard, On the calculation of most probable values of frequencyconstants for data arranged according to equidistant division of a scale, Proc. London Math. Soc. 29 (1898), 353–380.

[87] Peter J. Smith, A recursive formulation of the old problem of obtaining moments from cumulants and vice versa, Amer. Statist. 49 (1995), 217–218.

[88] Terry P. Speed, Cumulants and partition lattices, Austral. J. Statist. 25 (1983), 378–388.

[89] Roland Speicher, Multiplicative functions on the lattice of non-crossing partitions and free convolution, Math. Ann. 298 (1994), 611-628.

[90] Roland Speicher and Reza Woroudi, Boolean convolution, in: Free probability theory (Waterloo, ON, 1995), D. Voiculescu (ed.), Fields Inst. Commun., vol. 12, Amer. Math. Soc., Providence, R.I., 1997, pp. 267–279.

[91] Adriaan J. Stam, Polynomials of binomial type and compound Poisson processes, J. Math. Anal. Appl. 130 (1988), 493–508.

[92] Richard P. Stanley, Enumerative Combinatorics, vol. 1, Cambridge University Press, Cambridge, 1997.

[93] Benjamin Staude, Stefan Rotter and Sonja Grün, CuBIC: cumulant based inference of higher-order correlations in massively parallel spike trains, J. Comput. Neurosci. 29 (2010), 327–350.

[94] Alan Stuart and Keith Ord, Kendall’s Advanced Theory of Statistics, Vol. 1: Distribution Theory, sixth ed., Edward Arnold, London, U.K., 1994.

[95] Bernd Sturmfels and Piotr Zwiernik, Binary cumulant varieties, Ann. Combin. 17 (2013), 229–250.

[96] Bryan D. Taylor, Umbral presentations for polynomial sequences, Comput. Math. Appl. 41 (2001), 1085–1098.

[97] Thorvald N. Thiele, On an approximation formula, Tidsskr. Math. 3 (1873), 22–31.

[98] Thorvald N. Thiele, Elementär Jagttagelsesläre, Gyldendal, Copenhagen, 1897; translated into English: The theory of observations, London, 1903; reprinted in Ann. Math. Statist. 2 (1931), 165–308.

[99] David C. Torney, Binary cumulants, Adv. Appl. Math. 25 (2000), 34–40.

[100] Jacques Touchard, Nombres exponentiels et nombres de Bernoulli, Canad. J. Statist. 8 (1956), 305–320.

[101] John W. Tukey, Some sampling simplified, J. Amer. Statist. Assoc. 45 (1950), 501–519.

[102] Stephen B. Vardeman, Sheppard’s correction for variances and “the quantization noise model”, I.E.E.E. Trans. Instr. Meas. 54 (2005), 2117–2119.

[103] Dan Voiculescu, Symmetries of some reduced free product C^{∗}-algebras, in: Lecture Notes in Math., vol. 1132, Springer–Verlag, 1985, pp. 556–588.

[104] Dan Voiculescu, Lectures on free probability theory, in: Lectures on probability theory and statistics (Saint-Flour, 1998), P. Bernard (ed.), Lecture Notes in Math., vol. 1738, Springer–Verlag, 2000, pp. 281–349.

[105] Jan Vrbik, Population moments of sampling distributions, Comput. Statist. 20 (2005), 611–621.

[106] Herbert S. Wilf, generatingfunctionology, Academic Press, 2006.

[107] Christopher S. Withers, A simple expression for the multivariate Hermite polynomials, Statist. Probab. Lett. 47 (2000), 165–169.

[108] Piotr Zwiernik and Jim Q. Smith, Tree cumulants and the geometry of binary tree models, Bernoulli 18 (2012), 290–321.

[109] Piotr Zwiernik, L-cumulants, L-cumulant embeddings and algebraic statistics, J. Algebr. Statist. 3 (2012), 11–43.

Dept. of Mathematics, Computer Science and Economics, University of Basilicata, Campus Macchia Romana, Via dell’Ateneo Lucano 10, 85100 Potenza, Italy

E-mail address: elvira.dinardo@unibas.it