# Random First-Order Linear Discrete Models and Their Probabilistic Solution: A Comprehensive Study

Authors: M.-C. Casabán; J.-C. Cortés; J.-V. Romero; M.-D. Roselló
Journal: Abstract and Applied Analysis
 (2016)

---

## Abstract

This paper presents a complete stochastic solution represented by the first probability density function for random first-order linear difference equations. The study is based on Random Variable Transformation method. The obtained results are given in terms of the probability density functions of the data, namely, initial condition, forcing term, and diffusion coefficient. To conduct the study, all possible cases regarding statistical dependence of the random input parameters are considered. A complete collection of illustrative examples covering all the possible scenarios is provided.


---

## Body

## 1. Introduction and Motivation

The birth/death rates of species in biology, the volatility of assets in finance, the transmission rates of the spread of epidemics or social addictions in epidemiology, the diffusion and advection coefficients of mass transport processes in physics, and so forth are quantities that, in practice, involve uncertainty. Thus, their deterministic modelling is clearly limited. This motivates the search of mathematical models that consider randomness in their formulation. Deterministic differential and difference equations have been demonstrated to be useful mathematical representations for modelling numerous real problems. The consideration of randomness in these types of equations is a relatively recent research area whose main goal is to extend classical deterministic results to the random scenario. Regarding continuous models, most of the contributions have focussed on Itô-type stochastic differential equations. In this class of differential equations, uncertainty is considered through a Gaussian and stationary stochastic process (SP) called white noise, which is the derivative of the Wiener SP [1–3]. Autoregressive (AR) models can be interpreted as their discrete counterpart. These types of models are extensively used in time series analysis in statistics [4]. Some recent interesting models based on Itô-type stochastic differential equations include [5–7], for instance. Complementary to these approaches, uncertainty can be directly introduced in differential and difference equations by assuming that coefficients, source term, and/or initial/boundary conditions are random variables and/or stochastic processes. Under this approach, the probability distributions associated with RVs and SPs are not required to be Gaussian. This approach leads to the area usually referred to as random differential/difference equations [8], [9, p. 66]. Intensive studies on random differential and difference equations have been undertaken only over the last few decades. Currently, they are exerting a profound influence on the analysis of many problems in engineering and science [10, 11]. Most of these contributions are based on mean square calculus [8, 12–14].The solution of a random difference equation is a discrete SP, say{Zn:n≥0}. In dealing with random difference equations, the main goals are computing the solution SP and its statistical characteristics, such as the mean function, E[Zn], and the variance function, V[Zn]=E[(Zn)2]-(E[Zn])2. Despite being more complicated, the computation of the first probability density function (1-PDF), f1(z,n), associated with solution Zn is more convenient since from 1-PDF, besides determining the mean and the variance functions, one can also compute higher-order statistical moments of Zn:(1)EZnk=∫-∞∞zkf1z,ndz,n,k=0,1,2,….In the context of random differential equations, a number of contributions have dealt with the computation of the 1-PDF in specific problems appearing in physics [15–17] or in mathematics [18, 19]. A comprehensive study for the random first-order linear differential equation under general hypotheses has been recently published by the authors of [20]. The unifying element of all these contributions is the Random Variable Transformation (RVT) method. This technique permits, under certain hypotheses that will be specified later, the computation of the PDF of a random variable (RV) resulting after mapping another RV whose PDF is known [21–23]. Although RVT technique is a classical probability result, it must be pointed out that its application to study differential equations with randomness is likely due to Professor M. El-Tawil. He was the first author who had used RVT technique to present some approximate solutions of random differential equations. In [24], the RVT method together with different numerical schemes (finite difference and Runge-Kutta) is implemented to get the 1-PDF of the solution SP in solving both partial and/or ordinary (first- and/or second-order) random differential equations. His contribution includes the definition of probabilistic error and a formula for its computation as well.As continuation of the study initiated in [20], the aim of this paper is to develop a comprehensive study to determine the 1-PDF of the solution discrete SP of the following random initial value problem:(2)Zn+1=AZn+B,n=0,1,2,…,Z0=Γ0,taking advantage of RVT method. All the input parameters, Γ0,  B, and A, are assumed to be continuous RVs defined on a common probability space (Ω,F,P). Although RVT method is the unifying technique used to conduct the analysis of model (2) and also the one studied in [20], it is important to point out that both problems have distinctly different nature. Indeed, model (2) is discrete whereas problem faced in [20] is continuous counterpart. As we will see later, Proposition 1 constitutes the key result we have had to establish to conduct the analysis of problem (2). Its formulation is based on RVT method. Throughout the paper, significant differences between the new analytical expressions and graphical behaviour of the 1-PDF of the solution of (2) and its continuous counterpart will be also exhibited.As it also happens in the deterministic framework, in general, the study of random difference equations has been less prolific than of random differential equations. In [25], the authors study the mean square exponential stability of impulsive stochastic difference equations. In [26], one studies random matrix linear difference equations assuming that diffusion coefficient A in (2) is a deterministic matrix rather than a RV. In [26], the authors focus on the computation of the mean vector and variance-covariance matrix of the solution discrete SP instead of the 1-PDF. In this paper, we present a comprehensive study of model (2) assuming randomness in all the inputs, Γ0,  B, and A. This includes the general case where inputs are statistically dependent. From a statistical point of view, (2) is a generalization of autoregressive model of order 1, AR(1), where uncertainty is enclosed in the term B through white noise.For the sake of clarity in the presentation and in order to facilitate the comparison of the results obtained in this paper against the ones achieved in [20], we will keep the notation used in both contributions identical. Hence, the domain of the random inputs Γ0, B, and A will be denoted by(3)DΓ0=γ0=Γ0ω,ω∈Ω:γ0,1≤γ0≤γ0,2,DB=b=Bω,ω∈Ω:b1≤b≤b2,DA=a=Aω,ω∈Ω:a1≤a≤a2, respectively. From this point forward, we will omit the ω-notation when writing RVs. In this manner, for instance, we will write Γ0 rather than Γ0(ω). The same can be said for the notation of the PDFs that appear throughout this paper. For example, fΓ0(γ0) will denote the PDF of RV Γ0; fB,A(b,a) will denote the joint PDF of RVs B and A; we will write fΓ0,B,A(γ0,b,a) for the joint PDF of the random vector (Γ0,B,A), and so on. As usual, we will assume that any PDF is null outside its domain.Based on the same arguments exhibited in [20, Section 1], we will distinguish the thirteen cases listed in Table 1 to conduct our study. These casuistries consider whether the random difference equation (2) is homogeneous or nonhomogeneous as well as all possible cases regarding the random or deterministic nature of the input parameters Γ0, B, and A. Note that, by splitting the study in all these cases, the comparison of the results concerning the discrete problem (2) against its continuous counterpart is facilitated. Examples have also been devised with the same aim. Even more, in the majority of the examples, the same statistical distributions have been taken as in [20] to highlight better analogies and differences between both models.Table 1
List of the thirteen different cases considered to conduct the full study. This classification is made regarding whether the discrete initial value problem is homogeneous (H) or nonhomogeneous (NH) and the way that uncertainty is considered (one random variable or a random vector in two or three dimensions).
Type
Discrete initial value problem
Case
H
Z
n
+
1
=
A
Z
n
,
n
=
0,1
,
2
,
…
Z
0
=
Γ
0
(
I
)
I.1Γ0 is a random variable I.2  A is a random variable I.3  (Γ0,A) is a random vector
NH
Z
n
+
1
=
Z
n
+
B
,
n
=
0,1
,
2
,
…
Z
0
=
Γ
0
(
I
I
)
II.1Γ0 is a random variable  II.2  B is a random variable  II.3  (Γ0,B) is a random vector
P
ω
∈
Ω
:
A
ω
=
1
=
0  Zn+1=AZn+B,n=0,1,2,…Z0=Γ0(III)
III.1Γ0 is a random variable  III.2  B is a random variable  III.3  A is a random variable  III.4  (Γ0,B) is a random vector  III.5  (Γ0,A) is a random vector III.6  (B,A) is a random vector III.7  (Γ0,B,A) is a random vectorThe paper is organized as follows. Section2 is addressed to introduce the preliminaries related to RVT technique required to conduct our study. In this section, we establish a key result related to the PDF of the power transformation of RVs which will be crucial to deal with Case I.2 of Table 1, where uncertainty just enters in model (2) through the RV A. Section 3 is divided into three subsections where the 1-PDF of the discrete solution SP of (2) is determined for each one of Cases I, II, and III listed in Table 1. Illustrative examples covering the thirteen cases are provided throughout the paper. In the last section, we present our conclusions. Finally, we present an appendix where the main obtained results are collected in order to facilitate their practical use.


## 2. Preliminaries

As it has been pointed out in the previous section, the goal of this paper is to compute the1-PDF {f1(z,n):n≥0} of the solution SP of problem (2) in each one of the cases listed in Table 1. The key result to achieve this goal is the RVT method. This is a probabilistic technique that allows us to calculate the PDF of a random variable/vector which is obtained after mapping another random variable/vector whose PDF is known. Depending on the type of mapping as well as its dimension, several versions of RVT method can be established. Throughout this paper, the general scalar version and its specialization to the linear case, as well as the general multidimensional version, will be required. These results are stated in [20, Theorem 1, Proposition 2 and Theorem 4], respectively.Next, we will establish the following result concerning the PDF of a RV which is obtained after mapping another RV via a power transformation. This result will play a relevant role in the analysis of Case I.2. listed in Table1. It is important to underline the notion that power transformation is a distinctive feature to describe the solution SP of the discrete model (2) against the exponential transformation which appears when its continuous counterpart is dealt with. In this sense, the next result plays the same role as [20, Proposition 3] performed there.Proposition 1 (RVT technique: power transformation).
LetX be a continuous RV with domain DX=x:x1≤x≤x2 and PDF fX(x). Let one denote by δ(·) the Dirac delta function. Then, the PDF fY(y) of the power transformation Y=kXn, with n≥0 and k≠0, is given by the following:
(i) Ifn=0, (4)fYy=δy-k,-∞<y<+∞.
(ii) Ifn=1,(5)fYy=1kfXyk,kx1≤y≤kx2,ifk>0,kx2≤y≤kx1,ifk<0.
(iii) Ifn is even, then one has the following:
Case 1 (x1≥0):(6)fYy=1knyk1-n/nfX+ykn,kx1n≤y≤kx2n,ifk>0,kx2n≤y≤kx1n,ifk<0.
Case 2 (x2≤0):(7)fYy=1knyk1-n/nfX-ykn,kx2n≤y≤kx1n,ifk>0,kx1n≤y≤kx2n,ifk<0.
Case 3 (x1x2<0).
Case 3.1 (x2≥|x1|):(8)fYy=fY1y+fY2y, where(9)fY1y=1knyk1-n/nfX-ykn+fX+ykn, on the domain(10)0<y≤kx1n,ifk>0,kx1n≤y<0,ifk<0,fY2y=1knyk1-n/nfX+ykn,kx1n<y≤kx2n,ifk>0,kx2n≤y<kx1n,ifk<0.
Case 3.2 (x2<|x1|):(11)fYy=fY1y+fY2y, where(12)fY1y=1knyk1-n/nfX-ykn+fX+ykn, on the domain(13)0<y≤kx2n,ifk>0,kx2n≤y<0,ifk<0,fY2y=1knyk1-n/nfX-ykn,kx2n<y≤kx1n,ifk>0,kx1n≤y<kx2n,ifk<0.
(iv) Ifn≥3 and is odd, then one has the following:
Case 1 (x1>0orx2<0):(14)fYy=1knyk1-n/nfXykn,kx1n≤y≤kx2n,ifk>0,kx2n≤y≤kx1n,ifk<0.
Case 2 (x1x2≤0):(15)fYy=1knyk1-n/nfXykn,kx1n≤y<0,0<y≤kx2n,ifk>0,kx2n≤y<0,0<y≤kx1n,ifk<0.If k=0, for n≥0,(16)fYy=δy,-∞<y<+∞.Proof.
(i) Ifn=0, then Y=k w.p. 1 and its PDF is given by(17)fYy=δy-k,-∞<y<+∞.
(ii) Ifn=1, the mapping r is monotone on the whole domain of RV X; then the inverse function of r takes the form(18)x=sy=yk,whose derivative is given by(19)s′y=1k.Then, applying expression [20, Eq. (3)] and taking into account (18)-(19), one gets(20)fYy=1kfXyk,kx1≤y≤kx2,ifk>0,kx2≤y≤kx1,ifk<0.
(iii) Let us assumen is even. We distinguish three cases depending on the domain DX=x:x1≤x≤x2 of the RV X. In order to avoid any confusion, it is important to note that these cases are mutually exclusive.
Case 1  (
x
1
≥
0
). Assuming k>0, the mapping r is monotone in DX. Hence, the inverse function of r, denoted by s(y), takes the form(21)x=sy=+ykn,and its derivative, s′(y), is given by(22)s′y=+1knyk1-n/n.Then, applying [20, Eq. (3)] and taking into account (21)-(22), one gets(23)fYy=1knyk1-n/nfX+ykn,kx1n≤y≤kx2n.When k<0, the reasoning is analogous. Then, considering expression (23) when k>0, both cases for the sign of k can be expressed as follows:(24)fYy=1knyk1-n/nfX+ykn,on the domain(25)kx1n≤y≤kx2n,ifk>0,kx2n≤y≤kx1n,ifk<0.
Case 2  (
x
2
≤
0
). In this case, the mapping r is monotone in DX. Using an analogous development as we did in Case 1, one obtains(26)fYy=1knyk1-n/nfX-ykn, on the domain(27)kx2n≤y≤kx1n,ifk>0,kx1n≤y≤kx2n,ifk<0.
Case 3  (
x
1
x
2
<
0
). We will consider two subcases: x2≥|x1| and x2<|x1|. We split each subcase into appropriate subintervals in order to apply [20, Theorem 1] in order to compute the PDF.
Case 3.1  (
x
2
≥
|
x
1
|
). Let us consider the piece [x1,-x1]. On the subinterval [x1,0], the mapping r (denoted by r1, for the sake of clarity) is monotone and its inverse s1 is(28)x=s1y=-ykn,whose derivative, s1′(y), for y≠0, is given by(29)s1′y=-1knyk1-n/n.On the other hand, on the piece [0,-x1], its corresponding mapping r2 is monotone and its inverse s2 is(30)x=s2y=+ykn,and, for y≠0,(31)s2′y=+1knyk1-n/n.Notice that s1′(y)≠0 and s2′(y)≠0 if y≠0. Then, applying [23, Theorem 2.1.8] and taking into account (28)–(31), one gets(32)fY1y=fXs1ys1′y+fXs2ys2′y=1knyk1-n/nfX-ykn+fX+ykn,on the domain(33)0<y≤kx1n,ifk>0,kx1n≤y<0,ifk<0.As usual, we assume fY1(y)≡0 outside domain (33).
To complete the computation of PDFfY(y) on the whole domain, finally we consider the subinterval ]-x1,x2], where the RV X is positive. Hence, we are in Case 1 and according to (24) it follows that(34)fY2y=1knyk1-n/nfX+ykn,on the domain(35)k-x1n<y≤kx2n,ifk>0,kx2n≤y<k-x1n,ifk<0.Again, as usual, we assume fY2(y)≡0 outside domain (35). Notice that one satisfies (36)∫0kx1nfY1ydy+∫kx1nkx2nfY2ydy=1ifk>0,∫kx2nkx1nfY2ydy+∫kx1n0fY1ydy=1ifk<0. To summarize, from (32)–(35), the complete PDF of Y=kXn in this case is given by (37)fYy=fY1y+fY2y,0<y≤kx2n,ifk>0,kx2n≤y<0,ifk<0.
Case 3.2  (
x
2
<
|
x
1
|
). Let us consider the piece [-x2,x2]. Following analogous reasoning as in Case 3.1, according to (32), one obtains the piece of the PDF of the power transformation Y:(38)fY1y=1knyk1-n/nfX-ykn+fX+ykn,on the domain(39)0<y≤kx2n,ifk>0,kx2n≤y<0,ifk<0.We assume fY1(y)≡0 outside domain (39).
We complete the computation of PDFfY(y) on the whole domain considering the subinterval [x1,-x2[. In this subinterval, X is negative. As it was shown in Case 2, the PDF is given by(40)fY2y=1knyk1-n/nfX-ykn,on the domain(41)kx2n<y≤kx1n,ifk>0,kx1n≤y<kx2n,ifk<0.We assume fY2(y)≡0 outside domain (41).
To summarize, from (38)–(41), the complete PDF of Y=kXn in this case is(42)fYy=fY1y+fY2y,0<y≤kx1n,ifk>0,kx1n≤y<0,ifk<0.
(iv) Let us assume thatn≥3 and is odd. The mapping r is monotone on the whole domain of RV X; then the inverse function of r takes the form(43)x=sy=ykn,whose derivative, for y≠0, is given by(44)s′y=1knyk1-n/n.Notice that s′(y)≠0 if y≠0. Therefore, we distinguish two cases depending on the domain DX=x:x1≤x≤x2 of the RV X.
Case 1  (
x
1
>
0
or
x
2
<
0
). Applying [20, Theorem 1] and taking into account (43)-(44), one gets(45)fYy=1knyk1-n/nfXykn,kx1n≤y≤kx2n,ifk>0,kx2n≤y≤kx1n,ifk<0.
Case 2  (
x
1
x
2
≤
0
). Applying [20, Theorem 1] and taking into account (43)-(44), one gets(46)fYy=1knyk1-n/nfXykn,kx1n≤y<0,0<y≤kx2n,ifk>0,kx2n≤y<0,0<y≤kx1n,ifk<0.Finally, if k=0 and n≥0, then Y=0 with probability 1 (w.p. 1) and its PDF is given by(47)fYy=δy,-∞<y<+∞.


## 3. Case Study: Homogeneous Discrete Initial Value Problem (I)

This section is addressed to compute the1-PDF {f1(z,n):n≥0} of the solution discrete stochastic process {Zn:n≥0} of the homogeneous discrete initial value problem (I) in all different cases collected in Table 1. In this case, the solution {Zn:n≥0} can be expressed as follows:(48)Z0=Γ0,Zn=AnΓ0,n=1,2,….


### 3.1. Case I.1:Γ0 Is a Random Variable

For the sake of clarity in the presentation, we rewrite solution (48) by using the lowercase letter a in order to indicate the deterministic character of parameter A:(49)Z0=Γ0,Zn=anΓ0,n=1,2,….Letn≥0 be an arbitrary and fixed integer. Let us assume a≠0 and denote Z=Zn=anΓ0. By applying [20, Proposition 2] to(50)Y=Z,X=Γ0,α=an≠0,β=0, one gets the 1-PDF(51)f1z,n=1anfΓ0zan,n=0,1,2,…,z∈R.Note that ifa=0, from (49), it follows that Zn=0 w.p. 1 for each n≥1 and Z0=Γ0. So, the 1-PDF for the trivial case, a=0, can be written as(52)f1z,n=fΓ0z,n=0,δz,n=1,2,…,z∈R.In order to facilitate the comparison of the1-PDF of the solution of problem (49) against its continuous counterpart provided in [20], in the following example, Γ0 is assumed to be a Gaussian RV. Note that we are going to consider a standard distribution, although the method is also able to be applied to nonstandard distributions.Example 1.
Let us assumea≠0 and consider Γ0 a Gaussian distribution, Γ0~Nμ;σ2. Hence, applying (51), the 1-PDF of {Zn:n≥0} is given by(53)f1z,n=1an2πσ2e-z/an-μ2/2σ2,n=0,1,2,…,z∈R. It can be checked that f1(z,n) is a PDF for each n≥0. Figure 1 shows f1(z,n) for n∈0,…,10, in the particular case that Γ0~N(0;1): a=10/9 (a) and a=9/10 (b). Note the different behavior of 1-PDF depending on the modulus of the parameter a. This is in agreement with the expectation and variance of the solution which are given, respectively, by(54)EZn=0,VZn=a2nσ2,n=0,1,2,…. Indeed, in Figure 1, we observe that, for each n, the 1-PDF is symmetric about z=0, whereas, in the case that a=10/9>1 (a=9/10∈-1,1), it becomes flat (sharp) as n increases. This means that its variability around zero, that is, the variance, tends to infinity (zero).Figure 1
f
1
(
z
,
n
),  n∈0,1,…,10, in Example 1, where Γ0~Nμ=0;σ2=1;  a=10/9 (a) and a=9/10 (b).
(a)
(b)


### 3.2. Case I.2:A Is a Random Variable

In order to emphasize the deterministic nature of the initial conditionΓ0, we recast (48) by using the lowercase letter γ0:(55)Zn=Anγ0,n=0,1,2,….Letn≥0 be an arbitrary and fixed integer and denote Z=Zn=Anγ0. The RV Z represents power transformation of RV A; that is, Z can be written as Z=kAn. By applying Proposition 1 to Y=Z,  k=γ0, and X=A, one obtains the 1-PDF f1(z,n). For the sake of clarity, we do not provide the corresponding explicit expression for f1(z,n) since it just consists of substituting the previous identification. Below, we show an illustrative example.Example 2.
Let us assume thatA has a uniform distribution on the interval [-2,4], A~Un[-2,4], and γ0=1>0. Therefore, according to Proposition 1, the 1-PDF of {Zn:n≥0} is given by(56)
f
1
z
,
n
=
δ
z
-
1
,
n
=
0
,
z
∈
R
,
1
6
,
n
=
1
,
-
2
≤
z
≤
4
,
1
3
n
z
1
-
n
/
n
,
n
even
,
0
<
z
≤
2
n
,
1
6
n
z
1
-
n
/
n
,
n
even
,
2
n
<
z
≤
4
n
,
1
6
n
z
1
-
n
/
n
,
n
odd
,
n
≥
3
,
-
2
n
≤
z
<
0
,
0
<
z
≤
4
n
. It can be checked that f1(z,n) is a PDF for each n=0,1,2,…. Figure 2 shows f1(z,n) at different values of n.
This example exhibits a different behaviour of the1-PDF of the solution of (55) depending on whether n is odd or even.Figure 2
f
1
(
z
,
n
) in Example 2 at different values of n.


### 3.3. Case I.3:(Γ0,A) Is a Random Vector

Throughout this case, the joint PDF of the random vector(Γ0,A) will be denoted by fΓ0,A(γ0,a). Let n≥0 be an arbitrary and fixed integer and denote Z=Zn=AnΓ0. To compute the PDF of Z, first we will determine the joint PDF of the RVs Z and A by applying [20, Theorem 4] to the two-dimensional RV Y=r(X) with(57)X=Γ0A,Y=Y1Y2=r1Γ0,Ar2Γ0,A=AnΓ0A.From (57), the inverse transformation of r(X), X=r-1(Y)=s(Y), takes the form(58)
X
=
Γ
0
A
=
s
1
Y
1
,
Y
2
s
2
Y
1
,
Y
2
=
Y
1
Y
2
n
Y
2
. Taking into account the fact that ∂s2(y1,y2)/∂y1=0, the involved Jacobian simplifies to |J2|=1/|y2|n>0. Therefore, the joint PDF fY(y) is given by(59)fY1,Y2y1,y2=1y2nfΓ0,Ay1y2n,y2. Going back to the original RVs, that is, Z=AnΓ0=Y1 and A=Y2, one gets(60)fZ,Az,a=1anfΓ0,Azan,a,n=0,1,2,….Finally, considering the marginal density function of Z in (60), the 1-PDF of {Zn;n≥0} is given by(61)f1z,n=∫a1a21anfΓ0,Azan,ada,n=0,1,2,…,z∈R.Example 3.
Let(Γ0,A) be a random vector and let us assume that its PDF is given by(62)fΓ0,Aγ0,a=4γ0a,if0<γ0,a<1.By (61)-(62), the following 1-PDF of {Zn:n≥0} is obtained:(63)
f
1
z
,
n
=
4
z
∫
0
1
a
d
a
=
2
z
,
n
=
0
,
0
<
z
<
1
,
4
z
∫
z
1
1
a
d
a
=
-
4
z
ln
⁡
z
,
n
=
1
,
0
<
z
<
1
,
4
z
∫
z
n
1
1
a
2
n
-
1
d
a
=
2
z
1
-
z
2
1
-
n
n
1
-
n
,
n
=
2,3
,
…
,
0
<
z
<
1
. Figure 3 shows f1(z,n) for n∈0,1,2,…,10. From n≥2, we observe that the density of probability accumulates around z=0, which is in agreement with the asymptotic behaviour of the solution which tends to zero as n→∞.Figure 3
f
1
(
z
,
n
),  n∈{0,1,2,…,10}, in Example 3, where (Γ0,A) is a random vector, whose PDF is given by (62).


## 3.1. Case I.1:Γ0 Is a Random Variable

For the sake of clarity in the presentation, we rewrite solution (48) by using the lowercase letter a in order to indicate the deterministic character of parameter A:(49)Z0=Γ0,Zn=anΓ0,n=1,2,….Letn≥0 be an arbitrary and fixed integer. Let us assume a≠0 and denote Z=Zn=anΓ0. By applying [20, Proposition 2] to(50)Y=Z,X=Γ0,α=an≠0,β=0, one gets the 1-PDF(51)f1z,n=1anfΓ0zan,n=0,1,2,…,z∈R.Note that ifa=0, from (49), it follows that Zn=0 w.p. 1 for each n≥1 and Z0=Γ0. So, the 1-PDF for the trivial case, a=0, can be written as(52)f1z,n=fΓ0z,n=0,δz,n=1,2,…,z∈R.In order to facilitate the comparison of the1-PDF of the solution of problem (49) against its continuous counterpart provided in [20], in the following example, Γ0 is assumed to be a Gaussian RV. Note that we are going to consider a standard distribution, although the method is also able to be applied to nonstandard distributions.Example 1.
Let us assumea≠0 and consider Γ0 a Gaussian distribution, Γ0~Nμ;σ2. Hence, applying (51), the 1-PDF of {Zn:n≥0} is given by(53)f1z,n=1an2πσ2e-z/an-μ2/2σ2,n=0,1,2,…,z∈R. It can be checked that f1(z,n) is a PDF for each n≥0. Figure 1 shows f1(z,n) for n∈0,…,10, in the particular case that Γ0~N(0;1): a=10/9 (a) and a=9/10 (b). Note the different behavior of 1-PDF depending on the modulus of the parameter a. This is in agreement with the expectation and variance of the solution which are given, respectively, by(54)EZn=0,VZn=a2nσ2,n=0,1,2,…. Indeed, in Figure 1, we observe that, for each n, the 1-PDF is symmetric about z=0, whereas, in the case that a=10/9>1 (a=9/10∈-1,1), it becomes flat (sharp) as n increases. This means that its variability around zero, that is, the variance, tends to infinity (zero).Figure 1
f
1
(
z
,
n
),  n∈0,1,…,10, in Example 1, where Γ0~Nμ=0;σ2=1;  a=10/9 (a) and a=9/10 (b).
(a)
(b)


## 3.2. Case I.2:A Is a Random Variable

In order to emphasize the deterministic nature of the initial conditionΓ0, we recast (48) by using the lowercase letter γ0:(55)Zn=Anγ0,n=0,1,2,….Letn≥0 be an arbitrary and fixed integer and denote Z=Zn=Anγ0. The RV Z represents power transformation of RV A; that is, Z can be written as Z=kAn. By applying Proposition 1 to Y=Z,  k=γ0, and X=A, one obtains the 1-PDF f1(z,n). For the sake of clarity, we do not provide the corresponding explicit expression for f1(z,n) since it just consists of substituting the previous identification. Below, we show an illustrative example.Example 2.
Let us assume thatA has a uniform distribution on the interval [-2,4], A~Un[-2,4], and γ0=1>0. Therefore, according to Proposition 1, the 1-PDF of {Zn:n≥0} is given by(56)
f
1
z
,
n
=
δ
z
-
1
,
n
=
0
,
z
∈
R
,
1
6
,
n
=
1
,
-
2
≤
z
≤
4
,
1
3
n
z
1
-
n
/
n
,
n
even
,
0
<
z
≤
2
n
,
1
6
n
z
1
-
n
/
n
,
n
even
,
2
n
<
z
≤
4
n
,
1
6
n
z
1
-
n
/
n
,
n
odd
,
n
≥
3
,
-
2
n
≤
z
<
0
,
0
<
z
≤
4
n
. It can be checked that f1(z,n) is a PDF for each n=0,1,2,…. Figure 2 shows f1(z,n) at different values of n.
This example exhibits a different behaviour of the1-PDF of the solution of (55) depending on whether n is odd or even.Figure 2
f
1
(
z
,
n
) in Example 2 at different values of n.


## 3.3. Case I.3:(Γ0,A) Is a Random Vector

Throughout this case, the joint PDF of the random vector(Γ0,A) will be denoted by fΓ0,A(γ0,a). Let n≥0 be an arbitrary and fixed integer and denote Z=Zn=AnΓ0. To compute the PDF of Z, first we will determine the joint PDF of the RVs Z and A by applying [20, Theorem 4] to the two-dimensional RV Y=r(X) with(57)X=Γ0A,Y=Y1Y2=r1Γ0,Ar2Γ0,A=AnΓ0A.From (57), the inverse transformation of r(X), X=r-1(Y)=s(Y), takes the form(58)
X
=
Γ
0
A
=
s
1
Y
1
,
Y
2
s
2
Y
1
,
Y
2
=
Y
1
Y
2
n
Y
2
. Taking into account the fact that ∂s2(y1,y2)/∂y1=0, the involved Jacobian simplifies to |J2|=1/|y2|n>0. Therefore, the joint PDF fY(y) is given by(59)fY1,Y2y1,y2=1y2nfΓ0,Ay1y2n,y2. Going back to the original RVs, that is, Z=AnΓ0=Y1 and A=Y2, one gets(60)fZ,Az,a=1anfΓ0,Azan,a,n=0,1,2,….Finally, considering the marginal density function of Z in (60), the 1-PDF of {Zn;n≥0} is given by(61)f1z,n=∫a1a21anfΓ0,Azan,ada,n=0,1,2,…,z∈R.Example 3.
Let(Γ0,A) be a random vector and let us assume that its PDF is given by(62)fΓ0,Aγ0,a=4γ0a,if0<γ0,a<1.By (61)-(62), the following 1-PDF of {Zn:n≥0} is obtained:(63)
f
1
z
,
n
=
4
z
∫
0
1
a
d
a
=
2
z
,
n
=
0
,
0
<
z
<
1
,
4
z
∫
z
1
1
a
d
a
=
-
4
z
ln
⁡
z
,
n
=
1
,
0
<
z
<
1
,
4
z
∫
z
n
1
1
a
2
n
-
1
d
a
=
2
z
1
-
z
2
1
-
n
n
1
-
n
,
n
=
2,3
,
…
,
0
<
z
<
1
. Figure 3 shows f1(z,n) for n∈0,1,2,…,10. From n≥2, we observe that the density of probability accumulates around z=0, which is in agreement with the asymptotic behaviour of the solution which tends to zero as n→∞.Figure 3
f
1
(
z
,
n
),  n∈{0,1,2,…,10}, in Example 3, where (Γ0,A) is a random vector, whose PDF is given by (62).


## 4. Case Study: Nonhomogeneous Discrete Initial Value Problem (II)

In this section, we deal with the computation of the1-PDF {f1(z,n):n≥0} of the solution discrete SP {Zn:n≥0} of the nonhomogeneous discrete initial value problem (II). This will be done for every one of the cases considered in Table 1. Now, the solution {Zn:n≥0} has the following form:(64)Zn=Γ0+nB,n=0,1,2,….As we did in Section 3 and for the sake of clarity in the presentation, we will recast the input parameters in (64) by lowercase letters when they indicate deterministic quantities.


### 4.1. Case II.1:Γ0 Is a Random Variable

In this case, solution (64) takes the form(65)Zn=Γ0+nb,n=0,1,2,…. Let n≥0 be an arbitrary and fixed integer and denote Z=Zn=Γ0+nb. By applying [20, Proposition 2] to(66)Y=Z,X=Γ0,α=1≠0,β=nb,one obtains the 1-PDF(67)f1z,n=fΓ0z-nb,n=0,1,2,…,z∈R.Example 4.
LetΓ0 be a gamma RV of parameters α,β>0, Γ0~Gaα;β. Then, by (67), the 1-PDF of {Zn:n≥0} reads(68)f1z,n=1βαΓαz-nbα-1e-z-nb/β,n=0,1,2,…,nb≤z<+∞, where Γ(α) means the classical gamma function. Notice that, for each n=0,1,2,…, the domain of z follows from the corresponding domain of a gamma distribution. In Figure 4, the 1-PDF f1(z,n) is plotted at different values of n considering Γ0~Ga(2;1) and b=1/2. In this case, from the plot of f1(z,n), one observes that the expectation E[Zn] increases as n does. It is straightforward to check that the expectation lies on the straight line 1/2n+2, whereas the variance takes the constant value 2.Figure 4
f
1
(
z
,
n
),  n∈{0,1,2,…,10}, in Example 4, where Γ0~Ga(2;1) and b=1/2.


### 4.2. Case II.2:B Is a Random Variable

In this case, the solution discrete stochastic process (64) takes the form(69)Zn=γ0+nB,n=0,1,2,…. For n=0,  Z0=γ0 and the PDF is δ(z-γ0). If n≥1 is an arbitrary and fixed integer, denoting Z=Zn=γ0+nB, the PDF is obtained by applying [20, Proposition 2] to(70)Y=Z,X=B,α=n≠0,β=γ0.This leads to the 1-PDF(71)f1z,n=δz-γ0,n=01nfBz-γ0n,n=1,2,…,z∈R.Example 5.
LetB be a RV with χ2-distribution with ν degrees of freedom, B~χ2(ν), ν>0. Then, by (71), the 1-PDF of {Zn:n≥0} writes(72)f1z,n=δz-γ0,n=0,z∈R,2nν12ν/21Γν/2z-γ0ν/2-1e-1/2z-γ0/n2,n=1,2,…,γ0≤z<∞.For each n≥1, the domain of z has been determined considering the domain of χ2-distribution with ν degrees of freedom. For the sake of clarity, Figure 5 shows a 2D plot (a) and a 3D plot (b) of f1(z,n) at different values of n in the particular case that B~χ23 and γ0=1.Figure 5
f
1
(
z
,
n
) in Example 5, where B~χ23 and γ0=1, at different values of n. (a) 2D plot for n∈{1,2,4,10}. (b) 3D plot for n∈{1,2,…,10}.
(a)
(b)


### 4.3. Case II.3:Γ0,B Is a Random Vector

In accordance with the notation previously introduced,fΓ0,B(γ0,b) stands for the joint PDF of the random vector (Γ0,B). Let n≥0 be an arbitrary and fixed integer and denote Z=Zn=Γ0+nB. To compute the PDF of Z, first we will determine the joint PDF of the RVs Z and B by applying [20, Theorem 1] to the two-dimensional RV Y=r(X) with(73)X=Γ0B,Y=Y1Y2=r1Γ0,Br2Γ0,B=Γ0+nBB.From (73), the inverse transformation of r(X): X=r-1(Y)=s(Y) takes the form(74)X=Γ0B=s1Y1,Y2s2Y1,Y2=Y1-nY2Y2. By [20, Theorem 4] and taking into account the fact that |J2|=1≠0, the joint PDF fY(y) is given by(75)fY1,Y2y1,y2=fΓ0,By1-ny2,y2. Going back to the original RVs, that is, Z=Γ0+nB=Y1 and B=Y2, one gets(76)fZ,Bz,b=fΓ0,Bz-nb,b. Finally, considering the marginal density function of Z, one gets the 1-PDF of {Zn:n≥0}:(77)f1z,n=∫b1b2fΓ0,Bz-nb,bdb,n=0,1,2,…,z∈R.Example 6.
Let us consider the random vector(Γ0,B) whose joint PDF is defined by(78)fΓ0,Bγ0,b=14+14γ03b-14γ0b3,if-1<γ0,b<1.By (77), the 1-PDF of {Zn:n≥0} is(79)f1z,n=∫max⁡-1,z-1/nmin⁡1,z+1/nfΓ0,Bz-nb,bdb,–n+1≤z≤n+1,n=0,1,2,…. Computing the above integral, one gets(80)f1z,n=f1az,n,ifn=0,1,2,…,-n+1≤z≤-n+1,f1bz,n,ifn=0,1,2,…,-n+1≤z≤n-1,f1cz,n,ifn=0,1,2,…,n-1≤z≤n+1,
where(81)f1az,n=∫-1z+1/nfΓ0,Bz-nb,bdb=180n4-4n7-15n6z+20n31+z--4+z1+z4+n54-20z2+5n44+z-2z3+n2-4-5z+z5,f1bz,n=∫z-1/nz+1/nfΓ0,Bz-nb,bdb=1-n2+5n3+5z210n4,f1cz,n=∫z-1/n1fΓ0,Bz-nb,bdb=180n4-4n7-20n3-1+z+15n6z+-1+z44+z+n54-20z2+5n44-z+2z3-n24-5z+z5.
Figure6 shows two equivalent plots of f1(z,n) given by (80)-(81). It is straightforward to check that the expectation and variance of the solution Zn are given by(82)EZn=0,VZn=131+n2,n=0,1,2,…, respectively. Notice that the values of the expectation and variance obtained from expressions (82) agree with the plots shown in Figure 6, where one observes that the 1-PDF f1(z,n) is, for each n≥0, symmetric about z=0 and its support increases as n tends to infinity in such a way that the 1-PDF’s shape becomes flattened. Then, the variability about zero, in this case the variance, increases as n does.Figure 6
f
1
(
z
,
n
) in Example 6 at different values of n, where (Γ0,B) is a random vector whose PDF is given by (78). (a) 2D plot for n∈1,3,5,10. (b) 3D plot for n∈{1,2,…,10}.
(a)
(b)


## 4.1. Case II.1:Γ0 Is a Random Variable

In this case, solution (64) takes the form(65)Zn=Γ0+nb,n=0,1,2,…. Let n≥0 be an arbitrary and fixed integer and denote Z=Zn=Γ0+nb. By applying [20, Proposition 2] to(66)Y=Z,X=Γ0,α=1≠0,β=nb,one obtains the 1-PDF(67)f1z,n=fΓ0z-nb,n=0,1,2,…,z∈R.Example 4.
LetΓ0 be a gamma RV of parameters α,β>0, Γ0~Gaα;β. Then, by (67), the 1-PDF of {Zn:n≥0} reads(68)f1z,n=1βαΓαz-nbα-1e-z-nb/β,n=0,1,2,…,nb≤z<+∞, where Γ(α) means the classical gamma function. Notice that, for each n=0,1,2,…, the domain of z follows from the corresponding domain of a gamma distribution. In Figure 4, the 1-PDF f1(z,n) is plotted at different values of n considering Γ0~Ga(2;1) and b=1/2. In this case, from the plot of f1(z,n), one observes that the expectation E[Zn] increases as n does. It is straightforward to check that the expectation lies on the straight line 1/2n+2, whereas the variance takes the constant value 2.Figure 4
f
1
(
z
,
n
),  n∈{0,1,2,…,10}, in Example 4, where Γ0~Ga(2;1) and b=1/2.


## 4.2. Case II.2:B Is a Random Variable

In this case, the solution discrete stochastic process (64) takes the form(69)Zn=γ0+nB,n=0,1,2,…. For n=0,  Z0=γ0 and the PDF is δ(z-γ0). If n≥1 is an arbitrary and fixed integer, denoting Z=Zn=γ0+nB, the PDF is obtained by applying [20, Proposition 2] to(70)Y=Z,X=B,α=n≠0,β=γ0.This leads to the 1-PDF(71)f1z,n=δz-γ0,n=01nfBz-γ0n,n=1,2,…,z∈R.Example 5.
LetB be a RV with χ2-distribution with ν degrees of freedom, B~χ2(ν), ν>0. Then, by (71), the 1-PDF of {Zn:n≥0} writes(72)f1z,n=δz-γ0,n=0,z∈R,2nν12ν/21Γν/2z-γ0ν/2-1e-1/2z-γ0/n2,n=1,2,…,γ0≤z<∞.For each n≥1, the domain of z has been determined considering the domain of χ2-distribution with ν degrees of freedom. For the sake of clarity, Figure 5 shows a 2D plot (a) and a 3D plot (b) of f1(z,n) at different values of n in the particular case that B~χ23 and γ0=1.Figure 5
f
1
(
z
,
n
) in Example 5, where B~χ23 and γ0=1, at different values of n. (a) 2D plot for n∈{1,2,4,10}. (b) 3D plot for n∈{1,2,…,10}.
(a)
(b)


## 4.3. Case II.3:Γ0,B Is a Random Vector

In accordance with the notation previously introduced,fΓ0,B(γ0,b) stands for the joint PDF of the random vector (Γ0,B). Let n≥0 be an arbitrary and fixed integer and denote Z=Zn=Γ0+nB. To compute the PDF of Z, first we will determine the joint PDF of the RVs Z and B by applying [20, Theorem 1] to the two-dimensional RV Y=r(X) with(73)X=Γ0B,Y=Y1Y2=r1Γ0,Br2Γ0,B=Γ0+nBB.From (73), the inverse transformation of r(X): X=r-1(Y)=s(Y) takes the form(74)X=Γ0B=s1Y1,Y2s2Y1,Y2=Y1-nY2Y2. By [20, Theorem 4] and taking into account the fact that |J2|=1≠0, the joint PDF fY(y) is given by(75)fY1,Y2y1,y2=fΓ0,By1-ny2,y2. Going back to the original RVs, that is, Z=Γ0+nB=Y1 and B=Y2, one gets(76)fZ,Bz,b=fΓ0,Bz-nb,b. Finally, considering the marginal density function of Z, one gets the 1-PDF of {Zn:n≥0}:(77)f1z,n=∫b1b2fΓ0,Bz-nb,bdb,n=0,1,2,…,z∈R.Example 6.
Let us consider the random vector(Γ0,B) whose joint PDF is defined by(78)fΓ0,Bγ0,b=14+14γ03b-14γ0b3,if-1<γ0,b<1.By (77), the 1-PDF of {Zn:n≥0} is(79)f1z,n=∫max⁡-1,z-1/nmin⁡1,z+1/nfΓ0,Bz-nb,bdb,–n+1≤z≤n+1,n=0,1,2,…. Computing the above integral, one gets(80)f1z,n=f1az,n,ifn=0,1,2,…,-n+1≤z≤-n+1,f1bz,n,ifn=0,1,2,…,-n+1≤z≤n-1,f1cz,n,ifn=0,1,2,…,n-1≤z≤n+1,
where(81)f1az,n=∫-1z+1/nfΓ0,Bz-nb,bdb=180n4-4n7-15n6z+20n31+z--4+z1+z4+n54-20z2+5n44+z-2z3+n2-4-5z+z5,f1bz,n=∫z-1/nz+1/nfΓ0,Bz-nb,bdb=1-n2+5n3+5z210n4,f1cz,n=∫z-1/n1fΓ0,Bz-nb,bdb=180n4-4n7-20n3-1+z+15n6z+-1+z44+z+n54-20z2+5n44-z+2z3-n24-5z+z5.
Figure6 shows two equivalent plots of f1(z,n) given by (80)-(81). It is straightforward to check that the expectation and variance of the solution Zn are given by(82)EZn=0,VZn=131+n2,n=0,1,2,…, respectively. Notice that the values of the expectation and variance obtained from expressions (82) agree with the plots shown in Figure 6, where one observes that the 1-PDF f1(z,n) is, for each n≥0, symmetric about z=0 and its support increases as n tends to infinity in such a way that the 1-PDF’s shape becomes flattened. Then, the variability about zero, in this case the variance, increases as n does.Figure 6
f
1
(
z
,
n
) in Example 6 at different values of n, where (Γ0,B) is a random vector whose PDF is given by (78). (a) 2D plot for n∈1,3,5,10. (b) 3D plot for n∈{1,2,…,10}.
(a)
(b)


## 5. Case Study: Nonhomogeneous Discrete Initial Value Problem (III)

This section is addressed to determine the1-PDFs {f1(z,n):n≥0} of the solution SPs {Zn:n≥0} of problem (III) in Cases III.1–III.7 collected in Table 1. Now, the solution {Zn:n≥0} has the following form:(83)Z0=Γ0,Zn=Γ0-B1-AAn+B1-A,n=1,2,…,which is well defined due to the hypothesis Pω∈Ω:A(ω)=1=0.Analogously to the previous sections, for the sake of clarity in the presentation, we will rewrite each one of the involved parameters in (83) by lowercase letters when it denotes a deterministic quantity.


### 5.1. Case III.1:Γ0 Is a Random Variable

In this case, ifa≠0, solution (83) takes the form(84)Zn=anΓ0+b1-a1-an,n=0,1,2,….Let n≥0 be an arbitrary and fixed integer and denote Z=Zn=anΓ0+b(1-an)/(1-a). The application of [20, Proposition 2] to(85)Y=Z,X=Γ0,α=an≠0,β=b1-a1-anpermits computing the 1-PDF. This yields(86)f1z,n=1anfΓ0z1-a-b1-anan1-a,n=0,1,2,…,z∈R.For the trivial case a=0, solution (83) takes the form(87)Z0=Γ0,Zn=b,n=1,2,…, and hence the 1-PDF is given by(88)f1z,n=fΓ0z,n=0,δz-b,n=1,2,…,z∈R.Remark 7.
Notice that expression (51) obtained in Case I.1 is a particular case of (86) taking b=0. Similarly, if the parameter a tends to 1 in (86), one gets formula (67) of Case II.1.Example 7.
LetΓ0 be an exponential RV of parameter λ>0, Γ0~Exp(λ). Then, by (86), the 1-PDF of {Zn:n≥0} writes(89)
f
1
z
,
n
=
-
λ
a
n
e
-
λ
z
/
a
n
-
b
1
-
a
n
/
a
n
1
-
a
if
n
odd
,
z
<
b
1
-
a
n
1
-
a
,
λ
a
n
e
-
λ
z
/
a
n
-
b
1
-
a
n
/
a
n
1
-
a
if
n
even
,
z
>
b
1
-
a
n
1
-
a
,
a
<
0
,
(90)
f
1
z
,
n
=
λ
a
n
e
-
λ
z
/
a
n
-
b
1
-
a
n
/
a
n
1
-
a
,
z
>
b
1
-
a
n
1
-
a
,
a
>
0
, where the domain has been determined taking into account the domain of an exponential RV. Figure 7 shows f1(z,n) at different values of n depending on whether n is odd (a) or even (b) for λ=1,  a=-11/10, and b=1.Figure 7
f
1
(
z
,
n
) in Example 7 at different values of n depending on whether n is odd (a) or even (b), where Γ0~Exp(λ=1),  a=-11/10, and b=1.
(a)
(b)Below, we show an example with the aim of illustrating that once the1-PDF has been computed, important statistical moments of the solution SP, such as the expectation and variance, can be determined straightforwardly.Example 8.
Within the context of Example7 and assuming for illustrative purposes that, for instance,  a<0, the statistical moment of order n of Zn can be determined directly using expression (89) of f1(z,n) in the following way: (91)
m
Z
n
,
k
=
E
Z
n
k
=
∫
-
∞
b
1
-
a
n
/
1
-
a
z
k
f
1
z
,
n
d
z
,
if
n
is odd
,
∫
b
1
-
a
n
/
1
-
a
+
∞
z
k
f
1
z
,
n
d
z
,
if
n
is even
,
k
=
0,1
,
2
,
…
.
For example, takingλ=1, a=-11/10, and b=1, the mean and the variance of Zn are given by(92)
E
Z
n
=
m
Z
n
,
1
=
1
21
10
-
1
0
-
n
1
1
n
+
1
,
if
n
is odd
,
1
21
10
+
1
0
-
n
1
1
n
+
1
,
if
n
is even
,
V
Z
n
=
m
Z
n
,
2
-
m
Z
n
,
1
2
=
121
100
n
,
n
=
0,1
,
2
,
…
.Note that the expression of the variance does not depend on whether n is even or odd.
In Figure8, the expectation and variance of Zn have been plotted using expressions (92) to carry out computations.
In addition to computing the mean and the variance, further significant information related to the solution SP can be computed from the1-PDF, such as the probability of specific sets in which we could be interested. For example, the probability that the solution varies between the values v1=2 and v2=3 is given by (93)P2≤Zn≤3=∫23f1z,ndz=0ifnis odd,n≤11,e1/21-10-1n+11-e2/212n+45/11n+5-1nifnis odd,13≤n≤17,e1/21-1011-n10n+11ne310/11n-e5/11n2n+1ifnis odd,n≥19,e1/21-51910/11n+2e310/11n-e5/11n2n+1ifnis even.Figure 8
Mean (a) and variance (b) ofZn in Example 8, where Γ0~Exp(λ=1),  a=-11/10, and b=1.
(a)
(b)


### 5.2. Case III.2:B Is a Random Variable

Let us assumea≠-1. In this case, the discrete solution SP (83) takes the form(94)Z0=γ0,Zn=1-an1-aB+γ0an,n=1,2,…. Let n≥1 be an arbitrary and fixed integer and denote Z=Zn. By applying [20, Proposition 2] to(95)Z=Y,X=B,α=1-an1-a≠0,β=γ0an,and taking into account the fact that Z0=γ0, one obtains the 1-PDF(96)f1z,n=δz-γ0,n=0,1-a1-anfB1-a1-anz-γ0an,n=1,2,…,z∈R.Ifa=-1, then the discrete solution SP (83) is(97)Zn=1--1n2B+γ0-1n=γ0ifnis even,B-γ0ifnis odd.Then, forn≥0, an arbitrary but fixed integer, applying [20, Proposition 2] to (97), one obtains the 1-PDF of Zn defined by (97): (98)f1z,n=δz-γ0,neven,fBz+γ0,nodd,z∈R.Remark 9.
Notice that expression (71) in Case II.2 can be obtained assuming that parameter a tends to 1 in (96).Example 9.
LetB be a RV having a gamma distribution of parameters α,β>0; that is, B~Ga(α;β). Let us fix a≠-1. Then, by (96), the 1-PDF of {Zn:n≥0} reads(99)
f
1
z
,
n
=
δ
z
-
γ
0
if
n
=
0
,
z
∈
R
,
f
1
a
z
,
n
,
if
1
-
a
1
-
a
n
>
0
,
n
=
1,2
,
…
,
z
≥
γ
0
a
n
,
f
1
b
z
,
n
,
if
1
-
a
1
-
a
n
<
0
,
n
=
1,2
,
…
,
z
≤
γ
0
a
n
,where(100)f1az,n=1-a1-an1βαΓα1-a1-anz-γ0anα-1e1/βa-1/1-anz-γ0an,f1bz,n=a-11-an1βαΓα1-a1-anz-γ0anα-1e1/βa-1/1-anz-γ0an. Taking into account the fact that the domain of a gamma RV is (0,∞), one deduces the domain of z specified in (99) for each n≥1. Figure 9 shows, f1(z,n) at different values of n assuming that B~Ga(2;4), a=-1/2, and γ0=1. Note that, in accordance with (99), f1z,0=δz-1, -∞<z<∞.Figure 9
2D and 3D plots off1(z,n) in Example 9, where B~Ga(2;4), a=-1/2, and γ0=1. (a) n∈{1,2,3,4}. (b) n∈{1,2,…,10}.
(a)
(b)


### 5.3. Case III.3:A Is a Random Variable

So far, we have taken advantage of RVT method to compute the1-PDF of the solution of problem (2). The success of this approach has relied on the capability to find out an exact expression for the inverse transformation of the mapping that determines the solution in terms of the random inputs. However, in many practical cases, such expression can just be found in an approximate form rather than in an exact form. Under these circumstances, the RVT method can still be very useful. In fact, as it was shown in the analysis of Case III.3 of [20], the application of the Lagrange-Bürmann theorem [27] together with the RVT technique permits determining reliable approximations of the 1-PDF of the solution SP of the continuous counterpart of problem (2). Below, we provide an illustrative example dealing with a particular case of (83) where just the parameter A is assumed to be random. To avoid repetitions, we omit the theoretical development which can be found in [20].Example 10.
Throughout this example, we will use the notation introduced in [20]. Let A be a beta RV of parameters α=2, β=3: A~Be(2;3) and γ0=1 and b=1. In Figure 10, the approximation of f1(z,n) is shown at different values of n≥1. This plot has been made by considering expressions [20, eqs. (78)–(80)] with k=1, being A1=[0,1] because of monotony of r(A). To carry out computations, A1 has been divided into 7 subintervals in agreement with the process described in [20]. In each subinterval, an approximation of degree Nj=2 has been considered. If n=0, f1(z,0)=δ(z-γ0).Figure 10
2D and 3D plots off1(z,n) in Example 10, where A~Be(2;3),  γ0=1, and b=1. (a) n∈{0,1,2,3}. (b) n∈{0,1,2,…,10}.
(a)
(b)


### 5.4. Case III.4:Γ0,B Is a Random Vector

Let us denote byfΓ0,B(γ0,b) the joint PDF of the random vector (Γ0,B).Ifa≠0, the discrete solution SP (83) takes the form(101)Zn=1-an1-aB+anΓ0,n=0,1,2,….Let us fixn:n≥0 and denote Z=Zn=(1-an)/(1-a)B+anΓ0. To compute the 1-PDF of Z, first we will determine the joint PDF of the RVs Z and B by applying [20, Theorem 4] to the two-dimensional RV Y=r(X) being(102)
X
=
Γ
0
B
,
Y
=
Y
1
Y
2
=
r
1
Γ
0
,
B
r
2
Γ
0
,
B
=
1
-
a
n
1
-
a
B
+
a
n
Γ
0
B
.From (102), the inverse transformation of r(X),  X=r-1(Y)=s(Y), takes the form(103)
X
=
Γ
0
B
=
s
1
Y
1
,
Y
2
s
2
Y
1
,
Y
2
=
Y
1
-
1
-
a
n
1
-
a
Y
2
1
a
n
Y
2
. By [20, Theorem 4] and taking into account the fact that |J2|=1/|a|n≠0, the joint PDF fY(y) is given by(104)fY1,Y2y1,y2=1anfΓ0,By1-1-an1-ay21an,y2. Going back to the original RVs, that is, Z=1-an/1-aB+anΓ0=Y1 and B=Y2, one gets(105)fZ,Bz,b=1anfΓ0,Bz-1-an1-ab1an,b.Finally, considering the marginal PDF of Z in (105), one obtains the 1-PDF of Z for a≠0:(106)f1z,n=1an∫b1b2fΓ0,Bz-1-an1-ab1an,bdb,n=0,1,2,…,z∈R.For the trivial casea=0, solution (101) takes the form(107)Z0=Γ0Zn=B,n=1,2,…, and hence the 1-PDF of Z is given by the marginal PDFs of Γ0 and B:(108)
f
1
z
,
n
=
f
Γ
0
γ
0
=
∫
b
1
b
2
f
Γ
0
,
B
γ
0
,
b
d
b
,
n
=
0
,
f
B
b
=
∫
γ
0,1
γ
0,2
f
Γ
0
,
B
γ
0
,
b
d
γ
0
,
n
=
1,2
,
…
,
z
∈
R
.Remark 11.
Notice that expression (77) in Case II.3 can be obtained assuming that parameter a tends to 1 in (106).Example 11.
Let us assume that the random vectorη=(Γ0,B)T~N(μη,Ση) has a bidimensional Gaussian distribution whose mean vector and covariance matrix are given by(109)μη=μΓ0μB,Ση=σΓ02ρΓ0,BσΓ0σBρΓ0,BσZ0σBσB2,respectively. In (109), ρΓ0,B denotes the correlation coefficient of RVs Γ0 and B. Let us assume a≠0. By (106), the 1-PDF of {Zn:n≥0} reads as follows:(110)f1z,n=1an12πdet⁡Ση∫-∞+∞e-1/2ζ-μηTΣη-1ζ-μηdb,n=0,1,2,…,z∈R,being(111)
ζ
=
z
-
1
-
a
n
1
-
a
b
1
a
n
b
. To determine the domain of integration in (110), we have taken into account the fact that B is also a Gaussian RV. Figure 11 shows a graphical representation of f1(z,n) at different values of n assuming that μΓ0=1,  μB=0,  σΓ0=0.1,  σB=0.1,  ρΓ0,B=0.5, and a=-11/10.Figure 11
f
1
(
z
,
n
) in Example 11 at different values of n. The input parameters are assumed to be η=(Γ0,B)T~N(μη,Ση), where the mean vector and the variance-covariance matrix are defined by (109), with the values μΓ0=1,  μB=0,  σΓ0=0.1,  σB=0.1,  ρΓ0,B=0.5, and a=-11/10.


### 5.5. Case III.5:Γ0,A Is a Random Vector

In this case, solution (83) takes the form(112)Zn=AnΓ0+b1-A1-An,n=0,1,2,….Let us denote byfΓ0,A(γ0,a) the joint PDF of the random vector (Γ0,A). Let us fix n:n≥0 and denote Z=Zn=AnΓ0+b/(1-A)(1-An). To compute the 1-PDF of Z, first we will determine the joint PDF of the RVs Z and A by applying [20, Theorem 4] to the two-dimensional RV Y=r(X) with(113)
X
=
Γ
0
A
,
Y
=
Y
1
Y
2
=
r
1
Γ
0
,
A
r
2
Γ
0
,
A
=
A
n
Γ
0
+
b
1
-
A
1
-
A
n
A
.From (113), the inverse transformation of r(X), X=r-1(Y)=s(Y), takes the form(114)
X
=
Γ
0
A
=
s
1
Y
1
,
Y
2
s
2
Y
1
,
Y
2
=
Y
1
-
b
1
-
Y
2
1
-
Y
2
n
1
Y
2
n
Y
2
. Notice that Y2≠0 w.p. 1. By [20, Theorem 4] and taking into account the fact that |J2|=1/|y2|n≠0, the joint PDF fY(y) is given by(115)fY1,Y2y1,y2=1y2nfΓ0,Ay1-b1-y21-y2n1y2n,y2. Going back to the original RVs, that is, Z=AnΓ0+b/(1-A)(1-An)=Y1 and A=Y2, one gets(116)fZ,Az,a=1anfΓ0,Az-b1-a1-an1an,a.Finally, considering the marginal PDF of Z in (116), the 1-PDF of {Zn:n≥0} is given by(117)f1z,n=∫a1a21anfΓ0,Az-b1-a1-an1an,ada,n=0,1,2,…,z∈R.Example 12.
Let us consider the two-dimensional Gaussian vectorη=(Γ0,A)T~N(μη,Ση), whose mean and covariance matrix are given by(118)μη=μΓ0μA,Ση=σΓ02ρΓ0,AσΓ0σAρΓ0,AσΓ0σAσA2,respectively. In (118), ρΓ0,A denotes the correlation coefficient between RVs Γ0 and A. Notice that A has a Gaussian distribution; hence, a1=-∞ and a2=+∞. Then, according to (117), the 1-PDF of {Zn:n≥0} is given by(119)f1z,n=12πdet⁡Ση∫-∞+∞1ane-1/2ζ-μηTΣη-1ζ-μηda,n=0,1,2,…,z∈R, where(120)ζ=z-b1-a1-an1ana.
Figure12 shows a graphical representation of the 1-PDF f1(z,n) given by (117) at different values of n assuming that μΓ0=1,  μA=1.5,  σΓ0=0.1,  σA=0.1, ρΓ0,A=0.5, and b=1. For the sake of clarity, we have split the representation of f1(z,n) into two plots due to the significant differences in the vertical scales required depending on the values of n.Figure 12
f
1
(
z
,
n
) in Example 12 at different values of n. The input parameters are assumed to be η=(Γ0,A)T~N(μη,Ση), where the mean vector and the variance-covariance matrix are defined by (118), with the values μΓ0=1,  μA=1.5,  σΓ0=0.1,  σA=0.1,  ρΓ0,A=0.5, and b=1.


### 5.6. Case III.6:B,A Is a Random Vector

Solution (83) takes the form(121)Zn=Anγ0+1-An1-AB,n=0,1,2,….Let us denote byfB,A(b,a) the joint PDF of the random vector (B,A). Let us fix n:n≥1 and denote Z=Zn=Anγ0+(1-An)/(1-A)B. In order to compute the 1-PDF of Z, first we will determine the joint PDF of the RVs Z and A by applying [20, Theorem 4] to the two-dimensional RV Y=r(X) with(122)
X
=
B
A
,
Y
=
Y
1
Y
2
=
r
1
B
,
A
r
2
B
,
A
=
A
n
γ
0
+
1
-
A
n
1
-
A
B
A
.From (122), the inverse transformation of r(X), X=r-1(Y)=s(Y), takes the form(123)X=BA=s1Y1,Y2s2Y1,Y2=Y1-Y2nγ01-Y21-Y2nY2. By [20, Theorem 4] and taking into account the fact that |J2|=1-y2/1-(y2)n≠0, the joint PDF fY(y) is given by(124)fY1,Y2y1,y2=1-y21-y2nfB,Ay1-y2nγ01-y21-y2n,y2. Going back to the original RVs, that is, Z=Anγ0+1-An/1-AB=Y1 and A=Y2, one gets(125)fZ,Az,a=1-a1-anfB,Az-anγ01-a1-an,a,n=1,2,….Finally, considering the marginal PDF of Z in (125) and the case where n=0, which gives Z0=γ0, one gets the 1-PDF of {Zn:n≥0}:(126)f1z,n=δz-γ0,n=0,∫a1a21-a1-anfB,Az-anγ01-a1-an,ada,n=1,2,…,z∈R.Example 13.
So far, we have considered standard distributions in one or more dimensions to illustrate the obtained theoretical results. Now, we will assume that the joint PDF of the input parametersB and A is constructed by means of a copula transformation. Let us assume that B and A are uniform RVs defined on the interval ]0,1[; that is, B,A~Un(]0,1[). We transform these RVs by the Farlie-Gordon-Morgenstern copula [28], so that a two-dimensional RV (B,A) with joint PDF(127)fB,Ab,a=232-b-a+2ba,if0<b,a<1,is defined. This random vector satisfies the notion that the marginal distributions of fB,A(b,a) keep the one-dimensional distributions of A and B. Hereinafter, let us take γ0=0. Taking into account the fact that 0<a,b<1, by (126), for z and n previously fixed, one must calculate a such that (128)0<z1-a1-an<1.In general, the variable a cannot be determined in a closed form. Using Mathematica software to carry out computations numerically, in Figure 13, several representations of f1(z,n) have been plotted.Figure 13
2D and 3D plots off1(z,n) in Example 13 at different values of n, where γ0=0 and the PDF of the two-dimensional RV (B,A) is given by (127). (a) n∈{1,2,5,10}. (b) n∈{1,2,…,10}.
(a)
(b)


### 5.7. Case III.7:Γ0,B,A Is a Random Vector

In this last case, solution (83) takes the form(129)Zn=AnΓ0+B1-A1-An,n=0,1,2,….Let us denote byfΓ0,B,A(γ0,b,a) the joint PDF of the random vector (Γ0,B,A). Let us fix n:n≥0 and denote Z=Zn=AnΓ0+B/(1-A)(1-An). To determine the 1-PDF of Z, first we will calculate the joint PDF of the RVs Z, B, and A by applying [20, Theorem 4] to the three-dimensional RV Y=r(X) with(130)
X
=
Γ
0
B
A
,
Y
=
Y
1
Y
2
Y
3
=
r
1
Γ
0
,
B
,
A
r
2
Γ
0
,
B
,
A
r
3
Γ
0
,
B
,
A
=
A
n
Γ
0
+
B
1
-
A
1
-
A
n
B
A
.From (130), the inverse transformation of r(X),  X=r-1(Y)=s(Y), takes the form(131)X=Γ0BA=s1Y1,Y2,Y3s2Y1,Y2,Y3s3Y1,Y2,Y3=Y1-Y21-Y31-Y3n1Y3nY2Y3, where Y3≠0 w.p. 1. By [20, Theorem 4] and taking into account the fact that |J3|=1/|y3|n≠0, the joint PDF fY(y) is given by(132)fY1,Y2,Y3y1,y2,y3=1y3nfΓ0,B,Ay1-y21-y31-y3n1y3n,y2,y3.Going back to the original RVs, that is, Z=AnΓ0+B/1-A(1-An)=Y1,  B=Y2, and A=Y3, one gets(133)fZ,B,Az,b,a=1anfΓ0,B,Az-b1-a1-an1an,b,a.Finally, considering the marginal PDF of Z in (133), the 1-PDF of {Zn:n≥0} is given by(134)f1z,n=∫a1a2∫b1b21anfΓ0,B,Az-b1-a1-an1an,b,adbda,n=0,1,2,…,z∈R.Example 14.
Let us assume thatη=(Z0,B,A)T~N(μη,Ση); that is, η is a Gaussian random vector with mean μη=(μ1,μ2,μ3)T∈R3 and variance-covariance matrix Ση∈R3×3. By (134), the 1-PDF of Zn writes(135)f1z,n=12π2πdet⁡Ση∬-∞∞1ane-1/2ζ-μηTΣη-1ζ-μηdbda, where(136)ζ=z-b1-a1-an1anba. In order to illustrate the theoretical results previously established, let us fix the mean vector and the variance-covariance matrix as follows:(137)μη=112,Ση=110411141112.In this example, we do not make the 1-PDF, f1(z,n), explicit, since its expression is cumbersome. In Figure 14, we have plotted f1(z,n) at different values of n.Figure 14
f
1
(
z
,
n
),  n∈{0,1,2,3}, in Example 14. The input parameters (Γ0,B,andA) are assumed to be a Gaussian distribution with mean and variance-covariance matrix defined by (137).


## 5.1. Case III.1:Γ0 Is a Random Variable

In this case, ifa≠0, solution (83) takes the form(84)Zn=anΓ0+b1-a1-an,n=0,1,2,….Let n≥0 be an arbitrary and fixed integer and denote Z=Zn=anΓ0+b(1-an)/(1-a). The application of [20, Proposition 2] to(85)Y=Z,X=Γ0,α=an≠0,β=b1-a1-anpermits computing the 1-PDF. This yields(86)f1z,n=1anfΓ0z1-a-b1-anan1-a,n=0,1,2,…,z∈R.For the trivial case a=0, solution (83) takes the form(87)Z0=Γ0,Zn=b,n=1,2,…, and hence the 1-PDF is given by(88)f1z,n=fΓ0z,n=0,δz-b,n=1,2,…,z∈R.Remark 7.
Notice that expression (51) obtained in Case I.1 is a particular case of (86) taking b=0. Similarly, if the parameter a tends to 1 in (86), one gets formula (67) of Case II.1.Example 7.
LetΓ0 be an exponential RV of parameter λ>0, Γ0~Exp(λ). Then, by (86), the 1-PDF of {Zn:n≥0} writes(89)
f
1
z
,
n
=
-
λ
a
n
e
-
λ
z
/
a
n
-
b
1
-
a
n
/
a
n
1
-
a
if
n
odd
,
z
<
b
1
-
a
n
1
-
a
,
λ
a
n
e
-
λ
z
/
a
n
-
b
1
-
a
n
/
a
n
1
-
a
if
n
even
,
z
>
b
1
-
a
n
1
-
a
,
a
<
0
,
(90)
f
1
z
,
n
=
λ
a
n
e
-
λ
z
/
a
n
-
b
1
-
a
n
/
a
n
1
-
a
,
z
>
b
1
-
a
n
1
-
a
,
a
>
0
, where the domain has been determined taking into account the domain of an exponential RV. Figure 7 shows f1(z,n) at different values of n depending on whether n is odd (a) or even (b) for λ=1,  a=-11/10, and b=1.Figure 7
f
1
(
z
,
n
) in Example 7 at different values of n depending on whether n is odd (a) or even (b), where Γ0~Exp(λ=1),  a=-11/10, and b=1.
(a)
(b)Below, we show an example with the aim of illustrating that once the1-PDF has been computed, important statistical moments of the solution SP, such as the expectation and variance, can be determined straightforwardly.Example 8.
Within the context of Example7 and assuming for illustrative purposes that, for instance,  a<0, the statistical moment of order n of Zn can be determined directly using expression (89) of f1(z,n) in the following way: (91)
m
Z
n
,
k
=
E
Z
n
k
=
∫
-
∞
b
1
-
a
n
/
1
-
a
z
k
f
1
z
,
n
d
z
,
if
n
is odd
,
∫
b
1
-
a
n
/
1
-
a
+
∞
z
k
f
1
z
,
n
d
z
,
if
n
is even
,
k
=
0,1
,
2
,
…
.
For example, takingλ=1, a=-11/10, and b=1, the mean and the variance of Zn are given by(92)
E
Z
n
=
m
Z
n
,
1
=
1
21
10
-
1
0
-
n
1
1
n
+
1
,
if
n
is odd
,
1
21
10
+
1
0
-
n
1
1
n
+
1
,
if
n
is even
,
V
Z
n
=
m
Z
n
,
2
-
m
Z
n
,
1
2
=
121
100
n
,
n
=
0,1
,
2
,
…
.Note that the expression of the variance does not depend on whether n is even or odd.
In Figure8, the expectation and variance of Zn have been plotted using expressions (92) to carry out computations.
In addition to computing the mean and the variance, further significant information related to the solution SP can be computed from the1-PDF, such as the probability of specific sets in which we could be interested. For example, the probability that the solution varies between the values v1=2 and v2=3 is given by (93)P2≤Zn≤3=∫23f1z,ndz=0ifnis odd,n≤11,e1/21-10-1n+11-e2/212n+45/11n+5-1nifnis odd,13≤n≤17,e1/21-1011-n10n+11ne310/11n-e5/11n2n+1ifnis odd,n≥19,e1/21-51910/11n+2e310/11n-e5/11n2n+1ifnis even.Figure 8
Mean (a) and variance (b) ofZn in Example 8, where Γ0~Exp(λ=1),  a=-11/10, and b=1.
(a)
(b)


## 5.2. Case III.2:B Is a Random Variable

Let us assumea≠-1. In this case, the discrete solution SP (83) takes the form(94)Z0=γ0,Zn=1-an1-aB+γ0an,n=1,2,…. Let n≥1 be an arbitrary and fixed integer and denote Z=Zn. By applying [20, Proposition 2] to(95)Z=Y,X=B,α=1-an1-a≠0,β=γ0an,and taking into account the fact that Z0=γ0, one obtains the 1-PDF(96)f1z,n=δz-γ0,n=0,1-a1-anfB1-a1-anz-γ0an,n=1,2,…,z∈R.Ifa=-1, then the discrete solution SP (83) is(97)Zn=1--1n2B+γ0-1n=γ0ifnis even,B-γ0ifnis odd.Then, forn≥0, an arbitrary but fixed integer, applying [20, Proposition 2] to (97), one obtains the 1-PDF of Zn defined by (97): (98)f1z,n=δz-γ0,neven,fBz+γ0,nodd,z∈R.Remark 9.
Notice that expression (71) in Case II.2 can be obtained assuming that parameter a tends to 1 in (96).Example 9.
LetB be a RV having a gamma distribution of parameters α,β>0; that is, B~Ga(α;β). Let us fix a≠-1. Then, by (96), the 1-PDF of {Zn:n≥0} reads(99)
f
1
z
,
n
=
δ
z
-
γ
0
if
n
=
0
,
z
∈
R
,
f
1
a
z
,
n
,
if
1
-
a
1
-
a
n
>
0
,
n
=
1,2
,
…
,
z
≥
γ
0
a
n
,
f
1
b
z
,
n
,
if
1
-
a
1
-
a
n
<
0
,
n
=
1,2
,
…
,
z
≤
γ
0
a
n
,where(100)f1az,n=1-a1-an1βαΓα1-a1-anz-γ0anα-1e1/βa-1/1-anz-γ0an,f1bz,n=a-11-an1βαΓα1-a1-anz-γ0anα-1e1/βa-1/1-anz-γ0an. Taking into account the fact that the domain of a gamma RV is (0,∞), one deduces the domain of z specified in (99) for each n≥1. Figure 9 shows, f1(z,n) at different values of n assuming that B~Ga(2;4), a=-1/2, and γ0=1. Note that, in accordance with (99), f1z,0=δz-1, -∞<z<∞.Figure 9
2D and 3D plots off1(z,n) in Example 9, where B~Ga(2;4), a=-1/2, and γ0=1. (a) n∈{1,2,3,4}. (b) n∈{1,2,…,10}.
(a)
(b)


## 5.3. Case III.3:A Is a Random Variable

So far, we have taken advantage of RVT method to compute the1-PDF of the solution of problem (2). The success of this approach has relied on the capability to find out an exact expression for the inverse transformation of the mapping that determines the solution in terms of the random inputs. However, in many practical cases, such expression can just be found in an approximate form rather than in an exact form. Under these circumstances, the RVT method can still be very useful. In fact, as it was shown in the analysis of Case III.3 of [20], the application of the Lagrange-Bürmann theorem [27] together with the RVT technique permits determining reliable approximations of the 1-PDF of the solution SP of the continuous counterpart of problem (2). Below, we provide an illustrative example dealing with a particular case of (83) where just the parameter A is assumed to be random. To avoid repetitions, we omit the theoretical development which can be found in [20].Example 10.
Throughout this example, we will use the notation introduced in [20]. Let A be a beta RV of parameters α=2, β=3: A~Be(2;3) and γ0=1 and b=1. In Figure 10, the approximation of f1(z,n) is shown at different values of n≥1. This plot has been made by considering expressions [20, eqs. (78)–(80)] with k=1, being A1=[0,1] because of monotony of r(A). To carry out computations, A1 has been divided into 7 subintervals in agreement with the process described in [20]. In each subinterval, an approximation of degree Nj=2 has been considered. If n=0, f1(z,0)=δ(z-γ0).Figure 10
2D and 3D plots off1(z,n) in Example 10, where A~Be(2;3),  γ0=1, and b=1. (a) n∈{0,1,2,3}. (b) n∈{0,1,2,…,10}.
(a)
(b)


## 5.4. Case III.4:Γ0,B Is a Random Vector

Let us denote byfΓ0,B(γ0,b) the joint PDF of the random vector (Γ0,B).Ifa≠0, the discrete solution SP (83) takes the form(101)Zn=1-an1-aB+anΓ0,n=0,1,2,….Let us fixn:n≥0 and denote Z=Zn=(1-an)/(1-a)B+anΓ0. To compute the 1-PDF of Z, first we will determine the joint PDF of the RVs Z and B by applying [20, Theorem 4] to the two-dimensional RV Y=r(X) being(102)
X
=
Γ
0
B
,
Y
=
Y
1
Y
2
=
r
1
Γ
0
,
B
r
2
Γ
0
,
B
=
1
-
a
n
1
-
a
B
+
a
n
Γ
0
B
.From (102), the inverse transformation of r(X),  X=r-1(Y)=s(Y), takes the form(103)
X
=
Γ
0
B
=
s
1
Y
1
,
Y
2
s
2
Y
1
,
Y
2
=
Y
1
-
1
-
a
n
1
-
a
Y
2
1
a
n
Y
2
. By [20, Theorem 4] and taking into account the fact that |J2|=1/|a|n≠0, the joint PDF fY(y) is given by(104)fY1,Y2y1,y2=1anfΓ0,By1-1-an1-ay21an,y2. Going back to the original RVs, that is, Z=1-an/1-aB+anΓ0=Y1 and B=Y2, one gets(105)fZ,Bz,b=1anfΓ0,Bz-1-an1-ab1an,b.Finally, considering the marginal PDF of Z in (105), one obtains the 1-PDF of Z for a≠0:(106)f1z,n=1an∫b1b2fΓ0,Bz-1-an1-ab1an,bdb,n=0,1,2,…,z∈R.For the trivial casea=0, solution (101) takes the form(107)Z0=Γ0Zn=B,n=1,2,…, and hence the 1-PDF of Z is given by the marginal PDFs of Γ0 and B:(108)
f
1
z
,
n
=
f
Γ
0
γ
0
=
∫
b
1
b
2
f
Γ
0
,
B
γ
0
,
b
d
b
,
n
=
0
,
f
B
b
=
∫
γ
0,1
γ
0,2
f
Γ
0
,
B
γ
0
,
b
d
γ
0
,
n
=
1,2
,
…
,
z
∈
R
.Remark 11.
Notice that expression (77) in Case II.3 can be obtained assuming that parameter a tends to 1 in (106).Example 11.
Let us assume that the random vectorη=(Γ0,B)T~N(μη,Ση) has a bidimensional Gaussian distribution whose mean vector and covariance matrix are given by(109)μη=μΓ0μB,Ση=σΓ02ρΓ0,BσΓ0σBρΓ0,BσZ0σBσB2,respectively. In (109), ρΓ0,B denotes the correlation coefficient of RVs Γ0 and B. Let us assume a≠0. By (106), the 1-PDF of {Zn:n≥0} reads as follows:(110)f1z,n=1an12πdet⁡Ση∫-∞+∞e-1/2ζ-μηTΣη-1ζ-μηdb,n=0,1,2,…,z∈R,being(111)
ζ
=
z
-
1
-
a
n
1
-
a
b
1
a
n
b
. To determine the domain of integration in (110), we have taken into account the fact that B is also a Gaussian RV. Figure 11 shows a graphical representation of f1(z,n) at different values of n assuming that μΓ0=1,  μB=0,  σΓ0=0.1,  σB=0.1,  ρΓ0,B=0.5, and a=-11/10.Figure 11
f
1
(
z
,
n
) in Example 11 at different values of n. The input parameters are assumed to be η=(Γ0,B)T~N(μη,Ση), where the mean vector and the variance-covariance matrix are defined by (109), with the values μΓ0=1,  μB=0,  σΓ0=0.1,  σB=0.1,  ρΓ0,B=0.5, and a=-11/10.


## 5.5. Case III.5:Γ0,A Is a Random Vector

In this case, solution (83) takes the form(112)Zn=AnΓ0+b1-A1-An,n=0,1,2,….Let us denote byfΓ0,A(γ0,a) the joint PDF of the random vector (Γ0,A). Let us fix n:n≥0 and denote Z=Zn=AnΓ0+b/(1-A)(1-An). To compute the 1-PDF of Z, first we will determine the joint PDF of the RVs Z and A by applying [20, Theorem 4] to the two-dimensional RV Y=r(X) with(113)
X
=
Γ
0
A
,
Y
=
Y
1
Y
2
=
r
1
Γ
0
,
A
r
2
Γ
0
,
A
=
A
n
Γ
0
+
b
1
-
A
1
-
A
n
A
.From (113), the inverse transformation of r(X), X=r-1(Y)=s(Y), takes the form(114)
X
=
Γ
0
A
=
s
1
Y
1
,
Y
2
s
2
Y
1
,
Y
2
=
Y
1
-
b
1
-
Y
2
1
-
Y
2
n
1
Y
2
n
Y
2
. Notice that Y2≠0 w.p. 1. By [20, Theorem 4] and taking into account the fact that |J2|=1/|y2|n≠0, the joint PDF fY(y) is given by(115)fY1,Y2y1,y2=1y2nfΓ0,Ay1-b1-y21-y2n1y2n,y2. Going back to the original RVs, that is, Z=AnΓ0+b/(1-A)(1-An)=Y1 and A=Y2, one gets(116)fZ,Az,a=1anfΓ0,Az-b1-a1-an1an,a.Finally, considering the marginal PDF of Z in (116), the 1-PDF of {Zn:n≥0} is given by(117)f1z,n=∫a1a21anfΓ0,Az-b1-a1-an1an,ada,n=0,1,2,…,z∈R.Example 12.
Let us consider the two-dimensional Gaussian vectorη=(Γ0,A)T~N(μη,Ση), whose mean and covariance matrix are given by(118)μη=μΓ0μA,Ση=σΓ02ρΓ0,AσΓ0σAρΓ0,AσΓ0σAσA2,respectively. In (118), ρΓ0,A denotes the correlation coefficient between RVs Γ0 and A. Notice that A has a Gaussian distribution; hence, a1=-∞ and a2=+∞. Then, according to (117), the 1-PDF of {Zn:n≥0} is given by(119)f1z,n=12πdet⁡Ση∫-∞+∞1ane-1/2ζ-μηTΣη-1ζ-μηda,n=0,1,2,…,z∈R, where(120)ζ=z-b1-a1-an1ana.
Figure12 shows a graphical representation of the 1-PDF f1(z,n) given by (117) at different values of n assuming that μΓ0=1,  μA=1.5,  σΓ0=0.1,  σA=0.1, ρΓ0,A=0.5, and b=1. For the sake of clarity, we have split the representation of f1(z,n) into two plots due to the significant differences in the vertical scales required depending on the values of n.Figure 12
f
1
(
z
,
n
) in Example 12 at different values of n. The input parameters are assumed to be η=(Γ0,A)T~N(μη,Ση), where the mean vector and the variance-covariance matrix are defined by (118), with the values μΓ0=1,  μA=1.5,  σΓ0=0.1,  σA=0.1,  ρΓ0,A=0.5, and b=1.


## 5.6. Case III.6:B,A Is a Random Vector

Solution (83) takes the form(121)Zn=Anγ0+1-An1-AB,n=0,1,2,….Let us denote byfB,A(b,a) the joint PDF of the random vector (B,A). Let us fix n:n≥1 and denote Z=Zn=Anγ0+(1-An)/(1-A)B. In order to compute the 1-PDF of Z, first we will determine the joint PDF of the RVs Z and A by applying [20, Theorem 4] to the two-dimensional RV Y=r(X) with(122)
X
=
B
A
,
Y
=
Y
1
Y
2
=
r
1
B
,
A
r
2
B
,
A
=
A
n
γ
0
+
1
-
A
n
1
-
A
B
A
.From (122), the inverse transformation of r(X), X=r-1(Y)=s(Y), takes the form(123)X=BA=s1Y1,Y2s2Y1,Y2=Y1-Y2nγ01-Y21-Y2nY2. By [20, Theorem 4] and taking into account the fact that |J2|=1-y2/1-(y2)n≠0, the joint PDF fY(y) is given by(124)fY1,Y2y1,y2=1-y21-y2nfB,Ay1-y2nγ01-y21-y2n,y2. Going back to the original RVs, that is, Z=Anγ0+1-An/1-AB=Y1 and A=Y2, one gets(125)fZ,Az,a=1-a1-anfB,Az-anγ01-a1-an,a,n=1,2,….Finally, considering the marginal PDF of Z in (125) and the case where n=0, which gives Z0=γ0, one gets the 1-PDF of {Zn:n≥0}:(126)f1z,n=δz-γ0,n=0,∫a1a21-a1-anfB,Az-anγ01-a1-an,ada,n=1,2,…,z∈R.Example 13.
So far, we have considered standard distributions in one or more dimensions to illustrate the obtained theoretical results. Now, we will assume that the joint PDF of the input parametersB and A is constructed by means of a copula transformation. Let us assume that B and A are uniform RVs defined on the interval ]0,1[; that is, B,A~Un(]0,1[). We transform these RVs by the Farlie-Gordon-Morgenstern copula [28], so that a two-dimensional RV (B,A) with joint PDF(127)fB,Ab,a=232-b-a+2ba,if0<b,a<1,is defined. This random vector satisfies the notion that the marginal distributions of fB,A(b,a) keep the one-dimensional distributions of A and B. Hereinafter, let us take γ0=0. Taking into account the fact that 0<a,b<1, by (126), for z and n previously fixed, one must calculate a such that (128)0<z1-a1-an<1.In general, the variable a cannot be determined in a closed form. Using Mathematica software to carry out computations numerically, in Figure 13, several representations of f1(z,n) have been plotted.Figure 13
2D and 3D plots off1(z,n) in Example 13 at different values of n, where γ0=0 and the PDF of the two-dimensional RV (B,A) is given by (127). (a) n∈{1,2,5,10}. (b) n∈{1,2,…,10}.
(a)
(b)


## 5.7. Case III.7:Γ0,B,A Is a Random Vector

In this last case, solution (83) takes the form(129)Zn=AnΓ0+B1-A1-An,n=0,1,2,….Let us denote byfΓ0,B,A(γ0,b,a) the joint PDF of the random vector (Γ0,B,A). Let us fix n:n≥0 and denote Z=Zn=AnΓ0+B/(1-A)(1-An). To determine the 1-PDF of Z, first we will calculate the joint PDF of the RVs Z, B, and A by applying [20, Theorem 4] to the three-dimensional RV Y=r(X) with(130)
X
=
Γ
0
B
A
,
Y
=
Y
1
Y
2
Y
3
=
r
1
Γ
0
,
B
,
A
r
2
Γ
0
,
B
,
A
r
3
Γ
0
,
B
,
A
=
A
n
Γ
0
+
B
1
-
A
1
-
A
n
B
A
.From (130), the inverse transformation of r(X),  X=r-1(Y)=s(Y), takes the form(131)X=Γ0BA=s1Y1,Y2,Y3s2Y1,Y2,Y3s3Y1,Y2,Y3=Y1-Y21-Y31-Y3n1Y3nY2Y3, where Y3≠0 w.p. 1. By [20, Theorem 4] and taking into account the fact that |J3|=1/|y3|n≠0, the joint PDF fY(y) is given by(132)fY1,Y2,Y3y1,y2,y3=1y3nfΓ0,B,Ay1-y21-y31-y3n1y3n,y2,y3.Going back to the original RVs, that is, Z=AnΓ0+B/1-A(1-An)=Y1,  B=Y2, and A=Y3, one gets(133)fZ,B,Az,b,a=1anfΓ0,B,Az-b1-a1-an1an,b,a.Finally, considering the marginal PDF of Z in (133), the 1-PDF of {Zn:n≥0} is given by(134)f1z,n=∫a1a2∫b1b21anfΓ0,B,Az-b1-a1-an1an,b,adbda,n=0,1,2,…,z∈R.Example 14.
Let us assume thatη=(Z0,B,A)T~N(μη,Ση); that is, η is a Gaussian random vector with mean μη=(μ1,μ2,μ3)T∈R3 and variance-covariance matrix Ση∈R3×3. By (134), the 1-PDF of Zn writes(135)f1z,n=12π2πdet⁡Ση∬-∞∞1ane-1/2ζ-μηTΣη-1ζ-μηdbda, where(136)ζ=z-b1-a1-an1anba. In order to illustrate the theoretical results previously established, let us fix the mean vector and the variance-covariance matrix as follows:(137)μη=112,Ση=110411141112.In this example, we do not make the 1-PDF, f1(z,n), explicit, since its expression is cumbersome. In Figure 14, we have plotted f1(z,n) at different values of n.Figure 14
f
1
(
z
,
n
),  n∈{0,1,2,3}, in Example 14. The input parameters (Γ0,B,andA) are assumed to be a Gaussian distribution with mean and variance-covariance matrix defined by (137).


## 6. Conclusions

In this paper, we have provided general explicit formulae to compute the first probability density function (1-PDF) of the solution stochastic process to random first-order linear difference equations. It has been done in the general case where the involved random inputs are statistically dependent. The study has been based on the Random Variable Transformation technique. When solving random difference equations, most of the available studies focus on the computation of the solution stochastic process and its expectation and variance functions. However, the computation of explicit formulae to determine the 1-PDF is more advisable since it permits the computation of other higher-order moments and the probability of certain sets of interest as well. We have shown, through the theoretical development, that the study here presented generalizes its deterministic counterpart. In addition, all the theoretical results have been illustrated by a comprehensive list of examples. Finally, note that our analysis can be extended to determine the 1-PDF of the solution to random nonlinear first-order difference equations in future studies.



---

*Source: 6372108-2016-04-11.xml*