Analytically Computing the Moments of a Conic Combination of Independent Noncentral Chi-Square Random Variables and Its Application for the Extended Cox–Ingersoll–Ross Process with Time-Varying Dimension

Abstract

This paper focuses mainly on the problem of computing the ${\gamma }^{\mathrm{th}}$, $\gamma >0$, moment of a random variable $Y_n:={\sum }_{i=1}^n{\alpha }_i{X}_i$ in which the ${\alpha }_i$’s are positive real numbers and the $X_i$’s are independent and distributed according to noncentral chi-square distributions. Finding an analytical approach for solving such a problem has remained a challenge due to the lack of understanding of the probability distribution of $Y_n$, especially when not all ${\alpha }_i$’s are equal. We analytically solve this problem by showing that the ${\gamma }^{\mathrm{th}}$ moment of $Y_n$ can be expressed in terms of generalized hypergeometric functions. Additionally, we extend our result to computing the ${\gamma }^{\mathrm{th}}$ moment of $Y_n$ when $X_i$ is a combination of statistically independent $Z_i^2$ and $G_i$ in which the $Z_i$’s are distributed according to normal or Maxwell–Boltzmann distributions and the $G_i$’s are distributed according to gamma, Erlang, or exponential distributions. Our paper has an immediate application in interest rate modeling, where we can explicitly provide the exact transition probability density function of the extended Cox–Ingersoll–Ross (ECIR) process with time-varying dimension as well as the corresponding ${\gamma }^{\mathrm{th}}$ conditional moment. Finally, we conduct Monte Carlo simulations to demonstrate the accuracy and efficiency of our explicit formulas through several numerical tests.

1. Introduction

Consider a random variable $Y_n$ driven on a probability space $(\Omega ,\mathcal{F},\mathbb{P})$ defined by

$$Y_n:=\sum _{i=1}^n{\alpha }_i{X}_i$$

(1)

for an integer $n\ge 2$, where ${\alpha }_i>0$, and each $X_i$ is distributed according to a noncentral chi-square distribution with ${\nu }_i>0$ degrees of freedom and a noncentrality parameter ${\delta }_i\ge 0$, i.e., $X_i\sim {\chi }_{{\nu }_i}^2\left({\delta }_i\right)$, for all $i=1,\dots ,n$. We assume that the $X_i$’s are independent with respect to the $\sigma$-field $\mathcal{F}$ and probability measure $\mathbb{P}$.

This paper focuses mainly on the problem of computing the ${\gamma }^{\mathrm{th}}$ moment, $\gamma \in {\mathbb{R}}_{+}$, of $Y_n$ given by

$$\mathbb{E}\left[Y_n^{\gamma }\right]={\int }_0^{\infty }{y}^{\gamma }{f}_{Y_n}\left(y\right)dy$$

(2)

where $f_{Y_n}\left(y\right)$ is the probability density function (PDF) of $Y_n$, and $\mathbb{E}\left[X\right]$ denotes the expected value of a random variable X with respect to the probability measure $\mathbb{P}$. Utilizing the property of the noncentral chi-square random variables [1], we have that $\mathbb{E}\left[X_i^m\right]$ is finite for all non-negative integer m and $i=1,\dots ,n$. This result and the independence of the $X_i$’s ensure the integral on the RHS of (2) is always finite for all $\gamma \in {\mathbb{R}}_{+}$.

The random variable $Y_n$ is found in many statistical applications. In hypothesis testing, several test statistics converge in distribution toward a conic combination of independent noncentral chi-square random variables (see, e.g., [2,3,4,5]). Moreover, $f_{Y_n}\left(y\right)$ and $\mathbb{E}\left[Y_n^{\gamma }\right]$ play an interesting role in financial applications; see, e.g., [6,7,8,9,10,11]. Very recently, Rujivan and Rakwongwan [11], Chumpong et al. [6], and Rujivan [10] showed that the log-return realized variance when the underlying asset follows the extended Black–Scholes model can be expressed in terms of a conic combination of independent noncentral chi-square random variables. As a result, they derived the exact PDF of the log-return realized variance as well as an explicit formula for the ${\gamma }^{\mathrm{th}}$ moment of the log-return realized variance for $\gamma =\frac{1}{2},1$, yielding the first explicit pricing formulas for volatility swaps, volatility options, variance swaps, and variance options, respectively. Furthermore, Rujivan and Rakwongwan [11] utilized the approach proposed in Rujivan [12] for constructing an approximate formula for pricing volatility swaps for the Heston stochastic volatility model. On the other hand, Rujivan [10] proposed an approximate formula for pricing volatility swaps when the underlying asset evolves according to the constant elasticity of variance model.

We now return to the problem of computing the desired moments. Computing (2) is trivial when the values of the ${\alpha }_i$’s are equal and $\gamma =m$ is a non-negative integer. In other words, $Y_n$ reduces to a scaled noncentral chi-square random variable, the PDF of which is known, which in turn implies that $\mathbb{E}\left[Y_n^m\right]$ can be obtained in an explicit form since calculating the integral on the RHS of (2) can be worked out when $\gamma$ is a non-negative integer (see, e.g., [1,13]). On the other hand, it has been repeatedly shown in literature for several decades, see for example in [13,14,15,16,17,18,19,20,21,22,23,24,25,26], that finding an analytical approach for solving the nonlinear problem (2) is significantly more intricate since the PDF of $Y_n$ is not well-known, when the values of some ${\alpha }_i$’s are unequal, and it is not clear which existing representations of the PDF lend themselves to the calculation of the moments. This underlines the importance of our study.

Based on the above discussion, our paper has three aims which we now describe. The principal aim is to provide practitioners an accurate and efficient formula for computing $\mathbb{E}\left[Y_n^{\gamma }\right]$ for any integer $n\ge 2$ and $\gamma \in {\mathbb{R}}_{+}$, including both integers and nonintegers. The next aim is to illustrate further applications of $Y_n$ in interest rate modeling by adopting a Laguerre expansion for the PDF of $Y_n$ which is proposed in this paper to explicitly derive the transition probability density function (TPDF) of the extended Cox–Ingersoll–Ross (ECIR) process with time-varying dimension and which, to our knowledge, has never been found in explicit form until now. In addition, the ECIR process with time-varying dimension was intensively studied by Maghsoodi [27], where its TPDF was given in explicit form and used for pricing bond options for the case where the dimension is constant. The final aim is to utilize the explicit formula for $\mathbb{E}\left[Y_n^{\gamma }\right]$ obtained in this paper to find a novel formula for computing the ${\gamma }^{\mathrm{th}}$ conditional moment of the ECIR process which is accurate and more efficient in terms of computational complexity than existing methods, such as those in [12,28].

The rest of the paper is structured as follows. In Section 2, we derive a Laguerre expansion for the PDF of $Y_n$. In Section 3, by utilizing the Laguerre expansion, we write $\mathbb{E}\left[Y_n^{\gamma }\right]$, $n\in \mathbb{N}$, $\gamma \in {\mathbb{R}}_{+}$, in terms of generalized hypergeometric functions and analytically estimate its truncation errors. In that section, we also extend our result to computing $\mathbb{E}\left[Y_n^{\gamma }\right]$ when $X_i$ is a conic combination of $Z_i^2$ and $G_i$ in which the $Z_i$’s are distributed according to normal or Maxwell–Boltzmann distributions while the $G_i$’s are distributed according to gamma, Erlang, or exponential distributions, assuming that the $Z_i$’s and $G_i$’s are independent. Section 4 illustrates a usage of our result in analyzing the ECIR processes. This includes the first explicit formula for the TPDF of the ECIR process with time-varying dimension and novel explicit formula for the ${\gamma }^{\mathrm{th}}$ conditional moment of the ECIR process. In Section 5, all the explicit formulas proposed in this paper are validated with either Monte Carlo (MC) simulations or other formulas proposed in the literature through several numerical tests. The paper is concluded in Section 6. All proofs are provided in the appendices.

2. The PDF of Y_n

The PDF of $Y_n$ defined in (1) has been studied by many authors for several decades with various representations (see, for instance, in [16,17,21,24,25,26]). In this paper, we use the approach proposed in [16] to obtain a Laguerre expansion for the PDF of $Y_n$.

 Theorem 1.

The PDF of $Y_n$ given in (1) can be expressed as

$$f_{Y_n}\left(y\right)=f_{Y_n}^{\left(\beta \right)}\left(y\right):=\frac{e^{-\frac{y}{2\beta }}{y}^{\frac{\nu }{2}-1}}{{\left(2\beta \right)}^{\frac{\nu }{2}}}\sum _{k=0}^{\infty }\frac{k!}{\mathsf{\Gamma }(\frac{\nu }{2}+k)}{c}_k{\mathbf{L}}_k^{\left(\frac{\nu }{2}-1\right)}\left(\frac{y}{2\beta }\right)\forall y>0$$

(3)

where $\nu :={\sum }_{i=1}^n{\nu }_i$, $\beta >0$ can be arbitrarily chosen, $\mathsf{\Gamma }\left(x\right)$ is the gamma function, and ${\mathbf{L}}_k^{\left(\eta \right)}\left(x\right)$ is the generalized Laguerre function (see [29]). In addition, $c_k,k=0,1,2,\dots$, satisfy the recurrent relations:

$$c_0=1$$

(4)

and

$$c_k=\frac{1}{k}\sum _{j=0}^{k-1}{c}_j{d}_{k-j}\forall k\ge 1,$$

(5)

where

$$d_1=-\frac{1}{2\beta }\sum _{i=1}^n{\delta }_i{\alpha }_i+\frac{1}{2}\sum _{i=1}^n{\nu }_i\left(1-\frac{{\alpha }_i}{\beta }\right)$$

(6)

and

$$d_j=-\frac{j}{2}{\left(\frac{1}{\beta }\right)}^j\sum _{i=1}^n{\delta }_i{\alpha }_i{\left(\beta -{\alpha }_i\right)}^{j-1}+\frac{1}{2}\sum _{i=1}^n{\nu }_i{\left(1-\frac{{\alpha }_i}{\beta }\right)}^j\forall j\ge 2.$$

(7)

Proof. 

The proof is provided in Appendix A. □

A couple of remarks should be made about the free parameter $\beta$. First, we note that the impact of $\beta$ goes beyond Equation (3) as it also influences the $c_k$ coefficients through the recurrence relations (4)–(7). Second, the value of $\beta$ can also influence the convergence rate of (3). Indeed, if the $c_k$ coefficients diverge, then it would be more challenging to reliably approximate the infinite sum on the right-hand side of (3) by its truncated version. As a result of this, we follow the procedures of [16] to study and promote appropriate choices of $\beta$ in Section 3.

3. The ${\mathit{\gamma }}^{\mathrm{th}}$ Moment of Y_n

Computing $\mathbb{E}\left[Y_n^{\gamma }\right]$ as given in (2) can be achieved with any desired level of accuracy when the PDF of $Y_n$ is explicitly known. In this section, we use the Laguerre expansion (3) to obtain an explicit formula for $\mathbb{E}\left[Y_n^{\gamma }\right]$ as well as to showcase some interesting applications of the Laguerre expansion (3) in the following subsections.

3.1. Our Explicit Formula for $\mathbb{E}\left[Y_n^{\gamma }\right]$

From the Laguerre expansion (3) together with some properties of the generalized hypergeometric function

$${\displaystyle {}_2{\mathbf{F}}_1\left(a_1,a_2;b_1;z\right)=\sum _{k=0}^{\infty }\frac{{\left(a_1\right)}_k{\left(a_2\right)}_k}{{\left(b_1\right)}_k}\frac{z^k}{k!},}$$

where ${(·)}_k$ denotes the usual Pochhammer symbol, known from [30], and we derive a simple explicit formula for the ${\gamma }^{\mathrm{th}}$ moment of $Y_n$ for any $\gamma \in {\mathbb{R}}_{+}$.

 Theorem 2.

For any $\gamma \in {\mathbb{R}}_{+}$, we have

$$\mathbb{E}\left[Y_n^{\gamma }\right]={\left(2\beta \right)}^{\gamma }\sum _{k=0}^{\infty }{(-1)}^k\frac{\mathsf{\Gamma }\left(\gamma +k+\frac{\nu }{2}\right)}{\mathsf{\Gamma }\left(k+\frac{\nu }{2}\right)}{}_2{\mathit{F}}_1\left(-k,1-k-\frac{\nu }{2};1-k-\frac{\nu }{2}-\gamma ;1\right)c_k$$

(8)

where the coefficients $c_k,k=0,1,\dots$, are chosen according to (4)–(7), and the parameters $\nu ={\sum }_{i=1}^n{\nu }_i$ and $\beta >0$ can be arbitrarily chosen.

Proof. 

The proof is provided in Appendix B. □

Theorem 2 essentially expresses $\mathbb{E}\left[Y_n^{\gamma }\right]$ in terms of generalized hypergeometric functions. We remark that computing $\mathbb{E}\left[Y_n^{\gamma }\right]$ relies on the $c_k$ coefficients which can be obtained from the recursive Formulas (4)–(7), and we demonstrate later in our numerical study in Section 5 that implementing our Formula (8) for computing $\mathbb{E}\left[Y_n^{\gamma }\right]$ consumes significantly less time and effort than employing MC simulations.

In terms of applications, the result presented in the following corollary can be, for instance, applied to obtain an analytical formula for pricing volatility swaps in the discrete observation case based on the Black–Scholes model with time-varying risk-free interest rate and time-varying volatility as proposed in Theorem 3.1 of Rujivan [10].

 Corollary 1.

We have

$$\mathbb{E}\left[\sqrt{Y_n}\right]=\sqrt{2\beta }\frac{\mathsf{\Gamma }\left(\frac{\nu +1}{2}\right)}{\mathsf{\Gamma }\left(\frac{\nu }{2}\right)}\sum _{k=0}^{\infty }{}_2{\mathit{F}}_1\left(-k,\frac{\nu +1}{2};\frac{\nu }{2};1\right)c_k$$

(9)

where the coefficients $c_k,k=0,1,\dots$, are chosen according to (4)–(7) and the parameters $\nu ={\sum }_{i=1}^n{\nu }_i$ and $\beta >0$ can be arbitrarily chosen.

 Proof. 

The proof is provided in Appendix B. □

Another interesting special case of Theorem 1 is when $n=1$, that is, there is only one summand, say X. In this case, we can leverage Theorem 1 to compute a noninteger moment of any chi-square random variable as follows.

 Corollary 2.

For any $X\sim {\chi }_{\nu }^2\left(\delta \right)$ and $\gamma \in {\mathbb{R}}_{+}$, we have

$$\mathbb{E}\left[X^{\gamma }\right]=2^{\gamma }\sum _{k=0}^{\infty }\frac{\mathsf{\Gamma }\left(\gamma +k+\frac{\nu }{2}\right)}{\mathsf{\Gamma }\left(k+\frac{\nu }{2}\right)}{}_2{\mathit{F}}_1\left(-k,1-k-\frac{\nu }{2};1-k-\frac{\nu }{2}-\gamma ;1\right)\frac{1}{k!}{\left(\frac{\delta }{2}\right)}^k.$$

(10)

Proof. 

The proof is provided in Appendix B. □

3.2. Estimates for Truncation Errors of $\mathbb{E}\left[Y_n^{\gamma }\right]$

To implement $\mathbb{E}\left[Y_n^{\gamma }\right]$ on a computer, it is necessary to investigate truncation errors, that is, to quantify a loss due to replacing an infinite sum with a finite sum. This subsection derives an estimate for the truncation errors of $\mathbb{E}\left[Y_n^{\gamma }\right]$ by applying the results proposed in [16] as follows.

To begin with, based on the formula of $\mathbb{E}\left[Y_n^{\gamma }\right]$ in (8), we define

$$E_{k_1,k_2}^{\left(\gamma \right)}:={\left(2\beta \right)}^{\gamma }\sum _{k=k_1+1}^{k_2}{(-1)}^k\frac{\mathsf{\Gamma }\left(\gamma +k+\frac{\nu }{2}\right)}{\mathsf{\Gamma }\left(k+\frac{\nu }{2}\right)}{}_2{\mathbf{F}}_1\left(-k,1-k-\frac{\nu }{2};1-k-\frac{\nu }{2}-\gamma ;1\right)c_k$$

(11)

for any $k_1,k_2\in \mathbb{N}\cup \{0,\infty \}$ such that $k_1+1<k_2$. Therefore, $E_{K,\infty }^{\left(\gamma \right)}$ represents a truncation error of order K of $\mathbb{E}\left[Y_n^{\gamma }\right]$. To estimate this truncation error, we first derive valid bounds for $c_k$ for $k=1,2,\dots$, defined in (5)–(7).

 Lemma 1.

The coefficient $c_k$, $k\in \mathbb{N}$, satisfies

$$|c_k|\le e^{\frac{\delta }{2\zeta }}{\left(\frac{2k+\nu }{2k}\right)}^k{\left(\frac{2k+\nu }{\nu }\right)}^{\frac{\nu }{2}}{\zeta }^k,$$

(12)

where $\delta ={\sum }_{i=1}^n{\delta }_i$ and $\zeta =\frac{1}{\beta }{max}_{i\in \{1,\dots ,n\}}\left|\beta -{\alpha }_i\right|$. Moreover, if $\beta >\frac{1}{2}{max}_{i\in \{1,\dots ,n\}}{\alpha }_i$, then $0<\zeta <1$.

Proof. 

The proof is provided in Appendix B. □

From Lemma 1 above, we further define

$$B_{k_1,k_2}^{\left(\gamma \right)}\left(\zeta \right):={\left(2\beta \right)}^{\gamma }{e}^{\frac{\delta }{2\zeta }}\sum _{k=k_1+1}^{k_2}{b}_k(\gamma ,\nu ,\zeta )$$

(13)

for all $k_1,k_2\in \mathbb{N}\cup \{0,\infty \}$ such that $k_1+1<k_2$ and $\zeta >0$, where

$$b_k(\gamma ,\nu ,\zeta )=\frac{\mathsf{\Gamma }\left(\gamma +k+\frac{\nu }{2}\right)}{\mathsf{\Gamma }\left(k+\frac{\nu }{2}\right)}\left|{}_2{\mathbf{F}}_1\left(-k,1-k-\frac{\nu }{2};1-k-\frac{\nu }{2}-\gamma ;1\right)\right|{\left(\frac{2k+\nu }{2k}\right)}^k{\left(\frac{2k+\nu }{\nu }\right)}^{\frac{\nu }{2}}{\zeta }^k.$$

(14)

Utilizing Lemma 1, we obtain the following upper bound of a truncation error.

 Theorem 3.

Supposing that $\beta >\frac{1}{2}{max}_{i\in \{1,\dots ,n\}}{\alpha }_i$, then we have

$$\left|E_{K,\infty }^{\left(\gamma \right)}\right|\le B_{K,\infty }^{\left(\gamma \right)}\left(\zeta \right)\forall \gamma \in {\mathbb{R}}_{+},k\in \mathbb{N},$$

(15)

where $\zeta =\frac{1}{\beta }{max}_{i\in \{1,\dots ,n\}}\left|\beta -{\alpha }_i\right|$. Furthermore,

$$\underset{K\to \infty }{lim}{E}_{K,\infty }^{\left(\gamma \right)}=0.$$

(16)

Proof. 

The proof is provided in Appendix B. □

It should be mentioned from Theorem 3 that the inequality $\beta >\frac{1}{2}{max}_{i\in \{1,\dots ,n\}}{\alpha }_i$ ought to hold when we implement the Formula (8) for computing $\mathbb{E}\left[Y_n^{\gamma }\right]$ to ensure that the truncation occurring tends to zero when K approaches infinity.

Finally, we use Euler’s transformation [31] in order to show that (8) terminates when $\gamma =m$ is a non-negative integer.

 Theorem 4.

For any $m\in \mathbb{N}$, we have

$$\mathbb{E}\left[Y_n^m\right]={\left(2\beta \right)}^m\sum _{k=0}^m{(-1)}^k\frac{\mathsf{\Gamma }\left(m+k+\frac{\nu }{2}\right)}{\mathsf{\Gamma }\left(k+\frac{\nu }{2}\right)}{}_2{\mathit{F}}_1\left(-k,1-k-\frac{\nu }{2};1-k-\frac{\nu }{2}-m;1\right)c_k$$

(17)

where the coefficients $c_k,k=0,1,\dots$, are chosen according to (4)–(7) and the parameters $\nu ={\sum }_{i=1}^n{\nu }_i$ and $\beta >0$ can be arbitrarily chosen.

Proof. 

The proof is provided in Appendix B. □

Applying Corollary 2 and Theorem 4, an explicit formula for the $m^{\mathrm{th}}$ moment of noncentral chi-square random variables can be obtained as follows.

 Corollary 3.

For any $X\sim {\chi }_{\nu }^2\left(\delta \right)$ and $m\in \mathbb{N}$, we have

$$\mathbb{E}\left[X^m\right]=m!2^m\mathsf{\Gamma }\left(m+\frac{\nu }{2}\right)\sum _{k=0}^m\frac{{\left(\frac{\delta }{2}\right)}^k}{k!(m-k)!\mathsf{\Gamma }\left(k+\frac{\nu }{2}\right)}.$$

(18)

Proof. 

The proof is provided in Appendix B. □

3.3. Analytical Formulas for Other Conic Combinations of Independent Random Variables

In this subsection, we extend our previous results to various other types of random variables. First, instead of assuming that each component $X_i$ is a chi-square random variable, we assume instead that it is normally distributed with varying means and variances. Second, motivated by its application in queuing theory, we focus on the case where each component follows a gamma distribution with varying shape parameters. Note that gamma distributions are often used in queuing theory for modeling the distribution of certain types of waiting times, e.g., the excess water flow of a dam as explained in Mathai [32] and other problems in communication theory with respect to the performance of certain wireless transmission systems as described in Alouini et al. [33]. Third, we focus on the sum of independent Erlang distributed random variables, which lie at the core of many fields such as telecommunications, statistics, reliability theory, and risk analysis [34]. Last but not least, we focus on the sum of exponential random variables, which are often used in stochastic modeling thanks to its memoryless property, and the sum of squared Maxwell–Boltzmann random variables, which can be used to explain the molecular speed distribution of ideal gases [35], etc.

 Theorem 5.

Consider a random variable $Y_{(1,n)}:={\sum }_{i=1}^n{a}_{(1,i)}{Z}_i^2$, where each $a_{(1,i)}>0$ and each $Z_i$ is a normal random variable with mean ${\mu }_{(1,i)}\in \mathbb{R}$ and variance ${\sigma }_{(1,i)}^2$ for ${\sigma }_{(1,i)}>0$. Assuming that all summands are independent, then the PDF of $Y_{(1,n)}$, $\mathbb{E}\left[Y_{(1,n)}^{\gamma }\right]$, and $\mathbb{E}\left[Y_{(1,n)}^m\right]$ can be computed using (3), (8), and (17), respectively, for all $\gamma \in {\mathbb{R}}_{+}$ and integer $m\in \mathbb{N}$, by setting ${\alpha }_i=a_{(1,i)}{\sigma }_{(1,i)}^2$, ${\nu }_i=1$, and ${\delta }_i={\left(\frac{{\mu }_{(1,i)}}{{\sigma }_{(1,i)}}\right)}^2$ for all $i=1,\dots ,n$.

Proof. 

The proof is provided in Appendix B. □

 Theorem 6.

Consider a random variable $Y_{(2,n)}:={\sum }_{i=1}^n{a}_{(2,i)}{G}_i$, where each $a_{(2,i)}>0$ and each $G_i$ is distributed according to a gamma distribution with shape parameter ${\kappa }_{(2,i)}>0$ and scale parameter ${\theta }_{(2,i)}>0$. Assuming all summands are independent, then the PDF of $Y_{(2,n)}$, $\mathbb{E}\left[Y_{(2,n)}^{\gamma }\right]$, and $\mathbb{E}\left[Y_{(2,n)}^m\right]$ can be computed using (3), (8), and (17), respectively, for all $\gamma \in {\mathbb{R}}_{+}$ and integer $m\in \mathbb{N}$, by setting ${\alpha }_i=\frac{1}{2}{a}_{(2,i)}{\theta }_{(2,i)}$, ${\nu }_i=2{\kappa }_{(2,i)}$, and ${\delta }_i=0$ for all $i=1,\dots ,n$.

Proof. 

The proof is provided in Appendix B. □

 Theorem 7.

Consider a random variable $Y_{(3,n)}:={\sum }_{i=1}^n{a}_{(3,i)}{L}_i$, where each $a_{(3,i)}>0$ and each $L_i$ is distributed according to an Erlang distribution with shape parameter ${\kappa }_{(3,i)}\in \{1,2,\dots \}$ and rate parameter ${\lambda }_{(3,i)}>0$. Assuming that all summands $L_i$ are independent, then the PDF of $Y_{(3,n)}$, $\mathbb{E}\left[Y_{(3,n)}^{\gamma }\right]$, and $\mathbb{E}\left[Y_{(3,n)}^m\right]$ can be computed using (3), (8), and (17), respectively, for all $\gamma \in {\mathbb{R}}_{+}$ and integer $m\in \mathbb{N}$, by setting ${\alpha }_i=\frac{a_{(3,i)}}{2{\lambda }_{(3,i)}}$, ${\nu }_i=2{\kappa }_{(3,i)}$, and ${\delta }_i=0$ for all $i=1,\dots ,n$.

Proof. 

The proof is provided in Appendix B. □

 Theorem 8.

Consider a random variable $Y_{(4,n)}:={\sum }_{i=1}^n{a}_{(4,i)}{P}_i$, where each $a_{(4,i)}>0$ and each $P_i$ is distributed according to an exponential distribution with rate parameter ${\lambda }_{(4,i)}>0$. Assuming that all summands are independent, then the PDF of $Y_{(4,n)}$, $\mathbb{E}\left[Y_{(4,n)}^{\gamma }\right]$, and $\mathbb{E}\left[Y_{(4,n)}^m\right]$ can be computed using (3), (8), and (17), respectively, for all $\gamma \in {\mathbb{R}}_{+}$ and integer $m\in \mathbb{N}$, by setting ${\alpha }_i=\frac{a_{(4,i)}}{2{\lambda }_{(4,i)}}$, ${\nu }_i=2$, and ${\delta }_i=0$ for all $i=1,\dots ,n$.

Proof. 

The proof is provided in Appendix B. □

 Theorem 9.

Consider a random variable $Y_{(5,n)}:={\sum }_{i=1}^n{a}_{(5,i)}{W}_i^2$, where each $a_{(5,i)}>0$ and each $W_i$ is distributed according to a Maxwell–Boltzmann distribution with parameter ${\varphi }_{(5,i)}>0$. Assuming that all $W_i$’s are independent, then the PDF of $Y_{(5,n)}$, $\mathbb{E}\left[Y_{(5,n)}^{\gamma }\right]$, and $\mathbb{E}\left[Y_{(5,n)}^m\right]$ can be computed using (3), (8), and (17), respectively, for all $\gamma \in {\mathbb{R}}_{+}$ and integer $m\in \mathbb{N}$, by setting ${\alpha }_i=a_{(5,i)}{\varphi }_{(5,i)}^2$, ${\nu }_i=3$, and ${\delta }_i=0$ for all $i=1,\dots ,n$.

Proof. 

The proof is provided in Appendix B. □

4. Extensions to the ECIR Process with Time-Varying Dimension

The ECIR process is one of the most widely used processes to model interest rates and to price financial products such as zero-coupon bond, ex-coupon, moment swaps, options, and interest rate swaps. With time-dependent parameters, the ECIR process is capable of accounting for side information from potential political or economic events. Formally, according to Maghsoodi [27], the ECIR process, denoted by $V_t$, satisfies

$$d{V}_t=\kappa \left(t\right)(\theta \left(t\right)-V_t)dt+\sigma \left(t\right)\sqrt{V_t}d{W}_t$$

(19)

for $t\in (0,T]$ and $T>0$ with an initial value $V_0=v_0>0$, where the parameter functions $\theta \left(t\right)>0$, $\kappa \left(t\right)>0$, and $\sigma \left(t\right)>0$ are continuous on $[0,T]$, and $W_t$ is a standard Brownian motion under a probability space $(\Omega ,\mathcal{F},\mathbb{P})$ with a filtration ${\left({\mathcal{F}}_t\right)}_{0\le t\le T}.$ Note that $V_t$ reduces to a plain-vanilla Cox–Ingersoll–Ross process (CIR process) [36] provided that the relevant parameter functions are constants.

Focusing on the ECIR process (19), we define the dimension of $V_t$ as

$$\mathbf{d}\left(t\right):=\frac{4\kappa \left(t\right)\theta \left(t\right)}{{\sigma }^2\left(t\right)}$$

(20)

for $t\in [0,T]$, and this quantity plays an important role for deriving an expression for the distribution of $V_t$. Maghsoodi [27] discovered that when $\mathbf{d}\left(t\right)=d\ge 2$ for all $t\in [0,T]$ that included the CIR process, $V_t$ never hit zero almost surely and was in fact a scaled time-changed squared Bessel process; as a result, the TPDF of $V_t$ was explicitly given.

Moreover, Maghsoodi [27] showed that $V_t$ could be represented as a lognormal process through a stochastic time-change when $\mathbf{d}\left(t\right)\ge 2$ for all $t\in [0,T]$, but the TPDF of $V_t$ was not analytically derived. Consequently, it has been an open question until now how the TPDF of $V_t$ can be obtained in explicit form when $\mathbf{d}\left(t\right)$ is time-varying, based on the stochastic time-varying lognormal process representation.

To demonstrate our contribution in the current paper for solving this problem, we apply our previous results in Section 2 and Section 3 to explicitly derive the TPDF of $V_t$ as well as its ${\gamma }^{\mathrm{th}}$ conditional moment when $\mathbf{d}\left(t\right)$ is time-varying, provided that the following two assumptions hold:

 Assumption 1.

$\mathbf{d}\left(t\right)\ge 2$ for all $t\in [0,T]$.

 Assumption 2.

The derivative ${\mathbf{d}}^{\left(1\right)}\left(t\right)$ of $\mathbf{d}\left(t\right)$ with respect to t satisfies $0\le {\mathbf{d}}^{\left(1\right)}\left(t\right)<\infty$ for all $t\in [0,T]$.

4.1. The Exact TPDF of the ECIR Process with Time-Varying Dimension

To realize our objective, we firstly define a parameter function as follows:

$$\tau (t,s):=\frac{1}{4}{\int }_s^t{\sigma }^2\left(\zeta \right)e^{-{\int }_{\zeta }^t\kappa \left(u\right)du}d\zeta$$

(21)

for $0\le s\le t\le T$.

Peng and Schellhorn [37] showed that the ECIR process $V_t$ described by (19) could be represented as that of a convergent series of weighted independent noncentral chi-square and chi-square random variables. For the sake of completeness, we summarize their results in the following theorem.

 Theorem 10.

Supposing that Assumptions 1 and 2 hold, then the ECIR process $V_t$ described by (19) with an initial value $v_0>0$ can be expressed as

$$V_t\stackrel{law}{=}\underset{n\to \infty }{lim}\sum _{i=1}^n{\widehat{\alpha }}_i{\widehat{X}}_i$$

(22)

for any $t\in (0,T]$, where the random variables ${\widehat{X}}_i$ are independent and distributed according to noncentral chi-square and chi-square distributions as

$${\widehat{X}}_i\sim {\chi }_{{\widehat{\nu }}_i}^2\left({\widehat{\delta }}_i\right),$$

(23)

with the coefficients and parameters in (22) and (23) given by

$${\widehat{\alpha }}_i=\tau \left(t,(i-1)\frac{t}{n}\right)\forall i\in \{1,\dots ,n\},$$

(24)

$${\widehat{\nu }}_1=\mathbf{d}\left(0\right),$$

(25)

$${\widehat{\delta }}_1=\frac{v_0}{\tau (t,0)}{e}^{-{\int }_0^t\kappa \left(u\right)du},$$

(26)

and

$${\widehat{\nu }}_i={\mathbf{d}}^{\left(1\right)}\left((i-1)\frac{t}{n}\right)\frac{t}{n}\forall i\in \{2,\dots ,n\},$$

(27)

$${\widehat{\delta }}_i=0\forall i\in \{2,\dots ,n\}.$$

(28)

In particular, if $\mathbf{d}\left(s\right)=d\ge 2$ for all $s\in [0,t]$, then

$$V_t\sim \tau (t,0)·{\chi }_d^2\left(\frac{v_0}{\tau (t,0)}{e}^{-{\int }_0^t\kappa \left(u\right)du}\right)$$

(29)

Proof. 

See Theorem 3.1 in Peng and Schellhorn [37]. □

Peng and Schellhorn [37] also represented the TPDF of $V_t$ in terms of a limit of a sequence of convolutions of the PDFs of scaled noncentral chi-square and chi-square random variables.

Instead of utilizing the convolution property for independent random variables as shown by Peng and Schellhorn [37], we apply Theorem 1 to obtain the first explicit formula for the TPDF of $V_t$ with time-varying dimension $\mathbf{d}\left(t\right)$.

 Theorem 11.

The TPDF of $V_t$ defined by

$$f_{V_t}(v,t|v_0):=\mathbb{P}\left(V_t=v|V_0=v_0\right)$$

(30)

for $v,v_0>0$ and $t\in (0,T]$ can be expressed as

$$f_{V_t}(v,t|v_0)=\frac{e^{-\frac{v}{2\tau (t,0)}}{v}^{\frac{\mathbf{d}\left(t\right)}{2}-1}}{{\left(2\tau (t,0)\right)}^{\frac{\mathbf{d}\left(t\right)}{2}}}\sum _{k=0}^{\infty }\frac{k!}{\mathsf{\Gamma }(\frac{\mathbf{d}\left(t\right)}{2}+k)}{\widehat{c}}_k(t,v_0){\mathit{L}}_k^{\left(\frac{\mathbf{d}\left(t\right)}{2}-1\right)}\left(\frac{v}{2\tau (t,0)}\right)$$

(31)

where

$${\widehat{c}}_0(t,v_0)=1,$$

(32)

$${\widehat{c}}_k(t,v_0)=\frac{1}{k}\sum _{j=0}^{k-1}{\widehat{c}}_j(t,v_0){\widehat{d}}_{k-j}(t,v_0)\forall k\in \mathbb{N},$$

(33)

and

$${\widehat{d}}_1(t,v_0)=-\frac{1}{2\tau (t,0)}{v}_0{e}^{-{\int }_0^t\kappa \left(u\right)du}+\frac{1}{2}{\int }_0^t{\mathbf{d}}^{\left(1\right)}\left(s\right)\left(1-\frac{\tau (t,s)}{\tau (t,0)}\right)ds,$$

(34)

$${\widehat{d}}_j(t,v_0)=\frac{1}{2}{\int }_0^t{\mathbf{d}}^{\left(1\right)}\left(s\right){\left(1-\frac{\tau (t,s)}{\tau (t,0)}\right)}^{j}ds\forall j\in \mathbb{N}\setminus \left\{1\right\}.$$

(35)

In particular, if $\mathbf{d}\left(s\right)=d\ge 2$ for all $s\in [0,t]$, then

$${\widehat{c}}_k(t,v_0)={\left(-\frac{e^{-{\int }_0^t\kappa \left(u\right)du}}{2\tau (t,0)}\right)}^k\frac{v_0^k}{k!}\forall k\in \mathbb{N}\cup \left\{0\right\}.$$

(36)

Proof. 

The proof is provided in Appendix C. □

4.2. The ${\gamma }^{\mathrm{th}}$ Conditional Moment of the ECIR Process with Time-Varying Dimension

For $\gamma \in {\mathbb{R}}_{+}$ and a probability space $(\Omega ,\mathcal{F},\mathbb{P})$ with a filtration ${\left({\mathcal{F}}_t\right)}_{0\le t\le T}$, we define the ${\gamma }^{\mathrm{th}}$ conditional moment of the ECIR process $V_t$ as

$$U_E^{\left(\gamma \right)}\left(t\right|v_0):={\mathbb{E}}^{\mathbb{P}}\left[V_t^{\gamma }|{\mathcal{F}}_0\right]={\mathbb{E}}^{\mathbb{P}}\left[V_t^{\gamma }|V_0=v_0\right]={\int }_0^{\infty }{v}^{\gamma }{f}_{V_t}(v,t|v_0)dv$$

(37)

for $t\in (0,T]$ and $v_0>0$, where $f_{V_t}(v,t|v_0)$ is the TPDF of $V_t$ given in (31).

Rujivan [12] first presented a recursive formula for computing $U_E^{\left(\gamma \right)}\left(t\right|v_0)$ using a partial differential equation (PDE) approach. Alternatively, we apply Theorems 2, 4, and 11 to obtain a novel explicit formula for $U_E^{\left(\gamma \right)}\left(t\right|v_0)$ in the following theorem.

 Theorem 12.

Supposing that Assumptions 1 and 2 hold, then for any $\gamma \in {\mathbb{R}}_{+}$ we have

$$\begin{array}{cc}\hfill & U_E^{\left(\gamma \right)}\left(t\right|v_0)\hfill \\ & ={\left(2\tau (t,0)\right)}^{\gamma }\sum _{k=0}^{\infty }{(-1)}^k\frac{\mathsf{\Gamma }\left(\gamma +k+\frac{\mathbf{d}\left(t\right)}{2}\right)}{\mathsf{\Gamma }\left(k+\frac{\mathbf{d}\left(t\right)}{2}\right)}{}_2{\mathit{F}}_1\left(-k,1-k-\frac{\mathbf{d}\left(t\right)}{2};1-k-\frac{\mathbf{d}\left(t\right)}{2}-\gamma ;1\right){\widehat{c}}_k(t,v_0)\hfill \end{array}$$

(38)

for $t\in (0,T]$ and $v_0>0$, where ${\widehat{c}}_k(t,v_0),k=0,\dots ,$ are given in (32) and (33).

In particular, for any integer $m\in \mathbb{N}$,

$$\begin{array}{cc}\hfill & U_E^{\left(m\right)}\left(t\right|v_0)\hfill \\ & ={\left(2\tau (t,0)\right)}^m\sum _{k=0}^m{(-1)}^k\frac{\mathsf{\Gamma }\left(m+k+\frac{\mathbf{d}\left(t\right)}{2}\right)}{\mathsf{\Gamma }\left(k+\frac{\mathbf{d}\left(t\right)}{2}\right)}{}_2{\mathit{F}}_1\left(-k,1-k-\frac{\mathbf{d}\left(t\right)}{2};1-k-\frac{\mathbf{d}\left(t\right)}{2}-m;1\right){\widehat{c}}_k(t,v_0)\hfill \end{array}$$

(39)

for $t\in (0,T]$ and $v_0>0$.

Proof. 

The proof is provided in Appendix C. □

4.2.1. Comparison with Other Formulas

To date, the computation of the conditional moment has only been partially solved due to the unavailability of the transitional PDF. Indeed, the problem of computing the integral on the RHS of (37) with any stochastic differential equation (SDE) is typically addressed by the Feynman–Kac theorem, where the partial differential equation (PDE) is solved analytically, and some combinatorial techniques are used to simplify the system of recursive ordinary differential equations (ODEs) associated with the conditional moment; see, for instance, [38,39,40,41], for more details.

For a more concrete comparison, Rujivan [12] presented the first explicit formula for the $\gamma$ conditional moment of the ECIR process (19) with time-varying dimension as a power series. Rujivan demonstrated the effectiveness of this analytical approach over the other state-of-the-art techniques including the method by Dufresne [28] and MC simulations. This result has some similarities and dissimilarities to our work, which we shall explain. Our formula (38) expresses the $\gamma$ conditional moment of the ECIR process as an infinite series where the coefficients ${\widehat{c}}^k(t,v_0)$ can be computed analytically recursively. This offers a more efficient way than, say, using Equation (2.2) of Rujivan [12], which can achieve the same purpose of characterizing the conditional moment of the ECIR process. Therein, the parameters must also be computed recursively but the result of each iteration does not have a closed form. The method of Rujivan [12] is therefore more time-consuming and more prone to numerically accumulating errors.

5. Numerical Results and Discussions

As shown in Section 2, Section 3 and Section 4, our theoretical frameworks presented in this paper produce several new explicit formulas for computing the PDF of $Y_n$ defined in (1) and its moments, including the TPDF of the ECIR process (19) with time-varying dimension and its conditional moments. A natural question that may be raised by practitioners is whether these newly derived explicit formulas are accurate and efficient, especially considering that an infinite sum has to be truncated. Therefore, we intensively investigated the accuracies of our explicit formulas to confirm that there were no algebraic errors in the derivation processes as well as to demonstrate the efficiencies of our explicit formulas compared with either MC simulations or other explicit formulas proposed in the literature through a series of numerical examples which were coded in MATHEMATICA 11 and executed on a notebook with the following specifications: Intel(R) Core (TM) i5-6500, CPU @3.20GHz, 16GB RAM, Windows 10, 64 bit operating system.

5.1. The accuracy of Our Explicit Formula for $f_{Y_n}^{\left(\beta \right)}\left(y\right)$

In order to illustrate the accuracy of our explicit Formula (3), we introduced random variables as follows. For any $n\in \mathbb{N}$, we defined

$$Y_n^{\left(j\right)}:=\sum _{i=1}^n{\alpha }_i^{\left(j\right)}{X}_i^{\left(j\right)}$$

(40)

for $j=1,2,3$, where

$$X_i^{\left(j\right)}\sim {\chi }_{{\nu }_i^{\left(j\right)}}^2\left({\delta }_i^{\left(j\right)}\right)$$

(41)

for $i=1,\dots ,n$, and the parameters were set as

$${\alpha }_i^{\left(1\right)}={\alpha }_i^{\left(2\right)}={\alpha }_i^{\left(3\right)}=\frac{2i\alpha }{n(n+1)},$$

(42)

$${\nu }_i^{\left(1\right)}={\nu }_i^{\left(2\right)}={\nu }_i^{\left(3\right)}=\frac{i+3}{2},$$

(43)

$${\delta }_i^{\left(1\right)}=\frac{i}{10},{\delta }_i^{\left(2\right)}=0,{\delta }_i^{\left(3\right)}=\frac{(1-{(-1)}^i)i}{20},$$

(44)

for $\alpha >0$. Furthermore, for $n\ge 2$, we assumed that $X_i^{\left(j\right)},i=1,\dots ,n$ were independent for all $j=1,2,3$. By construction, we note that $\alpha$ represents the total conic coefficients, i.e., $\alpha ={\sum }_{i=1}^n{\alpha }_i^{\left(j\right)}$, and that each $Y_n^{\left(j\right)}$ constitutes a sum of independent noncentral chi-square random variables. Though, we note that through the transformations studied in Theorems 5–9, the distribution of $Y_n^{\left(j\right)}$ may be identical to those of other random sums. For instance,

$$Y_n^{\left(2\right)}\stackrel{law}{=}\sum _{i=1}^n\frac{2{\alpha }_i^{\left(2\right)}}{{\theta }_{(2,i)}}{G}_i\stackrel{law}{=}\sum _{i=1}^{n}2{\alpha }_i^{\left(2\right)}{\lambda }_{(3,i)}{L}_i,$$

where the $G_i$’s are independent random variables and are distributed according to a gamma distribution with the shape parameter ${\nu }_i^{\left(2\right)}/2$ and the scale parameter ${\theta }_{(2,i)}>0$, and similarly the $L_i$’s are independent random variables and are distributed according to an Erlang distribution with the shape parameter ${\nu }_i^{\left(2\right)}/2$ and the rate parameter ${\lambda }_{(3,i)}>0$.

Example 1.

We started by considering the PDF of $Y_6^{\left(1\right)}={\sum }_{i=1}^6{\alpha }_i^{\left(1\right)}{X}_i^{\left(1\right)}$ in which the values of parameters ${\alpha }_i^{\left(1\right)},{\nu }_i^{\left(1\right)}$, and ${\delta }_i^{\left(1\right)}$ for $i=1,\dots ,6$ are plotted in Figure 1a. The PDFs of ${\alpha }_i^{\left(1\right)}{X}_i^{\left(1\right)},i=1,\dots ,6$ which varied in range and shape, are shown in Figure 2a–f, respectively. The problem of analytically computing the PDF of $Y_6^{\left(1\right)}$ was investigated as follows.

To obtain the PDF of $Y_6^{\left(1\right)}$, denoted by $f_{Y_6^{\left(1\right)}}^{\left(\beta \right)}\left(y\right)$, we set $n=6,{\alpha }_i={\alpha }_i^{\left(1\right)},{\nu }_i={\nu }_i^{\left(1\right)},$ and ${\delta }_i={\delta }_i^{\left(1\right)}$ for $i=1,\dots ,6$ in (3). The $c_k$ coefficients of the Laguerre expansion (3) were computed by using (4)–(7) with $\beta >\frac{1}{2}{max}_i{\alpha }_i$. Then, the sequence of $c_k$’s was plotted (Figure 3a) showing that $c_k\to 0$ as $k\to \infty$. After inserting the values of the $c_k$’s into (3), a graph of $f_{Y_6^{\left(1\right)}}^{\left(\beta \right)}\left(y\right)$ was displayed, against the histogram of $Y_6^{\left(1\right)}$ obtained from MC simulations, as shown in Figure 4a. We clearly see from the figure that the histogram representing the PDF of $Y_6^{\left(1\right)}$ perfectly fit the graph of $f_{Y_6^{\left(1\right)}}^{\left(\beta \right)}\left(y\right)$ computed by using the Laguerre expansion (3).

To extend our study to various cases of the linear combination (40) as introduced in Section 3.3, with increasing values of n, we further considered the PDFs of $Y_{11}^{\left(1\right)}={\sum }_{i=1}^{11}{\alpha }_i^{\left(1\right)}{X}_i^{\left(1\right)},Y_{15}^{\left(2\right)}={\sum }_{i=1}^{15}{\alpha }_i^{\left(2\right)}{X}_i^{\left(2\right)}$, and $Y_{20}^{\left(3\right)}={\sum }_{i=1}^{20}{\alpha }_i^{\left(3\right)}{X}_i^{\left(3\right)}$. The values of the parameters ${\alpha }_i^{\left(j\right)},{\nu }_i^{\left(j\right)}$, and ${\delta }_i^{\left(j\right)}$ for $j=1,2,3$ computed by using (42)–(44) are plotted in Figure 1b–d. By following the procedure as previously described for determining the PDF of $Y_6^{\left(1\right)}$, we thus obtained the sequences of $c_k$’s as shown in Figure 3b–d along with the PDFs of $Y_{11}^{\left(1\right)},Y_{15}^{\left(2\right)}$, and $Y_{20}^{\left(3\right)}$, denoted by $f_{Y_{11}^{\left(1\right)}}^{\left(\beta \right)}\left(y\right)$, $f_{Y_{15}^{\left(2\right)}}^{\left(\beta \right)}\left(y\right)$, and $f_{Y_{20}^{\left(3\right)}}^{\left(\beta \right)}\left(y\right)$, as shown in Figure 4b–d, respectively. It is clearly seen from Figure 4b–d that the graphs of the PDFs markedly matched their corresponding histograms obtained from MC simulations.

Next, we illustrate the fitness as previously discussed by employing the Kolmogorov–Smirnov (K-S) test [42]. Figure 4a–d also display the p-values for the K-S tests computed from 1000 random samples generated by the PDFs $f_{Y_6^{\left(1\right)}}^{\left(\beta \right)}\left(y\right)$, $f_{Y_{11}^{\left(1\right)}}^{\left(\beta \right)}\left(y\right)$, $f_{Y_{15}^{\left(2\right)}}^{\left(\beta \right)}\left(y\right)$, and $f_{Y_{20}^{\left(3\right)}}^{\left(\beta \right)}\left(y\right)$ and their corresponding random variables $Y_6^{\left(1\right)}$, $Y_{11}^{\left(1\right)}$, $Y_{15}^{\left(2\right)}$, and $Y_{20}^{\left(3\right)}$, respectively. As shown in Figure 4a–d, the p-values obtained fell into the area of the acceptance region in which we set the significant levels for the K-S tests to be $0.1$; all the resulting null hypotheses of the K-S tests were accepted. Therefore, the PDFs $f_{Y_6^{\left(1\right)}}^{\left(\beta \right)}\left(y\right)$, $f_{Y_{11}^{\left(1\right)}}^{\left(\beta \right)}\left(y\right)$, $f_{Y_{15}^{\left(2\right)}}^{\left(\beta \right)}\left(y\right)$, and $f_{Y_{20}^{\left(3\right)}}^{\left(\beta \right)}\left(y\right)$ computed by using the Laguerre expansion (3) appeared consistent with their corresponding histogram of random samples generated from $Y_6^{\left(1\right)}$, $Y_{11}^{\left(1\right)}$, $Y_{15}^{\left(2\right)}$, and $Y_{20}^{\left(3\right)}$, respectively.

5.2. The Performance of Our Explicit Formula for $\mathbb{E}\left[Y_n^{\gamma }\right]$

Example 2.

In our next example, we demonstrate the performance of our explicit Formula (8) by selecting $Y_{11}^{\left(1\right)}$, $Y_{15}^{\left(2\right)}$, and $Y_{20}^{\left(3\right)}$ given in Example 1 to be our case study.

Firstly, we computed the values of $\mathbb{E}\left[{\left(Y_{11}^{\left(1\right)}\right)}^{\gamma }\right]$ for $\gamma \in [0,2]$ by using our explicit Formula (8) with $K=K_1=100$. Next, we plotted these values against the results obtained from MC simulations with $n_p=10,20,50,100,500,2000$, which were the numbers of sample paths used in the MC simulations, in Figure 5a–f, respectively. Next, we similarly applied this procedure to investigate the accuracy of our explicit Formula (8) for $\mathbb{E}\left[{\left(Y_{15}^{\left(2\right)}\right)}^{\gamma }\right]$ and $\mathbb{E}\left[{\left(Y_{20}^{\left(3\right)}\right)}^{\gamma }\right]$. The results obtained for $\mathbb{E}\left[{\left(Y_{15}^{\left(2\right)}\right)}^{\gamma }\right]$, where $\gamma =2,2.01,\dots ,3$, with $K=K_2=300$ and $\mathbb{E}\left[{\left(Y_{20}^{\left(3\right)}\right)}^{\gamma }\right]$, where $\gamma =3,3.01,\dots ,4$, with $K=K_3=400$, are displayed in Figure 6a,c,e and Figure 6b,d,f, respectively. Evidently, the variation of the approximate values from MC simulations decreased when $n_p$ increased for all chosen γ’s, demonstrating the convergence of the approximate values from MC simulations to the one computed by using our explicit Formula (8).

In order to verify our result presented in Theorem 3, we used (11) to compute the truncation errors of $\mathbb{E}\left[{\left(Y_{11}^{\left(1\right)}\right)}^{\gamma }\right]$, $\mathbb{E}\left[{\left(Y_{15}^{\left(2\right)}\right)}^{\gamma }\right]$, and $\mathbb{E}\left[{\left(Y_{20}^{\left(3\right)}\right)}^{\gamma }\right]$, denoted by $E_{k,\infty }^{(\gamma ,1)}$, $E_{k,\infty }^{(\gamma ,2)}$, and $E_{k,\infty }^{(\gamma ,3)}$, respectively, when $\gamma =\frac{1}{2}\notin \mathbb{N}$ and $\gamma =1,2\in \mathbb{N}$ for $k=0,\dots ,10$. The truncation errors obtained are tabulated in

Table 1. The truncation errors of $\mathbb{E}\left[{\left(Y_{11}^{\left(1\right)}\right)}^{\gamma }\right]$, $\mathbb{E}\left[{\left(Y_{15}^{\left(2\right)}\right)}^{\gamma }\right]$, and $\mathbb{E}\left[{\left(Y_{20}^{\left(3\right)}\right)}^{\gamma }\right]$, denoted by $E_{k,\infty }^{(\gamma ,1)}$, $E_{k,\infty }^{(\gamma ,2)}$, and $E_{k,\infty }^{(\gamma ,3)}$, respectively, computed in Example 2 by using (11) when $\gamma =\frac{1}{2}\notin \mathbb{N}$ and $\gamma =1,2\in \mathbb{N}$ for $k=0,\dots ,10$.

k	$\mathit{\gamma }=1/2$			$\mathit{\gamma }=1$			$\mathit{\gamma }=2$
	${\mathit{E}}_{\mathit{k},\infty }^{(\mathit{\gamma },\mathbf{1})}$	${\mathit{E}}_{\mathit{k},\infty }^{(\mathit{\gamma },\mathbf{2})}$	${\mathit{E}}_{\mathit{k},\infty }^{(\mathit{\gamma },\mathbf{3})}$	${\mathit{E}}_{\mathit{k},\infty }^{(\mathit{\gamma },\mathbf{1})}$	${\mathit{E}}_{\mathit{k},\infty }^{(\mathit{\gamma },\mathbf{2})}$	${\mathit{E}}_{\mathit{k},\infty }^{(\mathit{\gamma },\mathbf{3})}$	${\mathit{E}}_{\mathit{k},\infty }^{(\mathit{\gamma },\mathbf{1})}$	${\mathit{E}}_{\mathit{k},\infty }^{(\mathit{\gamma },\mathbf{2})}$	${\mathit{E}}_{\mathit{k},\infty }^{(\mathit{\gamma },\mathbf{3})}$
0	$6.2\times {10}^{-1}$	$1.2\times {10}^0$	$1.6\times {10}^0$	$3.5\times {10}^0$	$7.6\times {10}^0$	$1.2\times {10}^1$	$5.6\times {10}^1$	$1.6\times {10}^2$	$3.8\times {10}^2$
1	$6.5\times {10}^{-2}$	$1.9\times {10}^{-1}$	$2.8\times {10}^{-1}$	0	0	0	$1.3\times {10}^1$	$5.8\times {10}^1$	$1.5\times {10}^2$
2	$1.2\times {10}^{-2}$	$5.5\times {10}^{-2}$	$9.1\times {10}^{-2}$	0	0	0	0	0	0
3	$3.2\times {10}^{-3}$	$2.0\times {10}^{-2}$	$3.5\times {10}^{-2}$	0	0	0	0	0	0
4	$9.1\times {10}^{-4}$	$7.7\times {10}^{-3}$	$1.5\times {10}^{-2}$	0	0	0	0	0	0
5	$2.8\times {10}^{-4}$	$3.2\times {10}^{-3}$	$6.8\times {10}^{-3}$	0	0	0	0	0	0
6	$9.2\times {10}^{-5}$	$1.4\times {10}^{-3}$	$3.2\times {10}^{-3}$	0	0	0	0	0	0
7	$3.1\times {10}^{-5}$	$6.3\times {10}^{-4}$	$1.5\times {10}^{-3}$	0	0	0	0	0	0
8	$1.1\times {10}^{-5}$	$2.9\times {10}^{-4}$	$7.7\times {10}^{-4}$	0	0	0	0	0	0
9	$4.1\times {10}^{-6}$	$1.3\times {10}^{-4}$	$3.9\times {10}^{-4}$	0	0	0	0	0	0
10	$1.6\times {10}^{-6}$	$6.7\times {10}^{-5}$	$1.9\times {10}^{-4}$	0	0	0	0	0	0

. We clearly see from Table 1 that $E_{k,\infty }^{(\gamma ,j)}$ tended to zero when k increased for all $j=1,2,3$ and selected γ. In particular, $E_{k,\infty }^{(\gamma ,j)}=0$ for $k\ge \gamma$ when $\gamma =1,2$. This confirmed our result presented in Theorem 4 that our explicit Formula (17) could be used to compute $\mathbb{E}\left[Y_n^{\gamma }\right]$ without producing truncation errors when $\gamma \in \mathbb{N}$; it should be remarked from Table 1 that the truncation errors could be very large when $K+1$ was less than γ for $\gamma \in \mathbb{N}$. Although utilizing our explicit Formula (8) for computing $\mathbb{E}\left[Y_n^{\gamma }\right]$ when $\gamma \in {\mathbb{R}}_{+}$ and $\gamma \notin \mathbb{N}$ always produced a truncation error, our result presented in Theorem 3 ensured that the truncation error tended to zero as K increased.

We finish this example by illustrating the efficiency of our explicit Formula (8) over the MC simulations. As shown in Figure 5a–f and Figure 6a–f, we needed to increase the value of $n_p$ in the MC simulations in order to reduce the variations on the approximate values of $\mathbb{E}\left[{\left(Y_{11}^{\left(1\right)}\right)}^{\gamma }\right]$, $\mathbb{E}\left[{\left(Y_{15}^{\left(2\right)}\right)}^{\gamma }\right]$, and $\mathbb{E}\left[{\left(Y_{20}^{\left(3\right)}\right)}^{\gamma }\right]$, respectively, which could be time-consuming. For example, to yield an absolute difference between the exact value of $\mathbb{E}\left[{\left(Y_{20}^{\left(3\right)}\right)}^2\right]$ computed from our explicit Formula (8) and its approximate value obtained from MC simulations to be less than ${10}^{-3}$, our numerical experiment required MC simulations with $n_p={10}^7$ consuming 19 s, while implementing our explicit Formula (8) with $K=3$ took just $0.001$ s.

5.3. Extended Results for the ECIR Process with Time-Varying Dimension

We next shift our attention to our explicit formulas for the TPDF and the moments of ECIR processes with time-varying dimension. We specifically set the parameter functions of the ECIR process $V_t$ described by (19) with time-varying dimension $\mathbf{d}\left(t\right)$ defined in (20) as follows:

$$\kappa \left(t\right)=0.1+0.2t+0.3e^{{cos}^2(t+2)},$$

(45)

$$\theta \left(t\right)=0.1+0.5e^{2sin(t+2)},$$

(46)

and

$$\sigma \left(t\right)=\left(0.4+0.1t\right)e^{sin(t+2)}$$

(47)

for $t\in [0,3]$.

Figure 7a–f display variations of the three parameter functions (45)–(47) as well as $\mathbf{d}\left(t\right)$, ${\mathbf{d}}^{\left(1\right)}\left(t\right)$, and $\tau (t,0)$ defined in (21), respectively, as given in the figure. It should be noticed from Figure 7d that Assumption 1 was fulfilled on $[0,3]$, while Assumption 2 was violated, e.g., ${\mathbf{d}}^{\left(1\right)}\left(1.5\right)<0$, as shown in Figure 7e. Consequently, we set the time domain for this study to be ${\mathbf{D}}_1:=[0,1]$, ensuring that both Assumptions 1 and 2 were satisfied.

5.3.1. The Accuracy of Our Explicit Formula for the TPDF of the ECIR Process with Time-Varying Dimension

Example 3.

This example aimed to investigate the accuracy of our explicit Formula (31) for calculating the TPDF of $V_t$. We started by implementing the result presented in Theorem 10. Consider the convergence of the PDF of ${\widehat{Y}}_n(t,v_0):={\sum }_{i=1}^n{\widehat{\alpha }}_i{\widehat{X}}_i$ to the PDF of $V_t|v_0$ as n approaches infinity, written in (22). We set $v_0=1,2,$ and $t=0.1,0.5,1\in {\mathbf{D}}_1$. Then, we computed the approximate PDFs of ${\widehat{Y}}_n(t,v_0)$ for $n=1,2,3,$ based on the random samples drawn from populations distributed according to noncentral chi-square and chi-square distributions with a number of sample $n_p={10}^6$. On the other hand, as a benchmark, the approximate PDFs of $V_t|v_0$ for all $v_0=1,2,$ and $t=0.1,0.5,1$ were obtained by using the sample paths generated from (19) as follows.

Figure 8a–f show that the PDF of ${\widehat{Y}}_n(t,v_0)$ tended to align better with the histogram of $V_t|v_0$ representing the PDF of $V_t|v_0$ when n increased for all $v_0=1,2,$ and $t=0.1,0.5,1$, demonstrating ${\widehat{Y}}_n(t,v_0)$ converged in distribution to $V_t|v_0$ as n approached infinity. However, as shown in Figure 8a–f, the number of terms used in the summation ${\sum }_{i=1}^n{\widehat{\alpha }}_i{\widehat{X}}_i$ had to be increased in order to obtain a better approximation for the PDF of $V_t|v_0$. This is a major drawback of implementing (22), which requires the exact PDFs of ${\widehat{Y}}_n(t,v_0)$ and $V_t|v_0$ in order to estimate errors occurring for all n.

The problem mentioned above can completely be solved by employing our explicit Formula (31) for obtaining the exact PDF of $V_t|v_0$. Truncation errors occurring when implementing the infinite series in (31) can be estimated by applying Lemma 1. To demonstrate the accuracy of our explicit Formula (31), we computed $f_{V_t}(v,t|v_0)$ for $v_0=1,2,$$t=0.1,0.5,1$ and $v>0$ by setting $K=20$. Sequences of $\left|{\widehat{c}}_k(t,v_0)\right|$’s computed from (32) and (33) for $v_0=1,2$ and $t=0.1,0.5,1$ are displayed in Figure 9a,b, showing that the coefficients tended to zero when k increased. This ensured the truncation errors vanished when K approached infinity.

Figure 8a–f also display the graphs of $f_{V_t}(v,t|v_0)$ against the corresponding histograms of $V_t|v_0$ for $t=0.1,0.5,1\in {\mathbf{D}}_1$ and $v_0=1,2,$ obtained from the sample paths generated from (19). It is readily seen that the graphs of $f_{V_t}(v,t|v_0)$ gracefully matched the corresponding histograms. Following the K-S tests with the significance level of 10% employed in Example 1 to determine the equivalence between two distributions, the minimum of the p-values was $0.14$, implying that there was no significance difference between $f_{V_t}(v,t|v_0)$ obtained from our explicit Formula (31) and the corresponding histograms of the random samples generated from (19).

The validity of our explicit Formula (31) admittedly becomes questionable when Assumption 2 is violated. To investigate, we extended the time domain from ${\mathbf{D}}_1:=[0,1]$ to ${\mathbf{D}}_2:=[0,3]$ and considered $f_{V_t}(v,t|v_0)$ for $v>0$, $v_0=1,2$, and $t=0.03,0.75,1.5,2.25,3\in {\mathbf{D}}_2$. Figure 7e visually shows that at $t=1.5,3$, ${\mathbf{d}}^{\left(1\right)}\left(1.5\right)$ and ${\mathbf{d}}^{\left(1\right)}\left(3\right)$ were negative and hence Assumption 2 was violated, but $f_{V_t}(v,t|v_0)$ was finite for all $v>0$, $v_0=1,2$, and $t\in {\mathbf{D}}_2$, as shown in Figure 10a,b.Moreover, the graph of $f_{V_t}(v,t|v_0)$ appears to match the corresponding histograms of the random samples from (19) based on MC simulations because the minimum of the p-values from the K-S tests was $20\%$. The results obtained suggested that the condition in Assumption 2 could perhaps be replaced with ${\mathbf{d}}^{\left(1\right)}\left(t\right)<\infty$ for all $t\in [0,T]$ when our explicit Formula (31) was employed for computing the TPDF of the ECIR process (19). That being said, Assumption 2 remains an important ingredient of our analysis since the degrees of freedom ${\widehat{\nu }}_i$ given in (27) are not allowed to be negative from the definitions of noncentral chi-square and chi-square distributions.

5.3.2. The performance of Our Explicit Formula for $U_E^{\left(\gamma \right)}\left(t\right|v_0)$

Example 4.

In our last example, we illustrate the accuracy and efficiency of our explicit Formula (38) for computing the ${\gamma }^{\mathrm{th}}$ conditional moment of the ECIR process (19) with time-varying dimension. As previously discussed in Section 4.2.1, we demonstrate the advantages of using our explicit Formula (38) over the explicit formula presented in Theorem 2.1 of Rujivan [12] by setting the parameter functions of the ECIR process (19) to follow (45)–(47) as used in Example 3.

Firstly, we considered the accuracy of our explicit Formula (38) and the one written in Equation (2.2) of Rujivan [12] by setting $\gamma =0.5,1,1.5,2,2.5,3.$ Let $v_0=1$ and $t=1$. For each γ, we computed two sequences of $U_E^{(\gamma ,K)}\left(t\right|v_0)$ for $K=0,\dots ,5$, from our explicit formulas (38) and the one written in Equation (2.2) of Rujivan [12] with the number of terms, $K+1$, used in the infinite series on the RHS of (38) and Equation (2.2) of Rujivan [12]. Figure 11a–f display the graphs of the two sequences of $U_E^{(\gamma ,K)}\left(t\right|v_0)$ against the approximate value of $U_E^{\left(\gamma \right)}\left(t\right|v_0)$ from MC simulations based on the ECIR process (19), demonstrating that the two sequences of $U_E^{(\gamma ,K)}\left(t\right|v_0)$ converged to the corresponding approximate value of $U_E^{\left(\gamma \right)}\left(t\right|v_0)$ obtained from MC simulations when K increased.

It should also be pointed out from Figure 11a–f that the sequence of $U_E^{(\gamma ,K)}\left(t\right|v_0)$ computed from (38) converged to the corresponding approximate value of $U_E^{\left(\gamma \right)}\left(t\right|v_0)$ obtained from the MC simulations faster than the sequence of $U_E^{(\gamma ,K)}\left(t\right|v_0)$ computed from Equation (2.2) of Rujivan [12]. Moreover, as shown in Figure 11b,d,f, the two sequences of $U_E^{(\gamma ,K)}\left(t\right|v_0)$ coincided for $K+1\ge m$ when $\gamma =m$ was a positive integer, demonstrating the consistency of our closed-from Formula (39) and the explicit formula written in Equation (2.13) in Theorem 2.2 of Rujivan [12].

Secondly, we investigated the efficiency of our explicit Formula (38) over the one written in Equation (2.2) of Rujivan [12] for calculating $U_E^{\left(\gamma \right)}\left(t\right|v_0)$ by choosing $\gamma =0.5$ in our case study. Figure 12a,b illustrate the computational times (seconds) used to compute $U_E^{(\gamma ,K)}\left(t\right|v_0)$ from the two explicit formulas and the reduction (in folds) of computational time for $K=1,\dots ,10$, respectively. As expected, implementing our explicit Formula (38) consumed considerably less time and effort than implementing the explicit formula written in Equation (2.2) of Rujivan [12], in particular, the reduction was more than sixfold when $K\ge 10$.

6. Conclusions

In this paper, we presented the first explicit formula expressed in terms of generalized hypergeometric functions for computing the ${\gamma }^{\mathrm{th}}$ moment of a conic combination of n independent noncentral chi-square random variables defined in (1), when the conic coefficients were not all identical for any integer $n\ge 2$ and real number $\gamma \in {\mathbb{R}}_{+}$. Moreover, the truncation errors occurring by implementing our explicit formulas were determined analytically. We extended our result to various types of random variables which were independent and could be transformed to noncentral chi-square random variables. For validation purposes, several numerical examples were presented to show the performance of our explicit formulas compared with MC simulations. Furthermore, we highlighted an interesting application of our explicit formulas in interest rate modeling by expressing the exact TPDF of the ECIR process with time-varying dimension in terms of generalized Laguerre functions. As a result, a novel explicit formula for the ${\gamma }^{\mathrm{th}}$ conditional moment of the ECIR process was obtained and tested, and we concluded that the distinguished feature of our current analytical approach lay in its computational efficiency, which was superior to that of the other existing methods from the literature.