Journal of Probability and Statistics
Volume 2015 (2015), Article ID 369053, 24 pages
http://dx.doi.org/10.1155/2015/369053

Research Article

Polynomial Chaos Expansion Approach to Interest Rate Models

Luca Di Persio,¹ Gregorio Pellegrini,¹ and Michele Bonollo2,3

¹Department of Computer Science, University of Verona, Strada le Grazie 15, 37134 Verona, Italy
²Iason Ltd., Milan, Italy
³IMT Lucca, Piazza San Francesco 19, 55100 Lucca, Italy

Received 30 June 2015; Accepted 19 October 2015

Academic Editor: Z. D. Bai

Copyright © 2015 Luca Di Persio et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The Polynomial Chaos Expansion (PCE) technique allows us to recover a finite second-order random variable exploiting suitable linear combinations of orthogonal polynomials which are functions of a given stochastic quantity , hence acting as a kind of random basis. The PCE methodology has been developed as a mathematically rigorous Uncertainty Quantification (UQ) method which aims at providing reliable numerical estimates for some uncertain physical quantities defining the dynamic of certain engineering models and their related simulations. In the present paper, we use the PCE approach in order to analyze some equity and interest rate models. In particular, we take into consideration those models which are based on, for example, the Geometric Brownian Motion, the Vasicek model, and the CIR model. We present theoretical as well as related concrete numerical approximation results considering, without loss of generality, the one-dimensional case. We also provide both an efficiency study and an accuracy study of our approach by comparing its outputs with the ones obtained adopting the Monte Carlo approach, both in its standard and its enhanced version.

1. Introduction

In this paper, we use the Polynomial Chaos Expansion (PCE) technique to study the strong solution of the following Itô stochastic differential equation (SDE):

where is the initial value of the unknown stochastic process , and the terms and represent the drift and the volatility of the process, respectively, while is a standard Brownian motion and it is considered up to a certain finite time to which we refer, taking into account the financial setting, as the maturity time, for example, of some underlying. According to [1, Section 4.5], if the coefficients and are smooth enough, then we have the existence and uniqueness of the solution to (1).

The PCE approach to the aforementioned problem belongs to the so-called Uncertain Quantification (UQ) methods, which constitute a fundamental step towards the validation of numerical techniques to be used for critical decisions; see, for example, [2, 3] and the references therein. The term polynomial chaos was originally introduced by Wiener in his 1938 paper [4] in which he applies his generalized harmonic analysis (see, e.g., [5]). Then, Ghanem and Spanos in [6, Section 2.4, Subsection 3.3.6] combined the PCE with finite element method to compute the solution of PDEs in the presence of uncertain parameters. Then, such an approach has been applied in several frameworks. For example, for the simulation of probabilistic chemical reactions, see [7], for the stochastic optimal trajectories generation, see [8], for the sensitivity analysis, see [9]. In order to refine specific numerical methods, see [10], and also for a rather large set of problems arising in engineering and computational fluid dynamics (CFD), see, for example, [11, Chapter 6], [12, 13], and references therein.

Strictly speaking, PCE recovers a random variable in terms of a linear combination of functionals whose entries are known random variables called germs or basic variables. Several methods are available to compute these coefficients. In particular, we choose the Nonintrusive Spectral Projection (NISP) method which employs a set of deterministic realizations; see, for example, [11, Chapter 3] for further details.

In order to show the advantages of the PCE approach, the results of our PCE implementation have been compared versus the ones obtained using both the standard Monte Carlo (MC) and the quasi-Monte Carlo (QMC) techniques. The former is based on pseudorandom sampling and its convergence properties basically rely on the law of large numbers and the central limit theorem. The latter uses low-discrepancy sequences to simulate the process and it is theoretically based on the Koksma-Hlawka inequality; see [14], which provides error bounds for computations involving QMC method; see, for example, [15, Chapter 5].

Our analysis has been focused on three well known one-dimensional SDEs whose solutions are related to the Geometric Brownian Motion equity model; see, for example, [16, Chapter 11, Section 3], the Vasicek model, [17] and the Cox-Ingersoll-Ross (CIR) interest rate model, and [18]. We also refer to [19, Chapter 3, Section 3.2.3], for more details and references about such models.

We underline that the first two cases are meaningful because they have an analytical solution to the related SDE, so that we provide a rigorous study on convergence property of the PCE method applied to them. In the CIR case, as we do not have an analytical solution, the application of the PCE technique allows us also to derive some considerations on the convergence properties of the related approximating solutions.

The paper is organized as follows: in Section 2, the PCE machinery is described in a general setting, while Section 3 deals with the NISP approach. Section 4 points out how to apply PCE in order to decompose the solution of the considered SDE, while Section 5 is focused on the relation between the usual numerical method for solving (1) and the PCE technique. In Sections 6, 7, and 8, we consider numerical applications of the PCE approach to approximate the GBM equity and the Vasicek and the CIR model. In Section 9, we give an overview of the PCE results in terms of accuracy and computational time cost for everyone of the aforementioned considered models. Eventually, in the Appendix, we provide a table with the evaluation of the first six Hermite polynomials on a set of 40 equally spaced points in the interval .

2. Polynomial Chaos Expansion

Let be a probability space, where is the set of elementary events, is a -algebra of subsets of , and is a probability measure on .

We consider the Hilbert space of scalar real-valued random variables , whose generic element is a random variable defined on , such that

Notice that, as for Lebesgue spaces, the elements are equivalent classes of random variables. Note that is a Hilbert space endowed with the following scalar product:

with

being the usual norm. To shorten the notation, , and the related convergence has been always referred to as mean square convergence or strong convergence.

Among elements in , there is the class of basic random variables, which is used to decompose, as entries of functionals, the quantity of interest such as a random variable of interest, the solution of the SDE at time . We notice that not all the functions can be used to such a decomposition because they have to satisfy, for example, [20, Section 3], at least the following two properties:

(i)has finite raw moments of all orders;(ii)the distribution function of the basic random variable is continuous, with being its probability density function (pdf).

Because the increments of the Brownian motion are independent and normally distributed, from now on, we consider as basic random variable. The latter choice is motivated in Remark 1 in Section 3.1.

The elements in can be gathered in two groups: in the first one, we have the basic random variables; let us indicate them by , ruling the decomposition; for simplicity, let us call it , while the second is composed of generic elements; let us say , which we want to decompose using the elements of the first set, which is referred to as .

Let us denote by the -algebra generated by the basic random variable ; hence, . If we want to polynomially decompose the random variable in terms of , then has at least to be measurable with respect to the -algebra . Exploiting the Doob-Dynkin lemma, [21, Lemma 1.13], we have that is measurable, by detecting a Borel measurable function , such that . In what follows, without loss of generality, we restrict ourselves to consider the decomposition in . Moreover, the basic random variable determines the class of orthogonal polynomials indicated as the generalized polynomial chaos (gPC) basis. We underline that their orthogonality property is detected by means of the measure induced by in the image space , where denotes the Borel -algebra of . In particular, for each , we have

Because , the related set is represented by the family of Hermite polynomials defined on the whole real line; namely, , and

Figure 1 provides the graph of the first six orthonormal polynomials, achieved by scaling each in (6) by its norm in ; namely, , and is divided by . See also the Appendix where these polynomials are evaluated on equally spaced points in .

Figure 1: The Hermite polynomials up to degree .

Such polynomials are a maximal system in ; therefore, every finite second-order random variable can be approximated as follows:

for suitable coefficients depending on the random variable itself. We refer to (7) as the truncated PCE, at degree , of . Exploiting the previous definitions, taking , and considering the orthogonality property of the polynomials , we have

and as we also obtain

Moreover, converges in mean square sense into ; see, for example, [20, Section 3.1].

We underline that the PCE of approximates the -statistics using the coefficients appearing in (7); for example, the first two centered moments are determined by

3. Nonintrusive Spectral Projection

We consider the case in which a given functional governs the dynamic of a quantity we want to study. We assume that is characterized by a random output depending on the input , which is a random variable too. The goal is to describe the behavior of in dependent of the variations of and with respect to the action of the functional . In what follows, without loss of generality, we consider the 1-dimensional case, because the results we will obtain can be easily generalized to greater dimensions. In our setting, represents the solution functional to an SDE describing a specific interest rate model; see later Sections 6, 7, and 8 for further details, and we also allow the output to depend on a set of parameters; let us call it , with , which, in fact, characterizes the interest rate dynamic; namely,

In order to shorten the notation, the parameters are omitted in the definition. Our aim is to exploit the so-called Nonintrusive Spectral Projection (NISP) method to obtain the truncated PCE of the output , taking into account that the basic random variable inherits all the assumptions made in the previous section. Therefore, we have the following spectral projection:

where the coefficients are defined as in (8).

We recall that, for example, [11, Section 3.3], the NISP approach computes the scalar product in (8) by Gaussian quadrature formulae, in particular, , and we have the following Gauss-Hermite type result:

where and are the quadrature nodes and the weights of the one-dimensional Gaussian integration rule. Moreover, (14) requires the evaluation of the process at a well defined set of realizations for the input random variable . By the very definition of NISP, achieves spectral convergence with respect to for increasing degree ; see, for example, [11, Appendix B].

We note that the number of quadrature points is not linked to the degree of the truncation, a priori. Instead, the precision of the coefficients is related with ; thus a reasonable and standard choice which we have been using in what follows is to take , because the polynomials involved are at least of degree .

Moreover, in the following, we refer to as quadrature nodes as well as a set of independent realizations of the basic random variable.

3.1. Flowchart of the NISP Computation

In order to better explain the NISP approach, let us describe it using the associated flowchart in Figure 2.

Figure 2

We underline that the set of orthogonal polynomials, which constitutes a basis in , is computed off-line once the basic random variable has been given independently on the particular model we want to decompose. Analogously, the same happens for the set of realizations of and weights needed to apply the Gaussian quadrature method.

Different situation concerns the simulations of the process and related computations. In fact, in this case, we have to deal with the data of the type , which are sensible to the changes of the model parameters. This is why such a difference has been highlighted in Figure 2 using yellow and blue rectangles, which represent the off-line computations and the model correlated computations, respectively.

Remark 1. Let us focus our attention on the basic random variable. Other possibilities are available. For example, and , but the Gaussian choice is motivated by the fact that the Wiener increments are distributed as Gaussian variables. In particular, in our computations, we consider . Latter choice has been imposed by software restrictions. In particular, we have used the Scilab toolbox which provides a powerful environment to implement PCE decomposition, also giving to the user a very flexible and efficient tool for the concrete computations involved by the PCE method itself.
In particular, such toolbox implements only Hermite polynomials, which requires such input random variable in order to satisfy (5); see, for example, [22, 23], for further details.

Remark 2. The sampling of size of the PCE is achieved by detecting a unique sampling of size , of the basic random variable , and then we use (9), to each realization; namely,

and therefore the collection

is the required PCE sampling.

4. Polynomial Chaos Expansion and Stochastic Differential Equations

We consider a stochastic process which satisfies the following SDE:

Our aim is to write the PCE for the process at a given positive and finite time , namely the random variable which represents the process at that time. Without loss of generality, in what follows we consider the one-dimensional version of (17); hence, is a scalar, real-valued random variable. In particular, we focus our attention on three specific models of the type described by (17), namely, the Geometric Brownian Motion (GBM) and the Vasicek and CIR interest rate models. Because in each of the aforementioned models the random variable has finite variance, there exists a suitable probability space , where is well defined such that . It follows that the PCE method can be applied to approximate provided definition of the functional such that , for a suitable choice of the basic variable as seen in Section 3.1.

5. Usual Numerical Methods for SDE and PCE

Numerical methods for solving (1) exploit the independent increments of the Wiener process. For instance, by considering the Euler-Maruyama scheme, for example, [1, Section 10.2], for each , we have

where and is a set of increasing time steps . Moreover, is the value of the process at starting time and . To shorten the notation, let for each time step; hence, at each time step, the solution can be written as follows:

therefore adding a new independent random variable for every , . Hence, following the general approach discussed in Section 3.1, in order to apply the PCE method, we consider a set of i.i.d. Gaussian random variables . Then, the solution of (17) at time is

or equivalently

where is the aforementioned set of i.i.d. and we are in position to apply the PCE decomposition; see, for example, [20, Section 3.2] for further details concerning the multivariate decomposition. We underline that, even if we consider the Euler-Maruyama scheme, previous approach can be applied considering any other method which is based on independent increments of the Brownian motion, such as Milstein method; see, for example, [1, Section 10.2, Section ].

Let us recall that to get reliable results has to be great enough; usually , which implies a high computational cost as the complexity of the NISP method increases dramatically when more than inputs are involved. The latter problem is referred to as the curse of dimensionality; see [11].

6. Geometric Brownian Motion

In what follows, we analyze the Geometric Brownian Motion (GBM) defined as the solution of the following stochastic differential equation:

where and is the usual Wiener process. In particular, let us focus our attention on for .

6.1. PCE Approximation

The PCE method is based on the definition of a suitable process . We consider (22), taking into account the dependence of the process on the basic random variable .

Following the approach discussed in [24, Chapter 2] and [25, Chapter 4], instead of solving (22), we have studied a suitable transformation of it, in order to reduce the instability in numerical simulation. By defining and applying the Itō-Döebling lemma, we get

and once is computed, for any time , the original process is defined as

Because

the solution of (22) can be written through a functional , such that

Moreover, it embeds the transformed process (23) combined with (24). In particular, exploiting (25), it is natural to set as a Gaussian random variable. More specifically, we take due to the restrictions imposed by the software we have used to numerically implement the PCE method; see Section 3.1. We note that if the GBM model is used to describe the behaviour of the interest rate , for example, considering interest rate for equity markets, the related output at fixed time , that is, , depends not only on but also on other parameters, gathered by , which characterize its dynamics. To shorten the notation, in what follows, the dependence on is omitted.

Turning back to the analysis of , we apply the NISP method, exploiting (14), in order to obtain the truncated following PCE:

In particular, , where , is the evaluation of the process (26) at , defined by the Gaussian quadrature formula as stated in (14).

6.2. Numerical Application

In what follows, we give a PCE approximation of the solution to the SDE (22), for some particular values of the parameters involved; see Table 1.

Table 1: Parameters of GBM.

Parameters

Values
3%
15%
100
1

As (22) has an analytical solution, such a solution constitutes our benchmark; in particular, at time , we have

whose mean and variance are, respectively,

Figure 3 and Table 2 show the absolute errors of the mean and the variance; namely,

while Figure 4 and Table 2 show the absolute value of the relative error. Due to the logarithmic scale on 3, we consider their absolute values, allowing us to highlight the related order, rather than its numerical value. Therefore,

We note that such data have been computed for an increasing set of degrees .

Table 2: Absolute and relative errors of average and variance for PCE approximation of GBM at time , whose parameters are , , and starting value .

Degree of PCE
Average error
Variance error
Average relative error
Variance relative error

1
4.3214e − 03
6.2664e + 00
4.1937e − 05
2.5935e − 02

2
9.5314e − 06
6.2421e − 02
9.2497e − 08
2.5834e − 04

3
1.5686e − 07
4.3996e − 04
1.5223e − 09
1.8209e − 06

4
1.2075e − 07
2.5327e − 05
1.1718e − 09
1.0482e − 07

5
1.3828e − 07
2.8532e − 05
1.3420e − 09
1.1809e − 07

6
8.2988e − 08
2.4161e − 05
8.0536e − 10
9.9995e − 08

7
1.3581e − 07
2.7980e − 05
1.3179e − 09
1.1580e − 07

8
5.6821e − 08
2.4997e − 05
5.5142e − 10
1.0345e − 07

9
1.1640e − 07
2.7916e − 05
1.1296e − 09
1.1554e − 07

10
3.9095e − 08
2.5313e − 05
3.7939e − 10
1.0476e − 07

11
1.1656e − 07
2.7612e − 05
1.1312e − 09
1.1428e − 07

12
2.5694e − 08
2.5507e − 05
2.4934e − 10
1.0557e − 07

13
1.3620e − 07
2.7877e − 05
1.3217e − 09
1.1537e − 07

14
1.5107e − 08
2.5631e − 05
1.4661e − 10
1.0608e − 07

15
1.3547e − 07
2.7700e − 05
1.3146e − 09
1.1464e − 07

(a)

(b)

Figure 3: Semilogy scale plot of the absolute error of the mean (a) and the variance (b) computed via PCE for GBM at time , whose parameters are , , and starting value , for a set of degrees .

(a)

(b)

Figure 4: Semilogy scale plot of the absolute value of the relative error of the mean (a) and the variance (b) computed via PCE for GBM at time , whose parameters are , , and starting value , for a set of degrees .

Both average and variance errors are stationary, despite an expected spectral convergence for increasing degree . This can be motivated by the presence of an error, not due to PCE method that corrupts the expected spectral converge of . There are only two possible causes: the approximation of (22) with (26) and the error coming from numerical computations of used in (14).

Let us focus our attention on the latter: as displayed in Figures 5 and 6 and Table 3, by increasing the precision used to compute , the errors decrease. Hence, PCE errors are influenced by such approximation which is required to compute the coefficients. In Figure 5, we can see the spectral convergence of the error for the PCE approximation up to 4th degree for mean, respectively, up to the 6th degree for variance. Note that, despite the increased accuracy, the error displays again a stationary behavior, but actually no improvement on numerical methods can be implemented, which allows us to say that this is the only source of error.

Table 3: Absolute error of the average and the variance of PCE approximation of GBM at time (whose parameters are , , and starting value ) for higher precision at computing .

Degree of PCE
Average error
Variance error
Average relative error
Variance relative error

1
4.3212e − 03
6.2663e + 00
4.1935e − 05
2.5934e − 02

2
9.7228e − 06
6.2400e − 02
9.4355e − 08
2.5826e − 04

3
1.5941e − 08
4.1085e − 04
1.5470e − 10
1.7004e − 06

4
3.3636e − 10
2.0675e − 06
3.2642e − 12
8.5570e − 09

5
3.1876e − 10
9.3317e − 09
3.0934e − 12
3.8621e − 11

6
3.1841e − 10
9.8433e − 10
3.0900e − 12
4.0739e − 12

7
3.1777e − 10
8.6069e − 10
3.0838e − 12
3.5622e − 12

8
3.1750e − 10
8.8403e − 10
3.0812e − 12
3.6587e − 12

9
3.1645e − 10
8.2571e − 10
3.0709e − 12
3.4174e − 12

10
3.1707e − 10
8.9184e − 10
3.0770e − 12
3.6911e − 12

11
3.1784e − 10
8.7499e − 10
3.0845e − 12
3.6213e − 12

12
3.1901e − 10
9.1401e − 10
3.0958e − 12
3.7828e − 12

13
3.1812e − 10
8.6624e − 10
3.0872e − 12
3.5851e − 12

14
3.1837e − 10
8.9733e − 10
3.0896e − 12
3.7138e − 12

15
3.1686e − 10
8.6797e − 10
3.0749e − 12
3.5923e − 12

(a)

(b)

Figure 5: Semilogy scale plot for the absolute error of the mean (a) and variance (b) computed via PCE of GBM at time (whose parameters are , , and starting value ) with higher precision at computing .

(a)

(b)

Figure 6: Semilogy scale plot for the absolute value of the relative error of the mean (a) and variance (b) computed via PCE of GBM at time (whose parameters are , , and starting value ) with higher precision at computing .

Coming back to Table 1, let us compute a sampling of size of for , achieved with Latin Hypercube Sampling (LHS) technique; see [26, 27] and also [28] for further references, which is compared with the probability density function of the analytical solution to (22) at time , whose distribution is lognormal; namely,

and computed values shown in Figure 7 were the standard Monte Carlo sampling of shown.

(a)

(b)

Figure 7: Probability density function of GBM at time , whose parameters are , , and starting value (blue curve) and histogram of a sampling (size = ) of the for (a). (b) displays the analytical probability density function of GBM, for the same parameters values and a standard Monte Carlo sampling of size of the GBM-analytical solution.

Moreover, we have compared the two quantiles , where and , of the PCE approximation , writing their analytical values for the following set of degree: . Recall that the standard statistics for is the th sample quantile, namely, the th realization of the sampling of such that , where is the size of the sampling, and denotes the integer part of the real number within the brackets. In our analysis, we consider , while ; see, for example, [29] for further details.

In particular, we are aiming at computing the absolute error of with respect to the analytical values, which are

where represents the quantile of the normal random variable of zero mean and unitary variance. Moreover, the values of Table 1 are considered; thus the errors are defined as

for and . Obtained values are shown in Figure 8 and in Table 4.

Table 4: Absolute errors of the two quantiles and of the PCE approximation of GBM at time , whose parameters are , , and starting value .

Degree of PCE

1
5.7682e + 00
1.1964e + 01

2
4.3064e − 01
1.9530e + 00

3
8.0301e − 02
7.7120e − 01

4
7.9347e − 02
6.9096e − 01

5
8.0830e − 02
6.8884e − 01

6
8.0919e − 02
6.8898e − 01

7
8.0918e − 02
6.8900e − 01

8
8.0918e − 02
6.8900e − 01

9
8.0918e − 02
6.8900e − 01

10
8.0918e − 02
6.8900e − 01

11
8.0919e − 02
6.8900e − 01

12
8.0918e − 02
6.8900e − 01

13
8.0919e − 02
6.8900e − 01

14
8.0919e − 02
6.8900e − 01

15
8.0919e − 02
6.8899e − 01

Figure 8: Absolute errors of the two quantiles and of the PCE approximation of GBM at time , whose parameters are , , and starting value .

7. Vasicek Interest Rate Model

In what follows, we consider the Vasicek model; see [17]. We also refer to [30, Section 4.4.3], for further details, which is defined by

where , , and are positive, real-valued constants. The analytical solution to the SDE (35), at time , is

Therefore, is normally distributed and its mean and variance are given by

Note that we have used the notation in order to highlight that the Vasicek model (35) is used to describe the evolution of interest rates. We also underline that such a model belongs to the class of the so-called one factor short rate model which includes, for example, the Rendleman-Bartter model, the CIR model, the Hull-White model, and the Ho-Lee model, because depends only on one market risk factor; see, for example, [31, 32] and [19, Section 3.2.1].

It is worth mentioning that the Vasicek model solution allows for negative values for the interest rate evaluated at given time ; hence, its use in credit market framework is often criticized by practitioners.

7.1. PCE Approximation

Analogous to what we have done concerning the study of the GBM model, Section 6, we can rewrite the solution to (35) using a functional such that

for every fixed time , which is often referred to as maturity time. As in Section 6, we shorten the notation omitting the dependence on vector of parameters , and we consider the basic random variable to be Gaussian variable; namely, we take in order to respect the restrictions imposed by the software we have used to numerically implement the PCE method.

Applying the NISP procedure up to degree , the PCE decomposition of is given by

where the coefficients are detecting, using (14), and evaluating on a set of realizations of the basic random variable , namely, considering values . Moreover, in what follows, we set . The obtained simulated values for are indicated as follows:

7.2. Numerical Application

Let us implement the theoretical approach developed in previous section by a working example, considering the Vasicek model parameters as we have reported in Table 5.

Table 5: Parameters and initial value for the Vasicek interest rate model.

Parameters

Values
0.1
0.2
15%
110
1

Figure 9 displays the absolute error as well as the relative error of the variance, while Table 6 gathers the absolute error and the relative error of the mean and variance of , respectively, with respect to degree . In particular, , the absolute errors for the mean and for the variance are as follows:

while the relative errors are given by

Due to the logarithmic scale, we have used the semilogy scale to obtain the graph in Figure 9, and we consider their absolute values, highlighting the related order, rather than their numerical values.

Table 6: Absolute error of the average and the variance of at time , with respect to Vasicek parameters , , and .

Degree of PCE
Average error
Variance error
Average relative error
Variance relative error

0
2.8015e − 10
1.8544e − 02
3.1076e − 12
1.0000e + 00

1
2.8049e − 10
1.2490e − 03
3.1114e − 12
6.7349e − 02

2
2.8028e − 10
1.2490e − 03
3.1090e − 12
6.7349e − 02

3
2.8048e − 10
1.2490e − 03
3.1112e − 12
6.7349e − 02

4
2.8034e − 10
1.2490e − 03
3.1096e − 12
6.7349e − 02

5
2.8041e − 10
1.2490e − 03
3.1104e − 12
6.7349e − 02

6
2.8039e − 10
1.2490e − 03
3.1103e − 12
6.7349e − 02

7
2.8052e − 10
1.2490e − 03
3.1117e − 12
6.7349e − 02

8
2.8037e − 10
1.2490e − 03
3.1100e − 12
6.7349e − 02

9
2.8048e − 10
1.2490e − 03
3.1112e − 12
6.7349e − 02

10
2.8038e − 10
1.2490e − 03
3.1101e − 12
6.7349e − 02

11
2.8039e − 10
1.2490e − 03
3.1103e − 12
6.7349e − 02

12
2.8076e − 10
1.2490e − 03
3.1144e − 12
6.7349e − 02

13
2.8061e − 10
1.2490e − 03
3.1126e − 12
6.7349e − 02

14
2.8055e − 10
1.2490e − 03
3.1120e − 12
6.7349e − 02

15
2.8062e − 10
1.2490e − 03
3.1128e − 12
6.7349e − 02

(a)

(b)

Figure 9: Semilogy scale plot of the absolute error of the variance of evaluated at time (a) and the absolute relative error of the variance (b) of at time . In both cases, the parameters of the Vasicek interest rate model are , , and .

The absolute and relative errors of mean and variance are stationary; see Table 6; thus, there is an error that inhibits the spectral convergence of the PCE. The latter could be caused by the approximation of (35) with (38); indeed, the numerical approximations made during computations of produce error of order .

Let us first analyze the error of the variance by considering

for a fixed degree , and for a decreasing set of -s approaching zero. In particular, we consider , and then we obtain the results displayed in Figure 10, which points out how the error of the variance is influenced by the parameters; this agrees with the achieved value in previous computation (Table 6).

Figure 10: Log-Log scale plot of the variance absolute error of at for different compared with order 1. The other parameters of the Vasicek interest model are , , and initial value .

Coming back to the values in Table 5, let us compute a sampling of size of , with , whose histogram, Figure 11, is compared with the probability density function of the normal random variable . Moreover in Figure 11 we show a sampling obtained exploiting the standard Monte Carlo technique.

(a)

(b)

Figure 11: (a) Probability density function of Vasicek interest rate model, with parameters , , , and (blue curve) and histogram of a sampling (size = ) of the for . With the same parameters, in (b), we compare the analytical probability density function of the Vasicek model, with the one obtained using the Monte Carlo sampling of size , for , being .

Let us compute the two quantiles , where and , of the PCE approximation , where . Proceeding as in Section 6, we consider the th sample quantile, namely, the th realization of the sampling of such that , where is the size of the sampling and denotes the integer part of the real number within the brackets. In particular, we have , while , and we compute the absolute error of versus the analytical values . Therefore, we evaluate , where is the quantile of a normal random variable of mean and variance and, indicating its cumulative density function by , we have . The aforementioned errors for and , respectively, are shown in Figure 12 and in Table 7.

Table 7: Absolute errors of the two quantiles and of the PCE approximation of Vasicek interest rate model at time . The parameters are set as , , , and .

Degree of PCE

0
4.3141e − 02
5.7307e − 02

1
2.6231e − 01
3.4536e − 01

2
2.6231e − 01
3.4536e − 01

3
2.6231e − 01
3.4541e − 01

4
2.6231e − 01
3.4536e − 01

5
2.6231e − 01
3.4536e − 01

6
2.6231e − 01
3.4536e − 01

7
2.6231e − 01
3.4536e − 01

8
2.6231e − 01
3.4536e − 01

9
2.6231e − 01
3.4536e − 01

10
2.6231e − 01
3.4536e − 01

11
2.6231e − 01
3.4536e − 01

12
2.6231e − 01
3.4536e − 01

13
2.6231e − 01
3.4536e − 01

14
2.6231e − 01
3.4536e − 01

15
2.6231e − 01
3.4537e − 01

Figure 12: Absolute errors of the two quantiles and of the PCE approximation of Vasicek interest rate model at time . The parameters are set as , , , and .

8. CIR Interest Rate Model

In what follows, we consider the PCE method for the CIR interest rate model, introduced in 1985 by Cox et al.; see [18]. Consider

where , , and are positive and real-valued constants and is a standard Wiener process. Let us recall that (44) generalizes the Vasicek model seen in Section 7; in fact, even if the drift term is common between the two, the standard deviation factor is different and, in particular, it equals for the CIR model ensuring the fact of avoiding negative interest rates if and are positive. Moreover, if , we cannot have to be equal to zero. It is worth mentioning that, unlike the case of the Vasicek model, there is not an analytical solution to (44); nevertheless, the CIR process is ergodic and possesses a stationary distribution; see [19, Section 3.2.3], [16, Section 23.5], and [33] for further details.

The statistics of the CIR model for fixed are

8.1. PCE Approximation

In order to apply the PCE method to the process in (44) evaluated at a specific time , we consider the stochastic process

and then applying the Itō-Döebling lemma, satisfies the following SDE:

Therefore, we have canceled the interactions between the state of the original process and the increments of the Wiener process. The transformation defined by (46), for example, [24, Chapter 2] and [25, Chapter 4] for details, allows us to reduce simulations instability.

As seen from studying the GBM model, Section 6 and the Vasicek model and Section 7, exploiting the fact that the Wiener process in (46) is normally distributed, namely, , we can rewrite the solution of (44), through (46), using a suitable functional , in order to obtain

for every fixed end time . Moreover, proceeding as in previous Sections 6 and 7, we set the basic variable and we shorten the notation omitting the explicit dependence of the CIR model on the parameters vector , with , , and being positive constants. Applying the truncated, at degree , PCE to , we obtain

where the coefficients are detecting, using (14), and evaluating on a set of realizations of the basic random variable; namely, we compute at , taking and referring to these values as

8.2. Numerical Application

In order to develop a concrete example for the application of the PCE to the CIR model, we consider the parameters values shown in Table 8.

Table 8: Parameters involved in CIR interest model.

Parameters

Values
0.1
0.2
15%
110
1

Figure 13 and Table 9 show, respectively, the absolute errors of the mean of the variance of the following statistics:

while Figure 14 and Table 9 display the absolute value of the relative error, where, as we have assumed discussing the GBM and the Vasicek model, in Figure 4, in which we use a logarithmic scale, we report their absolute values, to highlight its order, rather than their numerical values. While the relative errors are given by

and the data are computed for an increasing set of degree .

Table 9: Error of the average and variance of the PCE approximation of CIR model at time , with parameters , , , and .

Degree of PCE
Average error
Variance error
Average relative error
Variance relative error

0
9.2721e − 03
2.7378e + 00
1.2547e − 04
1.0000

1
6.2749e − 04
1.8550e − 01
8.4911e − 06
0.0678

2
6.2739e − 04
1.8535e − 01
8.4898e − 06
0.0677

3
6.2748e − 04
1.8535e − 01
8.4909e − 06
0.0677

4
6.2741e − 04
1.8535e − 01
8.4900e − 06
0.0677

5
6.2748e − 04
1.8535e − 01
8.4909e − 06
0.0677

6
6.2742e − 04
1.8535e − 01
8.4901e − 06
0.0677

7
6.2748e − 04
1.8535e − 01
8.4909e − 06
0.0677

8
6.2743e − 04
1.8535e − 01
8.4902e − 06
0.0677

9
6.2748e − 04
1.8535e − 01
8.4909e − 06
0.0677

10
6.2744e − 04
1.8535e − 01
8.4903e − 06
0.0677

11
6.2748e − 04
1.8535e − 01
8.4909e − 06
0.0677

12
6.2744e − 04
1.8535e − 01
8.4904e − 06
0.0677

13
6.2749e − 04
1.8535e − 01
8.4910e − 06
0.0677

14
6.2744e − 04
1.8535e − 01
8.4904e − 06
0.0677

15
6.2749e − 04
1.8535e − 01
8.4910e − 06
0.0677

(a)

(b)

Figure 13: Semilogy scale plot of the absolute error of mean (a) and variance (b) computed via PCE of CIR model at time , with parameters , , , and .

(a)

(b)

Figure 14: Semilogy scale plot of the relative errors (a) and (b) computed for , for the CIR interest rate model at time , with parameters , , , and .

The errors in Figure 13 and in Table 9 are stationary; hence, there is not a spectral convergence when the degree of truncation increases, as we have already observed for the GBM model and for the Vasicek model. The latter can be implied by the model approximation we have chosen, namely, replacing (44) with (48), and/or by numerical approximations made during computations of , if we increase the precision in computing , and using these values to compute again the PCE for the same set of degree.

The related new values obtained for the mean, the variance, and error are shown in Table 10 and Figures 15 and 16. Although more precise values are computed, the stationary behavior still holds. Therefore, the lack of spectral convergence is due to having replaced (44) by means of (48).

Table 10: Error of the average and variance of the PCE approximation of the CIR interest rate model with higher precision in computing and parameters , , and .

Degree of PCE
Average error
Variance error
Average relative error
Variance relative error

0
9.2722e − 03
2.7378e + 00
1.2547e − 04
1.0000

1
6.2748e − 04
1.8550e − 01
8.4909e − 06
0.0678

2
6.2748e − 04
1.8535e − 01
8.4909e − 06
0.0677

3
6.2748e − 04
1.8535e − 01
8.4909e − 06
0.0677

4
6.2748e − 04
1.8535e − 01
8.4909e − 06
0.0677

5
6.2748e − 04
1.8535e − 01
8.4909e − 06
0.0677

6
6.2748e − 04
1.8535e − 01
8.4909e − 06
0.0677

7
6.2748e − 04
1.8535e − 01
8.4909e − 06
0.0677

8
6.2748e − 04
1.8535e − 01
8.4909e − 06
0.0677

9
6.2748e − 04
1.8535e − 01
8.4909e − 06
0.0677

10
6.2748e − 04
1.8535e − 01
8.4909e − 06
0.0677

11
6.2748e − 04
1.8535e − 01
8.4909e − 06
0.0677

12
6.2748e − 04
1.8535e − 01
8.4909e − 06
0.0677

13
6.2748e − 04
1.8535e − 01
8.4909e − 06
0.0677

14
6.2748e − 04
1.8535e − 01
8.4909e − 06
0.0677

15
6.2748e − 04
1.8535e − 01
8.4909e − 06
0.0677

(a)

(b)

Figure 15: Semilogy scale plot of the absolute error of mean and variance computed via PCE of the CIR interest rate model at time with higher precision at computing . Moreover, the parameters , , , and .

(a)

(b)

Figure 16: Semilogy scale plot of the relative errors (a) and (b) computed for with higher precision at computing , for . Moreover, the CIR interest rate parameters are , , , and .

Even if CIR model does not admit close solution, for example, [30, Section 4.4.4], it possesses a stationary distribution. Namely, by following [19, Section 3.2.3], features a noncentral chi-squared distribution. In particular, let us denote by the analytical probability density function of , for a positive end time , and then

where the subscript denotes the random variable considered. Moreover, the constants are defined as

The last two constants are the degree of freedom and noncentrality parameter, respectively. For further features, see [19, Section 3.2.3]. In order to achieve a natural number , let us set , where . The staring value is still , while the ending time is set to be .

Therefore, let us employ a LHS sampling of size of the PCE and a standard Monte Carlo sampling of the same size of , exploiting its noncentered chi-squared distribution. The results are displayed in Figure 17 and compared with the analytical probability density function .

(a)

(b)

Figure 17: (a) Probability density function of CIR interest rate model, whose parameters are , , , and (blue curve) and histogram of a sampling (size = ) of the for . (b) compares the analytical probability density function of CIR interest rate model, for the same parameters, with a standard Monte Carlo sampling of size of at time .

Let us compute the two quantiles , where and , of the PCE approximation , where . Proceeding as in Section 6, we consider the th sample quantile, namely, the th realization of the sampling of such that , where is the size of the sampling and denotes the integer part of the real number within the brackets. In particular, we have , while , and we compute the absolute error of versus the analytical values . Therefore, we evaluate , where is the quantile of a normal random variable of mean and variance and, indicating its cumulative density function by , we have . The aforementioned errors for and , are, respectively, shown in Figure 18 and Table 11.

Table 11: Absolute errors of the two quantiles and of the PCE approximation of CIR Interest rate model at time , , , , and .

Degree of PCE

0
3.9106e + 00
5.3240e + 00

1
1.9247e − 01
4.2524e − 01

2
1.5447e − 01
3.5282e − 01

3
1.5447e − 01
3.5282e − 01

4
1.5447e − 01
3.5282e − 01

5
1.5447e − 01
3.5282e − 01

6
1.5447e − 01
3.5282e − 01

7
1.5447e − 01
3.5282e − 01

8
1.5447e − 01
3.5282e − 01

9
1.5447e − 01
3.5282e − 01

10
1.5447e − 01
3.5282e − 01

11
1.5447e − 01
3.5282e − 01

12
1.5447e − 01
3.5282e − 01

13
1.5447e − 01
3.5282e − 01

14
1.5447e − 01
3.5282e − 01

15
1.5447e − 01
3.5282e − 01

Figure 18: Absolute errors of the two quantiles and of the PCE approximation of CIR interest rate model at time , , , , and .

Remark 3. In each of the above considered cases, the lack of spectral convergence arises. Nevertheless, such a point is not connected with the model nor with the related parameters; in fact, it is a general consequence due to the approximations implied by the concrete application of the PCE machinery used to compute the solution at a certain positive time .

9. Polynomial Chaos Expansion Compared with MC and QMC

In this section. we compare the PCE method with the standard Monte Carlo (MC) method and with the quasi-Monte Carlo (QMC) one. The accuracy in approximating the statistics of each model studied, namely, GBM, the Vasicek, and the CIR models, evaluated at time , and the computational time costs to detect their mean and variance are considered as criteria for comparison.

It is worth mentioning that both the MC and the QMC methods are usually applied to approximate analytical solutions; hence, from our interest models point of view, we can consider them only for the GBM and Vasicek cases, while, for the CIR model, we consider them in relation to the Euler-Maruyama scheme.

The usual Monte Carlo approximations require a set of independent realizations of the process . Then, the mean and variance are determined as

and then, exploiting the Law of large numbers, we have that both (55) and (56) converge to the correspondent true value, and hence their standard errors,

are more informative than the single execution error. The QMC method uses low-discrepancy sequence generated by Sobol algorithm that maximizes the uniformity of the sample points; see, for example, [15, Chapter 5]. Then, it uses such values to approximate the statistics we are interested in, or, generally speaking, to compute integrals of the following type:

where is the measure induced by the considered random variable which characterizes the related image space. By setting and , respectively, we can approximate the mean of the random variable and its variance, respectively. Moreover, in QMC computations, the parameters for accuracy are given by the absolute errors:

We underline that, for both MC and QMC, we consider only the effective time to compute the statistics, without taking into account the time consumed for the detection of the grid of samplings which is necessary to simulate the process. Moreover, we highlight that, in each of the cases we consider, namely, the GBM, the Vasicek, and the CIR models, all the computations concerning the PCE method have followed the schemes developed in Sections 6, 7, and 8, respectively, with accuracy estimated by means of the absolute error of the average and of the variance.

9.1. Geometric Brownian Motion

We consider the SDE in (11), where the parameters are taken as follows: , , , and . The MC and QMC methods approximate the statistics of the solution to the GBM model exploiting its analytical solution; see (28), for an increasing size of samplings , and we then compute, for each value of , the standard errors for the MC, respectively, and the absolute error for the QMC is computed. In order to compare the accuracy and computational time costs of the three aforementioned methods, we consider the plot in Figure 19 where the -axis represents the computational time costs of the methods, while the error of the statistics, namely, , , , , , and , are displayed along the -axis. Moreover, the -axis is normalized with respect to the highest data among MC, QMC, and PCE values. Actually on the -axis, we report related relative computational time costs shown, which are more effective than the absolute ones. We also note that each point of the plot represents the computations for different number of realization points, namely, for MC and QMC, and for the PCE method.

(a)

(b)

Figure 19: Semilogy scale plot which compares the accuracy and computational time costs for detecting the average (a) and variance (b) of the GBM by means of PCE, MC, and QMC approaches. The parameters of the GBM equity model are , , , and .

Figure 19 points out the high accuracy as well as the very low computational effort required by PCE to get the solution.

Let us display the values of computational time costs for the PCE; see Table 12. For MC and QMC methods, respectively, see Table 13. For the errors made using both MC and QMC methods, see Table 14.

Table 12: Time costs of computations for PCE, approximating GBM at for , , and .

Degree of PCE
Computational time costs

1
0.0170

2
0.0210

3
0.0340

4
0.0370

5
0.0530

6
0.0560

7
0.0710

8
0.0740

9
0.0930

10
0.0990

11
0.1090

12
0.1210

13
0.1370

14
0.1390

15
0.1480

Table 13: Computational time costs achieved by MC and QMC methods to detect the mean and variance of the process , at , with parameters , , and .

Size of sampling
Average
Variance

MC
QMC
MC
QMC

256
0.0020
0.0030
0.0020
0.0020

512
0.0020
0.0020
0.0040
0.0040

1024
0.0050
0.0050
0.0050
0.0060

2048
0.0090
0.0100
0.0110
0.0120

4096
0.0170
0.0190
0.0220
0.0200

8192
0.0350
0.0370
0.0430
0.0410

16384
0.0690
0.0740
0.0810
0.0830

32768
0.1410
0.2070
0.2040
0.1940

65536
0.3020
0.2820
0.3430
0.3310

Table 14: Errors of average and variance computed by MC and QMC methods for the GBM evaluated at time , with parameters , , and .

Size of sampling
Average
Variance

256
9.3149e − 01
1.7295e − 01
1.9672e + 01
6.5079e + 00

512
6.7405e − 01
9.4193e − 02
1.4553e + 01
3.6819e + 00

1024
4.8747e − 01
5.0744e − 02
1.0759e + 01
2.0760e + 00

2048
3.3705e − 01
2.7110e − 02
7.2722e + 00
1.1656e + 00

4096
2.3915e − 01
1.4388e − 02
5.1771e + 00
6.5155e − 01

8192
1.7213e − 01
7.5954e − 03
3.7926e + 00
3.6254e − 01

16384
1.2203e − 01
3.9921e − 03
2.6956e + 00
2.0082e − 01

32768
8.6400e − 02
2.0905e − 03
1.9111e + 00
1.1076e − 01

65536
6.0926e − 02
1.0913e − 03
1.3439e + 00
6.0845e − 02

For the computational time costs required for the detection of PCE’s coefficients, see (14), and its statistics, and see (10) and (11), which are irrelevant if compared with time required to evaluate the process, (26), at , which are defined by Gaussian quadrature formula as stated in (14).

This motivates the data in Table 12.

9.2. Vasicek Interest Rate Model

Proceeding as in Section 9.1, we consider the SDE (35) characterizing the Vasicek model for the following set of values for its parameters , , and , and we apply both the MC and QMC methods to approximate the analytical solution (36) to (35) and to compute the statistics of at time for an increasing size of samplings We use the same type of plot exploited in the previous section, in order to show the accuracy and the computational time costs and the errors of the PCE, the MC, the QMC method, respectively.

In particular, Figure 20 points out the high accuracy of PCE and the corresponding low computational effort to get these results: the key points are the low number of simulations required to get the solution of (35).

(a)

(b)

Figure 20: Semilogy scale plot of the error of mean (a) and variance (b) computed via MC, QMC, and PCE for the Vasicek model at time , with , , , and and . The displayed errors for the mean are , , and absolute error for PCE, while, for the variance, we have , , and absolute error for PCE.

In order to underline the efficiency and accuracy of the PCE method, let us show the related computational time costs, Table 15, compared to those of the MC and QMC methods, Table 16. We also highlight the numerical values concerning the absolute errors achieved using PCE method, Table 15, as well as MC and QMC approaches, Table 17.

Table 15: Time computational costs and absolute error of mean and variance of PCE, concerning Vasicek model at time . The parameters of the model are set as , , , and .

Degree of PCE
Simulation cost
Absolute mean error
Absolute variance error

0
0.0030
3.3322e − 10
2.2500e − 02

1
0.0420
3.4262e − 10
1.5000e − 07

2
0.0480
3.3637e − 10
1.5000e − 07

3
0.0840
3.4285e − 10
1.5000e − 07

4
0.0920
3.3759e − 10
1.5000e − 07

5
0.1340
3.4269e − 10
1.5000e − 07

6
0.1500
3.3836e − 10
1.5000e − 07

7
0.1870
3.4281e − 10
1.5000e − 07

8
0.1930
3.3869e − 10
1.5000e − 07

9
0.2280
3.4279e − 10
1.5000e − 07

10
0.2400
3.3906e − 10
1.5000e − 07

11
0.2810
3.4265e − 10
1.5000e − 07

12
0.2920
3.4008e − 10
1.5000e − 07

13
0.3360
3.4281e − 10
1.5000e − 07

14
0.3480
3.3984e − 10
1.5000e − 07

15
0.3900
3.4279e − 10
1.5000e − 07

Table 16: Computational time costs for determining the average and variance of the Vasicek interest rate , at time , achieved by MC and QMC approaches, with parameters , , , and .

Size of sampling
Average
Variance

MC
QMC
MC
QMC

1024
0.0030
0.0030
0.0040
0.0020

2048
0.0050
0.0040
0.0080
0.0050

4096
0.0070
0.0080
0.0090
0.0090

8192
0.0150
0.0150
0.0180
0.0190

16384
0.0290
0.0280
0.0370
0.0350

32768
0.0570
0.0570
0.0680
0.0710

65536
0.1090
0.1140
0.1350
0.1420

131072
0.2170
0.2220
0.3150
0.2780

262144
0.4360
0.4520
0.5360
0.5530

Table 17: Errors of the average and variance computed by MC and QMC for the Vasicek interest rate model at time , with parameters , , , and .

Size of sampling
Average
Variance

1024
4.7038e − 03
4.3575e − 04
1.0018e − 03
9.1509e − 05

2048
3.2845e − 03
2.3300e − 04
6.9060e − 04
4.6037e − 05

4096
2.3248e − 03
1.2366e − 04
4.8926e − 04
2.3142e − 05

8192
1.6647e − 03
6.5240e − 05
3.5473e − 04
1.1625e − 05

16384
1.1666e − 03
3.4249e − 05
2.4635e − 04
5.8363e − 06

32768
8.3019e − 04
1.7906e − 05
1.7644e − 04
2.9286e − 06

65536
5.8398e − 04
9.3287e − 06
1.2347e − 04
1.4689e − 06

131072
4.1560e − 04
4.8455e − 06
8.8436e − 05
7.3654e − 07

262144
2.9258e − 04
2.5103e − 06
6.1982e − 05
3.6919e − 07

As observed in the GBM case, the computational time costs required for the detection of PCE’s coefficients, (14), and its statistics, (10) and (11), are irrelevant if compared with time required to get the evaluation of the process in (38) at quadrature nodes , as we can see looking at data displayed in Table 15.

9.3. CIR Interest Rate Model

Even if the CIR model does not have a solution in closed form, we can still apply both the MC and the QMC methods, considering the numerical solution of the SDE (44).

As described in [24], the Euler-Maruyama scheme applied to the transformed process (47) is actually the Milstein scheme applied to the original SDE (46). Hence, let us integrate numerically the transformed SDE from up to using time steps. The mean and the variance are estimated, via MC and QMC, using a set of samplings of increasing size , setting the CIR model parameters as follows: , , and . We note that, while the MC method applies straightforward, the QMC implementation requires some further discussion which is, in fact, close to the analysis made in Section 5. Namely, the numerical scheme allows us to write the solution at as a functional of all the independent increments of the Wiener process, one for each time step.

As the dimensionality of the problem equals the number of time steps, that is, , then a multidimensional low-discrepancy sequence is computed, where the size represents the number of points belonging to this sequence. Once the the SDE is solved, the Euler-Maruyama scheme prescribes that each independent increment has to be replaced by an entry of a point belonging to such a multidimensional low-discrepancy sequence. Eventually, by applying (58) to the detected solution, we get the required statistics.

As we made in Sections 9.1 and 9.2, in Figure 21, we compare the computational time costs and its related accuracy reached by PCE, MC, and QMC methods, respectively. The data used to apply the PCE technique are computed following the approach developed in Section 8. We recall that the plot is characterized by setting along the -axis the computational time costs of the methods, while the errors of the statistics , , , and and the PCE absolute errors are placed along the -axis. Moreover, the time costs are normalized with respect to the highest time value among MC, QMC, and PCE; thus, we display their relative computational time costs, hence achieving more information than merely using the related absolute values.

(a)

(b)

Figure 21: Semilogy scale plot which compare the accuracy and computational time costs for detecting the average (a) and variance (b) of the CIR model at by means of PCE, MC, and QMC approaches, with parameters , , and .

Let us underline that Figure 21 points out the high accuracy as well as the very low computational effort required by the PCE method to get the solution. We can appreciate such an effectiveness property of the PCE approach by comparing the computational effort required by the PCE, Table 18, and by the MC and QMC method, Table 19. Moreover, the PCE errors suggest better accuracy for the latter, particularly if compared with MC and QMC; see Tables 19 and 20.

Table 18: Computational time costs and absolute error of mean and variance for PCE approximation of CIR interest rate model evaluated at , whose parameters are , , , and .

Degree of PCE
Simulation cost
Absolute mean error
Absolute variance error

0
0.0040
5.6192e − 03
2.4687e + 00

1
0.0170
2.6869e − 06
1.2408e − 03

2
0.0260
2.5466e − 06
1.1771e − 03

3
0.0340
2.6731e − 06
1.1776e − 03

4
0.0390
2.5744e − 06
1.1772e − 03

5
0.0570
2.6750e − 06
1.1775e − 03

6
0.0570
2.5896e − 06
1.1773e − 03

7
0.0700
2.6761e − 06
1.1775e − 03

8
0.0720
2.6001e − 06
1.1773e − 03

9
0.0920
2.6770e − 06
1.1775e − 03

10
0.0980
2.6101e − 06
1.1773e − 03

11
0.1140
2.6856e − 06
1.1775e − 03

12
0.1150
2.6160e − 06
1.1773e − 03

13
0.1300
2.6856e − 06
1.1775e − 03

14
0.1390
2.6194e − 06
1.1774e − 03

15
0.1890
2.6869e − 06
1.1775e − 03

Table 19: Computational cost of computation of statistics of CIR interest rate model at , determined by MC, QMC, with , , , and .

Size of sampling
Average
Variance

MC
QMC
MC
QMC

256
0.0070
0.0070
0.0060
0.0050

512
0.0080
0.0080
0.0080
0.0060

1024
0.0110
0.0110
0.0110
0.0100

2048
0.0190
0.0190
0.0180
0.0170

4096
0.0340
0.0370
0.0310
0.0310

8192
0.0680
0.0650
0.0600
0.0560

16384
0.1440
0.1240
0.1390
0.1110

32768
0.2700
0.2480
0.2400
0.2240

Table 20: Error of the average and variance MC, QMC approximations of CIR interest rate model at time for , , , and .

Size of sampling
Average
Variance

256
9.6603e − 02
1.1335e − 02
2.1158e − 01
1.7954e − 01

512
6.8815e − 02
1.1394e − 03
1.5168e − 01
1.3344e − 01

1024
5.0225e − 02
9.3599e − 04
1.1421e − 01
1.6659e − 01

2048
3.5422e − 02
2.2688e − 05
8.0321e − 02
8.1676e − 02

4096
2.4939e − 02
2.3517e − 04
5.6298e − 02
3.5725e − 02

8192
1.7372e − 02
3.4153e − 04
3.8630e − 02
4.4165e − 02

16384
1.2201e − 02
1.2216e − 04
2.6949e − 02
4.5768e − 03

32768
8.6404e − 03
6.2970e − 05
1.9112e − 02
2.2469e − 02

The computational time costs required for the detection of PCE’s coefficients, (14), and its statistics, (10) and (11), are irrelevant if compared with time required to get the evaluation of the process in (48) at the quadrature nodes . This motivates the time data in Table 18.

10. Conclusions

In this paper, we show how the Polynomial Chaos Expansion (PCE) method can be effectively used to approximate the solution of a given SDE, exploiting only some basic properties of the Brownian motion, which is assumed to drive the stochastic process we are interested in, and a suitable transformation of the original SDE we have taken into consideration. In particular, we have applied the PCE technique to approximate the solutions of some of the most relevant interest rate models, namely, the Geometric Brownian Motion (GBM) model, the Vasicek model, and the Cox-Ingersoll-Ross (CIR) model.

The provided numerical examples deeply discuss the convergence properties of the PCE method. Indeed, detailed analysis on error properties of the approximated statistics is given when analytical solution is available also providing a comparison with respect to the related analytical probability density function. The latter is not the case for the CIR model, conversely to what happens for both the GBM and the Vasicek models. Hence, in this latter case, the PCE approach is even more interesting also because it gives a general method to study those stochastic processes which are defined by SDEs which do not have solution in closed form.

We show by concrete examples how the PCE method overcomes the performances of the standard Monte Carlo (MC) method as well as those of the quasi-Monte Carlo (QMC) one. In particular, a close comparison between PCE, MC, and QMC techniques points out the advantage in using the PCE approach, especially in terms of computational costs. The latter result is due to the fact that the PCE method requires just few realizations of the random variable we want to approximate. In particular, even if the basic operations, required by the PCE, are more complicated than those implied by the MC or the QMC approach (weighted average versus quadrature formulas), the advantages in using the PCE technique are great and indubitable, as shown, for example, for the analysis of the CIR model, in Tables 18, 19, and 20.

Last but not least, it is worth mentioning that the PCE approach can be split in two distinct sections: the off-line computations and the real PCE approximation of the random variable we want to study. Thus, the basic random variable, the orthogonal polynomials, and some data used in the Gaussian quadrature formulas are computed once and for all, because they do not depend on the particular model, for example, the GBM, Vasicek, or the CIR models in our analysis. Hence, a clear computational effort saving is achieved when the model is sensible to time variation or if a recalibration of the involved parameters is required.

Appendix

Value of the Polynomials

In Table 21, the first six orthonormal Hermite polynomials are evaluated on a set of equally spaced points in the interval .

Table 21: Evaluation of first six orthonormal Hermite polynomials on the vector of equally spaced points from to .

−2.00000
1.00000
−2.82843
4.94975
4.94975
3.87836
0.25820

−1.89744
1.00000
−2.68338
4.38443
4.38443
2.37691
1.26342

−1.79487
1.00000
−2.53833
3.84887
3.84887
1.19518
1.83459

−1.69231
1.00000
−2.39328
3.34307
3.34307
0.29413
2.06903

−1.58974
1.00000
−2.24824
2.86701
2.86701
−0.36308
2.05174

−1.48718
1.00000
−2.10319
2.42071
2.42071
−0.81117
1.85611

−1.38462
1.00000
−1.95814
2.00417
2.00417
−1.08266
1.54464

−1.28205
1.00000
−1.81309
1.61737
1.61737
−1.20789
1.16962

−1.17949
1.00000
−1.66805
1.26033
1.26033
−1.21508
0.77387

−1.07692
1.00000
−1.52300
0.93305
0.93305
−1.13023
0.39138

−0.97436
1.00000
−1.37795
0.63551
0.63551
−0.97719
0.04808

−0.87179
1.00000
−1.23290
0.36773
0.36773
−0.77766
−0.23748

−0.76923
1.00000
−1.08786
0.12971
0.12971
−0.55115
−0.45346

−0.66667
1.00000
−0.94281
−0.07857
−0.07857
−0.31501
−0.59396

−0.56410
1.00000
−0.79776
−0.25709
−0.25709
−0.08441
−0.65840

−0.46154
1.00000
−0.65271
−0.40585
−0.40585
0.12764
−0.65073

−0.35897
1.00000
−0.50767
−0.52487
−0.52487
0.31028
−0.57879

−0.25641
1.00000
−0.36262
−0.61413
−0.61413
0.45486
−0.45358

−0.15385
1.00000
−0.21757
−0.67363
−0.67363
0.55485
−0.28856

−0.05128
1.00000
−0.07252
−0.70339
−0.70339
0.60594
−0.09896

0.05128
1.00000
0.07252
−0.70339
−0.70339
0.60594
0.09896

0.15385
1.00000
0.21757
−0.67363
−0.67363
0.55485
0.28856

0.25641
1.00000
0.36262
−0.61413
−0.61413
0.45486
0.45358

0.35897
1.00000
0.50767
−0.52487
−0.52487
0.31028
0.57879

0.46154
1.00000
0.65271
−0.40585
−0.40585
0.12764
0.65073

0.56410
1.00000
0.79776
−0.25709
−0.25709
−0.08441
0.65840

0.66667
1.00000
0.94281
−0.07857
−0.07857
−0.31501
0.59396

0.76923
1.00000
1.08786
0.12971
0.12971
−0.55115
0.45346

0.87179
1.00000
1.23290
0.36773
0.36773
−0.77766
0.23748

0.97436
1.00000
1.37795
0.63551
0.63551
−0.97719
−0.04808

1.07692
1.00000
1.52300
0.93305
0.93305
−1.13023
−0.39138

1.17949
1.00000
1.66805
1.26033
1.26033
−1.21508
−0.77387

1.28205
1.00000
1.81309
1.61737
1.61737
−1.20789
−1.16962

1.38462
1.00000
1.95814
2.00417
2.00417
−1.08266
−1.54464

1.48718
1.00000
2.10319
2.42071
2.42071
−0.81117
−1.85611

1.58974
1.00000
2.24824
2.86701
2.86701
−0.36308
−2.05174

1.69231
1.00000
2.39328
3.34307
3.34307
0.29413
−2.06903

1.79487
1.00000
2.53833
3.84887
3.84887
1.19518
−1.83459

1.89744
1.00000
2.68338
4.38443
4.38443
2.37691
−1.26342

2.00000
1.00000
2.82843
4.94975
4.94975
3.87836
−0.25820

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

P. E. Kloeden and E. Platen, Numerical Solution of Stochastic Differential Equations, Springer, New York, NY, USA, 1992.

S. Mahadevan, “Uncertainty quantification for decision-making in engineered systems,” in Proceedings of the International Symposium on Engineering under Uncertainty: Safety Assessment and Management (ISEUSAM-2012), pp. 97–117, Springer, Berlin, Germany, 2013.

J. E. Bickela and R. B. Bratvold, “From uncertainty quantification to decision making in the oil and gas industry,” Energy, Exploration & Exploitation, vol. 26, no. 5, pp. 311–325, 2008.

N. Wiener, “The homogeneous chaos,” American Journal of Mathematics, vol. 60, no. 4, pp. 897–936, 1938.

N. Wiener, “Generalized harmonic analysis,” Acta Mathematica, vol. 55, no. 1, pp. 117–258, 1930.

R. G. Ghanem and P. D. Spanos, Stochastic Finite Elements: a Spectral Approach, Springer, New York, NY, USA, 1991.

M. Villegas, F. Augustin, A. Gilg, A. Hmaidi, and U. Wever, “Application of the polynomial chaos expansion to the simulation of chemical reactors with uncertainties,” Mathematics and Computers in Simulation, vol. 82, no. 5, pp. 805–817, 2012.

J. Fisher and R. Bhattacharya, “Optimal trajectory generation with probabilistic system uncertainty using polynomial chaos,” Journal of Dynamic Systems, Measurement and Control, vol. 133, no. 1, Article ID 014501, 2011.

T. Crestaux, O. Le Maître, and J.-M. Martinez, “Polynomial chaos expansion for sensitivity analysis,” Reliability Engineering & System Safety, vol. 94, no. 7, pp. 1161–1172, 2009.

D. Xiu, Numerical Methods for Stochastic Computations, Princeton University Press, Princeton, NJ, USA, 2010.

O. P. Maître and O. M. Knio, Spectral Methods for Uncertainty Quantification, Springer, New York, NY, USA, 2010.

S. Oladyshkin and W. Nowak, “Data-driven uncertainty quantification using the arbitrary polynomial chaos expansion,” Reliability Engineering & System Safety, vol. 106, pp. 179–190, 2012.

G. Peccati and M. S. Taqqu, Wiener Chaos: Moments, Cumulants and Diagrams: A Survey with Computer Implementation, Springer, 2011.

F. G. Hickenell, “Koksma-Hlawka inequality,” in Encyclopedia of Statistical Sciences, John Wiley & Sons, 2006.

P. Glasserman, Monte Carlo Method in Financial Engineering, Springer, Berlin, Germany, 2003.

J. C. Hull, Options, Futures and Other Derivatives, Pearson/Prentice Hall, 2009.

O. Vasicek, “An equilibrium characterization of the term structure,” Journal of Financial Economics, vol. 5, no. 2, pp. 177–188, 1977.

J. C. Cox, J. E. Ingersoll, and S. A. Ross, “A theory of the term structure of interest rates,” Econometrica, vol. 53, no. 2, pp. 385–407, 1985.

D. Brigo and F. Mercurio, Interest Rate Models—Theory and Practice with Smile, Inflation and Credit, Springer, Berlin, Germany, 2006.

O. G. Ernst, A. Mugler, H.-J. Starkloff, and E. Ullmann, “On the convergence of generalized polynomial chaos expansions,” ESAIM: Mathematical Modelling and Numerical Analysis, vol. 46, no. 2, pp. 317–339, 2012.

O. Kallenberg, Foundations of Modern Probability, Springer, 1997.

M. Baudin, Scilab Enterprise, 2015, http://wiki.scilab.org/NISP%20Module

M. Baudin and J.-M. Martinez, NISP Toolbox Manual, Version 0.4, 2015, http://forge.scilab.org/index.php/p/nisp/downloads/get/nisp_toolbox_manual-v0.4.pdf.

S. M. Iacus, Simulation and Inference for Stochastic Differential Equation, vol. 1 of Springer Series in Statistics, Springer, New York, NY, USA, 2008.

P. Jackel, Monte Carlo Methods in Finance, John Wiley and Sons, Trowbridge, UK, 1999.

B. Tang, “Orthogonal array-based Latin hypercubes,” Journal of the American Statistical Association, vol. 88, no. 424, pp. 1392–1397, 1993.

A. B. Owen, “Orthogonal arrays for computer experiments, integration and visualization,” Statistica Sinica, vol. 2, no. 2, pp. 439–452, 1992.

R. L. Iman, “Latin hypercube sampling,” in Encyclopedia of Quantitative Risk Analysis and Assessment III, John Wiley & Sons, 2008.

H. A. David and H. N. Nagaraja, Order Statistics, Wiley Series in Probability and Statistics, John Wiley & Sons, New York, NY, USA, 3rd edition, 2003.

S. E. Shreve, Stochastic Calculus for Finance II, Springer, 2000.

R. Jarrow, “The term structure of interest rates,” Annual Review of Financial Economics, vol. 1, no. 1, pp. 69–96, 2009.

J. G. Andrew, “Cairns interest-rate models,” in Encyclopedia of Actuarial Science, John Wiley & Sons, 2004.

W. Sinkala, P. G. Leach, and J. G. O'Hara, “Zero-coupon bond prices in the Vasicek and CIR models: their computation as group-invariant solutions,” Mathematical Methods in the Applied Sciences, vol. 31, no. 6, pp. 665–678, 2008.