Manuscript submitted to
DCDS
Supplement Volume 2007

Website: www.aimSciences.org
pp. 1–10

arXiv:cs/0702163v1 [cs.CE] 28 Feb 2007

FIRST PASSAGE TIME FOR MULTIVARIATE JUMP-DIFFUSION
STOCHASTIC MODELS WITH APPLICATIONS IN FINANCE

Di Zhang and Roderick V.N. Melnik
Mathematical Modelling and Computational Sciences
Wilfrid Laurier University
Waterloo, ON, Canada N2L 3C5
Abstract. The “first passage-time” problem is an important problem with a wide
range of applications in mathematics, physics, biology and finance. Mathematically,
such a problem can be reduced to estimating the probability of a (stochastic) process
first to reach a critical level or threshold. While in other areas of applications the
FPT problem can often be solved analytically, in finance we usually have to resort
to the application of numerical procedures, in particular when we deal with jumpdiffusion stochastic processes (JDP). In this paper, we develop a Monte-Carlo-based
methodology for the solution of the FPT problem in the context of a multivariate
jump-diffusion stochastic process. The developed methodology is tested by using
different parameters, the simulation results indicate that the developed methodology
is much more efficient than the conventional Monte Carlo method, which establishes
itself as an efficient tool for further practical applications, such as the analysis of
default correlation and predicting barrier options in finance.

1. Introduction. In a jump-diffusion process (JDP), the dynamics of underlying
process have two random components: a continuous diffusion component and a discontinuous jump component [1], in which the jump component can be explained as
a sudden drop of process’s value. The first passage time (FPT) problem for jumpdiffusion processes has Attracted attention of researchers in such diverse fields as
queuing networks [2], computer vision [3], target recognition [4]. In the financial
world, many problems also require the information on the first passage time of
a stochastic process, for example, in modeling credit risk and valuing defaultable
securities [1], or in predicting barrier options [8]. Furthermore, it is now generally accepted that the geometric Brownian motion model for market behavior may
produce misleading results in [1, 9], such as mismatching credit spreads on corporate bonds, underlying derivative prices. Jump-diffusion processes have established
themselves as a sound alternative to the geometric Brownian motion model.
However, if we consider jumps in the process, except for very basic process types
where closed form solutions are available, when the jump sizes are doubly exponential or exponentially distributed [5], or when the jumps can have only nonnegative
values (assuming that the crossing boundary is below the process starting value)
[6], where closed form solutions are available, for most problems we can only resort
to the numerical procedures.
Monte Carlo simulation is a very promising candidate in dealing with FPT problems. However, in conventional Monte Carlo method, in order to avoid discretization
bias [7], we need to discretize the time horizon into small enough intervals, and to
2000 Mathematics Subject Classification. Primary: 60H35, 65C05; Secondary: 68U20.
Key words and phrases. Monte-Carlo simulations, first passage time, multivariate jumpdiffusion process.

1

2

DI ZHANG AND RODERICK V.N. MELNIK

evaluate the process at each discretized time that is very time-consuming. Recently,
Atiya and Metwally [8, 9] have developed a fast Monte Carlo-type numerical method
to solve the FPT problem for jump-diffusion process.
In many financial problems we have to deal with multiple processes in practice
and consider their jointly crossing the critical level, so it is very useful to develop
fast numerical procedure for multivariate jump-diffusion process. In this contribution, we extend the fast Monte Carlo-type numerical methodology to the more
general case that covers affine multivariate processes jump-diffusion. The developed
methodology can be easily extended to other financial applications and areas where
FPT problem arises.
The article is organized as follows: section 2 describes our mathematical model.
The algorithms are presented in section 3 and simulation results are given in section
4. Concluding remarks are given section 5.
2. Mathematical Model. In this section, first, we give brief discussion on the
affine jump-diffusion model, then we deduce the first passage time distribution
under multivariate jump-diffusion process. At last, we consider the issue of kernel
estimator.
2.1. Affine jump-diffusion. Affine jump-diffusion is a jump-diffusion process for
which the drift vector, “instantaneous” covariance matrix and jump intensities all
have affine dependence on the state vector.
Let us consider a complete probability space (Ω, F, P ) and an information filtration (Ft ), and suppose that X is a Markov process in some state space D ⊂ Rn ,
solving the stochastic differential equation [10]
dXt = µ(Xt )dt + σ(Xt )dWt + dZt ,

(1)

n

n

n×n

where W is an (Ft )-standard Brownian motion in R ; µ : D → R , σ : D → R
,
and Z is a pure jump process whose jumps have a fixed probability distribution ν
on Rn and arrive with intensity {λ(Xt ) : t ≥ 0}, for some λ : D → [0, ∞).
Definition 1. The above model is an affine model if [10]:
µ(Xt , t) = K0 + K1 Xt
(σ(Xt , t)σ(Xt , t)⊤ )ij = (H0 )ij + (H1 )ij Xj
λ(Xt ) = l0 + l1 · Xt ,
n

n×n

(2)
n×n

n×n×n

where K = (K0 , K1 ) ∈ R × R
, H = (H0 , H1 ) ∈ R
×R
, l = (l0 , l1 ) ∈
n
n×n
R
×R
and
the
“jump
transform”
(determines
the
jump-size
distribution)
ψ(c) =
R
n
exp(c,
z)dν(z),
for
c
∈
C
,
is
known
whenever
the
integral
is
well
defined.
The
n
R
“coefficients” (K, H, l, ψ) of X completely determine its distribution.
2.2. First passage time distribution. Now, we will consider the first passage
time distribution in the context of multiple processes. In order to obtain a computable multi-dimensional solutions of FPT distribution, we need to simplify Eq.
(1) and (2)based on the following assumptions:
1. Each Wt in Eq. (1) is independent;
2. K1 = 0, H1 = 0 and l1 = 0 that means the drift term, the diffusion process
(Brownian motion) and the arrival intensity are independent with state vector
X;
3. The distribution of jump-size Zt is also independent with X.

SIMULATIONS OF FIRST PASSAGE TIME

3

In this scenario, we can rewrite Eq. (1) as
dXt = µdt + σdWt + dZt ,

(3)

where
µ = K0 , σσ ⊤ = H0 , λ = l0 .
Atiya et al. [9] have deduced one-dimensional first passage time distribution in
time horizon [0, T ]. In order to judge multiple processes, from Eq. (3), we obtain
multi-dimensional formulas and simplify them into computable formulas. We will
highlight the main steps of our procedure below.
As defined, the multiple processes X can be written as X = [X1 , X2 , . . . ]⊤ . Let
us consider one of its components, a sub-process Xi that satisfies the following
stochastic differential equation:
X
σij dWj + dZi
dXi = µi dt +
j

= µi dt + σi dWi + dZi ,

(4)

where Wi is also a standard Brownian motion and σi is:
sX
2.
σij
σi =
j

We assume that in the interval [0, T ], Mi times of jumps happen for Xi . Let
the jump instants be T1 , T2 , · · · , TMi . Let T0 = 0 and TMi +1 = T . τj equals
interjump times, which is Tj − Tj−1 . Following the notation of [9], let Xi (Tj− ) be
the process value immediately before the jth jump, and Xi (Tj+ ) be the process
value immediately after the jth jump. The jump-size is then Xi (Tj+ ) − Xi (Tj− ),
and we can use this jump-size to generate Xi (Tj+ ) sequentially.
If we define Ai (t) as the event that process crossed the threshold Di (t) for the
first time in the interval [t, t + dt], we have
+
gij (t) = p(Ai (t) ∈ dt|Xi (Tj−1
), Xi (Tj− )).

(5)

If we only consider one interval [Tj−1 , Tj ], we can obtain [11, 12]
gij (t)

=

+
Xi (Tj−1
) − Di (t)
1
3
(t − Tj−1 )− 2 (Tj − t)− 2
2
2yi πσi
!
[Xi (Tj− ) − Di (t) − µi (Tj − t)]2
∗ exp −
2(Tj − t)σi2

+
[Xi (Tj−1
) − Di (t) + µi (t − Tj−1 )]2
∗ exp −
2(t − Tj−1 )σi2

where
+
[Xi (Tj−1
) − Xi (Tj− ) + µi τj ]2
1
exp −
yi = p
2τj σi2
σi 2πτj

!

!

,

(6)

.

After getting result in one interval, we combine the results to obtain the density
for the whole interval [0, T ]. Let B(s) be a Brownian bridge in the interval [Tj−1 , Tj ],

4

DI ZHANG AND RODERICK V.N. MELNIK

+
+
with B(Tj−1
) = Xi (Tj−1
), B(Tj− ) = Xi (Tj− ), the probability that the minimum of
B(si ) is always above the boundary level is [13]


+
+
Pij = P
inf
B(si ) > Di (t)|B(Tj−1
) = Xi (Tj−1
), B(Tj− ) = Xi (Tj− )
Tj−1 ≤si ≤Tj



+
2[Xi (Tj−1
)−Di (t)][Xi (Tj− )−Di (t)]

, if Xi (Tj− ) > Di (t),
1 − exp −
τj σi2
(7)
=

0,
otherwise.

Then B(si ) is below the threshold level, which means the default happens or
already happened, and its probability is 1 − Pij . Then let L(si ) ≡ Li denote the
index of the interjump period in which the time si falls in [TLi −1 , TLi ]. Also let Ii
represent the index of the first jump, which happened in simulated jump instant.
Ii

= min(j : Xi (Tk− ) > Di (t); k = 1, . . . , j, and
Xi (Tk+ ) > Di (t); k = 1, . . . , j − 1, and Xi (Tj+ ) ≤ Di (t)).

(8)

If no such Ii exists, we set Ii = 0.
Then we get the probability of the interval [0, T ]

=

+
P (Ai (si ) ∈ ds|Xi (Tj−1
), Xi (Tj− ), j = 1, . . . , Mi + 1)

QLi −1
P
if Li < Ii or Ii = 0,
 giLi (si ) k=1
QLi −1 ik QLi
Pik + k=1 Pik δ(si − TIi ) if Li = Ii ,
giLi (si ) k=1

0
if Li > Ii ,

(9)

where δ is the Dirac’s delta function.

2.3. The kernel estimator. For each Xi , after generating series of first passage
times si , we use a kernel density estimator with a Gaussian kernel to estimate
the first passage time density (FPTD) f . As described in [9], the kernel density
estimator is based on centering a kernel function of a bandwidth as follows:

where

N
1 X
K(h, t − si ),
fb =
N i=1



1
(t − si )2
.
K(h, t − si ) = p
exp − 2
h /2
π/2h

(10)

(11)

The optimal bandwidth in the kernel function K can be calculated as [15]:

−0.2
Z ∞
√
hopt = 2N π
(ft′′ )2 dt
,
(12)
−∞

where N is the number of generated points and ft is true density. Here we use the
approximation for the distribution as a gamma distribution as proposed in [9]:
ft =

αβ β−1
t
exp(−αt).
Γ(β)

In this case the functional becomes
Z ∞
5
X
Wi αi Γ(2β − i)
,
(ft′′ )2 dt =
2(2β−i) (Γ(β))2
0
i=1

(13)

(14)

SIMULATIONS OF FIRST PASSAGE TIME

5

where
W1 = A2 , W2 = 2AB, W3 = B 2 + 2AC, W4 = 2BC, W5 = C 2 ,
and
A = α2 , B = −2α(β − 1), C = (β − 1)(β − 2).
From Eq. (14), apparently, in order to get a nonzero bandwidth, we has constrained β to be at least equal to 3. Using this constraint, we can obtain the estib = 1 PN ti
mates of the parameters α and β via the method of moments: t = E(t)
i=1
N
q
PN 2
2
1
2
2
b
b
and E(t ) = N i=1 ti and the sample standard deviation is ζ = E(t ) − t . The
2
estimates are α
b = t2 and βb = t 2 ≥ 3.
ζ

ζ

The kernel estimator can be easily generalized to the multivariate case. Sup−
→
→
pose we consider X = [X1 , X2 , . . . , Xm ]⊤ , let t = [t1 , t2 , . . . , tm ], and −
si =
[s1i , s2i , . . . , smi ], sji (j = 1, 2, . . . , m) is the first passage time for Xj . Then, the
multivariate kernel density estimator with kernel K and window width h is defined
by [15]
N
1 X h −
→ →i
−
→
K h, ( t − −
si ) ,
(15)
fb( t ) =
N i=1
where



1 −
→⊤ −
→
−
→
2 −m/2
(16)
K( t ) = (2πh )
exp − 2 t t .
2h
And if we approximate the true density f as a unit m-variate normal density,
then the optimal bandwidth hopt is [15]

1/(m+4)
1
4
hopt = 1/(m+4)
.
(17)
2m + 1
N

3. Algorithms. In section 2, we have built up a multivariate jump-diffusion model
as describe in Eq. (3), and its first passage time distribution was also obtained in
section 2.2. In this section, we will discuss how to simulate the multivariate jumpdiffusion process efficiently via Monte Carlo method. In conventional Monte Carlo
method, the simulation is very straightforward, we divide the time horizon [0, T ]
into n small intervals [0, t1 ], [t1 , t2 ], · · · , [tn−1 , T ] and in each Monte Carlo run, we
need to calculate the value of Xi at each discretized time t. We should mention
that in order to exclude discretization bias, the number n must be large enough. It
is obvious that this conventional method is very time-consuming.
Recently, Atiya and Metwally [8, 9] have proposed two fast Monte Carlo type
methods, which are about 10-30 times faster than the conventional Monte Carlo approach. We called them uniform sampling (UNIF) method, which involves sampling
using uniform distribution, and inverse Gaussian density sampling (IG) method,
which uses inverse Gaussian density method for sampling.
In this article, we mainly focus on the uniform sampling (UNIF) method and
extend it to the multivariate jump-diffusion process. The major improvement of
UNIF method is that it only evaluate Xi at generated jump instants and between
each two jumps the process is a Brownian bridge, so we just consider the probability
of Xi crossing the boundary level in (Tj−1 , Tj ) instead of evaluating Xi at each
+
discretized time t. More exactly, we assume that the values of Xi (Tj−1
) and Xi (Tj− )
are known as two end points of Brownian bridge, we generate a variable si with
uniform distribution and by using Eq. (7) to see whether Xi (si ) is smaller than the

6

DI ZHANG AND RODERICK V.N. MELNIK

threshold level, if it defaults, then we have successfully generated a first passage
time si and can neglect the other intervals and perform another Monte Carlo run.
In [8, 9], the jump-diffusion process is involved in a univariate model. However,
our model is based on multivariate process in which Xi are correlated as described
in Eq. (3), so we need to consider several points as follows:
1. In this article, as a first step, we assume that the arrival rate λ for the Poisson
jump process and the distribution of (Tj − Tj−1 ) are the same for each Xi .
As for jump-size, we should generate them as given distribution, and it can
be different to reflect the different jump process for each Xi .
2. At present, we use exponential distribution (mean value µT ) for (Tj − Tj−1 )
and normal distribution (mean value µJ and standard deviation σJ ) for the
jump-size. Of course, we can use any distribution as desired.
3. If we consider m processes, i.e., X = [X1 , X2 , . . . , Xm ]⊤ , then we need an array IsDefault (whose size is m) to indicate whether process Xi has crossed the
threshold in this Monte Carlo run. If Xi has crossed, then we set IsDefault(i) =
1, and will not evaluate it during this Monte Carlo run.
Next, we will give a description of our algorithm, based on a multivariate extension of the algorithms proposed in [8, 9].
3.1. Uniform sampling method. Let us consider m processes in the given time
horizon [0, T ], as described above, we have generated the jump instant Tj by generating interjump times (Tj − Tj−1 ), besides we set IsDefault(i) = 0(i = 1, 2, · · · , m)
at first.
From Eq. (4), we can see that,
1. If jump doesn’t occur, the diffusion follows a standard Brownian motion,
+
Wi (T ) ∼ N (0, T ), so interjump size (Xi (Tj− ) − Xi (Tj−1
)) follows a normal
p
distribution of mean µi (Tj − Tj−1 ) and standard deviation σi Tj − Tj−1 .
After extend if necessary, we get
m
X
−
+
Xi (Tj ) ∼ Xi (Tj−1 ) + µi (Tj − Tj−1 ) +
σik N (0, Tj − Tj−1 ),
k=1

and the initial state is Xi (0) = Xi (T0+ ).
2. If jump occurs, the jump-size and direction of Zi (Tj ) are not fixed either. We
simulate the jump-size by a normal distribution, and of course we may generate it according to other distribution. Then we can compute the postjump
value:
Xi (Tj+ ) = Xi (Tj− ) + Zi (Tj ).

After generating beforejump and postjump value Xi (Tj− ) and Xi (Tj+ ) (j =
1, · · · , M , M is the total number of jumps for all the processes Xi ), we can compute
Pij according to Eq. (7). To recur the first passage time density (FPTD) fi (t), we
need to consider three conditions for each Xi that is still above the threshold:
+
1. First passage happens inside the interval. We know if Xi (Tj−1
) >
−
Di (t), Xi (Tj ) < Di (t), then the first passage happened in the time interval
[Tj−1 , Tj ]. To judge when the first passage happen, first we compute the
probability Pij of Xi always above the threshold according to Eq. (7), then
T −T
, and ui is a uniform
we generate si as si = bij ui + Tj−1 , where bij = j1−Pj−1
ij
random number in [0, 1]. If si also belongs to interval [Tj−1 , Tj ], then the first
passage time occurred in this interval, where si is the first passage time and

SIMULATIONS OF FIRST PASSAGE TIME

7

now we set IsDefault(i) = 1 to indicate process Xi has crossed the critical
level. In this condition, we can get conditional interjump first passage time
density of that specific interval
by Eq.(6).
To get the density of whole interval


Tj −Tj−1
b
[0, T ], we have fi,n (t) =
gij (si ) ∗ K(hopt , t − si ), where n is the
1−Pij

iteration number of Monte Carlo cycle.
2. First passage doesn’t happen in this interval. If si doesn’t belong to
interval [Tj−1 , Tj ], then the first passage time has not yet occur in this interval.
3. First passage happens in the right boundary of interval. If Xi (Tj+ ) <
Di (t), Xi (Tj− ) > Di (t), which follow the definition in Eq. (8), then obviously
TIi (set Ii = j) is the first passage time. Evaluate the density function using
kernel function fbi,n (t) = K(hopt , t − TIi ), and set IsDefault(i) = 1.
Then we increase j and examine the next interval and judge these three conditions
for each non-crossing process Xi again. For each Monte Carlo run, if we make a
rough assumption that the probability of Xi crossing the threshold

 is not correlated,
−
→
T −T
m
b
then we can obtain the multivariate FPTD as fn ( t ) = Πi=1 j1−Pj−1
gij (si ) ∗
ij
−
→ −
→
K(h , t − s ).
opt

After running N times of Monte Carlo cycle, we get the one-dimensional FPTD
PN b
PN b −
→
→
1
b−
of Xi as fbi (t) = N1
n=1 fi,n (t), and multivariate FPTD as f ( t ) = N
n=1 fn ( t )

4. Simulation results. In this section, as a demonstration, we will test the multivariate UNIF method on two-dimensional case. In order to check the efficient
and validity of the UNIF method, we use three examples with different arrival rate
λ = 1, 3, 8 for the Poisson jump process to judge the efficiency of our algorithms.
The parameters are as follows
X0 = [0, 0]⊤ , D(t) = [ln(0.9) − 0.002t, ln(0.95) − 0.012t]⊤


0.2 0.0
⊤
µ = [−0.002, −0.012] , σ =
,
0.0 0.2
µZ = [0, 0]⊤ , σZ = [0.2, 0.12]⊤.

where X0 is the starting value for the process, D( t) is the threshold, µ is the constant
instantaneous drift, σ represents the Brownian motion, and µZ and σZ are the mean
and standard deviations, respectively, of the jump sizes.
The simulation was carried out with total Monte Carlo runs N = 500, 000 in
horizon [0, 1]. Moreover, we have also carried out conventional Monte Carlo simulation with the same parameters, the estimated density functions are displayed in
Fig. 1-3. All the simulations were carried out on a 2.4 GHz AMD Opteron(tm)
Processor. The optimal bandwidth and CPU time are described in Table 1.
From Fig. 1-3, we can see that the multivariate UNIF method gives similar
density function as the conventional one, that check the validity of our algorithms.
An interesting phenomenon is that increasing λ will affect FPTD of X1 a lot whose
probability of crossing the threshold is low in the interval [0, T ].
Undoubtedly, from Table 1, one can easily realizes that the multivariate UNIF
approach is much more efficient than the conventional one, which establish it as a
fast methodology for practical applications.
5. Conclusion. In summary, we have studied the first passage time problem in the
context of multivariate jump-diffusion processes. We have extended the fast Monte
Carlo-type numerical method - the UNIF method – to multiple processes. From our

8

DI ZHANG AND RODERICK V.N. MELNIK

Figure 1. Example 1 (λ = 1): One-dimensional (top) and twodimensional density function (bottom) estimate using 100, 000 iterations for UNIF and conventional Monte Carlo approaches. The
discretization size of time horizon is ∆ = 0.0002 for conventional
Monte Carlo method.

Figure 2. Example 2 (λ = 3): One-dimensional (top) and twodimensional density function (bottom) estimate using 100, 000 iterations for UNIF and conventional Monte Carlo approaches. The
discretization size of time horizon is ∆ = 0.0002 for conventional
Monte Carlo method.

SIMULATIONS OF FIRST PASSAGE TIME

9

Figure 3. Example 3 (λ = 8): One-dimensional (top) and twodimensional density function (bottom) estimate using 100, 000 iterations for UNIF and conventional Monte Carlo approaches. The
discretization size of time horizon is ∆ = 0.0002 for conventional
Monte Carlo method.
Table 1. The optimal bandwidth hopt , and CPU time per Monte
Carlo run of the simulations. The first hopt is for X1 and the
second for X2 . All the simulations were performed with Monte
Carlo runs N = 100, 000, besides, for conventional Monte Carlo
(CMC) method, the discretization size of time horizon is ∆ =
0.0002.

Example 1
Example 2
Example 3

Optimal bandwidth CPU time
X1
X2
CMC 0.012664 0.006943 0.286642
UNIF 0.016030 0.013880 0.000527
CMC 0.011157 0.006582 0.284554
UNIF 0.012249 0.009443 0.000731
CMC 0.008894 0.005921 0.299156
UNIF 0.009117 0.006542 0.001222

simulation results, we can see that the multivariate UNIF approach is much more
efficient than the conventional Monte Carlo method, which illustrates that the developed methodology can provide an efficient tool for further practical applications,
such as the analysis of default correlation and predicting barrier options in finance.
Acknowledgements. We would like to thank NSERC for its support.
REFERENCES
[1] C. Zhou, The Term Structure of Credit Spreads with Jump Risk, Journal of Banking and
Finance 25 (2001), 2015–2040.

10

DI ZHANG AND RODERICK V.N. MELNIK

[2] G. Gorov, Ya. Kogan and N. Paradizov, Jump diffusion approximation in single-server systems with interruption of service and variable rate of arrival of calls, Avtomatika i Telemekhanika, 6 (1985), 44–51.
[3] S. C. Zhu, Stochastic jump-diffusion process for computing medial axes in Markov random
fields, IEEE Trans. Pattern Analysis and Machine Intelligence, 21 (1999), 1158–1169.
[4] M. Miller, U. Grenander, J. O’Sullivan and D. Snyder, Automatic target recognition organized
via jump-diffusion algorithm, IEEE Trans. Image Processing, 6 (1997), 157–174.
[5] S.G. Kou, H. Wang, First passage times of a jump diffusion process, Adv. Appl. Probab. 35
(2003), 504–531.
[6] I. Blake, W. Lindsey, Level-crossing problems for random processes, IEEE Trans. Inform.
Theory IT-19 (1973), 295–315.
[7] P.E. Kloeden, E. Platen, H. Schurz, “Numerical Solution of SDE Through Computer Experiments”, Third Revised Edition, Springer, Germany, 2003.
[8] S. Metwally, A. Atiya, Using brownian bridge for fast simulation of jump-diffusion processes
and barrier options, The Journal of Derivatives 10 (2002), 43–54.
[9] A.F. Atiya, S.A.K. Metwally, Efficient Estimation of First Passage Time Density Function
for Jump-Diffusion Processes, SIAM Journal on Scientific Computing 26 (2005), 1760–1775.
[10] D. Duffie, J. Pan, K. Singleton, Transform Analysis and Option Pricing for Affine JumpDiffusions, Econometrica 68 (2000), 1343–1376.
[11] W. Feller, “An Introduction to Probability Theory and Its Applications”, Vol. 1, Wiley, New
York, 1968.
[12] L.C.G. Rogers and D. Williams, “Diffusions, Markov Processes and Martingales”, Vol. 1,
Wiley, New York, 1994.
[13] I. Karatzas and S. Shreve, “Brownian Motion and Stochastic Calculus”, Springer-Verlag, New
York, 1991.
[14] O. Costa, M. Fragoso, Discrete time LQ-optimal control problems for infinite Markov jump
parameter systems, IEEE Trans. Automat. Control 40 (1995), 2076–2088.
[15] B.W. Silverman, “Density Estimation for Statistics and Data Analysis”, Chapman & Hall,
London, 1986.
[16] L. Devroye, “Non-uniform Random Variate Generation”, Springer-Verlag, New York, 1986.

Received Month Year; Sep.29, 2006.
E-mail address: rmelnik@wlu.ca

