A method to compute the transition probability density associated to a multifactor CoxIngersollRoss model of the term structure of interest rates with no drift term^{ *}
Lorella Fatone, Graziella Pacelli, Maria Cristina Recchioni, Francesco Zirilli
Warning: If your browser does not allow you to see the content of this website click here to download a postscript file that reproduces the material contained here.
We consider a ndimensional squareroot process and we obtain a formula involving series expansions for the associated transition probability density. The process mentioned previously (with no drift term) can be used to model interest rates as proposed, in the one dimensional case by Cox, Ingersoll and Ross [1.1] and in the multifactor case by Chen and Scott [1.2] and Nowman [1.3]. The use of the squareroot process to model the assets prices takes into account some spurious effects, for example the smile effect, that escapes the more usual lognormal models. The transition probability density associated to this process is needed in the evaluation of options on assets whose prices are described by the square root process. The formula that we propose for the transition probability density has been obtained using appropriately a perturbative expansion in the correlation coefficients of the square root process, the Fourier transform and the method of characteristics to solve first order hyperbolic partial differential equations. The computational effort needed to evaluate this formula is polynomial with respect to the dimension n of the space spanned by the squareroot process when the order where all the series involved in the transition probability density formula are truncated is fixed. This strategy gives an accuracy that some numerical tests show approximately constant for a wide range of values of n. This fact makes the proposed formula in its region of validity of practical interest when large values of n are considered (i.e. n> 4 or 5). The coefficients of the series expansions involved in the formula for the transition probability density are linear combinations of products of explicitely known one variable functions. The zeroth order term of the perturbative expansion is a well known formula. The coefficients of the first order terms of the expansion are given via explicit formulae and the coefficients of the higher order terms can be computed through elementary and involved calculations following a procedure described in [1.4]. The formula for the transition probability density truncated to first order in the series in the correlation coefficients is used as a tool to evaluate quantitatively and efficiently the price of some derivatives on several assets. In particular some examples of financial derivatives whose evaluation involves integrals in two, twenty and one hundred dimensions (i.e. n=2,20,100), that is derivatives on two, twenty and one hundred assets, where accurate results can be obtained are shown in [1.4]. An experiment shows that the formula derived for the transition probability density is well suited for parallel computing. This website contains some interactive tools that help the understanding of the companion paper [1.4] (see Section 4) and a portable software library (see Section 5 and Section 6) that makes possible to the user to exploit the formula that has been derived to evaluate the transition probability densities of its own models and the prices of the associated financial derivatives.
A detailed exposition of the material summarized in this website can be found in [1.4].
The main purpose of this website consists in providing to the user some tools to price European options on several underlyings (assets, interest rates,...). These tools have been used to price options on up to one hundred underlyings. The prices of the underlyings are assumed to be described by a ndimensional square root process (with no drift term) (see Section 3 for further details). The dynamics of the prices Eqs (3.1), (3.2) of Section 3 is characterized by the following parameters: the volatilities s_{i}, i=1,2,...,n and the correlation coefficients r_{i,j}, i=1,2,...,n, j=1,2,...,n. We note that since we have r_{i,i}=1, r_{i,j}= r_{j,i}, i=1,2,...,n, j=1,2,...,n the correlation coefficients that must be specified are n(n1)/2, so that the total number of parameters to be specified to determine the price dynamics is n+n(n1)/2.
In this website we use three different notions of time: the "physical" time t, the maturity time T, the time to maturity t=T t³0.
We remind that the price of an European option on assets whose prices satisfy (3.1), (3.2) at (physical) time t depends on the prices of the underlyings y_{1},y_{2},..., y_{n}, that we denote with y, at physical time t, the maturity time T, (T> t), and the payoff function P(x,T) associated to the option considered where x denotes the prices of the underlyings x_{1},x_{2},...,x_{n} at the maturity time T. Let t be the time to maturity, t=Tt, when the transition probability density p(x,t;y) associated to the square root process, i.e. the probability that starting from the prices x at the maturity time T one reaches the prices y at physical time t, is known the price of the European option at physical time t when the prices of the underlyings at time t is y is given by the following integral:
 (2.1) 
This website gives three tools to the user to exploit formula (2.1).
The first tool (see Section 4) is a demo that helps the user in the understanding of the technical content of this website and of the paper [1.4]. That is in Section 4 the user finds a Java Applet that computes the price of some financial products. In particular options on the minimum or on the maximum of two underlyings are considered. The user can choose the payoff function, the volatilities s_{i}, i=1,2, and the correlation coefficient r_{1,2} in a specified set of choices placed at his disposal. The price C(y,t) computed with formula (2.1) is visualized as an animation, that is is visualized through a sequence of graphs, parametrized by the time to maturity t, whose elements represent the surface described by C(y,t) when t is fixed and y ranges on a rectangular grid contained in the nonnegative quadrant of the y_{1}, y_{2} plane. The time ranges from t=0, that corresponds to the maturity time T, to t=1, that corresponds to the physical time t =T1.
The second tool (see Section 5) is the portable software library, UnivPM_MathFinSolver1, that computes the transition probability density p(x,t;y) according to the formula explained in Section 3. The library UnivPM_MathFinSolver1 consists of three modules, that is UnivPM_MathFinSolver1_1.dll, UnivPM_MathFinSolver1_2.dll, UnivPM_MathFinSolver1_3.dll. The modules of this library take as input parameters the number of the space variables n, the volatilities s_{i}, i=1,2,...,n , the correlation coefficients r_{i,j}, i=1,2,...,n1, j=i+1,i+2,...,n, the time to maturity t, the price vectors x, y, at the maturity time T and at the physical time t respectively, and return quantities that can be used to approximate the corresponding transition probability density p(x,t;y) associated with the ndimensional square root process defined by equations (3.1), (3.2) of Section 3.
The input parameters mentioned above must be chosen taking into account the following restrictions:
1<n<102
0.005<s_{i} <0.6 , i=1,2,...,n,
0.6<r_{i,j} <0.6 , i=1,2,...,n1, j=i+1,i+2,...,n,
x, y_{}Î [0.005,10]^{n},
0<t<1
Note that the choice n=1 is not allowed. In fact in this case there are no correlation coefficients and the procedure at the basis of our formula for the transition probability density, that is the power series expansion in the correlation coefficients has no sense. Moreover the case n=1 is a trivial case where an explicit formula is known for the transition probability density.
The three routines with exstension .dll (dynamic link library) of the library UnivPM_MathFinSolver1 can be called by a Fortran main code or by a C main code prepared by the user following some instructions that may be found in the file UnivPM_MathFinSolver1_readme.txt.
The third tool (see Section 6) consists in two main codes one implemented in C language and the other in Fortran language. These codes contain a Monte Carlo routine able to compute the integral appearing in (2.1). The Monte Carlo routine to evaluate the integral (2.1) links the routines UnivPM_MathFinSolver1_2.dll and UnivPM_MathFinSolver1_3.dll described in Section 5 to compute an appropriately rescaled version of the transition probability density p(x,t;x_{0}). A readme file Monte_Carlo_readme.txt explains how the user can insert his own data in these main codes. These data consist of the number of the space variables n, the volatilities s_{i}, i=1,2,...,n , the correlation coefficients r_{i,j}, i=1,2,...,n1, j=i+1,i+2,...,n, the time to maturity t, the vector x_{0} of the prices of the underlyings at the time to maturity t, the payoff function P(x,T), xÎ R^{n}_{+}. After being executed the C and Fortran main codes return the price of the desired option according to formula (2.1).
The reader not interested in mathematical details may jump to Section 4.
Let R be the real numbers, R^{n}, n ³ 1, be the ndimensional real Euclidean space and let x = (x_{1},x_{2},¼,x_{n})^{T} Î R^{n} be a generic vector where the superscript T means transposed.
Let t be a real variable that denotes the time ranging from 0 to +¥. Let X = (X_{1},X_{2},¼,X_{n})^{T} Î R^{n}, we consider the following stochastic differential equations:
 (3.1) 
 (3.2) 
 (3.3) 
Note that the stochastic process X(t)= (X_{1}(t), X_{2}(t),...,X_{n}(t))^{T} solution of (3.1), (3.2) belongs to R^{n}_{+} for any t > 0.
We denote with p(x,t;x_{0}), x Î R^{n}_{+}, t > 0, x_{0} Î R^{n}_{+} the transition probability density function associated to the stochastic process X(t) = (X_{1}(t),X_{2}(t),¼,X_{n}(t))^{T}, t > 0, solution of (3.1), (3.2). The function p(x,t;x_{0}), x Î R^{n}_{+}, t > 0 , x_{0}Î R^{n}_{+}, satisfies the FokkerPlanck equation:
 (3.4) 
 (3.5) 

 (3.7) 
 (3.8) 

x > 0, y > 0, t > 0, s > 0,  (3.9) 




 (3.14) 
We note that the coefficients of the power series in the time t in formula (3.8) depend on the time t. This dependence is a consequence of the approach used to derive the expansion of the transition probability density (see [1.4]). Moreover this coefficients are trascendentally small with respect to the powers of t when t goes to 0^{+}.
Let us remind what has already been written in Section 2 (formula (2.1)), that is let us consider a (European) financial derivative product on n assets whose prices satisfy the stochastic differential equation (3.1), (3.2) characterized by a maturity time T and a payoff function P(x,T). The price C(y,t) of the derivative at time t to maturity (the current time is t=Tt) when the prices of the assets at the current time t=Tt is y= ( y_{1},y_{2},...,y_{n})^{T} is given by:
 (3.15) 
Note that the software library provided in Section 5 computes the transition probability density p approximating it with the first two terms appearing on the right hand side of (3.6), and that the main codes of Section 6, using the library of Section 5, compute the price C(y,t) of a derivative chosen by the user computing the integral (3.15) (or equivalently (2.1)). The derivative is chosen by the user through the choice of the payoff function P(x,T), x Î R_{+}^{n}.
We propose a Java Applet that computes the price of some financial products. In particular options on the minimum or on the maximum of two underlyings are considered, that is options whose payoff function is given by:
P(x,T)=max(min(x_{1},x_{2})K,0) (4.1)
and by:P(x,T)=max(max(x_{1},x_{2})K,0) (4.2)
respectively. The quantity K is the strike price (K>0). In the Applet we choose K=1. The user can choose the payoff function among the two possibilities mentioned above, the volatilities s_{i}, i=1,2, among the two choices: s_{1}=0.4, s_{2}=0.5 or s_{1}=0.5, s_{2}=0.3 and the correlation coefficient r_{1,2} among three choices that is r_{1,2}=0.5, 0, 0.5. The price C(y,t) computed using formula (2.1) is visualized through an animation, that is through a sequence of graphs, parametrized by the time to maturity t, whose elements represent the surface described by C(y,t) when t is fixed and y=(y_{1},y_{2})^{T} ranges on a rectangular grid contained in the nonnegative quadrant of the y_{1}, y_{2} plane. In the sequence of graphs visualized the time ranges from t=0, that corresponds to the maturity time T, to t=1, that corresponds to the physical time t =T1.
The software library consists of three modules with extension .dll, i.e. UnivPM_MathFinSolver1_1.dll, UnivPM_MathFinSolver1_2.dll, UnivPM_MathFinSolver1_3.dll, that the user can download and run. We remind that the extension .dll means dynamic link library. The manual, i.e. UnivPM_MathFinSolver1_readme.txt helps the understanding of the procedures contained in the software library.
UnivPM_MathFinSolver1_1.dll computes the first order approximation (i.e.: the zeroth order approximation plus the first order corrections) of p(x,t;x_{0}) contained in (3.6)
UnivPM_MathFinSolver1_2.dll computes the zeroth order approximation of p(x,t;x_{0}) contained in (3.6) appropriately rescaled, that is computes the function p^{~}_{0}( x,t; x_{0}) defined by:
 (5.1) 
UnivPM_MathFinSolver1_3.dll for a given couple of indices i, j, i = 1,2,¼,n1, j = i+1,i+2,¼,n computes the first order correction r_{i,j}p_{1,i,j}( x,t; x_{0}) of p(x,t;x_{0}), contained in (3.6), appropriately rescaled, that is computes the function p^{~}_{1,i,j}( x,t; x_{0}) defined by:
 (5.2) 
The use of the scale factor appearing in formulae (5.1) and (5.2) is due to a change of variable that has been used to compute the integral appearing in formula (2.1), that is:
 (5.3) 
In fact the change of variable (5.3) and the use of a Monte Carlo routine with sample points distributed according with the standard Gaussian distribution gives a satisfactory approximation of the integral in (2.1).
The libraries UnivPM_MathFinSolver1_2.dll and UnivPM_MathFinSolver1_3.dll can be used profitably on a parallel machine. In fact, for example, we can compute in parallel the first order order correction terms p^{~}_{1,i,j}( x,t;x_{0}), i=1,2,...,n1, j=i+1,i+2,...,n, that is if we dispose of N_{p} processors, each processor can compute (n(n1))/(2N_{p}) first order correction terms.
You can download the Monte Carlo main codes that evaluate the integral (2.1) using the routines UnivPM_MathFinSolver1_2.dll and UnivPM_MathFinSolver1_3.dll. These two main codes are a C and a Fortran version of the implementation of the numerical method described previously that evaluates the integral (2.1) using the routines UnivPM_MathFinSolver1_2.dll and UnivPM_MathFinSolver1_3.dll.
Included with these main codes you will find a manual, i.e. Monte_Carlo_readme.txt, that helps the understanding of how to use the two versions of the main code and explains how to choose the input parameters in order to have a satisfactory performance of the routines UnivPM_MathFinSolver1_2.dll and UnivPM_MathFinSolver1_3.dll.
click here to download
We propose an experiment reported in [1.4] able to show the high parallel performance of the formulae implemented in the library UnivPM_MathFinSolver1 for the computation of p(x, t;x_{0}) through formula (3.6). For this experiment the numerical algorithm described previously has been coded in Fortran 77 language using MPI as message passing language and tested on the SP3 machine with 112 processors of the CASPUR  Roma  Italy computer center. Results similar to the ones shown here can be obtained with the routines UnivPM_MathFinSolver1_2.dll, UnivPM_MathFinSolver1_3.dll.
We consider n = 100 and we choose s_{1}=0.25 s_{2}=0.1 s_{i}=0.2, i = 3,4,¼, n, and r_{i,j} = 0.2, i = 1,2,¼,n1, j = i+1,i+2,¼,n and we evaluate the probability density function p(x,t;x_{0}) on the point x=x^{1}, x_{0} where x_{i}^{1} = 1.25, and x_{i,0} = 1.5, i = 1,2,¼,n at the time to maturity t = 0.5.
In order to approximate the series expansion (3.8) for i = 1,2,¼, n1, j = i+1,i+2,¼,n we consider:

 (7.2) 
Nt = 4  
processors  seconds 
5  21.10 
25  4.28 
50  2.15 
100  1.13 
Nt = 8  
processors  seconds 
5  559.88 
25  113.53 
50  56.85 
100  31.16 
We note that the speed up factors obtained are really satisfactory. In fact for example passing from 5 to 100 processors when Nt=4 we have a speed up factor (i.e. execution time of 5 processors divided by execution time of 100 processors divided by (100/5)) equal to 0.934 and when Nt=8 the same speed up factor is equal to 0.898. We note that passing from Nt= 4 to Nt = 8 the time required to make the evaluation of p_{0}(x,t;x_{0})+ p_{1,Nt}(x,t;x_{0}) increases dramatically. Table 2 shows the execution time required to evaluate p_{0}(x,t;x_{0})+ p_{1,Nt}(x,t;x_{0}) when the number of processors is fixed to be 5 and Nt = 2,4,8.
Nt  seconds 
2  1.61 
4  21.10 
8  559.88 
Note that when we pass from Nt= 2 to Nt= 4, that is, when we multiply Nt by a factor 2, the time required to compute the probability density is multiplied by a factor of about 13.10 . When we pass from Nt= 4 to Nt = 8 the same time is multiplied by a factor of about 26.53. To better understand this behavior, the observed times have been fitted with the formula c Nt ^{a} b^{Nt}, a>0, b>0, c>0. The result obtained with the fitting is a = 2.7 and b=1.421267, c=0.6132/(number of processors), where in the case considered in Table 2 the number of processors is 5. That is passing from Nt to 2Nt the time is multiplied by a factor of about 2^{2.7} (1.43527)^{Nt}. This behavior depends on the need of evaluating the functions L_{q}, q = 0,1,¼, defined in (3.11).
Table 3 shows that the time required to compute p_{0}( x,t;x_{0})+p_{1,Nt}( x,t;x_{0}) is O(n^{2}) when n® +¥. We choose x_{0} Î R^{n}_{+} and x Î R^{n}_{+} to be given respectively by x_{i,0} = 1, i = 1,2,¼,n, x_{i} = 1.1+0.001+0.05/n, i = 1,2,¼, n, and s_{1} = 0.5, s_{2} = 0.25, s_{i} = 0.65, i = 3,4,¼,n, n_{t} = 3, r_{i,j} = 0.1, i = 1,2,¼[n/2], j = i+1,i+2,¼,n, r_{i,j} = 0.2, i = [n/2]+1,[n/2]+2, ¼,n1, j = i+1,i+2,¼,n, for several values of the space dimension n, that is n = 16,32,64 and 128. Note that in Table 3 Nt is fixed, that is Nt is independent of n.
Table 3 has been obtained executing the computation on a personal computer Pentium IV 2.5 Ghz. The time shown in Table 3 has been measured using the Fortran 77 function secnds that allows to measure the time required by the computation in seconds.
n  seconds 
16  0.75 
32  2.85 
64  11.35 
128  45.35 
We note that when we go from n to 2n the time is multiplied by a factor of about 4. We note that our approach is very efficient, it is also accurate when the underlying factors (i.e assets or interest rates) are only slightly correlated. In fact in this case the first order approximation of p contained in (3.6) is a good approximation.
We invite the user to fill up the following
request of informations. We plan to use the informations acquired in this
way for the
following purposes:
1) improving the quality of this website and of the
software library UnivPM_MathFinSolver1,
2) enlarging our knowledge of mathematical finance through
the knowledge and understanding of
the problems of interest to the users of this website.

Acknowledgement: It is a pleasure to thank Mrs. Francesca Piersigilli of the Università di Camerino for the helpful assistance in the realization of the software library contained in this website.
JEL Classification Codes: G13, C63.