A Novel Analytical Formula for the Discounted Moments of the ECIR Process and Interest Rate Swaps Pricing

: This paper presents an explicit formula of conditional expectation for a product of polynomial functions and the discounted characteristic function based on the Cox–Ingersoll–Ross (CIR) process. We also propose an analytical formula as well as a very efﬁcient and accurate approach, based on the ﬁnite integration method with shifted Chebyshev polynomial, to evaluate this expectation under the Extended CIR (ECIR) process. The formulas are derived by solving the equivalent partial differential equations obtained by utilizing the Feynman–Kac representation. In addition, we extend our results to derive an analytical formula of conditional expectation of a product of mixed polynomial functions and the discounted characteristic function. The accuracy and efﬁciency of the proposed scheme are also numerically shown for various modeling parameters by comparing them with those obtained from Monte Carlo simulations. In addition, to illustrate applications of the obtained formulas in ﬁnance, analytical pricing formulas for arrears and vanilla interest rate swaps under the ECIR process are derived. The pricing formulas become explicit under the CIR process. Finally, the fractional ECIR process is also studied as an extended case of our main results.


Introduction
A conditional expectation has been widely used in many branches of science. It can be easily computed if the probability density function (PDF) is known. However, often, the densities are unknown or are in complicated forms. Let (Ω, F t , {F t } 0≤t≤T , Q) be a filtered probability space generated by an adapted stochastic process {r t } 0≤t≤T , where Ω is a sample space, Q is a risk-neutral measure and the family {F t } 0≤t≤T of σ-field on Ω parametrized by t ∈ [0, T] is a filtration. This paper focuses on the conditional expectation of a nonlinear function of the form: whose analytical formula has not been discovered, when α, λ, β, γ ∈ R and r t evolve according to the extended Cox-Ingersoll-Ross (ECIR) process [1] governed by the stochasic differential equation (SDE): where W t is a Wiener process. The parameter θ corresponds to the speed of adjustment to the mean of the invariant distribution µ and the parameter σ determines the state space of the diffusion and the shape of the invariant distribution. The Cox-Ingersoll-Ross (CIR) process, which was first introduced by Feller [2], is a special case of the ECIR process when The rest of the paper is organized as follows. Section 2 gives a brief overview of the ECIR process. The key methodology as well as main theorems are provided in Section 3. Section 4 introduces the numerical method for numerically evaluating the analytical formulas and provides numerical validations for the formulas. The formulas for pricing the arrears and vanilla swaps are presented in Section 5. Section 6 provides a study case of the fractional ECIR process adopt the idea from our main results given in Section 3. Section 7 presents the conclusions and discussions.

The Extended Cox-Ingersoll-Ross Process
For the ECIR process (2), in order to confirm that there exists a path-wise unique strong solution for the process r t in (2) and to avoid zero almost everywhere with respect to the probability measure Q for all t ∈ [0, T], the two following assumptions studied by Maghsoodi [5] are needed (see details in Theorems 2.1 and 2.4 of their work). Altogether, we need the following sufficient condition. Assumption 1 ([5]). The time dependent parameters θ(t), µ(t) and σ(t) in the ECIR process (2) are smooth and strictly positive, µ(t) σ 2 (t) is locally bounded on [0, T] and 2θ(t)µ(t) ≥ σ(t) 2 .
It is worth noting that the CIR process is a special case of the ECIR process in which the parameters are constants and the Assumption 1 still holds for this case.
To achieve our aim, one fundamental question arises: Why do we not use the CIR process's transition PDF directly? It is known that its transition PDF has the expression in terms of Bessel function of the first kind of order q, see [19,20], that is: where τ = T − t, c τ = 2θ σ 2 (1−e −θτ ) , u = c τ r t e −θτ , vs. = c τ r, q = 2θµ σ 2 − 1 and I q (·) is the Bessel function of the first kind of order q. As displayed above, since the CIR process's transition PDF is very complicated, the closed-form formulas for the conditional expectation (1) by applying this transition PDF are unavailable or complicated as well.
The situation becomes even more complicated for the ECIR process case; for instance, the ECIR(d) process presented by Egorov et al. [21]. Its dynamics are governed by the following time-inhomogeneous diffusion process: where θ, σ 0 are positive, σ 1 is a real and d is positive constant with the transition PDF Here, λ = r t v, G = e θτ v, v = 8σ 1 σ 2 0 e −θτ e 2σ 1 T − e 2σ 1 t , τ = T − t and again I q (·) is the Bessel function of the first kind. This transition PDF was first proposed by Maghsoodi [5]. To avoid using the transition PDF for calculating (1), we employ the Feynman-Kac representation. This representation offers a method of solving a conditional expectation of an Itô random process by deterministic methods. In fact, it provides a relation between an Itô random process and a PDE.

Main Results
In this section, we present the methodology used in this paper as well as the main results. Theorem 1 gives an analytical formula for the conditional expectation of a product of polynomial and exponential functions, E Q r γ T e −λr T − T t (αr s +β) ds | r t = r where r t follows the ECIR process and α, λ, β, γ ∈ R. Theorem 2 is a special case of Theorem 1 when γ is a non-negative integer. In such case, the infinite sum, which can potentially cause a truncation error in practice, can be reduced into a finite sum. It should be mentioned that our proposed formulas for the ECIR process are more general and cover the results given in [14][15][16]22]. The analytical formulas presented in Theorems 1 and 2 involve the computation of the solution of the Riccati differential equation. As the solution of such equation does not have an analytical form except for some special cases, such as when α = 0, we propose a very efficient numerical approach to solve the equation in Section 4. Theorem 1. Suppose that r t follows the ECIR process (2) with α, λ, β, γ ∈ R. Let 0 ≤ t ≤ T. Then, The function B is obtained by solving the Riccati differential equation Proof. Under the uniformly convergent assumption, the Feynman-Kac representation is applied to solve U := U γ E (r, τ) in (3) which satisfies the corresponding PDE [23], where the initial condition at τ = 0 is given by: Then, we compare the coefficients in (3) and (8) to obtain the conditions B(0) = −λ, A γ 0 (0) = 1 and A γ j (0) = 0 for j ∈ N. Next, we find the partial derivatives U τ , U r and U rr by using (3) which are: Now, these above partial derivatives and U in (3) are substituted into (7) by letting A j := A γ j (τ) and B := B(τ). Thus, (7) becomes: which can be simplified by employing (5) as: Consider (9) as a power series in r. Since r = 0 and the leading coefficient A 0 = 0, (9) holds under the following nonlinear differential equation: with the initial condition B(0) = −λ and the linear differential equations: with their initial conditions A 0 (0) = 1 and A j (0) = 0 for j ∈ N. This completes the proof and the solutions of those linear differential equations are given in (4).
Proof. Let γ = n ∈ N ∪ {0}. Consider (4) at the index j = n + 1, we get A n n+1 (τ) = 0 due to Q n+1 = 0 by (5). From the recurrence relation (4), we have A n j (τ) = 0 for all integers j ≥ n + 1. Consequently, the infinite sum (3) can be reduced to the finite sum (10). Theorem 3 shows that the formulas can be expressed in closed forms under the CIR process where all parameters θ(t) = θ, µ(t) = µ and σ(τ) = σ are constants. Theorem 4 extends the result (10) in Theorem 2 to derive an analytical formula for a product of two polynomial and one exponential functions E Q r n 1 t (αr u +β) du | r t = r for some constants α, β ∈ R and non-negative integers n 1 , n 2 .
From the result presented in (4) for j ∈ N, we obtain: Under the uniformly convergent assumption, the proof is completed.
Theorem 4. Suppose that r t follows the ECIR process (2) with β, γ ∈ R and n 1 , is a solution of the Riccati differential Equation (6) with the initial condition B(0) = x, and A γ j (τ; B) are coefficients described in (4).

Remark 1.
We can see that various statistics such as the first and second conditional moments, variances, mixed moments, covariances and correlations, can be derived as special cases of Theorem 4 where λ = α = β = 0. This corresponds to Rujivan's result [14]. As the maximum likelihood cannot be solved directly under ECIR, conditional moments are essential for parameter estimating methods such as the martingale estimating functions, generalized method of moments, and quasilikelihood method. Corollary 1. Suppose that r t follows the ECIR process (2), n ∈ N and τ = T − t ≥ 0, we have where V n E (r, τ) denotes U n E (r, τ) when λ = α = β = 0. Note that, without the exponential term, this corresponds to central moments and is a conditional variance for n = 2.
Proof. It is obtained directly by using the Binomial theorem together with Theorem 2.

Remark 2.
Results above can be extended to approximate the expectation of any real analytic functions f of a random variable on an open interval, where its Taylor expansion at x 0 converges pointwise to f . For example, consider f (x) = √ x which can be expressed using Taylor expansion as: By substituting x = r T and x 0 = E Q [r T | r t = r] where r t follows the ECIR process for 0 ≤ t ≤ T and taking the conditional expectation, we have that The above expansion can be applied to other functions, such as exponential, logarithmic or trigonometric functions. However, one needs to make sure that the approximation converges as there is no empirical analysis providing the convergence condition of the Taylor expansions in a general stochastic context. Moreover, there is empirical evidence indicating that increasing the number of terms in Taylor approximation does not necessarily yield more accurate estimation. The classic example is given by Zhu and Lian [24]. They showed that a better accuracy of the convexity correction approximation (CCA), used to estimate volatility swap prices, cannot be achieved by extending the second-order Taylor expansion to that of the third order and gave the condition when CCA provides decent estimates.
Proof. By using the definition of covariance, we have: Remark 3. The analytical form of E Q r T e − T s (αr u +β) du | r t = r can be easily obtained using the tower property and Theorem 2.

Numerical Procedures
The analytical Formulas (3) and (17) in Theorems 1 and 4, respectively, used for pricing swaps under the ECIR process in Section 3, require the evaluation of the two parameters B(τ) and A γ j (τ) which are the functions of time variable τ. As we know, under the ECIR process, (6) cannot be solved for analytical solution. Thus, the analytical form of the two parameters cannot be obtained. In this section, we propose a very efficient scheme to numerically evaluate the formulas. The approach consists of two parts. The first part is to numerically solve the Riccati differential Equation (6) using the FIM base on the Chebyshev polynomials [25]. Then, the algorithms to evaluate the formulas are illustrated.

FIM with Shifted Chebyshev Polynomial
This subsection introduces the major tool for solving the Riccati differential equation, that is, the FIM based on the shifted Chebyshev polynomial. This numerical scheme was proposed by Boonklurb et al. [25] in 2018, which provided a very accurate solution using less computational nodes and, thus, an inexpensive consuming time. However, we adjust this method slightly by employing the shifted Chebyshev polynomials, in order to be easily applied to our considering domain of the Riccati differential equation, which is [0, T] where T > 0.
Definition 1 describes the shifted Chebyshev polynomial. Its properties needed to construct the first and higher orders of the Chebyshev integration matrices, which are the main ingredients of the FIM, are detailed in Lemma 1.
(iii) The shifted Chebyshev matrix S at each node τ k defined by (18) is Then, its multiplicative inverse is Next, the first order Chebyshev integration matrix is constructed. Let n ∈ N be a number of computational nodes. Assume that f (τ) is an approximate solution of any function which is defined by the linear combination of the shifted Chebyshev polynomials S 0 (τ), S 1 (τ), S 2 (τ), . . . , S n (τ). Then, we have: where c i is an unknown coefficient to be considered. Let τ k be computational nodes that are generated by the zeros of the shifted Chebyshev polynomial S n+1 defined by (18) in ascending order. After that, each node τ k is plugged into (19). Then, we obtain the following matrix form: which is denoted by f = Sc. By Lemma 1(iii), we know that S is invertible. Thus, we get the coefficient c = S −1 f. Now, we consider a single-layer integration of f from 0 to τ k denoted by F(τ k ), we obtain: . . , n} are defined in Lemma 1(ii). As we substitute each nodal point τ k into the above equation, it can be written in the matrix form: matrix representation for the integral operator called the Chebyshev integration matrix. Thus, the single-layer integration F(τ k ) can be also expressed to another form as follows:

Numerical Procedure for Theorem 1
Next, we utilize the FIM with the shifted Chebyshev polynomial to seek the numerical solution for the Riccati differential Equation (6) which is the initial value problem. First of all, let us recall the Riccati differential equation again, with initial condition B(0) = −λ. To find its numerical formula, we discretize the domain [0, T] into n subintervals composed n + 1 nodes. These nodes are manufactured by the zeros of the shifted Chebyshev polynomial S n+1 defined in (18), that are τ 0 < τ 1 < τ 2 < · · · < τ n . By the idea of FIM, we have to transform the differential equation into the equivalent integral equation. Thus, to remove the derivative from (21), we take a single integral from 0 to τ k on both sides of (21). Then, we have: where C is an arbitrary constant emerged in the process of integration. Subsequently, we apply the FIM based on the shifted Chebyshev polynomials to approximate the integral terms contained in (22) by using (20). We obtain: As we vary the value of τ k for k ∈ {0, 1, 2, . . . , n} by the zeros (18) into the above equation, we have the system of nonlinear equations which can be expressed to the matrix form: It can be simplified and denoted by: where J = SS −1 is the Chebyshev integration matrix described in the previous section, the notation is the Hadamard product defined in [27] which means a product of elements in matrices at the same position and the other parameters contained in (23) are defined by the following: . . , τ n and e = 1, 1, 1, . . . , 1 has n + 1 components.
However, we see that (23) is in the matrix form of a nonlinear equation. Hence, the technique of linearization method is applied to (23) in order for convenient solving. Let m ∈ N, we then take the (m − 1)th and mth iterations into the variable B, which are, respectively, denoted by B m−1 and B m , for the nonlinear term in (23). The other terms are determined at the mth present iteration. Thus, (23) becomes which can be rearranged to where I is the (n + 1) × (n + 1) identity matrix and diag B m−1 is the diagonal matrix that diagonal entries are the elements of vector B m−1 . Moreover, we observe that C is the new unknown that occurred from the integration. Thus, we need to create one more equation by hiring the given initial condition B(0) = −λ. Accordingly, we use (19) to transform this supplementary condition into the matrix form. Thus, at the mth iteration, we obtain: where u = [S 0 (0), S 1 (0), S 2 (0), . . . , S n (0)] . Finally, we combine (24) and (25) to construct the iterative system of linear equations which has n + 2 unknowns as follows: Consequently, the iterative approximate solution B m can be sought by solving (26) in conjunction with an arbitrary initial guess of the iteration B 0 that makes the coefficient matrix to be invertible. Note that the stopping criterion for finding the mth iterative approximate solution B m is stopped when the mth error Euclidean norm m := B m − B m−1 of difference between the current and consecutively previous solutions is less than the given convergent tolerance TOL. Therefore, an approximate solution B(τ) at arbitrary point τ ∈ [0, T] can be estimated by: where and B m is the mth terminal solution obtained by solving (26). Now, we already have the approximate solutions B(τ) of the Riccati differential Equation (6) that will be used to estimate the coefficients A γ j (τ) for j ∈ N ∪ {0} contained in (4). Furthermore, we can also use the Chebyshev integration matrix J to represent the integral terms in (4). Thus, by applying (20) to (4), we have the estimated coefficients at the temporal variable τ k as follows: As the variables τ k for k ∈ {0, 1, 2, . . . , n} are varied as the zeros (18) into these above two equations for A γ 0 and A γ j , we obtain: where J k: is the kth row vector of the Chebyshev integration matrix J and the other parameters contained in the above expressions (28) are defined by: The elements of both P j and Q j can be directly calculated by (5). Note that the exponential vector exp {·} means the vector where its element is the exponential of each component in that position. In addition, we can use the obtained coefficient vectors A γ 0 and A γ j in (28) to approximate A γ j (τ) for j ∈ N ∪ {0} at any temporal variable τ ∈ [0, T] by applying the same idea as in (27). Then, we have: where v and S are defined as same as in (27) and A γ j for j ∈ N ∪ {0} can be directly found by computing (28). Now, we already have the approximate values of B(τ) and A γ j (τ) which are produced from solving (27) and (29), respectively. Finally, we can compute the numerical formula of (3) in Theorem 1 for approximating the arbitrage price by the following formula where M is the given maximum index. In practice, when this maximum index M is too large, it may result in a divergence of solution (30). Hence, from Theorem 2, we have proposed that if γ ∈ N ∪ {0}, the infinite sum can be reduced to a finite sum where M = γ. As a consequence, our numerical Formula (30) is well suitable to apply under this condition. Eventually, we validate our numerical Formula (30) and examine the accuracy by comparing it with the Monte Carlo (MC) simulations in the next section. For computational convenience, we summarize the above-mentioned procedure for finding the numerical solution of arbitrage price in the algorithm form as demonstrated by the flowchart provided in Figure 1.

Numerical Procedure for Theorem 4
Additionally, we can also adapt the concept of the numerical procedure using FIM based on Chebyshev polynomials in Section 4.1 to approximate the Formula (17) in Theorem 4. First, we start to solve the numerical solution B(τ; 0) denoted by B 2 (τ) from the following Riccati differential equation: with the initial condition B 2 (0) = 0. By using the FIM based on shifted Chebyshev polynomial, we can transform (31) into the matrix form (26) in which changes −λ to the initial value 0. Then, we obtain: where B2 m = [B 2 (τ 0 ), B 2 (τ 1 ), B 2 (τ 2 ), . . . , B 2 (τ n )] is the Riccati solution vector at mth iteration and the other parameters in (32) are defined in the same way as in Section 4.2. When we run the iterative system (32) until it converges, the solution vector B2 m is obtained. Thus, we can approximate B 2 (τ) at any τ ∈ [0, T 2 ] by using (27) where B2 m is the final iterative solution obtained by (32). Then, B(T 2 ; 0) = B 2 (T 2 ). Continuously, we use this obtained approximate solution at end point, or B 2 (T 2 ), to be the initial condition of B(τ; B 2 (T 2 )) denoted by B 1 (τ) for solving the Riccati differential equation: subject to the initial condition B 1 (0) = B 2 (T 2 ). (33) can be written in the matrix form: where B1 m = [B 1 (τ 0 ), B 1 (τ 1 ), B 1 (τ 2 ), . . . , B 1 (τ n )] is the Riccati solution vector at mth iteration and the other parameters in (34) are defined as in Section 4.2. When the linear system (34) is iterated until it converges, we then have the solution vector B1 m . Hence, the approximate solution B 1 (τ) can be sought at arbitrary τ ∈ [0, T 1 ] by employing (27). Then, we get B 1 (τ) ≈ v S −1 B1 m , where B1 m is the terminal solution after solving (34). Therefore, B(T 1 ; B 2 (T 2 ))) = B 1 (T 1 ). Next, we estimate the coefficients A n 2 j (τ; B 2 (T 2 )) and A n 1 +n 2 −j k (τ; B 1 (T 1 )) for any indices j, k ∈ N ∪ {0}. From Theorem 4, we have known that these coefficients depend on the obtained Riccati solutions B 2 (τ) and B 1 (τ), respectively. Therefore, we can approximate these coefficients by hiring the same idea as in (28). However, we need to slightly modify the variable P j (τ) in (5) to correspond with the Riccati solutions B 1 (τ) and B 2 (τ) denoted byP j (τ). Thus, the coefficients at arbitrary values j, γ ∈ N ∪ {0} can be estimated by: where P j = P j (τ 0 ),P j (τ 1 ),P j (τ 2 ), . . . ,P j (τ n ) in which its elements can be found bỹ Hence, by using (29), we obtain the approximate coefficients A γ j (τ) ≈ v S −1 A γ j for arbitrary value τ. Eventually, the numerical formula for arbitrage price in Theorem 4 can be computed by: For computational convenience, we provide the algorithm for approximating the numerical formula of Theorem 4 in form of the flowchart, as demonstrated in Figure 2.

Numerical Validation
In this section, the validation testing of the obtained formulas in Section 3 is given through comparison with the MC simulations based on the ECIR(d) process proposed by Egorov et al. [21], the transition PDF of which is a consequence studied by Maghsoodi [5], i.e., This process is recalled from Section 2. To compare the processes between (37) and (2), we have θ(t) = θ, µ(t) = dσ 2 (t) 4θ and σ(t) = σ 0 e σ 1 t , where θ, σ 0 and σ 1 are constant real numbers and d ∈ N such that d ≥ 2. These parameters make the functions θ(t), µ(t) and σ(t) satisfy the Assumption 1. Moreover, the process (37) also has the transition PDF provided in Section 2, which may be usable to validate the value of U γ E (r, τ). However, this transition density is not a suitable validation for a small τ because it behaves like the Dirac delta function which has an infinitely thin spike when τ → 0, see [21].
As a result, it also produces an inaccurate solution U γ E (r, τ). Therefore, for examining our formulas, we choose to employ the MC simulation to approximate them instead.
In these experiments, we also apply the technique of Euler-Maruyama (EM) discretization combined with the MC simulations to the process (37). This numerical scheme has been used to solve the ECIR process provided by Higham and Mao [28]. By EM method, we letr be a time-discretized approximation of r and discretize the time interval [0, T] into N steps that are 0 = t 0 < t 1 < t 2 < · · · < t N = T. Thus, the EM approximate is defined by: with the initial valuer t 0 = r t 0 , where ∆t = t i+1 − t i is the time step and Z i is a standard normal random variable. Next, we illustrate the validations of the formulas via two examples under the CIR and ECIR processes. According to (37), we select the parameters d = 2, θ = 1, σ 0 = 0.01 and σ 1 = 0, 1 that hold for Assumption 1. In addition, we can see that if σ 1 = 0 and σ 1 = 1, the process (37) becomes the CIR and ECIR processes, respectively, as shown in Examples 1 and 2.
Our MC simulations are based on the EM method, implemented by the basic packages in MATLAB, to obtain numerical solutions of the SDE (37) to calculate (1) by using the Trapezoid integration method in the integral term. In all our calculations, MATLAB R2019a and a PC with the following details were utilized: Intel(R) Core(TM) i7-5700HQ, CPU @2.70GHz, 16 For n = 1, we yield and for n = 2, we obtain for all r > 0 and τ = T − t ≥ 0, where the coefficients A n j (τ) and the solution B(τ) are given in Theorem 3. This example illustrates the formula based on the process (37) with σ 1 = 0 in the case of n ∈ N which actually gets the closed-form formula (11). Thus, validating the obtained closed-form formulas (39) and (40) under selecting the parameters d = 2, θ = 1, σ 0 = 0.01 and σ 1 = 0 for the process (37) with the MC simulations is shown in Figure 3. This implementation with MC method is simulated at each initial value r = 0.1, 0.2, 0.3, . . . , 1.6 to generate 10,000 sample paths of r t , where each path consists of 10,000 steps over the time lengths τ = 0.01, 0.1, 1, 2.
As displayed in Figure 3a, the results from MC simulations (circles) match completely with our formula (39) (solid lines) at every values r for n = 1. Similarly, the results from the MC simulations in Figure 3b also attach with the results obtained from (40). However, the essential disadvantage of MC simulation is that it consumes expensive computational time for approximating the value at each initial value of r. In contrast, our closed-form formulas produce the exact solution for all initial values r > 0 and also give a shorter time in computation. For n = 1, we yield and for γ = 2 to obtain for all r > 0 and τ = T − t ≥ 0, where the coefficients A n j (τ) and the solution B(τ) are provided in Theorem 1. In this example, we test the validation by using the parameters d = 2, θ = 1, σ 0 = 0.01 and σ 1 = 1. We see that these parameters produce the function σ(t) = 0.01e t for which the Riccati differential equation cannot be solved analytically. Thus, the numerical method needs to be applied that has created in Section 4, the FIM with shifted Chebyshev polynomial, to approximate the coefficients A n j (τ) and B(τ) for calculating the formulas (41) and (42). In order to measure the accuracy of the obtained approximate results (41) and (42) by comparison with the MC simulations, we use the mean absolute error defined by: where U n E and W n E are, respectively, the approximate solutions obtained from the FIM and MC methods, r i is the ith initial value and M is the number of initial values used to test. Hence, the comparisons of numerical results of the formulas (41) and (42) with the MC simulations measured by the above mean absolute errors are demonstrated in Table 1 using initial values r = 0.1, 0.2, 0.3, . . . , 1.6. For MC simulations, we vary the sample paths by 10,000, 20,000, 40,000 and 80,000 and each path is discretized into 10,000 steps. Table 1 shows that our approximate formulas for both n = 1, 2 quite match with the results of MC simulations. Especially, for a small length of time τ such as τ = 0.01 in Table 1, they are very close to each other. Moreover, we also observe that they are more coincide as the number of sample paths increases. Consequently, the MC simulations are most likely converge to our approximate formulas. This confirms that our proposed numerical scheme for approximating the formulas in Example 2 is very accurate. Although our method and the MC simulation obtain the results most closely, the average run time (ART) to find their results at each initial r for both is very different. The ARTs consumption for our method is relatively inexpensive, i.e., less than a second. In contrast, the MC simulation consumes tremendous computational times, especially with large path numbers, which can be seen in Table 1.

Interest Rate Swap Pricing
A fixed-to-floating interest rate swap is an agreement for two parties to exchange a series of cash flows where a buyer agrees to pay a floating interest rate to receive a fixed interest rate on a predetermined principle, called a notional principle, over a specified period of time [8]; see Figure 4. It has many potential uses such as in risk management, portfolio management, and speculation. For instance, a company borrowing 10 million dollars from a bank with an interest rate of 3% plus the London Interbank Offered Rate (LIBOR) can sell a fixed-to-floating interest rate swap to hedge against the exposure to fluctuations in interest rates. Because the company will be paying a fixed interest rate instead of LIBOR, it will be much easier for the company to come up with a plan to allocate enough capital in order to service the debt. An interest rate swap is the most traded overthe-counter (OTC) swap at present. According to the Bank for International Settlements report on OTC interest rate derivatives in 2019 [29], the outstanding notional amount for the interest rate swap is over 300 trillion dollars. We consider two types of interest rate swaps, namely the arrears swap and the vanilla swap. The difference between the two swaps is the time which the floating interest rate is fixed at the reset time. A floating payment for an arrears swap is based on an interest rate at a payment time, whereas for a vanilla swap, the interest rate used to calculate the floating payment is fixed at the time before the payment time, which is usually the time of the previous payment. Thus, for an arrears swap, the payment date and the reset time coincide. The formulas for these swaps under the CIR process were proposed by Mallier and Alobaidi [8] using the Green's function approach. Only the formula for the arrears swap is in a closed form. However, the formula is too complicated as it depends on the Kummer's and gamma functions. In addition, the closed-form formula for the vanilla swap has not been derived. To the best of our knowledge, analytical formulas for both swaps under the ECIR process have not been discovered in any literature yet. This section provides analytical formulas for pricing such swaps under the ECIR process from the perspective of a buyer, i.e., a person who pays a floating interest rate to receive a fixed one. Under the CIR process, a special case of the ECIR process, such formulas are explicit.

Arrears Swaps
As interest rate swaps are OTC products, they can be customized differently to suit a buyer's needs. We consider an interest rate swap contract such that fixed interest rate and floating interest rate payments are exchanged at every specified period of time, such as every three or six months over a predetermined time horizon such as two or ten years. At each payment date, an arrears swap buyer pays an interest on a notional principle determined by a floating interest rate at the time of the payment and receives a fixed interest. In other words, the reset time coincides with payment date.
Let P denotes a notional principle,r denotes a fixed rate, and r t denotes a floating rate at time t. Assuming that an arrears swap has a maturity T and N payment dates at T 1 , T 2 , T 3 , . . . , T N = T in an increment of ∆t, the payoff of such swap from a buyer's point of view at the ith payment date, V ar i , which is just the difference between the interest on a notional principle determined by the fixed interest rate and the floating interest rate; see Figure 5, can be expressed by V ar i = (r − r T i )∆tP. By the fundamental theorem of asset pricing [6], the no-arbitrage price for the arrears swap is the expectation of the sum of each payoff discounted to time t under the risk neutral measure. Mathematically, the noarbitrage price for the arrears swap, V ar , with an affine short rate discount, αr t + β where α, β ∈ R, can be expressed as: Theorem 5. Suppose that r t follows the ECIR process (2) with principle P > 0 and α, λ, β ∈ R. Let 0 ≤ t = T 0 < T 1 < T 2 < · · · < T N = T. Then, for all (r, τ i ) ∈ (0, ∞) 2 , where τ i = T i − t and ∆t = T i − T i−1 for i = 1, 2, 3, . . . , N.
Proof. The arrears swap price formula follows (43) with V ar i = (r − r T i )∆tP and then Applying Theorem 2 with n = 0 and n = 1 yields (44).

Corollary 3.
Suppose that r t follows the CIR process (2) with θ(t) = θ, µ(t) = µ and σ(t) = σ. According to Theorem 5, we have Proof. From the result in Theorem 5, we have where the coefficients A 0 0 , A 1 0 and A 1 1 are given in Theorem 3. With some algebraic manipulations, it can be easily checked that the right-hand-side of (46) is equal to the left-hand-side of (45).

Vanilla Swap
In contrast with an arrears swap, the interest rate determining the floating payment is set prior to the payment date, usually at the time of the previous payment. The cash flow of the vanilla swap is shown in Figure 6. For a vanilla swap buyer, the payoff at the ith payment date is V va i = (r − r T i−1 )∆tP. Similar to that of the arrears swap, the no-arbitrage price for a vanilla swap is In Theorem 6, we construct an analytical formula for pricing a vanilla swap by using Theorems 2 and 4. Corollary 4 shows that the formula is explicit under the CIR process. Theorem 6. Suppose that r t follows the ECIR process (2) with principle P > 0 and α, λ, β ∈ R.
for all (r, Proof. The vanilla swap price formula follows (47) with V ar i = (r − r T i−1 )∆tP and then Applying Theorem 4 to each term of the above conditional expectation; the first and second terms with n 1 = n 2 = 0 and the third term n 1 = 1, n 2 = 0, yields (48).
Clearly, our formulas for pricing IRSs are simpler and easier to use in practice compared with those in the literature. Under the CIR process, unlike the formulas by Mallier and Alobaidi' [8], which consist of infinite integrals, which are not necessarily integrable, and Kummer's function, our formulas are explicit. Under the ECIR process, these are the first analytical formulas for arrears and vanilla swaps pricing.

Examples
In this subsection, we compute the prices of both arrears and vanilla options under the modified ECIR(d) process which is to multiply the diffusion term of the ECIR(d) process proposed by Egorov et al. [21] with a parameter λ, i.e., Note that if λ = 1, the modified ECIR(d) process (50) becomes the original ECIR(d) process. The purpose of multiplying λ is to change the magnitude of the volatility without affecting the mean of the process. The parameter values used in this example are α = 1, β = 0, d = 5, θ = 0.5, σ = 0.15 and σ 1 = 0.001. The modeling parameters values are taken from Table 1 of Egorov et al. [21]. Figure 7 shows the prices of the arrears and vanilla swaps with 10-year maturity paid semi-annually computed using the formulas provided in Section 5 as functions of initial interest rates r for λ = 1, 2, 3, 4 where the notional principle P = 1 and the fixed ratē r = 0.05.
The prices are compared with those obtained with MC simulations with 10,000 sample paths where each path has 10,000 steps. We see that the results most likely agree. However, as mentioned before, the MC simulation is computationally expensive. For each point, MC uses 600 s on average, whereas the scheme proposed in Section 5 consumes less than 0.001 s.
As mentioned above, the parameter λ describes the process'volatility. This means if the value of λ increases, the volatility of the process also increases. In practice, the higher the volatility, the more expensive the price. This agrees with the results shown in Figure 7a,b in which the parameters λ are varied from 1 to 4.

Fractional ECIR Process
This paper also considers the fractional ECIR process which is a class of fractional Pearson diffusions [30,31] where the first derivative with respect to time variable in (7) is replaced by a Caputo fractional derivative [32] of the order α ∈ (0, 1). In real-world applications, time-fractional diffusion equations are useful in many branches of science, especially in finance. For instance, the time-fractional diffusion is used to model the delays between trades based on the continuous-time random walks, and they have been applied to extend the Black-Scholes formula in a subdiffusive regime [33,34].
Before embarking into the details of fractional ECIR process, we provide the basic definitions of fractional derivatives. The necessary notations and some important facts used throughout this section are also given. More details on definitions and basic results of fractional calculus can be found in [35]. and we also have Therefore, by comparing the coefficients in r, the above equation can be solved through a system of time-fractional differential Equations (55) and (56). Thus, the coefficient parameters A j (τ) can be calculated from (55) and (56). However, the solutions of (55) and (56) are always unavailable in closed form. Thus, a numerical method is needed. Next, to find the coefficient functions A j (τ), we investigate a numerical scheme to solve the time-fractional differential equations (55) and (56). Let us first uniformly discretize the time domain [0, τ] into m steps, i.e., τ i = i∆τ for i = 0, 1, 2, . . . , m, where ∆t is a time step. Then, we have: where A m j := A j (τ m ), θ m := θ(T − τ m ), µ m := µ(T − τ m ) and σ m := σ(T − τ m ). Then, we consider the time-fractional derivative terms of (58) and (59) by employing the Caputo fractional derivative in Definition 2. We obtain: After that, we use the first-order forward difference quotient to approximate the time derivative term in the above equation. For convenience, we also adjust the index in the last step of the following process by defining k := m − i − 1. Then, Note that for m = 1, the summation term is set to be zero. Therefore, the solutions of A j (τ m ) can be approximated by running the iterative system (61) and (62), respectively, staring with A 0 0 = 1 and A 0 j = 0 for j = 1, 2, 3, . . . , n. Finally, we can also obtain a numerical result of the time-fractional conditional moment (54) that is: Next, to test the efficiency of our proposed numerical scheme, we examine by varying the fractional order α ∈ (0, 1) of Examples 3 and 4 for both first-and second-moment of conditional expectations with the following parameters θ = 1.3, µ = 2 and σ = 0.1 at the terminal time T = 1.
By using the proposed numerical scheme described in this section, we examine the behavior of the solutions (54) for various values of α ∈ (0, 1) by discretizing the time domain [0, 1] into 100 steps distributed uniformly. Then, we can demonstrate the plotting solutions u 1 α (r, 1) and u 2 α (r, 1) for r ∈ (0, 1] of both examples by varying the fractional order α ∈ {0.5, 0.6, , 0.7, 0.8, 0.9}. We can obviously see that when α → 1, the behavior of the obtained approximate solutions, depicted in Figure 8, tends to the blue solid line, i.e., the integer order solution at α = 1 obtained by Formula (11) with parameters λ = α = β = 0.

Remark 4.
The idea presented in this section can be also applied to the generalised geometric Brownian motion (GBM), which the standard and subdiffusive GBM arise as special cases, see [39] for more details.

Conclusions and Discussions
In this paper, we first study the conditional expectation of the product of polynomial and exponential functions, E Q r γ T e −λr T − T t (αr s +β) ds | r t = r where α, λ, β, γ ∈ R and r t follows ECIR process. By solving the PDE (7) derived from the Feynman-Kac representation, we archive an analytical formula of the conditional expectation given in (1) for the ECIR process in terms of the infinite sum of analytical expressions as shown in Theorem 1. Interestingly, the infinite sum can be reduced to a finite sum if γ ∈ N ∪ {0}, as shown in Theorem 2. We also show that under the CIR process, which is a special case of the ECIR process, as the parameters given in (4) can be integrable. The Riccati differential Equation (5) can be solved analytically. The formula of the expectation can be expressed in the closed form as in Theorem 3. The formula is then extended to derive the expectation of a product of two polynomials and exponential functions, E Q r n 1 s r n 2 T e − T t (αr u +β) du | r t = r where n 1 , n 2 ∈ N ∪ {0} and 0 ≤ t < s < T, see Theorem 4.
In addition, in Section 4, we propose the numerical scheme constructed using the FIM with shifted Chebyshev polynomials to evaluate the expectation (1) when the Riccati differential equation cannot be solved analytically and show its efficiency as well as its accuracy by comparing it with the MC simulations.
To illustrate an application of such formulas in finance, we apply the Formula (10) in Theorem 2 and the formula (17) in Theorem 4 to derive the analytical pricing formulas for two interest rate swaps, namely the arrears swap and vanilla swaps under the ECIR process. It is shown that the formulas for such swaps are explicit under the CIR process. The numerical pricing examples for both swaps under the ECIR process are also given. Finally, in Section 6, we provide a study case of the fractional ECIR process and adopt the idea from our main results.

Acknowledgments:
We are grateful for a variety of valuable suggestions from the anonymous referees which have substantially improved the quality and presentation of the results. All errors are the author's own responsibility. This research project is supported by Second Century Fund (C2F), Chulalongkorn University.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: Over-the-counter PDE Partial differential equation PDF Probability density function SDE Stochastic differential equation