Some new Stein operators for product distributions

TL;DR

This paper develops a general method for deriving Stein operators for the product of two independent random variables, extending previous work and applying to various distributions including normals and gamma, with insights into their complexity.

Contribution

It introduces a unified framework for Stein operators of product distributions, covering non-centered normals, gamma, and variance-gamma distributions, expanding prior results.

Findings

Derived Stein operators for products of independent variables.

Provided a simple derivation of the characteristic function for product of normals.

Explained the increased complexity of the PDF in non-centered cases.

Abstract

We provide a general result for finding Stein operators for the product of two independent random variables whose Stein operators satisfy a certain assumption, extending a recent result of Gaunt, Mijoule and Swan \cite{gms18}. This framework applies to non-centered normal and non-centered gamma random variables, as well as a general sub-family of the variance-gamma distributions. Curiously, there is an increase in complexity in the Stein operators for products of independent normals as one moves, for example, from centered to non-centered normals. As applications, we give a simple derivation of the characteristic function of the product of independent normals, and provide insight into why the probability density function of this distribution is much more complicated in the non-centered case than the centered case.

Tables

Table 1. Table 1: Stein operators for products of normal random variables.

Product $P$	Stein operator $A_{P}f\left(x\right)$ (here we set $\sigma :={\sigma }_X{\sigma }_Y$)
$N(0,{\sigma }_X^2)\times N(0,{\sigma }_Y^2)$	${\sigma }^2\left(x{f}^{\prime \prime }\left(x\right)+xf\left(x\right)\right)-xf\left(x\right)$
$N(\mu ,{\sigma }_X^2)\times N(\mu ,{\sigma }_Y^2)$	${\sigma }^{3}x{f}^{\left(3\right)}\left(x\right)+{\sigma }^2\left(\sigma -x\right)f^{\prime \prime }\left(x\right)-\sigma \left(x+\left(\sigma +{\mu }^2\right)\right)f^{\prime }\left(x\right)$
	$+\left(x-{\mu }^2\right)f\left(x\right)$
$N({\mu }_X,{\sigma }_X^2)\times N({\mu }_Y,{\sigma }_Y^2)$	${\sigma }_X^4{\sigma }_Y^{4}x{f}^{\left(4\right)}\left(x\right)+{\sigma }_X^4{\sigma }_Y^4{f}^{\left(3\right)}\left(x\right)-{\sigma }_X^2{\sigma }_Y^2\left(2x+{\mu }_X{\mu }_Y\right)f^{\prime \prime }\left(x\right)$
	$-\left({\sigma }_X^2{\sigma }_Y^2+{\mu }_X^2{\sigma }_Y^2+{\mu }_Y^2{\sigma }_X^2\right)f^{\prime }\left(x\right)+\left(x-{\mu }_X{\mu }_Y\right)f\left(x\right)$

Equations

$$A{f}_h(x)=h(x)-\mathrm{I}\mathrm{E}h(X),$$

$$P(MD)M=MP(MD+I),\text{and}DP(MD)=P(MD+I)D$$

$${\tau }_{a}P(MD)=P(MD){\tau }_a.$$

$$A=M-Q(MD)-P(MD)D$$

$$A_Z=R_1(MD)D^2+R_2(MD)D+R_3(MD)+M{R}_4(MD),$$

$$R_1(U)$$

$$R_2(U)$$

$$R_3(U)$$

$$R_4(U)$$

$$\mathrm{I}\mathrm{E}[Xf(Z)]$$

$$=\mathrm{I}\mathrm{E}[Q(U){\tau }_{Y}f(X)+P(U)D{\tau }_{Y}f(X)]$$

$$=\mathrm{I}\mathrm{E}[{\tau }_{Y}Q(U)f(X)+Y{\tau }_{Y}P(U)Df(X)]$$

$$\mathrm{I}\mathrm{E}[Xf(Z)]$$

$$\mathrm{I}\mathrm{E}[Yf(Z)]=\mathrm{I}\mathrm{E}[Q(U)f(Z)+XP(U)Df(Z)].$$

$$\mathrm{I}\mathrm{E}[X(I-P(U)DP(U)D)f(Z)]=\mathrm{I}\mathrm{E}[(Q(U)+Q(U)P(U)D)f(Z)],$$

$$\mathrm{I}\mathrm{E}[X(I-P(U+I)P(U)D^2)f(Z)]=\mathrm{I}\mathrm{E}[(Q(U)+Q(U)P(U)D)f(Z)].$$

$$\mathrm{I}\mathrm{E}[Zf(Z)]$$

$$\mathrm{I}\mathrm{E}[X^{2}f(Z)]$$

$$=\mathrm{I}\mathrm{E}[Q(U)M{\tau }_{Y}f(X)+P(U)DM{\tau }_{Y}f(X)]$$

$$=\mathrm{I}\mathrm{E}[M{\tau }_{Y}Q(U+I)f(X)+{\tau }_{Y}P(U)(U+I)f(X)]$$

$$\mathrm{I}\mathrm{E}[X^{2}f(Z)]$$

$$\mathrm{I}\mathrm{E}[(M-P(U)P(U)(U+I)D)f(Z)]=\mathrm{I}\mathrm{E}[X\left(Q(U)+Q(U+I)P(U)D\right)f(Z)].$$

$$\mathrm{I}\mathrm{E}[X(Q(U)P(U)D+Q(U-I))f(Z)]$$

$$=\mathrm{I}\mathrm{E}[(UP(U-I)-(P(U){)}^{2}P(U+I)(U+I)D^2$$

$$+(Q(U)+Q(U)P(U)D)Q(U-I))f(Z)].$$

$$A=M-\alpha T_a-\beta T_{b}D.$$

$$A_Z=(M-{\alpha }^2{T}_a^2-{\beta }^2{T}_b^2{T}_{1}D)(T_{a-1}-\beta T_b{T}_{a+1}D)-2{\alpha }^2\beta T_a^2{T}_b{T}_{a+1}D.$$

$$A_Z$$

$$=M{D}^3+(I-M)D^2-(M+(1+{\mu }^2)I)D+M-{\mu }^{2}I.$$

$$A_{Z}f(x)$$

$$=x(f^{(3)}(x)-f^{\prime \prime }(x))+(f^{\prime \prime }(x)-f^{\prime }(x))+x(f^{\prime }(x)-f(x)).$$

$$A_Z=(M-T_{r+\mu }^2-{\mu }^2{T}_{1}D)(T_{r+\mu -1}+\mu T_{r+\mu +1}D)+2\mu T_{r+\mu }^2{T}_{r+\mu +1}D.$$

$$f(x)=\frac{1}{\sigma\sqrt{\pi}\Gamma(\frac{r}{2})}\mathrm{e}^{\frac{\theta}{\sigma^{2}}(x-\mu)}\bigg{(}\frac{|x-\mu|}{2\sqrt{\theta^{2}+\sigma^{2}}}\bigg{)}^{\frac{r-1}{2}}K_{\frac{r-1}{2}}\bigg{(}\frac{\sqrt{\theta^{2}+\sigma^{2}}}{\sigma^{2}}|x-\mu|\bigg{)},$$

$$A=(M-4{\theta }^2{T}_{r/2}^2-{\sigma }^4{T}_r^2{T}_{1}D)(T_{r/2-1}-{\sigma }^4{T}_r{T}_{r/2+1}D)-8{\theta }^2{\sigma }^2{T}_{r/2}^2{T}_r{T}_{r/2+1}D.$$

$$Af(x)=(M-{\sigma }^4{T}_r^2{T}_{1}D)(T_{r/2-1}-{\sigma }^4{T}_r{T}_{r/2+1}D)f(x).$$

$$Ag(x)$$

$${\sigma }^2(M-\mu )D^2+(r{\sigma }^2+2\theta (M-\mu ))D+(r\theta -(M-\mu ))I,$$

Full text

60E15,

62E15,

Stein’s method,

Stein operators,

product distributions,

product of independent normal random variables,

and

a]The University of Manchester b]INRIA Paris c]Université libre de Bruxelles

1 Introduction

In 1972, Charles Stein [25] introduced a powerful technique for deriving explicit bounds in normal approximations. Shortly after, in 1975, Louis Chen [5] adapted the method to the Poisson distribution, and since then Stein’s method has been extended to a wide variety of distributional approximations. For a given target distribution $p$, the first step in the general procedure is to find a suitable operator $A$ acting on a class of functions $\mathcal{F}$ such that ${\rm I\kern-1.60004ptE}[Af(X)]=0$ for all $f\in\mathcal{F}$, where the random variable $X$ has distribution $p$. The operator $A$ is called the Stein operator, and for continuous distributions is typically a differential operator; for the $N(\mu,\sigma^{2})$ distribution, the classical operator is $Af(x)=\sigma^{2}f^{\prime}(x)-(x-\mu)f(x)$. This leads to the Stein equation

[TABLE]

where $h$ is a real-valued function. If $A$ is well chosen then, for a given $h$, the Stein equation (1.1) can be solved for $f_{h}$, and the problem of estimating the proximity of the distribution of a random variable $W$ of interest to the distribution of the target random variable $X$, as measured by $|{\rm I\kern-1.60004ptE}h(W)-{\rm I\kern-1.60004ptE}h(X)|$, reduces to one of bounding $|{\rm I\kern-1.60004ptE}[Af_{h}(W)]|$. For a detailed account of the method we refer the reader to the monograph Stein [26].

In addition to the normal and Poisson distributions, Stein’s method has been adapted to many classical distributions, such as the exponential (Chatterjee, Fulman and Röllin [4]), gamma (Luk [18]) and Laplace (Pike and Ren [22]), as well as quite general families of distributions, such as the Pearson family (Schoutens [24]), variance-gamma distributions (Gaunt [8]) and a wide class of distributions satisfying a certain diffusive assumption (Döbler [7], Kusuoka and Tudor [15]); for an overview see Ley, Reinert and Swan [16]. As such, over the years, a number of techniques have been developed for finding Stein operators for a variety of distributions. These include the density method (Stein [26], Ley, Reinert and Swan [16], Ley and Swan [17], Mijoule, Reinert and Swan [19]), the generator method (Barbour [3], Götze [14]), the differential equation duality approach (Gaunt [10], Ley, Reinert and Swan [16]), and probability generating function and characteristic function based approaches of Upadhye, Čekanavičius and Vellaisamy [27] and Arras et al. [2]. The corpus of literature concerning Stein operators and their applications is now vast, and it continues growing at a steady pace. Stein operators provide handles on target distributions which are in some sense just as important and natural characteristics of a probability distribution as its moments, its moment generating function, its p.d.f., c.d.f. or even its characteristic function. Finding tractable Stein operators is thus, naturally, an important question.

In this paper we pursue the work begun in Gaunt [9] and Gaunt [11] concerning the following question : “given two independent random variables $X$ and $Y$ with Stein operators $A_{X}$ and $A_{Y}$, can one find a Stein operator for $Z=XY$?” More specifically, the present paper is a complement (sequel) to our paper Gaunt, Mijoule and Swan [13] where we developed an algebraic technique for finding Stein operators for products of independent random variables with polynomial Stein operators satisfying a technical condition. Let $M(f)=(x\mapsto xf(x))$, $D(f)=(x\mapsto f^{\prime}(x))$ and $I$ be the identity operator. We say that the absolutely continuous variates $X$ and $Y$ have polynomial Stein operators if they allow Stein operators of the form $A=\sum_{i,j}a_{ij}M^{i}D^{j}$ for $a_{ij}$ some real numbers. The highest value of $j$ such that $a_{ij}\neq 0$ is called the order of the operator. In Gaunt, Mijoule and Swan [13] we provided a method for deriving operators under the technical assumption that $\#\left\{j-i\,|\,a_{ij}\neq 0\right\}\leq 2$ (see Assumption 3 and Lemma 2.6 of Gaunt, Mijoule and Swan [13] for more details on this condition). For such random variables, Proposition 2.12 of Gaunt, Mijoule and Swan [13] gives a polynomial Stein operator for the product $XY$. A number of classical random variables have Stein operators which satisfy this assumption, such as the $N(0,\sigma^{2})$ distribution with Stein operator $\sigma^{2}D-M$, with others including the gamma, beta, and even some more exotic distributions such as the zero-mean symmetric variance-gamma distribution and PRR distribution of Pëkoz, Röllin and Ross [21]. However, some very natural densities do not satisfy the assumption. In fact, even the non-centered normal distribution does not satisfy this assumption, as its Stein operator $\sigma^{2}D+\mu I-M$ instead satisfies $\#\left\{j-i\,|\,a_{ij}\neq 0\right\}=3$. In Proposition 2.1, we shall address the natural problem of extending the result of Gaunt, Mijoule and Swan [13] to treat the product of two independent random variables satisfying this new assumption. Here we have only added one level of complexity in the operator; nevertheless, as we will see later on, it is sufficient to include the classical cases of non-centered normal and non-centered gamma, and a more general sub-family of the variance-gamma distributions. Also, as noted in Remark 2.3, the proof technique is novel and seems to be a useful addition to the toolkit for finding Stein operators.

The Stein operators for the products of independent normal random variables are particularly theoretically interesting, and we devote Section 3 to exploring some of their properties. For the case of two independent centered normals a second order Stein operator was obtained by Gaunt [9], whereas, rather curiously, we find a third order operator for the product of two i.i.d. normals, and a fourth order operator for the product of two independent general normals; see Table 1. It is an important and natural question to ask whether our operators have minimal order amongst all Stein operators with polynomial coefficients. We believe this is the case but are unable to prove it. However, in Section 3.1, we are able to provide a brute force approach for verifying this assertion for polynomial coefficients up to a particular order. This brute force approach is very general and in principle can be applied to any polynomial Stein operators. In Section 3.2, we prove that our Stein operators for products of independent normals characterise the distribution. We do this by appealing to a more general result, Proposition 3.2, which treats distributions that are determined by their moments.

For the Stein operator of Gaunt [9] for the product of two independent standard normal random variables, it was possible to solve the corresponding Stein equation and bound the derivatives of the solution. As a result, Gaunt [9] was able to derive explicit bounds for product normal approximations. However, it seems to be beyond the scope of existing techniques in the Stein’s method literature to solve and then bound the derivatives of the solution to our more complicated third and fourth order Stein equations for products of non-centered normals. It should be noted, though, that there is still great utility to Stein equations even when it is not possible to obtain bounds for the solution. For example, as has been demonstrated in several papers such as Nourdin, Peccati and Swan [20], Arras et al. [1] and Arras et al. [2], Stein operators can be used for comparison of probability distributions directly without solving Stein equations. We also stress that Stein operators are also of use in applications beyond proving approximation theorems; for example, in obtaining distributional properties (Gaunt [9], Gaunt [11], Gaunt, Mijoule and Swan [13]). Indeed, in Section 3.3, we use our Stein operators to obtain a simple derivation of the characteristic function of two independent normals, and also provide valuable insight into why there is a dramatic increase in complexity in the probability density function from the centered to non-centered case.

2 New Stein operators for product distributions

2.1 A general result

Throughout this paper, we shall make the following assumptions, which were also made in Gaunt, Mijoule and Swan [13]; we refer the reader to that paper for some remarks on these assumptions.

**Assumption ** 1). $X$ admits a smooth density $p$ with respect to the Lebesgue measure on ${\rm I\kern-1.60004ptR}$; this density is defined and non-vanishing on some (possibly unbounded) interval $J\subseteq{\rm I\kern-1.60004ptR}$. 2). $X$ admits an operator $A$ acting on $\mathcal{F}$ which contains the set of smooth functions with compact support $\mathcal{C}_{0}^{\infty}({\rm I\kern-1.60004ptR})$.

Let $P$ be a real polynomial. Then it is easily proved (by checking it when $P$ is a monomial, then by linearity) that

[TABLE]

(recall the notations $M(f)=(x\mapsto xf(x))$, $D(f)=(x\mapsto f^{\prime}(x))$ and $I$ the identity operator from the introduction). Now, for $a\in{\rm I\kern-1.60004ptR}\setminus\{0\}$, let $\tau_{a}(f)=\left(x\mapsto f(ax)\right)$. Simple computations show that (see Gaunt, Mijoule and Swan [13, Lemma 2.5]) $\tau_{a}M=aM\tau_{a}$ and $D\tau_{a}=a\tau_{a}D.$ This implies that for any real polynomial $P$,

[TABLE]

Proposition 2.1**.**

Let $X$ and $Y$ be i.i.d. with common Stein operator of the form

[TABLE]

for $P,Q$ two real polynomials. Then, a Stein operator for $Z=XY$ is

[TABLE]

where, for $U=MD$,

[TABLE]

Proof.

Let $Z=XY$ and $f\in\mathcal{F}$. Denote $U=MD$. We have

[TABLE]

Similarly,

[TABLE]

Replace $f$ with $P(U)Df$ in (2.4) and add up to (2.3) to get

[TABLE]

which is also, using (2.1),

[TABLE]

Now using (2.4) and conditioning, we can compute

[TABLE]

We also have

[TABLE]

Thus we obtain by (2.6)

[TABLE]

Apply (2.7) to $P(U)Df$ and add up to (2.5) applied to $Q(U-I)f$ to obtain

[TABLE]

Apply the preceding equation to $Q(U)f$ and subtract to (2.7) applied to $Q(U-I)f$ to get the result. ∎

The case that $P$ and $Q$ are polynomials of degree one is important, as it is applicable to non-centered normal and non-centered gamma random variables, as well as a general sub-family of the variance-gamma distributions. To this end, let us define the operator $T_{r}:=MD+rI$. We note that the limit of $T_{r}$ as $r\rightarrow\infty$ is ill-defined, but we do have $\lim_{r\rightarrow\infty}r^{-1}T_{r}=I$ (see Gaunt, Mijoule and Swan [13, Remark 2.3]).

Corollary 2.2**.**

Let $\alpha,\beta\in{\rm I\kern-1.60004ptR}$ and $a,b\in{\rm I\kern-1.60004ptR}\cup\{\infty\}$ (if either $a$ or $b$ are set to $+\infty$, then we proceed as described above). Let $X,Y$ be i.i.d. with common Stein operator

[TABLE]

Then, a Stein operator for $Z=XY$ is

[TABLE]

Proof.

Set $Q(U)=\alpha(U+a)$ and $P(U)=\beta(U+b)$ in (2.2). A calculation then verifies that (2.8) and (2.2) are equivalent operators in this case (up to a factor $\alpha$). ∎

Remark 2.3**.**

The proof of Proposition 2.1 involves applying certain equations to test functions of the form $Lf$, where $L$ is a linear differential operator. This allowed us to cancel terms to obtain (2.2). We consider this technique to be a useful addition to the toolkit for finding Stein operators. Indeed, this approach was recently used by Gaunt [12] to find Stein operators for the $H_{3}(Z)$ and $H_{4}(Z)$, where $H_{n}$ is the $n$-th Hermite polynomial and $Z\sim N(0,1)$. In Section 2.2.4, we also use the technique to derive a Stein operator for the product of independent non-centered normals with different means.

Remark 2.4**.**

We attempted to generalise Proposition 2.1 so that $X$ and $Y$ are no longer identically distributed, for which $X$ and $Y$ have Stein operators of the form $A_{X}=M-Q_{X}(MD)-P_{X}(MD)D$ and $A_{Y}=M-Q_{Y}(MD)-P_{Y}(MD)D$. We were only able to find a Stein operator for the product $XY$ under the very restrictive condition that $P_{Y}(U)Q_{X}(U)Q_{X}(U+I)=P_{X}(U)Q_{Y}(U)Q_{Y}(U+I)$. This Stein operator had the unusual feature of not being symmetric in $X$ and $Y$. In certain simple cases, we can, however, apply the proof technique of Proposition 2.1 to derive a Stein operator for the product of two non-identically distributed random variables; see Section 2.2.4.

Remark 2.5**.**

Note that, whilst the Stein operator for $X$ and $Y$ in Proposition 2.1 satisfies the condition $\#\left\{j-i\,|\,a_{ij}\neq 0\right\}=3$, the Stein operator (2.2) for their product satisfies $\#\left\{j-i\,|\,a_{ij}\neq 0\right\}=4$. Thus, it is not possible to iterate Proposition 2.1 to find a Stein operator for product of three i.i.d. random variables. This is in contrast to the work of Gaunt, Mijoule and Swan [13] which was carried out under the assumption $\#\left\{j-i\,|\,a_{ij}\neq 0\right\}=2$.

2.2 Examples

2.2.1 Product of non-centered

normals

Assume $X$ and $Y$ are independent standard normal random variables. A Stein operator for $X+\mu$ (or $Y+\mu$) is $A=D-M+\mu I$. Applying Corollary 2.2 with $\alpha=\mu$, $\beta=1$ and $a=b=\infty$ gives the following Stein operator for $Z=(X+\mu)(Y+\mu)$:

[TABLE]

(Here, and for the rest of this paper, we consider the unit variance case; the extension to general case follows from a straightforward rescaling and the resulting Stein operator for the product is given in Table 1.) Note that when $\mu=0$ the operator becomes

[TABLE]

Taking $g(x)=f^{\prime}(x)-f(x)$ then yields $A_{Z}f(x)=\tilde{A}_{Z}g(x)=xg^{\prime\prime}(x)+g^{\prime}(x)-xg(x)$, which we recognise as the product normal Stein operator that was obtained by Gaunt [9].

2.2.2 Product of non-centered gammas

Assume $X$ and $Y$ are distributed as a $\Gamma(r,1)$, with p.d.f. $p(x)=\frac{1}{\Gamma(r)}x^{r-1}\mathrm{e}^{-x}$, $x>0$, and let $\mu\in{\rm I\kern-1.60004ptR}$. A Stein operator for $X+\mu$ (or $Y+\mu$) is $A=T_{r+\mu}-\mu D-M.$ Corollary 2.2 applied with $\alpha=1$, $\beta=-\mu$, $a=r+\mu$, $b=\infty$ yields the following fourth-order Stein operator for $Z=(X+\mu)(Y+\mu)$:

[TABLE]

Note also that when $\mu=0$, this operator reduces to $(M-T_{r}^{2})T_{r-1}$, which is the product gamma Stein operator of Gaunt [11] applied to $T_{r-1}f$ instead of $f$.

2.2.3 Product of variance-gamma random variables

The variance-gamma distribution with parameters $r>0$, $\theta\in{\rm I\kern-1.60004ptR}$, $\sigma>0$, $\mu\in{\rm I\kern-1.60004ptR}$ has p.d.f.

[TABLE]

$x\in{\rm I\kern-1.60004ptR}$, where $K_{\nu}(x)=\int_{0}^{\infty}\mathrm{e}^{-x\cosh(t)}\cosh(\nu t)\,\mathrm{d}t$, $x>0$, is the modified Bessel function of the second kind. If a random variable $W$ has density (2.10) then we write $W\sim\mathrm{VG}(r,\theta,\sigma,\mu)$. A $\mathrm{VG}(r,\theta,\sigma,0)$ Stein operator is given by $\sigma^{2}T_{r}D+2\theta T_{r/2}-M$ (see Gaunt [8]). Applying Corollary 2.2 with $\alpha=2\theta,\beta=\sigma^{2},a=r/2,b=r$, we get the following Stein operator for the product of two independent $\mathrm{VG}(r,\theta,\sigma,0)$ random variables:

[TABLE]

Note that when $\theta=0$ we have

[TABLE]

Defining $g:{\rm I\kern-1.60004ptR}\rightarrow{\rm I\kern-1.60004ptR}$ by $xg(x)=-(T_{r/2-1}-\sigma^{4}T_{r}T_{r/2+1}D)f(x)$ gives

[TABLE]

which is in agreement with the product variance-gamma Stein operator given in Section 3.2 of Gaunt, Mijoule and Swan [13]. Lastly, we note that the $\mathrm{VG}(r,\theta,\sigma,\mu)$ Stein operator of Gaunt [8], as given by

[TABLE]

satisfies $\#\left\{j-i\,|\,a_{ij}\neq 0\right\}=4$ when $\mu\not=0$, and therefore one cannot apply Proposition 2.1 or Corollary 2.2 to find a Stein operator for the product of two such variates.

2.2.4 Product of non-identically distributed non-central normals

By working on a case-by-case basis it is possible to use the proof technique of Proposition 2.1 to find Stein operators for the product of two non-identically distributed random variables, whose Stein operators satisfy the assumptions of the proposition. We find that a Stein operator for the product of independent normals $N(\mu_{X},1)$ and $N(\mu_{Y},1)$ is

[TABLE]

Let us now provide a derivation of this Stein operator. Let $X$ and $Y$ be independent standard normal random variables and define $Z=(X+\mu_{X})(Y+\mu_{Y})$. We will use repeatedly the fact that ${\rm I\kern-1.60004ptE}[Wg(W)]={\rm I\kern-1.60004ptE}[g^{\prime}(W)]$ for $W\sim N(0,1)$, as well as conditioning arguments, and we let ${\rm I\kern-1.60004ptE}_{W}[]$ stand for the expectation conditioned on $W$. Let $f:{\rm I\kern-1.60004ptR}\rightarrow{\rm I\kern-1.60004ptR}$ be four times differentiable and such that ${\rm I\kern-1.60004ptE}|Zf^{(i)}(Z)|<\infty$ for $i=0,1,\ldots 4$ and ${\rm I\kern-1.60004ptE}|f^{(i)}(Z)|<\infty$ for $i=0,1,2,3$, where $f^{(0)}\equiv f$. Then

[TABLE]

By again applying a conditioning argument we obtain

[TABLE]

(and the same applies to ${\rm I\kern-1.60004ptE}[Yf^{\prime}(Z)]$). Hence

[TABLE]

Isolating the expressions depending on $Y$ from (2.12) and (2.13), we obtain two different equations:

[TABLE]

and

[TABLE]

Substract (2.15) to (2.2.4) to get

[TABLE]

from which we deduce that (2.11) is a Stein operator for $Z$. $\Box$

Lastly, we note that applying the operator (2.9) to $f(x)=g^{\prime}(x)+g(x)$ yields

[TABLE]

which we recognise as the Stein operator (2.11) in the special case $\mu_{X}=\mu_{Y}=\mu$.

2.2.5 Sums of products of normals

Let us begin by noting a simple result, that has perhaps surprisingly not previously been stated explicitly in the literature. Suppose $X,X_{1},\ldots,X_{n}$ are i.i.d., with Stein operator $A_{X}f(x)=\sum_{k=0}^{m}(a_{k}x+b_{k})f^{(k)}(x)$, where $m\geq 1$ and the $a_{k}$ and $b_{k}$ are real-valued constants. Let $W=\sum_{j=1}^{n}X_{j}.$ Then, by conditioning,

[TABLE]

Thus, a Stein operator for $W$ is given by

[TABLE]

Remark 2.6**.**

Identity (2.16) actually generalises similar observations for score functions and Stein kernels, for which such an additive stability is well-known, see Nourdin, Peccati and Swan [20].

Since the coefficients in the Stein operators (2.9) and (2.11) are linear, we can use (2.16) to write down a Stein operator for the sum $W=\sum_{i=1}^{r}X_{i}Y_{i}$, where $(X_{i})_{1\leq i\leq r}\sim N(\mu_{X},1)$ and $(Y_{i})_{1\leq i\leq r}\sim N(\mu_{Y},1)$ are independent. When $\mu_{X}=\mu_{Y}=\mu$, we have

[TABLE]

and when $\mu_{X}$ and $\mu_{Y}$ are not necessarily equal, we have

[TABLE]

When $\mu_{X}=\mu_{Y}=0$, the random variable $W$ follows the $\mathrm{VG}(r,0,1,0)$ distribution (see Gaunt [8], Proposition 1.3). Taking $g=f^{\prime}-f$ in (2.17) gives $A_{W}g(x)=xg^{\prime\prime}(x)+rg^{\prime}(x)-xg(x)$, which we recognise as the $\mathrm{VG}(r,0,1,0)$ Stein operator that was obtained in Gaunt [8].

3 Some results concerning the Stein operators for products of independent normal random variables

3.1 On the minimality of the operators

The operator (2.8) is at most a seventh order differential operator. However, for particular cases, such as the product of two i.i.d. non-centered normals, the operator reduces to one of lower order, see Section 2.2.1. Whilst we believe that the third order operator (2.9) is a minimal order polynomial operator, we have no proof of this claim (nor do we have much intuition as to whether the seventh order operator (2.8) is of minimal order). We believe this question of minimality to be important:

Conjecture 3.1**.**

There exists no second order Stein operator (acting on smooth functions with compact support) with polynomial coefficients for the product of two independent non-centered normal random variables.

One can use a brute force approach to verify the conjecture for polynomials of fixed order (if the conjecture holds). Such results would be worthwhile in practice, because a third order Stein operator with linear coefficients may be easier to use in applications than one of second order with polynomial coefficients of degree greater than one.

Let us now us the brute force approach to prove that there is no second order Stein operator with linear coefficients for the product of two independent non-centered normals (generalisations are obvious). Let $X$ and $Y$ be independent $N(1,1)$ random variables and let $Z=XY$. Suppose that there was such a Stein operator for $Z$, then it would be of the form $A_{Z}f(x)=\sum_{j=0}^{2}(a_{0,j}+a_{1,j}x)f^{(j)}(x)$, where $f^{(0)}\equiv f$. Now, if $A_{Z}$ was a Stein operator for $Z$, we would have ${\rm I\kern-1.60004ptE}[A_{Z}f(Z)]=0$ for all $f$ in some class $\mathcal{F}$ that contains the monomials $\{x^{k}\,:\,k\geq 1\}$. Taking $f(x)=x^{k}$, $k=0,1,\ldots,5$, we obtain six equations for six unknowns. Letting $\mu_{k}$ denote ${\rm I\kern-1.60004ptE}Z^{k}$, we have $\mu_{1}=1$, $\mu_{2}=4$, $\mu_{3}=16$, $\mu_{4}=100$, $\mu_{5}=676$ and $\mu_{6}=5776$. This leads to the system of equations

[TABLE]

We used Mathematica to compute that the determinant of the matrix corresponding to this system of equations is $783\,360\not=0$. Therefore, there is a unique solution, which is clearly $a_{1,2}=\cdots=a_{0,0}=0$. Thus, there does not exist a second order Stein operator with linear coefficients for $Z$.

Similarly, one can show that there is no third order Stein operator with linear coefficients for the product of two independent normals with different means. Here we took $X\sim N(1,1)$ and $Y\sim N(2,1)$, and sought a Stein operator of the form $A_{Z}f(x)=\sum_{j=0}^{3}(a_{0,j}+a_{1,j}x)f^{(j)}(x)$. We then used the monomials $f(x)=x^{k}$, $k=0,1,\ldots,7$, to generate eight linear equations in eight unknowns, and found the determinate of the matrix corresponding to this system of equations to be $10\,157\,222\,707\,200\not=0$.

3.2 Characterisation by the operators

We begin with a simple general result, which perhaps surprisingly has not previously been given in the literature. The proof technique has, however, appeared in the literature; see the proof of Lemma 5.2 of Ross [23] for case of the exponential distribution.

Proposition 3.2**.**

Suppose that the law of the random variable $X$, supported on $I\subset{\rm I\kern-1.60004ptR}$, is determined by its moments. Let the operator $A_{X}=\sum_{i=1}^{n}\sum_{j=1}^{p}a_{i,j}M^{j}D^{i}$, where $a_{i,j}\in{\rm I\kern-1.60004ptR}$, act on a class of functions $\mathcal{F}$ which contains all polynomial functions. Suppose $A_{X}$ is a Stein operator for $X$: that is, for all $f\in\mathcal{F}$,

[TABLE]

Now, let $m=\max_{i,j}(j-i)-\min_{i,j}(j-i)-1$, where the maxima and minima are taken over all $i,j$ such that $a_{i,j}\not=0$. Suppose that the first $m$ moments of $Y$ are equal to those of $X$ and that

[TABLE]

for all $f\in\mathcal{F}$. Then $Y$ has the same law as $X$.

Proof.

We prove that all moments of $Y$ are equal to those of $X$. As the moments of $X$ determine its law, verifying this proves the Proposition. The monomials $\{x^{k}\,:\,k\geq 1\}$ are contained in the class $\mathcal{F}$, so applying $f(x)=x^{k}$, $k\geq m$, to (3.2) yields the recurrence

[TABLE]

where $C_{k}=k(k-1)\cdots(k-i+1)$ if $k-i+1>0$ and $C_{k}=0$ otherwise. We have that ${\rm I\kern-1.60004ptE}Y^{0}=1$ and we are given that ${\rm I\kern-1.60004ptE}Y^{k}={\rm I\kern-1.60004ptE}X^{k}$ for $k=1,\ldots,m$. We can then use forward substitution in (3.3) to (uniquely) obtain all moments of $Y$. Due to (3.1), ${\rm I\kern-1.60004ptE}[A_{X}f(X)]=0$ for all $f\in\mathcal{F}$, and so it follows by the above reasoning that

[TABLE]

But this is same recurrence relation as (3.3) and, since ${\rm I\kern-1.60004ptE}Y^{k}={\rm I\kern-1.60004ptE}X^{k}$ for $k=1,\ldots,m$, it follows that ${\rm I\kern-1.60004ptE}Y^{k}={\rm I\kern-1.60004ptE}X^{k}$ for all $k\geq m$ as well. ∎

If we have obtained a Stein operator $A_{X}$ for a random variable $X$, then Proposition 3.2 tells us that the operator characterises the law of $X$ if $X$ is determined by its moments. This characterisation is weaker than those typically found in Stein’s method literature, as it involves moment conditions on the random variable $Y$. This is perhaps not surprising, because the characterisations given in the literature have mostly been found on a case-by-case basis, whereas ours applies to a wide class of distributions.

The distribution of the product of two independent normal distributions is determined by its moments, which can be seen from the existence of its moment generating function $M(s)$ for all $|s|<1$; see Section 3.3.1. The following full characterisation of the distribution is thus immediate from Proposition 3.2.

Proposition 3.3**.**

(i) Let $W$ be a real-valued random variable whose first three moments are equal to that of the random variable $Z=XY$, where $X\sim N(\mu_{X},1)$ and $Y\sim N(\mu_{Y},1)$ are independent. Then $W$ is equal in law to $Z$ if and only if

[TABLE]

for all functions $f\in C^{4}({\rm I\kern-1.60004ptR})$ such that ${\rm I\kern-1.60004ptE}|Zf^{(j)}(Z)|<\infty$ for $0\leq j\leq 4$, and ${\rm I\kern-1.60004ptE}|f^{(k)}(Z)|<\infty$ for $0\leq k\leq 3$, where $f^{(0)}\equiv f$.

(ii) Now suppose that $\mu_{X}=\mu_{Y}=\mu$, and that the first two moments of $W$ are equal to those of $Z$. Then $W$ is equal in law to $Z$ if and only if

[TABLE]

for all $f\in C^{3}({\rm I\kern-1.60004ptR})$ such that ${\rm I\kern-1.60004ptE}|Zf^{(j)}(Z)|<\infty$, $0\leq j\leq 3$, and ${\rm I\kern-1.60004ptE}|f^{(k)}(Z)|<\infty$, $0\leq k\leq 2$.

Proposition 3.2 can be used to prove that some other Stein operators given in the literature fully characterise the distribution. For example, the Stein operator for the product of $n$ independent Beta random variables of Gaunt [11] is characterising, since this product is supported on $(0,1)$ and thus the distribution is determined by its moments.

3.3 Applications of the operators

3.3.1 Characteristic function

As the Stein operator (2.11) has linear coefficients, it turns out to be straightforward to use the characterising equation (3.3) to find a formula for the characteristic function of the random variable $Z=XY$, where $X\sim N(\mu_{X},1)$ and $Y\sim N(\mu_{Y},1)$ are independent.

On taking $f(x)=\mathrm{e}^{\mathrm{i}tx}$ in the characterising equation (3.3) and setting $\phi(t)={\rm I\kern-1.60004ptE}[\mathrm{e}^{\mathrm{i}tZ}]$, we deduce that $\phi(t)$ satisfies the differential equation

[TABLE]

It should be noted that $f(x)=\mathrm{e}^{\mathrm{i}tx}$ is a complex-valued function; here we have applied the characterising equation to the real and imaginary parts of $f$, which are themselves real-valued functions. Solving (3.5) subject to the condition that $\phi(0)=1$ then gives that

[TABLE]

Setting $s=\mathrm{i}t$ yields a formula for the moment generating function $M(s)={\rm I\kern-1.60004ptE}[\mathrm{e}^{sZ}]$, which is well-defined for $|s|<1$. We doubt these formulas are new, but it is interesting to note that we were able to obtain such a simple proof via the Stein characterisation.

3.3.2 Probability density function

Let $X\sim N(\mu_{X},1)$ and $Y\sim N(\mu_{Y},1)$ be independent, and let $Z=XY$. For $\mu_{X}=\mu_{Y}=0$, it is a well-known and easy to prove result that the p.d.f. is given by $p_{Z}(x)=\frac{1}{\pi}K_{0}(|x|)$, $x\in{\rm I\kern-1.60004ptR}$. However, in general, the p.d.f. takes a much more complicated form (see Cui et al. [6])): for $x\in{\rm I\kern-1.60004ptR}$,

[TABLE]

It is possible to use the Stein operators for the product $Z$ to gain insight into why there is such a dramatic increase in complexity from the zero mean case to non-zero mean case. To see this, we recall a duality result given in Remark 2.7 of Gaunt, Mijoule and Swan [13] (see also Section 4 of that paper for further details). If $V$ admits a smooth density $p$, which solves the differential equation $Bp=0$ with $B=\sum_{i,j}b_{ij}M^{j}D^{i}$, then a Stein operator for $V$ is given by $A=\sum_{i,j}(-1)^{i}b_{ij}D^{i}M^{j}$, and similarly given a Stein operator for $V$ one can write down a differential equation satisfied by $p$. In this manner, we can write down differential equations satisfied by the density $p_{W}$ of the random variable $W=\sum_{i=1}^{r}X_{i}Y_{i}$, where the $X_{i}$ and $Y_{i}$ are independent copies of $X$ and $Y$ respectively, using the Stein operators (2.17) and (2.18) for this distribution. When $\mu_{X}=\mu_{Y}=\mu$, we have

[TABLE]

and in general

[TABLE]

In the special case $\mu_{X}=\mu_{Y}=0$, the density of $Z$ satisfies the modified Bessel differential equation $xp_{Z}^{\prime\prime}(x)+p_{Z}^{\prime}(x)-xp_{Z}(x)=0$.

From Section 3.1 and the duality result of Gaunt, Mijoule and Swan [13], we know there do not exist differential equations for $p_{Z}$ with linear coefficients with a lower degree than (3.8) and (3.9). Moreover, we were unable to transform (3.8) or (3.9) into a well-understood class, such as the Meijer $G$-function differential equation. Therefore, the increase in complexity in the p.d.f. $p_{Z}$ of $Z$ from the zero mean to non-zero mean case can be understood from the increase in complexity of the differential equation satisfied by $p_{Z}$. Also, due to the above reasoning, it seems plausible that formula (3.7) cannot be simplified further.

Finally, we note that there is not a severe increase in complexity in the differential equations satisfied by $W$ from the $r=1$ case to the general case. To the best of our knowledge, a formula for general $r\geq 1$ has not been obtained in the literature, and even if the differential equations (3.8) and (3.9) are not ultimately used to derive such a formula, they do indicate that the formula should be at a similar level of complexity to that of (3.7), and thus provide motivation for obtaining such a formula. We note that such a result would be of interest due to the occurrence of such random variables in, for example, electrical engineering applications, see Ware and Lad [28].

Acknowledgements

RG acknowledges support from EPSRC grant EP/K032402/1 and is currently supported by a Dame Kathleen Ollerenshaw Research Fellowship. RG is grateful to Université de Liège, FNRS and EPSRC for funding a visit to University de Liège, where some of the details of this project were worked out. YS acknowledges support by the Fonds de la Recherche Scientifique - FNRS under Grant MIS F.4539.16. Part of GM’s research was supported by a Welcome Grant from Université de Liège. We would like to thank the referee for their careful reading of our paper and their helpful comments.