A method to compute the transition probability density associated to a multifactor Cox-Ingersoll-Ross model of the term structure of interest rates with no drift term *

Lorella Fatone, Graziella Pacelli, Maria Cristina Recchioni, Francesco Zirilli

 

 

  1. Abstract
  2. Introduction
  3. Notation and Mathematical Preliminaries
  4. An example of interactive evaluation of some financial derivatives  
  5. UnivPM_MathFinSolver1: a software library that computes the transition probability density of the n-dimensional square root process
  6. A Monte Carlo routine to evaluate financial derivatives using the routines UnivPM_MathFinSolver1_2.dll and UnivPM_MathFinSolver1_3.dll  
  7. The speed up factors obtainable with UnivPM_MathFinSolver1 using parallel computing  
  8. Informations request
  9. Contact us

 

 

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.

 

 

 

 

1. Abstract

We consider a n-dimensional square-root 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 square-root process to model the assets prices takes into account some spurious effects, for example the smile effect, that escapes the more usual log-normal 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 square-root 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].

References

 

[1.1]
COX, J.C., INGERSOLL, J.E. Jr. and ROSS, S.A.,"A theory of the term structure of interest rates", Econometrica 53, 385-408 (1985).

 

[1.2]
CHEN, R. and SCOTT, L., "Pricing interest rate options in a two-factor Cox-Ingersoll-Ross model of the term structure ", The Review of Financial Studies 5, 613-636 (1992).

 

[1.3]
NOWMAN, K.B. "Gaussian estimation and forecasting of multi-factor term structure models with an application to Japan and the United Kingdom", Asia-Pacific Financial Markets 8, 23-34 (2001).

 

[1.4]
FATONE, L., PACELLI, G., RECCHIONI, M.C. and ZIRILLI, F., "A method to compute the transition probability density associated to a multifactor Cox-Ingersoll-Ross model of the term structure of interest rates with no drift term", Journal of Nonlinear Analysis: Hybrid Systems 2, 144-183 (2008).

 

 

 

 

 

 

 

2. Introduction

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 n-dimensional 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 si, i=1,2,...,n and the correlation coefficients ri,j, i=1,2,...,n, j=1,2,...,n. We note that since we have ri,i=1, ri,j= rj,i, i=1,2,...,n, j=1,2,...,n the correlation coefficients that must be specified are n(n-1)/2, so that the total number of parameters to be specified to determine the price dynamics is n+n(n-1)/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 y1,y2,..., yn, 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 x1,x2,...,xn at the maturity time T. Let t be the time to maturity, t=T-t, 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:

C(y,t) =  
ò
 R+n  
 p(x,t;y)P(x,T) dx,   y Î R+n,   t>0,
(2.1)
where t=T-t and Rn+ = { x = (x1,x2,¼,xn)T Î Rn | xi > 0,  i = 1,2,¼,n} (see Section 3). Note that p(x,t;y) depends on the n+n(n-1)/2 parameters that determine the price dynamics (3.1), (3.2).

 

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 si, i=1,2, and the correlation coefficient r1,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 y1, y2 plane. The time ranges from t=0, that corresponds to the maturity time T, to t=1, that corresponds to the physical time t =T-1.

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 si, i=1,2,...,n , the correlation coefficients ri,j, i=1,2,...,n-1, 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 n-dimensional 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<si <0.6 , i=1,2,...,n,

-0.6<ri,j <0.6 , i=1,2,...,n-1, 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;x0). 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 si, i=1,2,...,n , the correlation coefficients ri,j, i=1,2,...,n-1, j=i+1,i+2,...,n, the time to maturity t, the vector x0 of the prices of the underlyings at the time to maturity t, the payoff function P(x,T), xÎ Rn+. 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.

 

 

 

3. Notation and Mathematical Preliminaries

Let R be the real numbers, Rn, n ³ 1, be the n-dimensional real Euclidean space and let x = (x1,x2,¼,xn)T Î Rn 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 = (X1,X2,¼,Xn)T Î Rn, we consider the following stochastic differential equations:

dXi = si|Xi|1/2 dWi(t), t > 0, i = 1,2,¼,n,
(3.1)

X(0) = x0,
(3.2)
where x0 = (x1,0,x2,0,¼,xn,0)T Î Rn+, si > 0, i = 1,2,¼,n are the volatilities, dWi, i = 1,2,¼,n are the stochastic differentials of n standard Wiener processes Wi(t), t > 0, i = 1,2,¼,n such that:

< dWi dWj > = ri,jdt, i,j = 1,2,¼,n,
(3.3)
where < · > denotes the expected value of ·   . Due to the fact that < dWi > = 0, and ri,i=1, i = 1,2,¼,n, we have that ri,j, i,j = 1,2,¼,n are the correlation coefficients of the variables dWi, i=1,2,...,n. The quantities ri,j, i,j = 1,2,¼,n are assumed to be constants. We have ri,j = rj,i, i,j = 1,2,¼,n. Let G = ((ri,j)), i,j = 1,2,¼,n be the correlation matrix, we assume that G is a positive definite matrix.

Note that the stochastic process X(t)= (X1(t), X2(t),...,Xn(t))T solution of (3.1), (3.2) belongs to Rn+ for any t > 0.

We denote with p(x,t;x0), x Î Rn+, t > 0, x0 Î Rn+ the transition probability density function associated to the stochastic process X(t) = (X1(t),X2(t),¼,Xn(t))T, t > 0, solution of (3.1), (3.2). The function p(x,t;x0),  x Î Rn+, t > 0 , x0Î Rn+, satisfies the Fokker-Planck equation:

p
t
 = n
å
 i,j = 1 
 1
2
 2
xixj
 æ
è 		sisjri,j   __
Öxi
 
   __
Öxj
 
 p ö
ø 		,  x Î Rn+, t > 0,
(3.4)
with the following initial condition:

p(x,0;x0) = d (x-x0), x Î Rn+,
(3.5)
where d( x- x0) is the Dirac's delta concentrated in x = x0. In [1.4] we obtain the following formula for p( x,t;x0), x Î Rn+, t > 0, x0Î Rn+, solution of (3.4), (3.5):

p( x,t;x0)
=
p0( x,t;x0)+ n-1
å
 i = 1 
 n
å
 j = i+1 
 ri,jp1,i,j( x,t;x0)+o æ
è 		é
ë 		n-1
å
 i = 1 
 n
å
 j = i+1 
 ri,j2 ù
û 		1/2
 
 ö
ø
x Î Rn+,  t > 0,  x0Î Rn+,   ri,j® 0, i = 1,2,¼, n-1,  j = i+1,i+2,¼,n,
(3.6)
where o(·) is the Landau symbol, and the functions p0(x,t;x0) and p1,i,j(x,t;x0), i=1,2,...,n-1, j=i+1,i+2,..., n are given by:

p0( x,t;x0) = n
Õ
 i = 1 
 Y0(xi,t;xi,0,si),    x Î Rn+,  t > 0, x0Î Rn+,
(3.7)

p1,i,j( x,t;x0) = n
Õ
 m = 1,m ¹ i,m ¹ j  
 Y0(xm,t;xm,0,sm) +¥
å
 p = 0 
 tp+1
p+1
 p
å
 m = 0 
 Lp-m(xi,t;xi,0,si)Lm(xj,t;xj,0,sj),
x Î Rn+,  t > 0, x0Î Rn+,
(3.8)
where the function Y 0 is given by following formula:

Y0(x,t;y,s) = 2
s2  t
 æ
ç
è 		y
x
 ö
÷
ø 		1/2

 
 exp(-2(x1/2-y1/2 )2/(s2 t))   K1 æ
ç
è 		4(xy)1/2
s2t
 ö
÷
ø
 x > 0, y > 0, t > 0, s > 0, (3.9)
and:

Ym(x,t;y,s) = æ
ç
è 		2
s2t
 ö
÷
ø 		m+1

 
 exp(-2(x1/2-y1/2 )2/(s2 t)) m
å
 j = 0 
 æ
ç
è
m
j
 ö
÷
ø 		(-1)m-j æ
ç
è 		x
y
 ö
÷
ø 		(m-1-j)/2

 
 Km-1-j æ
ç
è 		4 (xy)1/2
s2t
 ö
÷
ø 		,
       x > 0 y > 0,  t > 0,  s > 0,  m = 1,2,¼,(3.10)
where Kj(x) = e-xIj(x), and Ij , j = ¼,-2,-1,0,1,2,¼, are the modified Bessel functions of integer order (see [3.1]), pag.374-378). Finally the function Lq, q=0,1,... is given by the following formula:

 
Lq(x,t;y,s) = s y1/2 æ
ç
è 		s2
2
 ö
÷
ø 		q

 
 q
å
 m = 0 
 2m
å
 s = 0 
 y-s/2
2m-s
å
 n = [(2m-s+1)/2] 
 ym-n-s/2
(2m-s-n)!
 vs,m,n n
å
 m = [(2n-2m+s)/2] 
 h2m-s,n,mym (-1)m
m!
 Yq+m-s-n+m+1(x,t;y x,s),
 x > 0, y > 0, t > 0, s > 0, q=0,1,...   ,
(3.11)
where

h2m-s,n,m
=
m
å
 j = [(2n-(2m-s))/2] 
 (-1)m-j æ
ç
è
m
j
 ö
÷
ø 		æ
ç
è
2j
2n-(2m-s)
 ö
÷
ø 		(2m-s)-n-1
Õ
 l = 0 
 (j+l),
m=0,1,..., s=0,1,...,2m, n=[(2m-s+1)/2],...,2m-s, m=[(2n-2m+s)/2],..., n,
(3.12)
vs,m,n = æ
ç
ç
ç
ç
è
1
2
s
 ö
÷
÷
÷
÷
ø 		(2n-2m+2s-1)!!
2n-m+s
 s
+å
 m = 2  
 æ
ç
ç
ç
ç
è
1

2
s-m
 ö
÷
÷
÷
÷
ø 		[m/2]
å
 q= 1
 (2n-2m+2s-2q-1)!!

2n-m+s-q    

cq(y)

y-q/2

 æ
ç
è
m-q-1
m-2q
 ö
÷
ø 		(-1)m-2q æ
ç
è 		s2
2
 ö
÷
ø 		-q
,
 
            m=0,1,..., s=0,1,...,2m, n=[(2m-s+1)/2],...,2m-s,     (3.13)

cq(y) = (-s2)q
q! 32qyq/2
 q-1
Õ
 m = 0 
 (4-(2m+1)2),   y>0,   q = 1,2,¼,
(3.14)
where [·] denotes the integer part of ·, åkm = 0 when m < k, 0! = 1, (-1)!! = 1 and (2m-1)!! = Õmj=1 (2j-1), m = 1,2,¼.

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=T-t) when the prices of the assets at the current time t=T-t is y= ( y1,y2,...,yn)T is given by:

C(y,t) =    
ò
 R+n  
 p(x,t;y)P(x,T) dx,    yÎ R+n,   t>0.
(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.

Reference

[3.1]
ABRAMOWITZ, M. and STEGUN, I.A., Handbook of mathematical functions, Dover, New York, 1970.

 

 

 

4. An example of interactive evaluation of some financial derivatives

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(x1,x2)-K,0)      (4.1)

and by:

P(x,T)=max(max(x1,x2)-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 si, i=1,2, among the two choices: s1=0.4, s2=0.5 or s1=0.5, s2=0.3 and the correlation coefficient r1,2 among three choices that is r1,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=(y1,y2)T ranges on a rectangular grid contained in the nonnegative quadrant of the y1, y2 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 =T-1.

 

 

 

5. UnivPM_MathFinSolver1: a software library that computes the transition probability density of the n-dimensional square root process

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;x0) contained in (3.6)

click here to download

UnivPM_MathFinSolver1_2.dll computes the zeroth order approximation of p(x,t;x0) contained in (3.6) appropriately rescaled, that is computes the function   p~0( x,t; x0) defined by:

p~0( x,t; x0) = (2pt)n/2 æ
è 		n
Õ
 m = 1 
 sm xm1/2 exp(2(xm1/2-xm,01/2 )2/(s2m t)) ö
ø		 p0(x,t;x0),  x Î Rn+, t > 0,  x0 Î Rn+.
(5.1)
click here to download

UnivPM_MathFinSolver1_3.dll for a given couple of indices i, j, i = 1,2,¼,n-1, j = i+1,i+2,¼,n computes the first order correction ri,jp1,i,j( x,t; x0) of p(x,t;x0), contained in (3.6), appropriately rescaled, that is computes the function p~1,i,j( x,t; x0) defined by:

p~1,i,j( x,t; x0) = (2pt)n/2 æ
è 		n
Õ
 m = 1 
 sm xm1/2 exp(2(xm1/2-xm,01/2 )2/(s2m t)) ö
ø 		ri,jp1,i,j( x,t; x0),  x Î Rn+, t > 0,  x0 Î Rn+.
(5.2)
click here to download

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:

zi = 2
s t1/2
 (xi1/2-xi,01/2 ),  i = 1,2,¼,n.
(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;x0), i=1,2,...,n-1, j=i+1,i+2,...,n, that is if we dispose of Np processors, each processor can compute (n(n-1))/(2Np) first order correction terms.

 

 

 

 

6. A Monte Carlo routine to evaluate financial derivatives using the routines UnivPM_MathFinSolver1_2.dll and UnivPM_MathFinSolver1_3.dll

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

 

 

 

 

7. The speed up factors obtainable with UnivPM_MathFinSolver1 using parallel computing

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;x0) 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 s1=0.25 s2=0.1 si=0.2, i = 3,4,¼, n, and ri,j = 0.2, i = 1,2,¼,n-1, j = i+1,i+2,¼,n and we evaluate the probability density function p(x,t;x0) on the point x=x1, x0 where xi1 = 1.25, and xi,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,¼, n-1, j = i+1,i+2,¼,n we consider:
p1,i,jNt( x,t; x0) = n
Õ
 m = 1,m ¹ i,m ¹ j  
 Y0(xm,t;xm,0,sm) Nt
å
 p = 0 
 tp+1
p+1
 p
å
 l = 0 
 Lp-l(xi,t;xi,0,si)Ll(xj,t;xj,0,sj),
x Î Rn+,   x0 Î Rn+,  t > 0,  x0 Î Rn+,
(7.1)
that is the sum of the first Nt+1 terms of the series appearing in (3.8). In particular let

p1,Nt( x,t;x0) = n-1
å
 i = 1 
 n
å
 j = i+1 
 ri,jp1,i,jNt( x,t;x0),
(7.2)
Table 1 shows the execution time required to compute p0( x,t;x0)+p1,Nt( x,t;x0) when Nt = 4 and Nt = 8 and the number of processors used is 5, 25, 50 or 100. The time is measured using the function MPI-WTIME().

Table 1: Time versus number of processors, n=100

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 p0(x,t;x0)+ p1,Nt(x,t;x0) increases dramatically. Table 2 shows the execution time required to evaluate p0(x,t;x0)+ p1,Nt(x,t;x0) when the number of processors is fixed to be 5 and Nt = 2,4,8.

Table 2: Time versus Nt using 5 processors, n=100

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 bNt, 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 22.7 (1.43527)Nt. This behavior depends on the need of evaluating the functions Lq, q = 0,1,¼, defined in (3.11).

Table 3 shows that the time required to compute p0( x,t;x0)+p1,Nt( x,t;x0) is O(n2) when n® +¥. We choose x0 Î Rn+ and x Î Rn+ to be given respectively by xi,0 = 1, i = 1,2,¼,n, xi = 1.1+0.001+0.05/n, i = 1,2,¼, n, and s1 = 0.5, s2 = 0.25, si = 0.65, i = 3,4,¼,n, nt = 3, ri,j = 0.1, i = 1,2,¼[n/2], j = i+1,i+2,¼,n, ri,j = -0.2, i = [n/2]+1,[n/2]+2, ¼,n-1, 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.

Table 3: Time versus n using 1 processor, Nt=3

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.

 

 

 

8. Informations request

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.

First Name
Last Name
Job Title
Company
Address 1
Address 2
City
State
ZIP
Country
Phone
Fax
Cell
E-mail
Please specify
your professional
interests
Please give your
evaluation of this website and of
the library UnivPM_MathFinSolver1

 

 

 

9. Contact us

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.

Entry n.