# A Multivariate Stochastic Hybrid Model with Switching Coefficients and Jumps: Solution and Distribution

Authors: D. P. Siu; G. S. Ladde
Journal: Journal of Probability and Statistics
 (2011)

---

## Abstract

In this work, a class of multidimensional stochastic hybrid dynamic models is studied. The system under investigation is a first-order linear nonhomogeneous system of Itô-Doob type stochastic differential equations with switching coefficients. The switching of the system is governed by a discrete dynamic which is monitored by a non-homogeneous Poisson process. Closed-form solutions of the systems are obtained. Furthermore, the major part of the work is devoted to finding closed-form probability density functions of the solution processes of linear homogeneous and Ornstein-Uhlenbeck type systems with jumps.


---

## Body

## 1. Introduction

The study of stochastic hybrid systems exhibiting both continuous and discrete dynamics has been an area of great interest over the years. The properties of various types of stochastic hybrid systems have been studied extensively. Davis [1, 2] introduced a piecewise-deterministic Markov process, where transitions between discrete modes are triggered by random events and deterministic conditions for hitting the boundary, while the continuous-state process between jumps for the model is governed by a deterministic differential equation. Hespanha [3] proposed a model where transitions between modes are triggered by stochastic events much like transitions between states of a continuous-time Markov chains. Hu et al. [4] proposed a stochastic hybrid system where the deterministic differential equations for the evolution of the continuous-state process are replaced by Itô-Doob type stochastic differential equations [5, 6]. However, in this proposed model, the transitions are only triggered by hitting the boundaries. Siu and Ladde [7] studied a stochastic hybrid dynamic process where the transitions of its discrete time state are governed by either a non-homogeneous Poisson process or triggered by hitting the boundaries, while the continuous state is governed by a stochastic diffusion. A study of a stochastic hybrid system whose continuous time component is stochastic and altered by transitions of a finite-state Markov chain can be found in Chandra and Ladde [8], Korzeniowski and Ladde [9], and Ladde [10]. Yin et al. [11] provided an algorithm of the numerical solutions for a class of jump diffusions with regime switching. Mao and Yuan [12] and Yin and Zhu [13] summarize a wide range of properties of stochastic differential equations with Markovian switching coefficients. In view of applications, stochastic hybrid systems have been employed to diverse fields of studies, such as communication networks [3, 14], air traffic management [15, 16], and insurance risk models [17].In this work, we attempt to solve two fundamental problems in the stochastic modeling of dynamic processes described by Itô-Doob type stochastic differential equations with jumps. First, we find closed-form solutions of the stochastic hybrid systems. By using the closed-form solutions, we determine the closed-form probability density functions of solution processes of special cases of the general systems. The presented method provides an accessible way of obtaining the probability density functions without solving the Fokker-Planck partial differential equations [18] or approximating their solutions.The rest of the paper is organized as follows. In Section2, closed-form solution processes of the multidimensional systems are obtained through utilizing the result of A. G. Ladde and G. S. Ladde [19] piecewisely on the intervals between jumps.  The problem of finding the closed-form probability density functions is investigated in a systematic and coherent way. In Section 3 the probability density function of the solution process of one-dimensional linear homogeneous system of Itô-Doob type stochastic differential equations is derived. This is an extension of the geometric Brownian motion processes [20]. Then, by using the concept of modal matrix, the probability density function of the solution process of n-dimensional linear homogeneous systems is obtained in Section 4. The probability distribution of the solution process of the system with continuous dynamic consisting of only drift part and additive noise, namely, a Ornstein-Uhlenbeck system [18], is extended to hybrid system in Section 5. Some concluding remarks are given in Section 6.


## 2. Model Formulation

In this section, we develop a conceptual stochastic model for dynamic processes in chemical, biological, engineering, medical, physical, and social science [19, 21, 22] that are under the influence of discrete time events. The continuous time dynamic of the stochastic model between jumps follows a first-order linear non-homogeneous system of Itô-Doob type stochastic differential equations. At random times, governed by a non-homogeneous Poisson process, the coefficients of the continuous time dynamic are switched, and the process is multiplied by a random factor which results in a discontinuous jump.Letx(t) be a real n-dimensional process. Ak and Bkj are n×n matrices for any k∈ℕ∪{0} and j=1,2,…,q. Let Ckr be n-dimensional vectors for any k∈ℕ∪{0} and r=1,2,…,p. Let the continuous dynamic of the process x(t) be determined by the following system of stochastic differential equations:(2.1)dx(t)=AN(t)x(t)dt+∑j=1qBN(t)jx(t)dwj(t)+∑r=1pCN(t)rdw¯r(t),t≥t0,x(t0)=x0>0,
where w(t)=(w1(t),…,wq(t)) and w¯(t)=(w¯1(t),…,w¯p(t)) are independent q-dimensional and p-dimensional standard Wiener processes, and N(t) is a non-homogeneous Poisson process with intensity λ(t). Here, we denote x=(x1,x2,…,xn)>0 as xi>0 for all i=1,2,…,n. When Ckr=0 for all k and r, system (2.1) reduces to a first-order linear homogeneous system of Itô-Doob type stochastic differential equations, given below
(2.2)dx(t)=AN(t)x(t)dt+∑j=1qBN(t)jx(t)dwj(t),t≥t0,x(t0)=x0>0.
Here, Ak and Bkj are such that solution process x(t) of (2.2) is nonnegative.To obtain solution process of system (2.1), we first consider the solution process of the initial-value system when A, Bj, and Cr’s are fixed over time, that is, the solution process on a subinterval between jumps. Under the condition that the matrices A,B1,B2,…,Bq pairwise commute, the solution can be explicitly obtained; see A. G. Ladde and G. S. Ladde [19] or Movellan [23]. We state the result in the following lemma.Lemma 2.1.
Letx(t)≡x(t,t0,x0) be the solution of the following initial value problem (IVP):
(2.3)dx(t)=Ax(t)dt+∑j=1qBjx(t)dwj(t)+∑r=1pCrdw¯r(t),t≥t0,x(t0)=x0,
then the x(t) is expressed by
(2.4)x(t)≡x(t,t0,x0)=exp⁡[(A-12∑j=1q(Bj)2)(t-t0)+∑j=1qBj(wj(t)-wj(t0))]x0+∑r=1p∫t0texp⁡[(A-12∑j=1q(Bj)2)(t-s)+∑j=1qBj(wj(t)-wj(s))]Crdw¯r(s)
provided that ABj=BjA and BjBj′=Bj′Bj for all j,j′=1,2,…,q.Letx(t,Tk,xk) be the solution to system (2.3) with t0∶=Tk, x0∶=xk, A∶=Ak, Bj∶=Bkj, and Cr∶=Ckr. Now, we consider the following system of two interconnected stochastic dynamics:(2.5)dx(t)=Ak-1x(t)dt+∑k=1qBk-1jx(t)dwj(t)+∑r=1pCk-1rdw¯r(t),Tk-1≤t<Tk,x(Tk-1)=xk-1,xk=zkx(Tk-,Tk-1,xk-1),
where zk, k=1,2,3,…, are iid positive random variables with z0=1, and x(Tk-,Tk-1,xk-1)=limt→Tk-x(t,Tk-1,xk-1). Here, we assume that N(t), w(t), w¯(t), and zk are independent.By applying Lemma2.1, piecewisely, on each interval between jumps, we then obtain the solution process of system (2.5). The result is given in the following proposition.Proposition 2.2.
IfAkBkj=BkjAk and BkjBkj′=Bkj′Bij for all k∈ℕ∪{0} and j,j′=1,2,…,q, then the solution to the system (2.5) is given by
(2.6)x(t)=(∏k=1N(t)zk)exp⁡[(AN(t)-12∑j=1q(BN(t)j)2)(t-TN(t))+∑j=1qBN(t)j(wj(t)-wj(TN(t)))]×∏k=1N(t)exp⁡[(Ak-1-12∑j=1q(Bk-1j)2)(Tk-Tk-1)+∑j=1qBk-1j(wj(Tk)-wj(Tk-1))]x0+∑i=1N(t)(∏k=iN(t)zk)exp⁡[(AN(t)-12∑j=1q(BN(t)j)2)(t-TN(t))+∑j=1qBN(t)j(wj(t)-wj(TN(t)))]×∏k=i+1N(t)exp⁡[(Ak-1-12∑j=1q(Bk-1j)2)(Tk-Tk-1)+∑j=1qBk-1j(wj(Tk)-wj(Tk-1))]×[∑r=1p∫Ti-1Tiexp⁡[(Ai-1-12∑j=1q(Bi-1j)2)(Ti-s)+∑j=1qBi-1j(wj(Ti)-wj(s))]Ci-1rdw¯r(s)]+∑r=1p∫TN(t)texp⁡[(AN(t)-12∑j=1q(BN(t)j)2)(t-s)+∑j=1qBN(t)j(wj(t)-wj(s))]CN(t)rdw¯r(s).Proof.
By applying the result of (2.4) on the subintervals [Tk-1,Tk), for k=1,…,N(t), and [TN(t),t), we have the solution to system (2.5) as the following piecewise function:
(2.7)x(s)={x(s,t0,x0)x(t0)=x0,t0≤s<T1,x(s,T1,x1)x(T1)=x1,T1≤s<T2,⋯x(s,Ti,xi)x(Ti)=xk,Tk≤s<Tk+1,⋯x(s,TN(t),xN(t))x(TN(t))=xN(t),TN(t)≤s<t,
where x(s,Tk,xk) is a solution process in (2.4) with t0∶=Tk, x0∶=xk, A∶=Ak, Bj∶=Bkj, and Cr∶=Ckr, then, we have
(2.8)x(t)=x(t,TN(t),xN(t))=exp⁡[(AN(t)-12∑j=1q(BN(t)j)2)(t-TN(t))+∑j=1qBN(t)j(wj(t)-wj(TN(t)))]xN(t)+∑r=1p∫TN(t)texp⁡[(AN(t)-12∑j=1q(BN(t)j)2)(t-s)+∑j=1qBN(t)j(wj(t)-wj(s))]CN(t)rdw¯r(s)=zN(t)exp⁡[(AN(t)-12∑j=1q(BN(t)j)2)(t-TN(t))+∑j=1qBN(t)j(wj(t)-wj(TN(t)))]×[exp⁡⁡[(AN(t)-1-12∑j=1q(BN(t)-1j)2)(TN(t)-TN(t)-1)+∑j=1qBN(t)-1j(wj(TN(t))-wj(TN(t)-1))]xN(t)-1+∑r=1p∫TN(t)-1TN(t)exp⁡[(AN(t)-1-12∑j=1q(BN(t)-1j)2)(TN(t)-s)+∑j=1qBN(t)-1j(wj(TN(t))-wj(s))]CN(t)-1rdw¯r(s)]+∑r=1p∫TN(t)texp⁡[(AN(t)-12∑j=1q(BN(t)j)2)(t-s)+∑j=1qBN(t)j(wj(t)-wj(s))]CN(t)rdw¯r(s).
Next, substitute  xN(t)-1=zN(t)-1x(TN(t)-1-,TN(t)-2,xN(t)-2)=zN(t)-1x(TN(t)-1,TN(t)-2,xN(t)-2). The term x(TN(t)-1-,TN(t)-2,xN(t)-2) can be replaced by x(TN(t)-1,TN(t)-2,xN(t)-2) because the solution process is continuous between jumps. Repeating the substitution gives the desired result.In the following, we present two important particular by-products of Proposition2.2. When Ckr=0 for all k and r, system (2.5) reduces to the following first-order linear homogeneous system of Itô-Doob type stochastic differential equations with jumps:(2.9)dx(t)=Ak-1x(t)dt+∑k=1qBk-1jx(t)dwj(t),Tk-1≤t<Tk,x(Tk-1)=xk-1,xk=zkx(Tk-,Tk-1,xk-1).The solution of the above system is given in the following corollary. The result follows from Proposition2.2 by letting Ckr be zero for all k and r.Corollary 2.3.
IfAkBkj=BkjAk and BkjBkj′=Bkj′Bij for all k∈ℕ∪{0} and j,j′=1,2,…,q, then the solution of system (2.9) is given by
(2.10)x(t)=(∏k=1N(t)zk)exp⁡[(AN(t)-12∑j=1q(BN(t)j)2)(t-TN(t))+∑j=1qBN(t)j(wj(t)-wj(TN(t)))]×∏k=1N(t)exp⁡[(Ak-1-12∑j=1q(Bk-1j)2)(Tk-Tk-1)+∑j=1qBk-1j(wj(Tk)-wj(Tk-1))]x0.In the case whenBkj are zeros for all j and k, system (2.5) becomes a linear system with additive noise. The continuous dynamics between jumps are now governed by Ornstein-Uhlenbeck equations as follows:(2.11)dx(t)=Ak-1x(t)dt+∑r=1pCk-1rdw¯r(t),Tk-1≤t<Tk,x(Tk-1)=xk-1,xk=zkx(Tk-,Tk-1,xk-1).
Denote Ck=(Ck1,Ck2,…,Ckp) for all k, then the above system can be rewritten as(2.12)dx(t)=Ak-1x(t)dt+Ck-1dw¯(t),Tk-1≤t<Tk,x(Tk-1)=xk-1,xk=zkx(Tk-,Tk-1,xk-1),
where Ck's are n×p matrices, and w¯(t) is a p-dimensional standard Wiener process.Corollary 2.4.
The solution of system (2.11) is given by
(2.13)x(t)=(∏k=1N(t)zk)exp⁡[AN(t)(t-TN(t))]∏k=1N(t)exp⁡[Ak-1(Tk-Tk-1)]x0+∑i=1N(t)(∏k=iN(t)zk)exp⁡[AN(t)(t-TN(t))]∏k=i+1N(t)exp⁡[Ak-1(Tk-Tk-1)]×[∑r=1p∫Ti-1Tiexp⁡[Ai-1(Ti-s)]Ci-1rdw¯r(s)]+∑r=1p∫TN(t)texp⁡[AN(t)(t-s)]CN(t)rdw¯r(s).


## 3. Probability Distribution of One-Dimensional Linear Homogeneous Models

In this section, we will derive the probability density function of the solution process of the scalar version of system (2.9). Now, x(t) takes values in ℝ+, and in this case, Ak and Bk are scalars for all k in the system (2.9) and the solution process (2.10). Some auxiliary results are presented below. The following lemma provides the joint density function of the jump times given the number of jumps due to the non-homogeneous Poisson process N(t). The proof of this result can be found in many textbooks, for example, see Cox and Lewis [24] or Crowder et al. [25].Lemma 3.1.
For a non-homogeneous Poisson processN(t) with intensity λ(t), the joint density function of the jump times T1,T2,…,Tl conditioned on N(t)=l is given by
(3.1)fT1,T2,…,Tl∣N(t)=l(t1,t2,…,tl)=l!∏k=1lλ(tk)(∫0tλ(u)du)l.Next lemma gives the conditional probability density function ofx(t) given the number of jumps and the jump times.Lemma 3.2.
Given thatN(t)=l, and T1=t1,…,Tl=tl, x(t) has a probability density function as
(3.2)fx(t)∣N(t)=l,t1,…,tl(x)=1x∫-∞∞h*l(ln⁡x-s)ϕ(s;μ,σ)ds,x>0,
where h*l is the lth convolution of the common probability density function h of lnzk, for k=1,2,…,l, and ϕ(·;μ,σ) denotes the normal density function with mean μ and variance σ.Proof.
Given thatN(t)=l, and T1=t1,…, Tl=tl, from (2.10), we have
(3.3)ln⁡x(t)=∑k=1lln⁡zk+ln⁡x0+[(Al-Bl22)(t-tl)+Bl(w(t)-w(tl))]+∑k=1l[(Ak-1-Bk-122)(tk-tk-1)+Bk-1(w(tk)-w(tk-1))]=V+S,
where we denote V=∑k=1llnzk and S as the sum of last three terms in (3.3).
Sinceh is the common probability density function of lnzk, V as the sum of l iid random variables has the probability density function as the lth convolution h*l(v). We further note that S is the sum of l+1 independent normal variables due to the independent increment property of Wiener process. Then S is normally distributed with mean and variance as(3.4)μ≡μ(t,t1,…,tl)=ln⁡x0+(Al-Bl22)(t-tl)+∑k=1l(Ak-1-Bk-122)(tk-tk-1),σ≡σ(t,t1,…,tl)=Bl2(t-tl)+∑k=1lBk-12(tk-tk-1).
Since zk and w(t) are independent, then V and S are also independent. Using transformation method [26] on lnx(t)=V+S, we can obtain the conditional probability density function of lnx(t) as
(3.5)fln⁡x(t)∣N(t)=l,t1,…,tl(x)=∫-∞∞h*l(x-s)ϕ(s;μ,σ)ds.
If follows that
(3.6)fx(t)∣N(t)=l,t1,…,tl(x)=fln⁡x(t)∣N(t)=l,t1,…,tl(ln⁡x)1x=1x∫-∞∞h*l(ln⁡x-s)ϕ(s;μ,σ)ds,x>0.Having obtained the conditional probability density function ofx(t), we can derive the marginal probability distribution of the solution process in the one-dimensional case as follows.Proposition 3.3.
The probability density function of the scalar version of the solution processx(t) to the system (2.5) is given by
(3.7)fx(t)(x)=∑l=0∞[∫0t⋯∫0t3∫0t2∫-∞∞h*l(ln⁡x-s)ϕ(s;μ,σ)ds∏k=1lλ(tk)dt1dt2⋯dtl×1xexp⁡[-∫0tλ(u)du]].Proof.
From (3.1) and (3.2), the joint density function of x(t), T1,…, Tl given the condition N(t)=l is given by
(3.8)fx(t),T1,…,Tl∣N(t)=l(x,t1,…,tl)=fx(t)∣N(t)=l,t1,…,tl(x)fT1,T2,…,Tl∣N(t)=l(t1,t2,…,tl)=1x∫-∞∞h*l(ln⁡x-s)ϕ(s;μ,σ)dsl!∏k=1lλ(tk)(∫0tλ(u)du)l.
Integrating with respect to t1,t2,…, and tl gives the probability density function of x(t) given that N(t)=l as
(3.9)fx(t)∣N(t)=l(x)=∫0t⋯∫0t3∫0t2∫-∞∞h*l(ln⁡x-s)ϕ(s;μ,σ)ds∏k=1lλ(tk)dt1dt2⋯dtll!x(∫0tλ(u)du)l.
By multiplying the above density by the probability of l jumps, namely,
(3.10)Pr(N(t)=l)=exp⁡[-∫0tλ(u)du](∫0tλ(u)du)ll!,
and then taking the summation over l, the marginal probability density function of x(t) given in (3.7) is established.


## 4. Probability Distribution of Multivariate Linear Homogeneous Models

In this section, we will derive the probability density function of the solution processx(t) for the n-dimensional stochastic system (2.9) under the following assumptions. (i)
The drift and diffusion coefficients,Ak and Bkj for all k∈ℕ∪{0} and j=1,2,…,q, are diagonalizable.(ii)
The coefficients in each regime pairwise commute, that is,AkBkj=BkjAk and BkjBkl=BklBkj for all k∈ℕ∪{0} and j,l=1,2,…,q.(iii)
Ak∈𝒞 for k∈ℕ∪{0}, where 𝒞 denotes the set of n×n diagonalizable matrices whose eigenvector matrix, M, is such that M-1x>0 for any x>0.We first consider the stochastic system on the interval between jumps. Given thatN(t)=l and T1=t1, T2=t2,…,Tl=tl, consider the following SDE on [tk,tk+1), for some k=0,1,…,l:(4.1)dx(s)=Akx(s)dt+∑j=1qBkjx(s)dwj(s),tk≤s<tk+1,x(tk)=xk.
In the following, we provide the necessary background material that will be used, subsequently. The following result provides a way to find a modal matrix that can diagonalize the coefficients in the above system.Theorem 4.1 (see [27]).
A set of diagonalizable matrices commutes if and only if the set is simultaneously diagonalizable, that is, there exists an invertible matrix that can diagonalize all the matrices simultaneously.Remark 4.2.
In fact, the set of diagonalizable and commuting matrices shares the same set of independent eigenvectors. The eigenvector matrix is the one that simultaneously diagonalizes all the matrices in this set, and the resulting diagonal elements are the eigenvalues of each matrix.We recall [26] that a random vector x is said to have a n-dimensional multivariate normal distribution with mean and covariance matrix, μx and Σx, if its probability density function is given by(4.2)f(x)=(2π)-n/2|Σx|-1/2exp⁡[-12(x-μx)TΣx-1(x-μx)].
The following lemma gives the useful fact that the linear transformation of a multivariate normal random vector is again multivariate normally distributed.Theorem 4.3 (see [26]).
Ifx has a multivariate normal distribution with mean and covariance matrix, μx and Σx, then y=Px+c as a linear transformation of x follows also multivariate normal distribution with mean and covariance matrix, μy=Pμx+c and Σy=PΣxPT.Proof.
Sinceyi=∑j=1nPijxj+ci and every linear combination of normal random variables is still normal, then y follows a multivariate normal distribution. By the linearity of expectation we have
(4.3)μy=E[y]=E[Pμx+c]=PE[x]+c=Pμx+c,Σy=E[(y-μy)(y-μy)T]=E[[P(x-μx)][P(x-μx)]T]=E[P(x-μx)(x-μx)TPT]=PE[(x-μx)(x-μx)T]PT=PΣxPT.To find the probability density function ofx(s) satisfying the SDE (4.1), we need to introduce some notations and definitions that will be used, subsequently. First note that according to Theorem 4.1 and Remark 4.2, Aks and Bkjs have the same eigenvector matrix, denoted by Mk. Moreover, Ãk≡Mk-1AkMk=diag(ãk1,…,ãkn) and B̃kj≡Mk-1BkjMk=diag(b̃k,j1,…,b̃k,jn) where {ãk1,…,ãkn} and {b̃k,j1,…,b̃k,jn} are the sets of eigenvalues of Ak and Bkj, respectively, for all k,j. Next, we define a linear transformation y=Mk-1x, then x=Mky. Define lny=(lny1,lny2,…,lnyn) [28].The below proposition gives the probability distribution of solution processx(t) on the interval [tk,tk+1) over which the coefficients are constant.Lemma 4.4.
Under assumptions (i)–(iii), the processx(s) satisfying the SDE (4.1) has a probability density function given by
(4.4)fx(s)(x)=(2π)-n/2|Σk(s)|-1/2|Mk-1|∏i=1n1(Mk-1x)i×exp⁡[-12(ln⁡(Mk-1x)-μk(s))T(Σk(s))-1(ln⁡(Mk-1x)-μk(s))],fors∈[tk,tk+1),
where
(4.5)μk(s)=(ln⁡y1(tk)+(ãk1-12∑j=1q(b̃k,j1)2)(s-tk)ln⁡y2(tk)+(ãk2-12∑j=1q(b̃k,j2)2)(s-tk)⋮ln⁡yn(tk)+(ãkn-12∑j=1q(b̃k,jn)2)(s-tk)),Bk*=(b̃k,11b̃k,21⋯b̃k,q1b̃k,12b̃k,22⋯b̃k,q2⋮⋮⋮⋮b̃k,1nb̃k,2n⋯b̃k,qn),Σk(s)=(s-tk)Bk*(Bk*)T.Proof.
Fors∈[tk,tk+1), we have y(s)=Mk-1x(s), and x(s)=Mky(s). By multiplying Mk-1 on both sides of the SDE (4.1), we obtain the transformed SDE as follows:
(4.6)Mk-1dx(s)=Mk-1Akx(s)ds+∑j=1qMk-1Bkjx(s)dwj(s),⟹dMk-1x(s)=Mk-1AkMky(s)ds+∑j=1qMk-1BkjMky(s)dwj(s),⟹dy(s)=Ãky(s)ds+∑j=1qB̃kjy(s)dwj(s),
where Ãk and B̃kj are diagonal matrices as defined before.
From the application of Lemma2.1 with Ckr=0, the solution process of the transformed system (4.6) is(4.7)y(s)=exp⁡[(Ã-12∑j=1q(B̃kj)2)(s-tk)+∑j=1qB̃kj(wj(s)-wj(tk))]yk,
for s∈[tk,tk+1).
It follows from assumption (iii) thaty(s)=Mk-1x(s)>0 since x(s)>0, then, we can rewrite system (4.7) in the following form:(4.8)ln⁡y(s)=μk(s)+Bk*(w(s)-w(tk)),
where μk(s) and Bk* are defined above.
Sincew(t) is a standard Wiener process, then w(s)-w(tk) has a multivariate normal distribution with mean zero and covariance matrix (s-tk)In, where In is the n×n identity matrix. Then, by Theorem 4.3, lny(s) as a linear transformation of w(s)-w(tk) is also multivariate normally distributed with mean μk(s) and covariance matrix Σk(s), where(4.9)Σk(s)=Bk*[(s-tk)In](Bk*)T=(s-tk)Bk*(Bk*)T,
and the (u,v)th element of Σk(s) is (s-tk)∑j=1qb̃k,jub̃k,jv, for u,v=1,2,…,n. The probability density function of lny(s) is
(4.10)fln⁡y(s)(ỹ)=(2π)-n/2|Σk(s)|-1/2exp⁡[-12(ỹ-μk(s))T(Σk(s))-1(ỹ-μk(s))].
We now apply the method of transformation from yln(s) to y(s), then,
(4.11)fy(s)(y)=fln⁡y(s)(yln⁡)|det⁡(J1)|,
where lny=(lny1,lny2,…,lnyn), and J1 is the Jacobian matrix. The Jacobian determinant can be computed as
(4.12)det⁡(J1)=|∂ln⁡y1∂y1∂ln⁡y1∂y2⋯∂ln⁡y1∂yn∂ln⁡y2∂y1∂ln⁡y2∂y2⋯∂ln⁡y2∂yn⋮⋮⋮⋮∂ln⁡yn∂y1∂ln⁡yn∂y2⋯∂ln⁡yn∂yn|=|1y10⋯001y2⋯0⋮⋮⋱⋮00⋯1yn|=∏i=1n1yi.
Then, from (4.10) and (4.11), we have the probability density function of y(s) as
(4.13)fy(s)(y)=(2π)-n/2|Σk(s)|-1/2∏i=1n1yi×exp⁡[-12(ln⁡y-μk(s))T(Σk(s))-1(ln⁡y-μk(s))].
Next step is to find the probability density function of x(s) by using method of transformation from y(s) to x(s). Since x(s)=Mky(s), or y(s)=Mk-1x(s), we have
(4.14)fx(s)(x)=fy(s)(Mk-1x)|det⁡(J2)|,
where
(4.15)det⁡(J2)=det⁡(∂∂xv∑j=1n(Mk-1)ujxj)u,v=det⁡((Mk-1)uv)u,v=det⁡(Mk-1).
The result follows from combining (4.13) and (4.14).Now, we are to develop the conditional probability density function of solution processx(t), which is a result parallel to the one-dimensional case in Lemma 3.2.Lemma 4.5.
Under assumptions (i)–(iii) and given thatN(t)=l, and T1=t1, T2=t2,…,Tl=tl, the solution process x(t) has a conditional probability density function as
(4.16)fx(t)∣N(t)=l,t1,…,tl(x)=∫R+n⋯∫R+nfx(t),x(tl),…,x(t1)∣N(t)=l,t1,…,tl(x,xl,…,x1)dxl⋯dx1=∫R+n⋯∫R+n(2π)-n/2|Σl(t)|-1/2|Ml-1|∏i=1n1(Ml-1x)i×exp⁡[-12(ln⁡(Ml-1x)-μl(t))T(Σl(t))-1(ln⁡(Ml-1x)-μl(t))]×∏k=0l-1[(ln⁡(Mk-1(xk+1zk+1))-μk(tk+1))∫0∞(2π)-n/2|Σk(tk+1)|-1/2|Mk-1|∏i=1n1(Mk-1(xk+1/zk+1))i×exp⁡[(ln⁡(Mk-1(xk+1zk+1))-μk(tk+1))-12(ln⁡(Mk-1(xk+1zk+1))-μk(tk+1))T(Σk(tk+1))-1×(ln⁡(Mk-1(xk+1zk+1))-μk(tk+1))(ln⁡(Mk-1(xk+1zk+1))-μk(tk+1))T]g(zk+1)1zk+1ndzk+1]dxl⋯dx1,
where g is the common probability density function of zk, k=1,2,…,l.Proof.
We will apply the result in Lemma4.4 piecewisely to the system (2.9) under the conditions N(t)=l and T1=t1,T2=t2,…,Tl=tl. First, we note that the joint probability density function of (x(t),x(tl),…,x(t1)) can be expressed as
(4.17)fx(t),x(tl),…,x(t1)∣N(t)=l,t1,…,tl(x,xl,…,x1)=fx(t)∣x(tl),N(t)=l,t1,…,tl(x)fx(tl)∣x(tl-1),N(t)=l,t1,…,tl(xl)⋯fx(t1)∣N(t)=l,t1,…,tl(x1),
then, for k=0,1,…,l-1, consider that x(tk+1)=x(tk+1-)zk+1 as a product of two random variables where the first one has the probability density function given in (4.4), and zk+1 is the random jump factor at time tk+1. By the independence of x(tk+1-) and zk+1, we then have
(4.18)fx(tk+1)|x(tk),N(t)=l,t1,…,tl(xk+1)=∫0∞fx(tk+1-)|x(tk),N(t)=l,t1,…,tl(xk+1zk+1)g(zk+1)1zk+1ndzk+1=∫0∞(2π)-n/2|Σk(tk+1)|-1/2|Mk-1|∏i=1n1(Mk-1(xk+1/zk+1))i×exp⁡[-12(ln⁡(Mk-1(xk+1zk+1))-μk(tk+1))T(Σk(tk+1))-1×(ln⁡(Mk-1(xk+1zk+1)))-μk(tk+1)(ln⁡(Mk-1(xk+1zk+1))-μk(tk+1))(ln⁡(Mk-1(xk+1zk+1))-μk(tk+1))T]g(zk+1)1zk+1ndzk+1,
then (4.17) can be written as
(4.19)fx(t),x(tl),…,x(t1)∣N(t)=l,t1,…,tl(x,xl,…,x1)=(2π)-n/2|Σl(t)|-1/2|Ml-1|∏i=1n1(Ml-1x)i×exp⁡[-12(ln⁡(Ml-1x)-μl(t))T(Σl(t))-1(ln⁡(Ml-1x)-μl(t))]×∏k=0l-1[12(ln⁡(Mk-1(xk+1zk+1))-μk(tk+1))T∫0∞(2π)-n/2|Σk(tk+1)|-1/2|Mk-1|∏i=1n1(Mk-1(xk+1/zk+1))i×exp⁡[-12(ln⁡(Mk-1(xk+1zk+1))-μk(tk+1))T(Σk(tk+1))-1×(ln⁡(Mk-1(xk+1zk+1))-μk(tk+1))12(ln⁡(Mk-1(xk+1zk+1))-μk(tk+1))T]g(zk+1)1zk+1ndzk+1∫0∞].
The conditional probability density function of x(t) given in (4.16) is obtained by integrating (4.19) with respect to x1,x2,…, and xl.Now, we can derive the unconditional probability distribution of the solution process of then-dimensional system given in the following proposition.Proposition 4.6.
Under assumptions (i)–(iii), the probability density function of the solution processx(t) to the n-dimensional system (2.9) is given by
(4.20)fx(t)(x)=∑l=0∞[(ln⁡(Mk-1(xk+1zk+1))-μk(tk+1))∫0t⋯∫0t3∫0t2×[(ln⁡(Mk-1(xk+1zk+1))-μk(tk+1))∫R+n⋯∫R+n(2π)-n/2|Σl(t)|-1/2|Ml-1|∏i=1n1(Ml-1x)i×exp⁡[-12(ln⁡(Ml-1x)-μl(t))T(Σl(t))-1(ln⁡(Ml-1x)-μl(t))]×∏k=0l-1[(ln⁡(Mk-1(xk+1zk+1))-μk(tk+1))∫0∞(2π)-n/2|Σk(tk+1)|-1/2|Mk-1|∏i=1n1(Mk-1(xk+1/zk+1))i×exp⁡[(ln⁡(Mk-1(xk+1zk+1))-μk(tk+1))-12(ln⁡(Mk-1(xk+1zk+1))-μk(tk+1))T(Σk(tk+1))-1×(ln⁡(Mk-1(xk+1zk+1))-μk(tk+1))(ln⁡(Mk-1(xk+1zk+1))-μk(tk+1))T]×g(zk+1)1zk+1ndzk+1(ln⁡(Mk-1(xk+1zk+1))-μk(tk+1))∏i=1n]dxl⋯dx1∏i=1n1(Ml-1x)i]×∏i=1lλ(ti)dt1dt2⋯dtlexp⁡[-∫0tλ(u)du]].Proof.
The proof follows the argument in Proposition3.3 by incorporating the random jumps.Remark 4.7.
It is obvious that the result in Proposition4.6 yields the one-dimensional result as a special case. As a result of this, the proof for Proposition 4.6 is considered to be an alternative proof of the one-dimensional result in Proposition 3.3.


## 5. Probability Distribution of an Ornstein-Uhlenbeck Model with Jumps

In this section, we will derive the probability distribution of an Ornstein-Uhlenbeck model with jumps described by system (2.12). To obtain the desired result, we need the following lemma which gives the probability distribution of an Ornstein-Uhlenbeck process which is the continuous dynamic between jumps in system (2.12).Lemma 5.1 (see [18, 29]).
The solution processx(t) of the Ornstein-Uhlenbeck equation
(5.1)dx(t)=Ax(t)dt+Cdw¯(t),t≥t0,x(t0)=x0
follows a multivariate normal distribution with mean μ(t) and covariance matrix Σ(t), where
(5.2)μ(t)=exp⁡A(t-t0)x0,Σ(t)=∫t0teA(t-u)VeAT(t-u)du,
and V=CCT.Suppose that the number of jumps and the jump times are given, then, by applying the above lemma piecewisely, the following result gives the conditional probability density function ofx(t).Lemma 5.2.
Under the conditionsN(t)=l, and T1=t1, T2=t2,…,Tl=tl, the solution process x(t) for system (2.12) has a conditional probability density function as
(5.3)fx(t)∣N(t)=l,t1,…,tl(x)=∫Rn⋯∫Rnfx(t),x(tl),…,x(t1)∣N(t)=l,t1,…,tl(x,xl,…,x1)dxl⋯dx1=∫Rn⋯∫Rn(2π)-n/2|Σl(t)|-1/2exp⁡[-12(x-μl(t))T(Σl(t))-1(x-μl(t))]×∏k=0l-1[(xk+1zk+1-μk(tk+1))T∫0∞(2π)-n/2|Σk(tk+1)|-1/2×exp⁡[-12(xk+1zk+1-μk(tk+1))T(Σk(tk+1))-1(xk+1zk+1-μk(tk+1))]×g(zk+1)1zk+1ndzk+1(xk+1zk+1-μk(tk+1))T]dxl⋯dx1,
where g is the common probability density function of zk, k=1,2,…,l.Proof.
As we noted before that the joint probability density function of(x(t),x(tl),…,x(t1)), under the conditions N(t)=l and T1=t1, T2=t2,…,Tl=tl, can be expressed as
(5.4)fx(t),x(tl),…,x(t1)∣N(t)=l,t1,…,tl(x,xl,…,x1)=fx(t)∣x(tl),N(t)=l,t1,…,tl(x)fx(tl)∣x(tl-1),N(t)=l,t1,…,tl(xl)⋯fx(t1)∣N(t)=l,t1,…,tl(x1),
then, by applying the result in Lemma 5.1 on each interval between jumps [tk,tk+1), we have, for s∈[tk,tk+1),
(5.5)fx(s)∣x(tk),N(t)=l,t1,…,tl(x)=ϕ(x;μk(s),Σk(s))=(2π)-n/2|Σk(s)|-1/2exp⁡[-12(x-μk(s))T(Σk(s))-1(x-μk(s))],
where
(5.6)μk(s)=exp⁡Ak(s-tk)xk,Σk(s)=∫tkseAk(s-u)VkeAkT(s-u)du,
and Vk=CkCkT. Then, for k=0,1,…,l-1, consider that x(tk+1)=x(tk+1-)zk+1 as a product of two random variables where the first one has the probability density function given in (5.5), and zk+1 is the random jump factor at time tk+1. By the independence of x(tk+1-) and zk+1, then
(5.7)fx(tk+1)∣N(t)=l,t1,…,tl(xk+1)=∫0∞fx(tk+1-)∣N(t)=l,t1,…,tl(xk+1zk+1)g(zk+1)1zk+1ndzk+1=∫0∞(2π)-n/2|Σk(tk+1)|-1/2×exp⁡[-12(xk+1zk+1-μk(tk+1))T(Σk(tk+1))-1(xk+1zk+1-μk(tk+1))]×g(zk+1)1zk+1ndzk+1.
Then the conditional joint probability density function (5.4) can be written as
(5.8)fx(t),x(tl),…,x(t1)∣N(t)=l,t1,…,tl(x,xl,…,x1)=(2π)-n/2|Σl(t)|-1/2exp⁡[-12(x-μl(t))T(Σl(t))-1(x-μl(t))]×∏k=0l-1[(xk+1zk+1-μk(tk+1))T∫0∞(2π)-n/2|Σk(tk+1)|-1/2×exp⁡[-12(xk+1zk+1-μk(tk+1))T(Σk(tk+1))-1(xk+1zk+1-μk(tk+1))]×g(zk+1)1zk+1ndzk+1(xk+1zk+1-μk(tk+1))T].
The conditional probability density function of x(t) given in (5.3) is obtained by integrating (5.8) with respect to x1,x2,…, and xl.Finally, following the same argument in Proposition3.3, we present the unconditional probability distribution of the solution process of the Ornstein-Uhlenbeck models with jumps given in the following proposition.Proposition 5.3.
The probability density function of the solution processx(t) of system (2.12) is given by
(5.9)fx(t)(x)=∑l=0∞[∫0t⋯∫0t3∫0t2[(xk+1zk+1-μk(tk+1))T∫Rn⋯∫Rn(2π)-n/2|Σl(t)|-1/2×exp⁡[-12(x-μl(t))T(Σl(t))-1(x-μl(t))]×∏k=0l-1[∫0∞(2π)-n/2|Σk(tk+1)|-1/2×exp⁡[-12(xk+1zk+1-μk(tk+1))T(Σk(tk+1))-1(xk+1zk+1-μk(tk+1))]×g(zk+1)1zk+1ndzk+1(xk+1zk+1-μk(tk+1))]dxl⋯dx1]∏i=1lλ(ti)dt1dt2⋯dtl×exp⁡[-∫0tλ(u)du](xk+1zk+1-μk(tk+1))T].


## 6. Concluding Remarks

In this work, the closed-form solutions of general linear non-homogeneous stochastic hybrid systems are obtained. The methods of finding probability density functions of closed-form solutions are initiated for the linear homogeneous system, and for system with drift and additive noise (Ornstein-Uhlenbeck system). This approach provides a procedure of finding the probability density functions without solving the Fokker-Planck equations. In fact, the Fokker-Planck equation corresponding to system (2.9) has state-dependent coefficients. For example, for n=1 the Fokker-Planck equation corresponding to system (2.9) even in the absence of discrete time interventions is given by(6.1)∂f∂t=-∂∂x[Axf]+∂2∂x2[B2x2f]=(2B2-A)f+(4B2-A)x∂∂x+B2x2∂2∂x2.
As a result of this, equations of this type are not easily solvable in closed-form solutions. In future, we attempt to find the probability distribution of the solution of the general linear non-homogeneous system. In addition, by employing nonlinear transformation, we hope to develop probability distributions for nonlinear systems.



---

*Source: 720614-2011-08-28.xml*