Next Article in Journal
Using Machine Learning to Profile Asymmetry between Spiral Galaxies with Opposite Spin Directions
Next Article in Special Issue
Construction of General Implicit-Block Method with Three-Points for Solving Seventh-Order Ordinary Differential Equations
Previous Article in Journal
Spin Sum Rule of the Nucleon in the QCD Instanton Vacuum
Previous Article in Special Issue
On a Nonlocal Coupled System of Hilfer Generalized Proportional Fractional Differential Equations
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Analytically Pricing Formula for Contingent Claim with Polynomial Payoff under ECIR Process

Department of Mathematics and Computer Science, Faculty of Science, Chulalongkorn University, Bangkok 10330, Thailand
*
Author to whom correspondence should be addressed.
Symmetry 2022, 14(5), 933; https://doi.org/10.3390/sym14050933
Submission received: 3 April 2022 / Revised: 26 April 2022 / Accepted: 2 May 2022 / Published: 4 May 2022
(This article belongs to the Special Issue Differential Equations and Applied Mathematics)

Abstract

:
Contingent claims, such as bonds, swaps, and options, are financial derivatives whose payoffs depend on uncertain future real values of underlying assets which emphasize various real-world applications. In general, valuations for contingent claims can be derived from the conditional expectations of underlying assets. For a simple process, the moments are usually directly obtained from its transition probability density function (PDF). However, if the transition PDF is unavailable in simple form, the derivations of the moments and the contingent claim prices will not be accessible in closed forms. This paper provides a closed-form formula for pricing contingent claims with nonlinear payoff under a no-arbitrage principle when underlying assets follow the extended Cox–Ingersoll–Ross (ECIR) process with the symmetry properties of the Brownian motion. The formula proposed here is a consequence of successfully solving an explicit solution for a system of recurrence partial differential equations in which its solution subtly depends on the conditional moments. For the particular CIR process, we obtain simple closed-form formulas by solving the Riccati differential equation. Furthermore, we carry out a complete investigation of the convergent case for those formulas. In case such as that of the unsolvable Riccati differential equation, ECIR case, a numerical method for numerically evaluating the mentioned analytical formulas and numerical validations for the formulas are examined. The validity and efficiency of the formulas are numerically demonstrated by comparison with results from Monte Carlo simulations for various modeling parameters. Finally, the proposed formula is applied to the value contingent claims such as coupon bonds, interest rate swaps, and arrears swaps.

1. Introduction

Conditional expectations are useful statistical values applied in many branches in science, especially in describing behaviors of observed data, and are often studied from a probabilistic viewpoint based on transition probability density functions (PDFs) of the data. Many applications in finance and economics require the knowledge of conditional expectations, for example, in the valuation of financial products, e.g., contingent claims such as coupon bonds, swaps, discount factor, etc., as can be seen in the works of Ben-Ameur et al. [1] and Grasselli [2] for more details.
To value financial products, the no-arbitrage principle is essential; if arbitrage exists, it is guaranteed that the investor can make profit from nothing, which provides an investment opportunity with infinite return. To satisfy the arbitrage-free property, symmetric information is required [3]. The valuation of a contingent claim is usually investigated based on the conditional expectation (1) under a filtered risk-neutral probability space ( Ω , F t , { F t } 0 t T , Q ) in the form
E Q e t T α x s + β d s f T + t T e t s α x u + β d u g s d s | F t ,
where α , β R and { x t } 0 t T is an adapted stochastic process describing the underlying asset, and for some real value functions f T and g s with 0 t < s T .
In this work, we consider the conditional expectations (1) in the form
E Q x T γ e t T ( α x s + β ) d s | F t = E Q x T γ e t T ( α x s + β ) d s | x t = x ,
where γ R and to be more specific, the process x t is considered as a short rate described by the extended Cox–Ingersoll–Ross (ECIR) process [4]. The case that γ = 1 in (2) was also studied in Duffie et al.’s work [5] for the class of affine jump-diffusion processes with generalized payoff, including some case studies such as those of Ornstein–Uhlenbeck (OU) and the Cox–Ingersoll–Ross (CIR) or squared processes. To be more general, the Wiener process W t in the ECIR process can be generalized by using the mixed fractional Brownian motion, i.e., a linear combination of the Wiener process W t and the fractional Brownian motion W t H with the Hurst parameter, H ( 0 , 1 ) . The mixed fractional Brownian motion has some useful characteristic, e.g., it is arbitrage free and gains more attention for applications in finance. An example of an application based on the mixed fractional Brownian motion is described in [6].
Note that the contingent claim prices in the forms (1) and (2) have many useful applications. In finance, a non-arbitrage price at time t of financial derivatives is considered on a conditional expectation under a risk-neutral measure of their discounted payoff, more details of which can be found in [7]. Therefore, the valuations of financial derivatives, such as the coupon bonds, variance swaps, interest rate swaps, and options, always involve calculating the forms of conditional expectations (1) and (2).
Since many processes describing financial assets have no transition PDFs in simple form, the conditional expectations, the valuation of many financial products involving (1) and (2) usually are not accessible in closed form; thus, alternative methods are required. In practice, when analytical formulas for the expectations are not known in concise form, a practical method such as the Monte Carlo simulation is required, which has disadvantages in term of computational time. This paper aims to propose a closed-form formula for the conditional expectation (1) based on the solution of a partial differential equation (PDE) according to the Feynman–Kac representation, without requiring the knowledge of the transition PDF. In some applications, the related PDEs may have no closed form solution and numerical methods are required in order to obtain the results, for example, Ahmadian and Ballestra [8] proposed the finite element method to solve ruin-related problems, and Liang and Zou [9] studied the valuation of credit contingent interest rate swap with credit rating migration using the alternating direction implicit method.
The CIR process is a diffusion process satisfying the Pearson Equation [10] and involving a wide variety of issues in many branches—more details on this can be found in [11]. This process was initially introduced by Feller [12] as a population growth stochastic model and becomes popular in finance when Cox et al. [13] applied it to describe the evolution behavior of short-term interest rates. Even though the CIR process is very useful in terms of pricing financial derivatives, especially short-term interest rates, the process has a limitation on its constant parameters, which are not suitable for modeling time-varying observed data. A lot of strong empirical evidence has found that extreme movements in finance-based practices tend to be assumed in function of the time, more details on which can be found in [4,14,15]. In 1990, Hull and White [4] proposed a novel SDE such that the dynamics of the CIR process can be governed by time-depending parameters, which is called the ECIR process. The ECIR process is so attractive as a practical model to price the European bond option. In 2003, Egorov et al. [16] presented the transition PDF of the ECIR process in a complicated form of the modified Bessel function of the first kind and proposed a method to receive a closed-form approximation of the transition PDF through the Hermite approximation. This is one of the most practically used empirical evidences to confirm that it is not easy to obtain the conditional expectation (2) by using the transition PDF of the ECIR process.
For this work, we assume that x s is governed by the ECIR process and R ( t , r ) is a bounded continuous discount rate function; the value of the asset at initial time t can be rewritten as V t = E Q e t T R ( s , x s ) d s f ( x T ) F t , as can also be seen in [17]. With R 0 and f ( x ) = x γ where γ R , under the probability measure P, Dufresne [18] derived a closed-form formula for the conditional moments, E P x T γ F t , for some sufficient conditions on γ and the parameters in the CIR process x t . In 2007, under a risk-neutral probability Q based on the CIR process, Ben-Ameur et al. [1] estimated an ex-coupon holding value at time T, E Q e t T x s d s f ( x T ) F t , where f denotes the value of the bond at time T. Their result is applicable but not in closed-form solutions. Moreover, they needed to find the joint distribution of the random vector x t + δ , t t + δ x t d t where δ > 0 . This is characterized by its Laplace transform which can be described by the conditional expectation (2) when γ = 0 and β = 0 . Recently, under the CIR process, Grasselli [2] directly determined a mathematical expression of the conditional expectation (2). However, their expression was expressed in terms of a product of the confluent hypergeometric function and the gamma functions. This may be hard to work with in some cases.
We move our focus onto the ECIR process. In 2016, under probability measures, an analytical formula was proposed by Rujivan [19] which was extended from Dufresne’s approach [18] to the ECIR process for the case γ R . In 2018, an explicit formula for the conditional expectations of a product of polynomial and exponential function, in the form E P x T γ e λ x T F t , was analytically derived by Sutthimat et al. [20] for the case γ , λ R . Their results cover the results in such formulas of the Rujivan’s present [19] in the case of λ = 0 . Indeed, both works on the ECIR process have a limitation. A major concern for their formulas in Theorems 1 and 2 of their works [19,20] is that the coefficients A γ k ( τ ) may not be integrable to receive the exact integrations. Some numerical methods for integrations are required to manipulate those integral terms in this very reason. However, both results presented in [19,20] did not provide any methods to overcome this issue. Both results are not ready for practical applications. In our analysis, we also present a numerical method to deal with this challenge.
The useful applications of (2) under the CIR and ECIR processes need to be mentioned. To price interest rate swaps (IRSs), a financial contingent claim, which is a financial derivative whose payoff depends on the uncertain future real value of other underlying assets, is assumed together to follow the CIR process. The IRS is one of the common types of contingent claim derivatives as a modified version of swaps. Normally, the cash flows of IRSs on the payment dates are the same as the forward rate agreements (FRAs) which are the contact that the forward rates can be fixed by an investor. In brief, an IRS is a form of series of FRAs. We give some interesting works under the assumption that the discount rate is continuously compounded, which have been well studied in the literature and can apply our result of (2) to those works. In 2004, Mallier and Alobaidi [21] supposed that the risk-neutral interest rates follow the CIR process. By utilizing the Green’s function approach, they provided analytical expressions of the swap values for two well-known types of IRSs, which are the arrears and vanilla swaps. Their analytical expressions, a sum of values of the FRAs, which was in a closed form for an arrears swap but very complicated because it depends on the gamma and the Kummer’s functions. However, for a vanilla swap, their result was not in closed form and much more complicated than the results of those arrears swaps. Unlike the results of Mallier and Alobaidi, Moreno and Platania [22] provided a mathematical formula for the FRA values in 2015 for a special case of the ECIR process, namely the cyclical square-root model; more details on this can be found in their Proposition 8. Thus, an interest rate swap valuation was straightforwardly obtained as a consequence of this proposition. In fact, Proposition 8 in their work consists of the Mathieu cosine and sine functions [23] and the parameter A ( τ ) given in this proposition may not be exactly integrable. Thamrongrat and Rujivan [24] recently published an analytical formula for pricing IRSs in terms of bond prices based on the ECIR process which was performed under a discrete discount rate.
This paper successfully worked out an analytical formula of the conditional expectations (2) for the ECIR process in terms of analytical expression. Furthermore, their consequences were investigated without requiring the transition PDF of the ECIR process. Additionally, for the CIR process, the formulas were reduced to concise forms which give a greater advantage than the other approaches in the literature. Furthermore, the ECIR process-facilitated valuation of financial derivatives is provided by using our proposed results. Under ECIR process, this paper further suggests a numerical algorithm of the conditional expectations (2) in case the Riccati differential equation may not be exactly solved.
This paper is organized as follows. A brief overview of the CIR process as well as the ECIR process are provided in Section 2. The key methodology is mentioned in Section 3 to address the main relevant concept for our main result, which is an analytical formula of the conditional expectations (2) of the ECIR process. Section 4 gives a numerical method to work with the generalized Riccati differential equation and one major concerned limitation of our formula is discussed here. Section 5 validates our formulas and discusses the analytical formulas’ advantages compared with the Monte Carlo (MC) simulations. In Section 6, some financial applications are demonstrated based on our proposed formula. The aim of this study is recapitulated and concluded in Section 7.

2. The Extended Cox–Ingersoll–Ross Process

In this paper, we assume that the interest rate x t follows the ECIR process under a risk-neutral probability measure Q, which is a diffusion model whose solution satisfies the following SDE [4],
d x t = θ ( t ) μ ( t ) x t d t + σ ( t ) x t d W t , 0 t T .
The well-known W t is a Wiener process or Brownian motion whose increments are generated by the symmetry of mean zero Gaussian distribution. Sometimes, the parameters in (3) are referred to as follows: θ is the speed of adjustment to the long-term mean μ , while σ indicates to the state space of the diffusion. The two assumptions explored by Maghsoodi [15] are required to demonstrate that there is a path-wise unique strong solution for the ECIR process x t and to avoid zero a.e. with regard to the probability measure P for a specified time t during a time period [ 0 , T ] ; more details on this can be found in Theorems 2.1 and 2.4 of [15]. We thus require the following sufficient condition.
Assumption 1.
Time parameters θ ( t ) , μ ( t ) and σ ( t ) in (3) are smooth and strictly positive. The time function μ ( t ) σ 2 ( t ) is locally bounded and 2 θ ( t ) μ ( t ) σ ( t ) 2 on [ 0 , T ] .
To achieve our aim, a common question arises: why not directly use the transition PDF of the CIR process? It is known that its transition PDF has an expression in a form of Gamma density function and Laguerre polynomials; more details on this can be found in [25,26]. The transition PDF can be written in an explicit form as
p x , T x t , t = c τ e ( u + v ) v u q / 2 I q 2 u v ,
where τ = T t , c τ = 2 θ σ 2 1 e θ τ , u = c τ x t e θ τ , v s . = c τ x , q = 2 θ μ σ 2 1 and I q ( · ) is the ordered q Bessel function of the first kind,
I q ( x ) = k = 0 x 2 2 k + q 1 Γ ( k + 1 ) Γ ( k + q + 1 ) .
Since the transition PDF is complicated, as shown above, solving the closed-form formulas for (2) by applying the transition PDF is more complicated.
It becomes even more difficult in the ECIR process, for example, as in the ECIR ( d ) process observed by Egorov et al. in 2003 [16]. Its dynamics are followed by a time-inhomogeneous diffusion process as
d x t = θ σ 0 2 d 4 θ e 2 σ 1 t x t d t + σ 0 e σ 1 t x t d W t ,
where θ , σ 0 are positive, σ 1 is real and d is positive. Its transition PDF was first proposed by Maghsoodi [15],
p x , T x t , t = 1 2 G e λ + G r 2 G r λ d 2 4 I d 2 1 λ G r
with λ = x t v , G = e θ τ v , v = 8 σ 1 σ 0 2 e θ τ e 2 σ 1 T e 2 σ 1 t , τ = T t and again I q ( · ) is the Bessel function of the first kind. To avoid using those of the transition PDFs for solving (2), this paper applies Feynman–Kac representation which offers a method for solving a conditional expectation of an Itô random process by deterministic implementations, more details on which can be found in [27,28,29,30].

3. Main Results

3.1. Closed-Form Formula: Conditional Expectation

The first contribution of this section is to provide an integral form formula for the pricing of the contingent claim with the polynomial payoff (2) by solving the PDE according to the Feynman–Kac representation.
Theorem 1.
Let 0 t T and x t follow the ECIR process (3) with α , β , γ R . Then,
E Q x T γ e t T ( α x s + β ) d s x t = x = e B ( τ ) x j = 0 A j γ ( τ ) x γ j = : U E γ ( x , τ ) ,
for all ( x , τ ) D E ( 0 , ) × [ 0 , ) , where τ = T t 0 . Under the assumption that the infinite series in (4) uniformly converges on D E , the coefficients in (4) can be expressed as
A 0 γ ( τ ) = e 0 τ P 0 ( u ) d u , A j γ ( τ ) = e 0 τ P j ( u ) d u 0 τ e 0 u P j ( s ) d s Q j ( u ) A j 1 γ ( u ) d u ,
for j N , where
P j ( τ ) = θ ( T τ ) μ ( T τ ) B ( τ ) θ ( T τ ) ( γ j ) + σ 2 ( T τ ) ( γ j ) B ( τ ) β , Q j ( τ ) = ( γ j + 1 ) θ ( T τ ) μ ( T τ ) + 1 2 σ 2 ( T τ ) ( γ j ) ,
and B ( τ ) can be obtained by solving the following Riccati differential equation
B ( τ ) = 1 2 σ 2 ( T τ ) B 2 ( τ ) θ ( T τ ) B ( τ ) α , B ( 0 ) = 0 .
Proof. 
The Feynman–Kac representation is used to solve U : = U E γ ( x , τ ) , which is defined as a series in (4) that satisfies the appropriate PDE under the uniformly convergent assumption. Thus, we have
U τ + θ ( T τ ) μ ( T τ ) x U x + 1 2 x σ 2 ( T τ ) U x x α x + β U = 0
with the initial condition at τ = 0
U E γ ( x , 0 ) = E Q x T γ e T T ( α x s + β ) d s | x T = x = x γ .
Comparing the coefficients at τ = 0 in (4) and (9) to receive the initial conditions B ( 0 ) = 0 , A 0 γ ( 0 ) = 1 and A j γ ( 0 ) = 0 for j N . Then, we compute (8) using (4) to find the partial derivatives U τ , U x and U x x . Consequently, we have
0 = e B ( τ ) x j = 0 A j γ ( τ ) x γ j + x B ( τ ) j = 0 A j γ ( τ ) x γ j + θ ( T τ ) μ ( T τ ) x e B ( τ ) x j = 0 ( n j ) A j γ ( τ ) x γ j 1 + B ( τ ) j = 0 A j γ ( τ ) x γ j + 1 2 x σ 2 ( T τ ) e B ( τ ) x j = 0 ( n j ) ( n j 1 ) A j γ ( τ ) x γ j 2 + 2 B ( τ ) j = 0 ( n j ) A j γ ( τ ) x γ j 1 + B 2 ( τ ) j = 0 A j γ ( τ ) x γ j α x + β e B ( τ ) x j = 0 A j γ ( τ ) x γ j ,
which can be simplified as
0 = A 0 γ ( τ ) B ( τ ) + θ ( T τ ) B ( τ ) 1 2 σ 2 ( T τ ) B 2 ( τ ) + α x γ + 1 + A 0 γ ( τ ) + P 0 ( τ ) A 0 γ ( τ ) A 1 γ ( τ ) B ( τ ) + θ ( T τ ) B ( τ ) 1 2 σ 2 ( T τ ) B 2 ( τ ) + α x γ + j = 1 A j γ ( τ ) + P j ( τ ) A j γ ( τ ) + Q j ( τ ) A j 1 γ ( τ ) A j + 1 γ ( τ ) B ( τ ) + θ ( T τ ) B ( τ ) 1 2 σ 2 ( T τ ) B 2 ( τ ) + α x γ j
Considering (11) as a series solution in x, we receive the system of ODEs as follows,
0 = A 0 γ ( τ ) P 0 ( τ ) A 0 γ ( τ ) , 0 = A j γ ( τ ) P j ( τ ) A j γ ( τ ) Q j ( τ ) A j 1 γ ( τ ) , 0 = B ( τ ) + θ ( T τ ) B ( τ ) 1 2 σ 2 ( T τ ) B 2 ( τ ) + α ,
with their initial conditions B ( 0 ) = 0 , A 0 γ ( 0 ) = 1 and A j γ ( 0 ) = 0 for j N . Thus, the solutions of the ODEs are given in (5) and (7) as required. □
Notice that the function B is the solution of a generalized Riccati differential equation. It is well known that its solution is not available and can be only treated in some cases. To overcome this problem, one can employ the numerical method discussed in the following section.
By examining (4) when γ = n N 0 and γ = m 2 θ ( τ ) μ ( τ ) σ 2 ( τ ) for m N , the infinite sum in (4) is terminated at order n and can be written as in the following theorems.
Corollary 1.
According to Theorem 1 with γ = n N 0 , we have
U E n ( x , τ ) = E Q x T n e t T ( α x s + β ) d s | x t = x = e B ( τ ) x j = 0 n A j n ( τ ) x n j ,
for all ( x , τ ) D E ( 0 , ) × [ 0 , ) , τ = T t 0 , where the coefficients A j n ( τ ) are defined by (5) and (6), and B ( τ ) is the solution of (7).
Proof. 
Considering (5) when j = n + 1 obtains Q j ( τ ) = Q n + 1 ( τ ) = 0 , it implies that the coefficient A j n ( τ ) = A n + 1 n ( τ ) = 0 . Since the coefficient A j n ( τ ) is a type of recurrence problem involving the initial condition A n + 1 n ( τ ) = 0 , A j n ( τ ) = 0 for all j n + 1 . Thus, the infinite sum (4) can be reduced to the finite sum as shown in (12). □
Corollary 2.
According to Theorem 1 with γ = m 2 θ ( τ ) μ ( τ ) σ 2 ( τ ) for all τ 0 , m N , we have
U E γ ( x , τ ) = E Q x T γ e t T ( α x s + β ) d s | x t = x = e B ( τ ) x j = 0 m A j γ ( τ ) x γ j ,
for ( x , τ ) D E ( 0 , ) × [ 0 , ) , τ = T t 0 , where the coefficients A j γ ( τ ) are defined by (5) and (6), and B ( τ ) is the solution of (7).
Proof. 
The proof can be shown similarly as in Corollary (1) by considering (5) and (6). □
It is worth mentioning that our proposed formulas for the ECIR process generalize some results in the literature. Furthermore, this section provides concise formulas for the CIR process reduced from the previous theorems where the parameters depending on time θ ( t ) = θ , μ ( t ) = μ and σ ( τ ) = σ are constant functions. In this case, Riccati differential Equation (7) is solvable; thus, integral functions can be exactly solved as shown in the following theorems.
Corollary 3.
Suppose that x t follows the CIR process with α θ 2 2 σ 2 , β , γ R . Let 0 t T . Then,
U C γ ( x , τ ) : = E Q x T γ e t T ( α x s + β ) d s | x t = x = e B ( τ ) x j = 0 A j γ ( τ ) x γ j ,
for all ( x , τ ) D C ( 0 , ) × [ 0 , ) , τ = T t 0 . Under the assumption that the infinite series in (14) uniformly converges on D C , the coefficients in (14) can be expressed as
A 0 γ ( τ ) = H 0 ( τ ) , A j γ ( τ ) = H j ( τ ) k = 1 j 2 Q k k e ρ τ 1 ( ρ θ ) + e ρ τ ( ρ + θ ) j ,
for j N , where ρ = θ 2 + 2 α σ 2 and
Q j = ( γ j + 1 ) θ μ + 1 2 σ 2 ( γ j ) , H j ( τ ) = exp θ 2 μ σ 2 β + ( γ j ) + θ μ σ 2 ρ τ 2 ρ ( ρ θ ) + e ρ τ ( ρ + θ ) 2 ( γ j ) + θ μ σ 2 .
In addition, function B given in (7) can be solved as
B ( τ ) = 2 α ( e ρ τ 1 ) ρ ( e ρ τ + 1 ) + θ ( e ρ τ 1 ) .
Proof. 
Determining (7) with constant parameters θ ( t ) = θ , μ ( t ) = μ and σ ( τ ) = σ , the explicit solution for the Riccati differential Equation (7) is
B ( τ ) = 2 α ( e ρ τ 1 ) ρ ( e ρ τ + 1 ) + θ ( e ρ τ 1 ) ,
where ρ = θ 2 + 2 α σ 2 , and its integration is
0 τ B ( u ) d u = 2 σ 2 ln 2 ρ e ( ρ + θ 2 ) τ ( ρ θ ) + e ρ τ ( ρ + θ ) .
Thus, we have
e 0 τ P j ( u ) d u = exp 0 τ θ μ B ( u ) θ ( γ j ) + σ 2 ( γ j ) B ( u ) β d u = exp θ 2 μ σ 2 β + ( γ j ) + θ μ σ 2 ρ τ 2 ρ ( ρ θ ) + e ρ τ ( ρ + θ ) 2 ( γ j ) + θ μ σ 2 ,
which are defined as H j ( τ ) for j N 0 . From Theorem 1, A 0 γ ( τ ) = e 0 τ P 0 ( u ) d u = H 0 ( τ ) and
A 1 γ ( τ ) = H 1 ( τ ) 0 τ 1 H 1 ( u ) Q 1 A 0 γ ( u ) d u = 2 H 1 ( τ ) Q 1 e ρ τ 1 ( ρ θ ) + e ρ τ ( ρ + θ ) .
From the result presented in (5) for j N , we obtain
A j γ ( τ ) = H j ( τ ) 0 τ 1 H j ( u ) Q j A j 1 γ ( u ) d u = H j ( τ ) 0 τ 1 H j ( u ) Q j H j 1 ( u ) k = 1 j 1 2 Q k k e ρ u 1 ( ρ θ ) + e ρ u ( ρ + θ ) j 1 d u = H j ( τ ) k = 1 j 1 2 Q k k 0 τ e ρ u 2 ρ ( ρ θ ) + e ρ u ( ρ + θ ) 2 e ρ u 1 ( ρ θ ) + e ρ u ( ρ + θ ) j 1 d u = H j ( τ ) k = 1 j 2 Q k k e ρ τ 1 ( ρ θ ) + e ρ τ ( ρ + θ ) j .
Under the uniformly convergent assumption, this completes the proof. □
The case in which γ = n N 0 , U C γ ( x , τ ) in (14) can be expressed as a finite term of a power series in x as follows.
Corollary 4.
According to Corollary 3 with γ = n N 0 , we have
U C n ( x , τ ) = E Q x T n e t T ( α x s + β ) d s | x t = x = e B ( τ ) x j = 0 n A j n ( τ ) x n j ,
for ( x , τ ) D C ( 0 , ) × [ 0 , ) , τ = T t 0 , where the coefficients A j n ( τ ) and B ( τ ) are defined by (15), (16) and (17).
Proof. 
The proof is rather trivial by combining Corollary 1 with Corollary 3. □
Furthermore, in this case and in that where γ = m 2 θ μ σ 2 for m N , the following theorem can readily be reduced from Corollary 3.
Corollary 5.
According to Corollary 3 with γ = m 2 θ μ σ 2 for m N , we have
U C γ ( x , τ ) = E Q x T γ e t T ( α x s + β ) d s | x t = x = e B ( τ ) x j = 0 m A j γ ( τ ) x γ j ,
for ( x , τ ) D C ( 0 , ) × [ 0 , ) , τ = T t 0 , where coefficients A j γ ( τ ) and B ( τ ) are defined by (15), (16) and (17).
Proof. 
The proof is relatively trivial by combining Corollary 2 with Corollary 3. □
Remark 1.
Consequently, our approach formulas can be extended to calculate the mixed polynomial payoff. By applying the tower property for 0 t < s T where τ 1 = s t and τ 2 = T s , the price of the T-claim with mixed polynomial payoff for ECIR process (3) can be expressed as
E Q x s x T e t T ( α x u + β ) d u x t = x = E Q x s e t s ( α x u + β ) d u E Q x T e s T ( α x u + β ) d u x s x t = x .
Remark 2.
The benefits of this work, when α = β = 0 , can be performed by special cases of our proposed theorems, such as the first and the second conditional moments, variance, central moment, mixed moment, covariance, and correlation.

3.2. Closed-Form Formula: Unconditional Expectation

Under Assumption 1, this section proposes two corollaries which are deduced from the conditional formula for the CIR process presented in Corollary 4 to the unconditional formula whereas τ . It should be noted that the following formulas are no longer dependent on the given value of x.
Corollary 6.
Suppose that x t follows the CIR process and α , β > 0 . The price of the T-claim with polynomial payoff at equilibrium, for all n N 0 , x > 0 and τ = T t 0 , is
lim τ U C n ( x , τ ) = lim T E Q x T n e t T ( α x s + β ) d s | x t = x = 0 .
Proof. 
By considering H j ( τ ) in (16), for all j = 0 , 1 , 2 , , n , it can be rewritten as
H j ( τ ) = 2 ρ exp ρ 2 + θ 2 μ σ 2 β 2 ( n j ) + θ μ σ 2 τ ( ρ θ ) + e ρ τ ( ρ + θ ) 2 ( n j ) + θ μ σ 2 .
It is not difficult to see that ρ 2 + θ 2 μ σ 2 β 2 ( n j ) + θ μ σ 2 < ρ for all j. So, lim τ H j ( τ ) = 0 and then
lim τ A j ( τ ) = k = 1 j 2 Q k k lim τ H j ( τ ) e ρ τ 1 ( ρ θ ) + e ρ τ ( ρ + θ ) j = 0 ,
for all j = 0 , 1 , 2 , , n . Thus, lim τ U C n ( x , τ ) = 0 . □
Corollary 7.
Suppose that x t follows the CIR process and α , β = 0 . The price of the T-claim with polynomial payoff at equilibrium, for all n N 0 , x > 0 and τ = T t 0 , is
U n : = lim τ U C n ( x , τ ) = lim T E Q x T n | x t = x = k = 1 n 2 θ μ + σ 2 ( n k ) 2 θ .
Proof. 
According to the (20) in Corollary 4 with α , β = 0 , it is then obtained ρ = θ . As τ , we note that lim τ H n ( τ ) = 1 and lim τ H j ( τ ) = 0 for all j = 0 , 1 , , n 1 . In other words, the coefficient terms A j ( τ ) of x n j approach 0, for k = 0 , 1 , , n 1 . The remaining term, j = n , is only determined. From (17) in Corollary 3, since α = 0 , B ( τ ) = 0 and yields
U n = lim τ A n ( τ ) = lim τ H n ( τ ) k = 1 n 2 Q k k e θ τ 1 2 θ e θ τ n = k = 1 n 2 θ μ + σ 2 ( n k ) 2 θ ,
as required. □
Remark 3.
After some algebraic manipulations, we can show that Rujivan’s formula [19] reduces to our Formula (22).
Remark 4.
It is obvious to what extent Corollary 4 applies to unconditional cases when α , β 0 as displayed in Corollaries 6 and 7. In particular, it is easily checked that, when β < 0 , lim τ U C n ( x , τ ) converges towards 0 if and only if θ 2 μ σ 2 β θ μ < ρ .

3.3. Analysis of Convergence

The aim of this section is to investigate Corollary 3 and whether the infinite sum given in (14) converges. Notice that series (14) converges whenever A j γ = 0 for some j N 0 . To see that the factor of A j γ in (15) only depends on Q j which can be zero. In other words, the series (14) converges if and only if Q j = 0 for some j N 0 . There are only two cases that Q j = 0 , i.e., γ = n N 0 and γ = m 2 θ μ σ 2 for m N . The convergence cases of (14) are already provided in Corollaries 4 and 5. However, the case of Q j 0 for any j N 0 , and the series (14) diverges as follows.
lim n A n + 1 γ ( τ ) x γ n 1 A n γ ( τ ) x γ n = lim n 2 n + 1 ( n + 1 ) ! H n + 1 ( τ ) k = 1 n + 1 Q k e ρ τ 1 ( ρ θ ) + e ρ τ ( ρ + θ ) n + 1 x γ n 1 2 n n ! H n ( τ ) k = 1 n Q k e ρ τ 1 ( ρ θ ) + e ρ τ ( ρ + θ ) n x γ n = lim n 2 e ρ τ Q n + 1 x ( n + 1 ) e ρ τ 1 ( ρ θ ) + e ρ τ ( ρ + θ ) 2 ρ ( ρ θ ) + e ρ τ ( ρ + θ ) 2 .
Since Q n + 1 is a second-degree polynomial in n, the above expression is O ( n ) , thus, by ratio test, the (14) diverges.

4. Numerical Procedures

The valuation of the contingent claim with polynomial payoff based on the ECIR process through Theorem 1, the Formula (4), is an infinite sum of coefficients in (5). These coefficients are defined as in the integral forms and depend on many parameters. Under certain circumstances, i.e., when parameters are complicated, the integral cannot be precisely evaluated, or when the Riccati differential Equation (7) cannot be solved directly. Thus, numerical methods are required to approximate the coefficients. In this section, we numerically investigate the coefficients in (5) by utilizing numerical schemes based on the symmetry concept to approximate the Formula (4).
Let us first consider the Riccati differential Equation (7). From Corollary 3, if the Riccati differential equation has constant coefficients, it has the exact solution, as shown in (18). However, if it has variable coefficients, the analytical solution is not easily obtained. In this case, one needs to approximate the solution by a numerical method; for example, in this work, we use the fourth-order Runge–Kutta (RK4) method [31]. Thus, we are concerned with the following initial value problem:
B ( s ) = 1 2 σ 2 ( T s ) B 2 ( s ) θ ( T s ) B ( s ) α , B ( 0 ) = 0
for s [ 0 , τ ] . We uniformly divide [ 0 , τ ] into m subintervals generated by s i = i h , i = 0 , 1 , , m , where h = τ m is the step size. Then, we denote (23) by
f ( s , B ) : = 1 2 σ 2 ( T s ) B 2 θ ( T s ) B α .
Let B i = B ( s i ) ; then B 0 = 0 . By employing the RK4 method, we have four increments as follows:
k 1 = h f s i , B i , k 2 = h f s i + h 2 , B i + k 1 2 , k 3 = h f s i + h 2 , B i + k 2 2 , k 4 = h f s i + h , B i + k 3 ,
thus, we obtain that
B i + 1 = B i + 1 6 k 1 + 2 k 2 + 2 k 3 + k 4 .
Now, we have the approximate solutions B i of the Riccati differential equation at each nodal point s i [ 0 , τ ] , i = 0 , 1 , , m . We denote B = [ B 0 , B 1 , , B m ] . Afterwards, this vector solution B is used to estimate the coefficients in (5). Since (5) is in integral form, in this work, we construct matrix representation for integration based on the concept of trapezoidal rule. By considering an integral function from the initial point s 0 to each point s i , i = 0 , 1 , , m , it is approximated by the trapezoidal rule. We obtain:
F ( s 0 ) : = s 0 s 0 f ( ξ ) d ξ 0 , F ( s 1 ) : = s 0 s 1 f ( ξ ) d ξ h 2 f ( s 0 ) + f ( s 1 ) , F ( s 2 ) : = s 0 s 2 f ( ξ ) d ξ h 2 f ( s 0 ) + 2 f ( s 1 ) + f ( s 2 ) , F ( s m ) : = s 0 s m f ( ξ ) d ξ h 2 f ( s 0 ) + 2 f ( s 1 ) + + 2 f ( s m 1 ) + f ( s m ) .
From these integrations, we can construct the integration matrix by
F ( s 0 ) F ( s 1 ) F ( s 2 ) F ( s m ) = 0 h 2 h 2 h 2 h h 2 h 2 h h h 2 f ( s 0 ) f ( s 1 ) f ( s 2 ) f ( s m )
and denote this by F = J f . This J is called the integration matrix, which is easily computed. We will then approximate the integral terms of A j γ ( τ ) for j N 0 = : N { 0 } in (5) using the integration matrix J , and we have
A 0 γ = e J P 0 , A j γ = e J P j J e J P j Q j A j 1 γ ,
where A j γ = A j γ ( s 0 ) , A j γ ( s 1 ) , A j γ ( s m ) , P j = P j ( s 0 ) , P j ( s 1 ) , , P j ( s m ) and Q j = Q j ( s 0 ) , Q j ( s 1 ) , , Q j ( s m ) ; the elements of P j and Q j can be directly calculated by (6). The notation ⊙ is the Hadamard product defined in [32] as the product of element-wise at the same positions in matrices. In this work, we use the exponential function of a matrix to denote the matrix whose element is the exponential of the element in that component.
Finally, we obtain the numerical formula for the pricing of the T-claim with the polynomial payoff (4) by
U E γ ( x , τ ) e B ( s m ) x j = 0 A j γ ( s m ) x γ j ,
where B ( s m ) and A j γ ( s m ) are the last components of vector solutions B and A j γ described above, respectively. Moreover, we can reduce the number of computational points m, but still preserve the accuracy by using other numerical integration approaches such as Simpson’s rule, Newton–Cotes, quadrature formula, etc., as can be seen in [33] for more details and references.

5. Experimental Validations

In this section, we verify the formula proposed in Section 3 by comparing with the Monte Carlo (MC) simulation based on the following ECIR process,
d x t = θ σ 0 2 d e 2 σ 1 t 4 θ x t d t + σ 0 e σ 1 t x t d W t .
Comparing with (3), we use θ ( t ) = θ , μ ( t ) = σ 0 2 d e 2 σ 1 t 4 θ and σ ( t ) = σ 0 e σ 1 t , where θ , σ 0 are positive numbers, σ 1 is a real number and with an integer d 2 ; thus, the Assumption 1 is satisfied. In general, the expectation for the ECIR process (24) may be directly computed using the transition density of x t presented by Egorov et al. [16]. However, this approach will not produce an accurate value of U E γ ( x , τ ) for a small τ , and to overcome this problem, the MC simulation is employed to approximate the value of U E γ ( x , τ ) .
The MC simulations presented in this paper are based on the EM scheme which are implemented by the MATLAB software to receive numerical solutions of (24) for evaluating (2). In this case, we use the trapezoidal integration for the integral term. MATLAB R2020a and a laptop with the following specifications were used in all of our calculations: Windows 10 Education, 64-bit Operating System, Intel(R) Core(TM) i7-8550U, CPU @1.80GHz, 8.0 GB RAM.

5.1. Closed-Form Formulas for CIR and ECIR Processes with MC Simulations

In this experiment, Euler–Maruyama (EM) discretization was applied for the ECIR process (24). Higham and Mao proved the accuracy of approximations by the EM scheme [34]. According to (24), we used θ = 1 , d = 2 , σ 0 = 1 , and in particular, we set σ 1 = 0 , 1 for the CIR and ECIR processes, respectively.
Example 1
(CIR case, σ 1 = 0 ). Formula (20) with α = β = 0.01 and for γ = n = 1 , 2 :
This example demonstrates the closed-form Formula (20) based on (24) in the case of γ = n N 0 and σ 1 = 0 , CIR process. To validate Formula (20) in Corollary 4, we use parameters θ = 1 , σ 0 = 0.1 and d = 2 in the process (24), and employ MC simulations with different initial values x = 0.1 , 0.2 , , 1 to generate 10 , 000 sample paths of x t , where each path consists of 10 , 000 time steps over two different time intervals 0 , 0.1 and 0 , 1 . The validations are performed through the comparisons between the Formula (20) and MC simulations based on γ = 1 , 2 .
As presented in Figure 1, the numerical results from MC simulations (colored circles) match completely with the results from Formula (20) (solid lines) for each x = 0.1 , 0.2 , , 1 . Thus, the agreement of the results has validated the accuracy of Formula (20).
Example 2
(ECIR case, σ 1 = 1 ). Formula (12) with α = β = 1 and for γ = n = 1 , 2 :
In this example, we validate the Formula (12) of U E n ( x , τ ) based on the process (24) in the case of γ = n N 0 with parameters θ = 1 , σ 0 = σ 1 = 1 and d = 2 by comparing with the MC simulation. Since the parameters are no longer constant, the solution B ( τ ) of (7) cannot be solved analytically; thus, the numerical scheme RK4 (see Section 4) is applied to approximate the solution B ( τ ) with 500 step sizes and the integral terms of the coefficients in (5) are approximated by the trapezoidal rule using 500 subintervals. MC simulations are performed using a number of sample paths, including 5000, 10,000, 20,000 and 400,000, where each path consists of 10,000 steps over different time intervals [ 0 , τ ] for τ = 0.01 , 0.1 , 1 , 2 .
Table 1 shows the comparisons between the numerical Formula (12) of U E n ( x , τ ) and MC simulation for γ = 1 , 2 at each initial value x = 0.1 , 0.2 , , 1.0 . The accuracy is measured by the mean absolute differences (MADs) over initial values x. Table 1 shows that the obtained MADs are very small and become smaller as the number of sample paths increases for all cases of γ and T. This result confirms the accuracy of the Formula (12) as compared with the MC simulations. The average run times (ARTs) of the MC simulations for different numbers of paths are displayed in Table 1. Obviously, the ARTs of MC simulations are much more than that from the formula, which is approximately 0.3 s, especially when using a large number of paths.

5.2. Numerical Approximation of the Proposed Formulas with a Finite Sum

According to Section 3.3, when Q j are non-zero for all j N 0 , the series (14) diverges. In this section, we study the level of accuracy of the formulas in Theorems 1 and Corollary 3, when the values are estimated using their partial sum. We denote U E γ , K for the approximation of U E γ by a partial sum of (4) up to order γ K . In order to find a suitable number K, before comparing with MC simulations, we need to measure the significant difference of the value of U E γ , K at each K N , which is defined by a sequence of absolute relative differences (ARDs),
D E γ , K ( x , τ ) : = U E γ , K ( x , τ ) U E γ , K 1 ( x , τ ) U E γ , K ( x , τ ) ,
for all ( x , τ ) ( 0 , ) × [ 0 , ) . Furthermore, the accuracy of U E γ , K ( x , τ ) compared with MC for fixing the suitable K is measured via the absolute relative errors (AREs) defined by
E E γ , K ( x , τ ) : = U E γ , K ( x , τ ) U E γ , M ( x , τ ) U E γ , K ( x , τ ) ,
for all ( R , τ ) ( 0 , ) × [ 0 , ) , where U E γ , M is the result of E Q x T γ e t T ( α x s + β ) d s x t = x received from MC simulations.
The sequences of ARDs of D E γ ( x , 0.01 ) are displayed in Table 2 for K = 5 , 10 , 15 , 20 with parameters θ = 1 , σ 0 = σ 1 = 1 and d = 2 , except x = 0.01 , 1, 5 for γ = 1.5 , 0.5 , 0.5 , 1.5 . In the case of infinite sum of U E γ ( x , 0.01 ) , we consider the parameters γ = 1.5 , 0.5 , 0.5 , 1.5 . According to Table 1, the received ARDs are likely improved when K increases up to K = 20 , showing that U E γ , K for these K can already produce good approximations to U E γ . To verify this claim, the results of U E γ , K with K = 10 are constructed to compare with MC simulations, and the results are shown in the next example.
Example 3.
The partial sum of (4) up to order γ K with α = β = 1 for γ = 1.5 , 0.5 , 0.5 , 1.5 .
The comparison results between the formulas U E γ , 10 ( x , 0.01 ) and MC simulations are shown in Table 3. MC simulations are performed by 5000, 10,000, 20,000, and 40,000 sample paths using 10,000 discretized steps. Table 3 shows the results of MC simulations that closely match with the approximate Formula (4) with better approximations (smaller AREs) when the number of sample paths increases. This confirms that the finite partial sum approximation of (4) is very accurate as compared by MC simulation.

6. Contingent Claims Pricing

In the context of pricing an option, assume that the underlying asset is set up to follow the ECIR process (3); we first define the following process
V t : = E Q e t T α x s + β d s f T + t T e t s α x u + β d u g s d s | F t , 0 t T ,
where f T and g are nonnegative functions. In particular, according to Karatzas and Shreve’s exercise 8.13 in [35], the process V t in (25) gives the unique wealth process with the initial wealth x; more details on this can be found in [35]. This is also called the valuation process of a contingent claim ( g , f T ) , where f T is the terminal payoff at maturity and g = g t , F t | 0 t T is the payoff rate.
This section illustrates an application for valuing the contingent claim with a date of maturity T, which depends on the underlying asset x t following the ECIR or CIR process. The analytical formulas for a contingent claim are provided in the following theorems.
Proposition 1.
Let x t follow the ECIR process (3) with α , β R and n 1 , n 2 N 0 . Suppose that f T x T n 1 and g s x s n 2 for 0 t s T , then
V t = U E n 1 ( x , τ ) + 0 τ U E n 2 ( x , v ) d v ,
for ( x , τ ) D E ( 0 , ) × [ 0 , ) , τ = T t , and U E n 1 and U E n 2 are given in Corollary 1.
Proof. 
Applying Fubini’s theorem and Corollary 1 yields
V t = E Q x T n 1 e t T α x s + β d s | x t = x + t T E Q x s n 2 e t s α x u + β d u | x t = x d s = U E n 1 ( x , τ ) + t T U E n 2 ( x , s t ) d s .
Setting v = s t obtains (26) as required. □
Remark 5.
Suppose that f T k = 0 n 1 a k x T k and g s k = 0 n 2 b k x s k , where x t follows the ECIR process with α , β R and n 1 , n 2 N 0 , for some sequences of real numbers ( a 0 , a 1 , , a n 1 ) , ( b 0 , b 1 , , b n 2 ) in which a n 1 and b n 2 are not zero. According to Proposition 1, we have
V t = E Q e t T α x s + β d s k = 0 n 1 a k x T k + t T e t s α x u + β d u k = 0 n 2 b k x s k d s | x t = x = k = 0 n 1 a k U E k ( x , τ ) + k = 0 n 2 0 τ b k U E k ( x , v ) d v .
Furthermore, for a CIR process, the above equation can readily be reduced to the following form
V t = k = 0 n 1 a k U C k ( x , τ ) + k = 0 n 2 b k j = 0 k 0 τ e B ( v ) x A j ( v ) d v x k j .
Corollary 8.
Suppose that x t follows the CIR process. According to Proposition 1 with α = 0 , β = r > 0 (also called fixed rate), and n 1 , n 2 N 0 , we have
V t = U C n 1 ( x , τ ) + 1 e ( r + n 2 θ ) τ r + n 2 θ x n 2 + j = 1 n 2 k = 1 j Q k k θ k = 0 j ( 1 ) j k + 1 j k e r + n 2 k θ τ 1 r + n 2 k θ x n 2 j ,
for ( x , τ ) D C ( 0 , ) × [ 0 , ) , τ = T t , and Q k is given in Corollary 3.
Proof. 
From (20) in Corollary 4 with α = 0 , then ρ = θ , B ( τ ) = 0 , and H j ( τ ) = e ( r + ( γ 2 j ) θ ) τ for all τ 0 and j N 0 . Recalling Remark 5, we have
V t = U C n 1 ( x , τ ) + j = 0 n 2 0 τ A j ( v ) d v x n 2 j .
First, considering 0 τ A j ( v ) d v for only j = 0 yields
0 τ A 0 ( v ) d v = 0 τ H 0 ( v ) d v = 1 e ( r + n 2 θ ) τ r + n 2 θ ,
and the remaining terms, for j = 1 , 2 , , n 2 ,
0 τ A j ( v ) d v = 0 τ k = 1 j 2 Q k k H j ( v ) e θ v 1 2 θ e θ v j d v = k = 1 j Q k k θ 0 τ e ( r + ( n 2 j ) θ ) v e θ v 1 e θ v j d v = k = 1 j Q k k θ 0 τ e ( r + n 2 θ ) v ( e θ v 1 ) j d v = k = 1 j Q k k θ k = 0 j ( 1 ) j k j k 0 τ e ( r + ( n 2 k ) θ ) v d v = k = 1 j Q k k θ k = 0 j ( 1 ) j k + 1 j k e ( r + n 2 k θ ) τ 1 r + n 2 k θ .
The benefits of these theorems to some well-known pricing instruments are shown in the following examples.
Example 4.
Zero-coupon bond.
The valuation of a zero coupon bond at time t with expiration date T, p ( t , T ) is given by the expression
p ( t , T ) = E Q e t T α x s + β d s | x t = x
where x t follows the ECIR process. Applying Corollary 1 by setting γ = 0 , we obtain the formula for the price of the zero coupon bond
p ( t , T ) = A 0 0 ( τ ) e B ( τ ) x .
In the case that x t follows the CIR process, Corollary 4 is used to produce the closed-form formula for valuing the zero coupon bond
p ( t , T ) = H 0 ( τ ) e B ( τ ) x
where
H 0 ( τ ) = exp θ 2 μ σ 2 β + θ μ ρ σ 2 τ 2 ρ ( ρ θ ) + e ρ τ ( ρ + θ ) 2 θ μ σ 2 , B ( τ ) = 2 α ( e ρ τ 1 ) ρ ( e ρ τ + 1 ) + θ ( e ρ τ 1 ) ,
and ρ = θ 2 + 2 α σ 2 .
Remark 6.
If we set α = 1 and β = 0 for the CIR process, we obtain the identical formula for the zero-coupon bond which appears in many pieces of literature.
Example 5.
Two bonds interest rate swap.
In this example, we apply the Corollary 1 for pricing the value of fixed rate for a floating swap, in which one company agrees to pay a fixed interest rate and receives in exchange a floating rate, see [21]. We consider the interest swap as the difference between the two bonds. From the point of view of the fixed ratepayer, the value of the interest rate swap, denoted by P s w a p , is P s w a p : = B f l o a t B f i x where B f l o a t is the value of floating rate bond, and B f i x is the value of fixed rate bond; see [36] for more details.
Suppose that the value of the swap is zero at the initial time t and the London Interbank Offered Rate (LIBOR), then zero rates are used as discount rates, denoted by x t , which follows the ECIR process. Then
B f i x = i = 1 N k i e t T i ( α x s + β ) d s + L e t T ( α x s + β ) d s B f l o a t = ( L + k 0 ) e t T 1 ( α x s + β ) d s
for some integer N 2 , where t is the initial time, T i is the time until the i th payment is exchanged; k t is the fixed payment made at time t; x t is the LIBOR zero rates corresponding to maturity t; and L is the notional principal in swap agreement. Thus, the value of the interest rate swap at time T is
E Q P s w a p | x t = x = E Q ( L + k 0 ) e t T 1 ( α x s + β ) d s i = 1 N k i e t T i ( α x s + β ) d s + L e t T ( α x s + β ) d s | x t = x
To calculate (28), Corollary 1 can be applied by setting γ = 0 .
Example 6.
Arrears swap
An arrears swap, also known as a delayed reset swap, is one of the traded instruments in the over-the-counter market, in which two companies or financial institutes decide to exchange periodic payments with another. In this interest rate fixed for floating swap, the floating rate paid on a payment date is based on the interest rate observed at the end of the reset period, as can be seen in [36] for more details.
Let x f i x be a fixed rate, x t be a floating rate at time t, and P be a notional principle. Suppose that an arrears swap has an expiration date T with N payment dates at t = T 0 < T 1 < < T N = T in an increment of Δ t = T i T i 1 , i = 1 , 2 , , N . The payoff of such a swap from a floating rate payer’s point of view at the i th payment date, V i a r , is the difference between interest in a notional principle considered by the fixed and floating interest rates, which can be expressed in the form V i a r = x f i x x T i Δ t P . By the fundamental theorem of asset pricing [3], a no-arbitrage price at any time t of the arrears swap, V ar , is the conditional expectation of the sum of each payoff discounted to the initial time t = 0 , which is
V a r = E Q i = 1 N V i a r e t T i α x s + β d s x t = x = Δ t P x f i x i = 1 N E Q e t T i α x s + β d s x t = x i = 1 N E Q x T i e t T i α x s + β d s x t = x .
By applying Corollary 1 and setting γ = 0 and γ = 1 , the value of the arrears swap (29) can be obtained as an analytical form. It should be noted that the fair value for paying the fixed rate is
x f i x = i = 1 N E Q x T i e t T i α x s + β d s x t = x i = 1 N E Q e t T i α x s + β d s x t = x = i = 1 N U E 1 ( x , i Δ t ) i = 1 N U E 0 ( x , i Δ t ) .

7. Conclusions

In this work, we proposed the analytical formula for a contingent claim with the polynomial payoff under the ECIR process represented as the conditional expectation of the product of polynomial and exponential functions, E Q x T γ e t T ( α x s + β ) d s x t = x where α , β , γ R and x t follow the ECIR process (3). By solving (8) from the Feynman–Kac representation, the analytical formula of (2) for the ECIR process was constructed in terms of an infinite sum of analytical expressions in Theorem 1. Interestingly, the infinite sum is reduced to a finite sum if γ N 0 in Corollaries 1 and 2. Under the CIR process, the parameters in (5) can be evaluated and the Riccati differential Equation (6) can be analytically solved; thus, the formula in Theorem 1 for (2) can be rapidly deduced to the formula in Corollary 3. In addition, as shown in Corollaries 4 and 5, the formulas can be stated in closed form. The formula for the unconditional expectations are also obtained by taking T , as can be seen in Corollaries 6 and 7.
The proposed closed-form formulas validated the validity and efficiency by comparing the results with Monte Carlo simulations as illustrated in Section 5. The applications of the proposed formulas for contingent claims are described in Section 6 for financial products such as zero-coupon bond, two bonds interest rate swap, an arrears swap. For these applications under the ECIR process (3), the formula of (2) is extended for a more general form, E Q e t T α x s + β d s f T + t T e t s α x u + β d u g s d s | F t , for some real value functions f T and g s where 0 t < s T . The proposed formulas can also be considered as generalized results of other formulas which appeared in the literature.

Author Contributions

Conceptualization, F.N. and K.M.; methodology, F.N. and K.M.; software, F.N.; validation, K.M.; formal analysis, K.M.; investigation, F.N.; writing—original draft preparation, F.N.; writing—review and editing, K.M.; visualization, F.N.; supervision, K.M.; project administration, K.M. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

Authors thank the Development and Promotion of Science and Technology Talents (DPST), which is a project of the Institute for the Promotion of Teaching Science and Technology (IPST) for their kind support.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
ARDabsolute relative difference
ARTaverage run time
CIRCox–Ingersoll–Ross
ECIRextended Cox–Ingersoll–Ross
EMEuler–Maruyama
FRAforward rate agreement
IRSinterest rate swap
LIBORLondon Interbank Offered Rate
MADmean absolute difference
MCMonte Carlo
ODEordinary differential equation
OUOrnstein–Uhlenbeck
PDEpartial differential equation
PDFprobability density function
RK4fourth-order Runge–Kutta
SDEstochastic differential equation

References

  1. Ben-Ameur, H.; Breton, M.; Karoui, L.; L’Ecuyer, P. A dynamic programming approach for pricing options embedded in bonds. J. Econ. Dyn. Control 2007, 31, 2212–2233. [Google Scholar] [CrossRef] [Green Version]
  2. Grasselli, M. The 4/2 stochastic volatility model: A unified approach for the Heston and the 3/2 model. Math. Financ. 2017, 27, 1013–1034. [Google Scholar] [CrossRef]
  3. Eide, I.B. Arbitrage and asymmetric information. Preprint Ser. Pure Math. 2008, 24, 1–22. [Google Scholar]
  4. Hull, J.; White, A. Pricing interest-rate-derivative securities. Rev. Financ. Stud. 1990, 3, 573–592. [Google Scholar] [CrossRef]
  5. Duffie, D.; Pan, J.; Singleton, K. Transform analysis and asset pricing for affine jump-diffusions. Econometrica 2000, 68, 1343–1376. [Google Scholar] [CrossRef] [Green Version]
  6. Ahmadian, D.; Ballestra, L. Pricing geometric Asian rainbow options under the mixed fractional Brownian motion. Phys. Stat. Mech. Its Appl. 2020, 555, 124458. [Google Scholar] [CrossRef]
  7. Delbaen, F.; Schachermayer, W. A general version of the fundamental theorem of asset pricing. Math. Ann. 1994, 300, 463–520. [Google Scholar] [CrossRef]
  8. Ahmadian, D.; Ballestra, L.V. The finite element method: A high-performing approach for computing the probability of ruin and solving other ruin-related problems. Math. Methods Appl. Sci. 2021, 44, 12640–12653. [Google Scholar] [CrossRef]
  9. Liang, J.; Zou, H. Valuation of credit contingent interest rate swap with credit rating migration. Int. J. Comput. Math. 2020, 97, 2546–2560. [Google Scholar] [CrossRef]
  10. Pearson, K. Tables for Statisticians and Biometricians; University Press: Cambridge, UK, 1914. [Google Scholar]
  11. Ditlevsen, S.; Rubio, A.C.; Lansky, P. Transient dynamics of Pearson diffusions facilitates estimation of rate parameters. Commun. Nonlinear Sci. Numer. Simul. 2020, 82, 105034. [Google Scholar] [CrossRef]
  12. Feller, W. Two singular diffusion problems. Ann. Math. 1951, 54, 173–182. [Google Scholar] [CrossRef]
  13. Cox, J.C.; Ingersoll, J.E., Jr.; Ross, S.A. An intertemporal general equilibrium model of asset prices. Econom. J. Econom. Soc. 1985, 53, 363–384. [Google Scholar] [CrossRef] [Green Version]
  14. Hansen, L.P.; Scheinkman, J.A. Back to the future: Generating moment implications for continuous-time Markov processes. Econom. J. Econom. Soc. 1993, 63, 767–804. [Google Scholar]
  15. Maghsoodi, Y. Solution of the extended CIR term structure and bond option valuation. Math. Financ. 1996, 6, 89–109. [Google Scholar] [CrossRef]
  16. Egorov, A.V.; Li, H.; Xu, Y. Maximum likelihood estimation of time-inhomogeneous diffusions. J. Econom. 2003, 114, 107–139. [Google Scholar] [CrossRef]
  17. Lamberton, D.; Lapeyre, B. Introduction to Stochastic Calculus Applied to Finance; CRC Press: Boca Raton, FL, USA, 2007. [Google Scholar]
  18. Dufresne, D. The Integrated Square-Root Process; Research Paper No. 90; University of Melbourne: Melbourne, Australia, 2001; pp. 1–34. [Google Scholar]
  19. Rujivan, S. A closed-form formula for the conditional moments of the extended CIR process. J. Comput. Appl. Math. 2016, 297, 75–84. [Google Scholar] [CrossRef]
  20. Sutthimat, P.; Mekchay, K.; Rujivan, S. Explicit Formula for Conditional Expectations of Product of Polynomial and Exponential Function of Affine Transform of Extended Cox-Ingersoll-Ross Process. J. Physics Conf. Ser. 2018, 1132, 012083. [Google Scholar] [CrossRef]
  21. Mallier, R.; Alobaidi, G. Interest rate swaps under CIR. J. Comput. Appl. Math. 2004, 164, 543–554. [Google Scholar] [CrossRef] [Green Version]
  22. Moreno, M.; Platania, F. A cyclical square-root model for the term structure of interest rates. Eur. J. Oper. Res. 2015, 241, 109–121. [Google Scholar] [CrossRef]
  23. Frenkel, D.; Portugal, R. Algebraic methods to compute Mathieu functions. J. Phys. Math. Gen. 2001, 34, 3541. [Google Scholar] [CrossRef]
  24. Thamrongrat, N.; Rujivan, S. An analytical formula for pricing interest rate swaps in terms of bond prices under the extended Cox-Ingersoll-Ross model. Songklanakarin J. Sci. Technol. 2020, 43, 987–992. [Google Scholar]
  25. Chihara, T.S. An Introduction to Orthogonal Polynomials; Courier Corporation: New York, NY, USA, 2011. [Google Scholar]
  26. Leonenko, G.M.; Phillips, T.N. High-order approximation of Pearson diffusion processes. J. Comput. Appl. Math. 2012, 236, 2853–2868. [Google Scholar] [CrossRef] [Green Version]
  27. Chumpong, K.; Mekchay, K.; Thamrongrat, N. Analytical formulas for pricing discretely-sampled skewness and kurtosis swaps based on Schwartz’s one-factor model. Songklanakarin J. Sci. Technol 2021, 43, 1–6. [Google Scholar]
  28. Sutthimat, P.; Rujivan, S.; Mekchay, K.; Rakwongwan, U. Analytical formula for conditional expectations of path-dependent product of polynomial and exponential functions of extended Cox–Ingersoll–Ross process. Res. Math. Sci. 2022, 9, 1–17. [Google Scholar] [CrossRef]
  29. Sutthimat, P.; Mekchay, K. Closed-form formulas for conditional moments of inhomogeneous Pearson diffusion processes. Commun. Nonlinear Sci. Numer. Simul. 2022, 106, 106095. [Google Scholar] [CrossRef]
  30. Boonklurb, R.; Duangpan, A.; Rakwongwan, U.; Sutthimat, P. A Novel Analytical Formula for the Discounted Moments of the ECIR Process and Interest Rate Swaps Pricing. Fractal Fract. 2022, 6, 58. [Google Scholar] [CrossRef]
  31. Burden, R.L.; Faires, J.D.; Burden, A.M. Numerical Analysis; Cengage Learning: Boston, MA, USA, 2015. [Google Scholar]
  32. Caro-Lopera, F.J.; Leiva, V.; Balakrishnan, N. Connection between the Hadamard and matrix products with an application to matrix-variate Birnbaum-Saunders distributions. J. Multivar. Anal. 2012, 104, 126–139. [Google Scholar] [CrossRef] [Green Version]
  33. Boonklurb, R.; Duangpan, A.; Treeyaprasert, T. Modified finite integration method using Chebyshev polynomial for solving linear differential equations. J. Numer. Ind. Appl. Math. 2018, 12, 1–19. [Google Scholar]
  34. Higham, D.J.; Mao, X. Convergence of Monte Carlo simulations involving the mean-reverting square root process. J. Comput. Financ. 2005, 8, 35–61. [Google Scholar] [CrossRef] [Green Version]
  35. Karatzas, I.; Shreve, S. Brownian Motion and Stochastic Calculus; Springer: Berlin/Heidelberg, Germany, 2014; Volume 113. [Google Scholar]
  36. Hull, J. Option, Future, & Other Derivatives; Prentice Hall: Hoboken, NJ, USA, 2002. [Google Scholar]
Figure 1. The validations of U C n ( x , τ ) for τ = 0.1 , 0.5 , 1 , 1.5 and x = 0.1 , 0.2 , , 1.0 : (a) the first conditional moments; and (b) the second conditional moments.
Figure 1. The validations of U C n ( x , τ ) for τ = 0.1 , 0.5 , 1 , 1.5 and x = 0.1 , 0.2 , , 1.0 : (a) the first conditional moments; and (b) the second conditional moments.
Symmetry 14 00933 g001
Table 1. The MADs between estimated solutions of Formula (12) and MC simulations.
Table 1. The MADs between estimated solutions of Formula (12) and MC simulations.
γ No. of Paths τ ART (s)
0.01 0.1 12
150008.756 × 10 4 2.428 × 10 3 1.884 × 10 3 4.578 × 10 4 18.26
10,0005.149 × 10 4 1.386 × 10 3 1.013 × 10 3 3.746 × 10 4 39.23
20,0003.129 × 10 4 8.440 × 10 4 7.492 × 10 3 2.363 × 10 4 73.19
40,0001.791 × 10 4 7.463 × 10 4 6.305 × 10 4 1.436 × 10 4 134.27
250001.599 × 10 3 3.162 × 10 3 9.097 × 10 3 6.633 × 10 3 18.34
10,0006.147 × 10 4 2.182 × 10 3 5.019 × 10 3 3.873 × 10 3 39.35
20,0002.980 × 10 4 1.643 × 10 3 3.390 × 10 3 2.865 × 10 3 73.63
40,0002.773 × 10 4 7.568 × 10 4 2.459 × 10 3 2.305 × 10 3 135.69
Table 2. The ARDs D E γ , K ( x , 0.01 ) .
Table 2. The ARDs D E γ , K ( x , 0.01 ) .
x K γ
−1.5−0.50.51.5
0.152.605 × 10 4 2.416 × 10 6 2.985 × 10 8 4.980 × 10 9
104.959 × 10 6 1.262 × 10 8 3.498 × 10 11 9.897 × 10 13
158.822 × 10 7 1.030 × 10 9 1.226 × 10 12 1.375 × 10 14
207.289 × 10 7 4.866 × 10 10 3.202 × 10 13 1.912 × 10 15
152.929 × 10 9 2.445 × 10 11 3.019 × 10 13 5.490 × 10 14
105.576 × 10 16 1.277 × 10 18 3.539 × 10 21 1.091 × 10 22
159.921 × 10 22 1.043 × 10 24 1.240 × 10 27 1.515 × 10 29
208.196 × 10 27 4.926 × 10 30 3.238 × 10 33 2.107 × 10 35
559.460 × 10 13 7.834 × 10 15 9.672 × 10 17 1.772 × 10 17
105.763 × 10 23 1.309 × 10 25 3.627 × 10 28 1.127 × 10 29
153.281 × 10 32 3.421 × 10 35 4.068 × 10 38 5.012 × 10 40
208.674 × 10 41 5.170 × 10 44 3.399 × 10 47 2.230 × 10 49
Table 3. AREs E E γ , 10 ( x , 0.01 ) of approximations U E γ , 10 ( x , 0.01 ) and MC simulations.
Table 3. AREs E E γ , 10 ( x , 0.01 ) of approximations U E γ , 10 ( x , 0.01 ) and MC simulations.
x No. of Paths γ
1 . 5 0 . 5 0 . 5 1 . 5
0.150008.879 × 10 3 1.724 × 10 3 1.904 × 10 3 5.297 × 10 3
10,0004.499 × 10 3 1.439 × 10 3 1.209 × 10 3 2.196 × 10 3
20,0001.089 × 10 3 6.156 × 10 4 5.260 × 10 4 1.532 × 10 3
40,0009.658 × 10 4 5.739 × 10 4 1.281 × 10 4 8.521 × 10 4
150004.498 × 10 3 7.423 × 10 4 7.291 × 10 4 2.107 × 10 3
10,0001.681 × 10 3 7.044 × 10 4 6.195 × 10 4 9.662 × 10 4
20,0005.316 × 10 4 5.866 × 10 4 3.168 × 10 4 9.356 × 10 4
40,0001.468 × 10 4 1.600 × 10 4 2.496 × 10 4 7.086 × 10 4
550001.260 × 10 3 1.531 × 10 4 2.136 × 10 4 4.115 × 10 4
10,0008.906 × 10 4 1.327 × 10 4 1.785 × 10 4 3.865 × 10 4
20,0008.220 × 10 4 7.713 × 10 5 1.128 × 10 4 1.216 × 10 4
40,0009.445 × 10 5 1.588 × 10 5 3.812 × 10 5 8.439 × 10 5
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Nualsri, F.; Mekchay, K. Analytically Pricing Formula for Contingent Claim with Polynomial Payoff under ECIR Process. Symmetry 2022, 14, 933. https://doi.org/10.3390/sym14050933

AMA Style

Nualsri F, Mekchay K. Analytically Pricing Formula for Contingent Claim with Polynomial Payoff under ECIR Process. Symmetry. 2022; 14(5):933. https://doi.org/10.3390/sym14050933

Chicago/Turabian Style

Nualsri, Fukiat, and Khamron Mekchay. 2022. "Analytically Pricing Formula for Contingent Claim with Polynomial Payoff under ECIR Process" Symmetry 14, no. 5: 933. https://doi.org/10.3390/sym14050933

APA Style

Nualsri, F., & Mekchay, K. (2022). Analytically Pricing Formula for Contingent Claim with Polynomial Payoff under ECIR Process. Symmetry, 14(5), 933. https://doi.org/10.3390/sym14050933

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop