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

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.

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 Section 2, 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 𝑛-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.

Let 𝑥(𝑡) be a real 𝑛-dimensional process. 𝐴_𝑘 and 𝐵^𝑗_𝑘 are 𝑛×𝑛 matrices for any 𝑘∈ℕ∪{0} and 𝑗=1,2,…,𝑞. Let 𝐶^𝑟_𝑘 be 𝑛-dimensional vectors for any 𝑘∈ℕ∪{0} and 𝑟=1,2,…,𝑝. Let the continuous dynamic of the process 𝑥(𝑡) be determined by the following system of stochastic differential equations:𝑑𝑥(𝑡)=𝐴_𝑁(𝑡)𝑥(𝑡)𝑑𝑡+^𝑞_𝑗=1𝐵^𝑗_𝑁(𝑡)𝑥(𝑡)𝑑𝑤_𝑗(𝑡)+^𝑝_𝑟=1𝐶^𝑟_𝑁(𝑡)𝑑𝑤_𝑟(𝑡),𝑡≥𝑡₀,𝑥𝑡₀=𝑥₀>0, (2.1) where 𝑤(𝑡)=(𝑤₁(𝑡),…,𝑤_𝑞(𝑡)) and 𝑤(𝑡)=(𝑤₁(𝑡),…,𝑤_𝑝(𝑡)) are independent 𝑞-dimensional and 𝑝-dimensional standard Wiener processes, and 𝑁(𝑡) is a non-homogeneous Poisson process with intensity 𝜆(𝑡). Here, we denote 𝑥=(𝑥¹,𝑥²,…,𝑥^𝑛)>0 as 𝑥^𝑖>0 for all 𝑖=1,2,…,𝑛. When 𝐶^𝑟_𝑘=0 for all 𝑘 and 𝑟, system (2.1) reduces to a first-order linear homogeneous system of Itô-Doob type stochastic differential equations, given below 𝑑𝑥(𝑡)=𝐴_𝑁(𝑡)𝑥(𝑡)𝑑𝑡+^𝑞_𝑗=1𝐵^𝑗_𝑁(𝑡)𝑥(𝑡)𝑑𝑤_𝑗(𝑡),𝑡≥𝑡₀,𝑥𝑡₀=𝑥₀>0. (2.2) Here, 𝐴_𝑘 and 𝐵^𝑗_𝑘 are such that solution process 𝑥(𝑡) 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 𝐴, 𝐵^𝑗, and 𝐶^𝑟’s are fixed over time, that is, the solution process on a subinterval between jumps. Under the condition that the matrices 𝐴,𝐵¹,𝐵²,…,𝐵^𝑞 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. Let 𝑥(𝑡)≡𝑥(𝑡,𝑡₀,𝑥₀) be the solution of the following initial value problem (IVP): 𝑑𝑥(𝑡)=𝐴𝑥(𝑡)𝑑𝑡+^𝑞_𝑗=1𝐵^𝑗𝑥(𝑡)𝑑𝑤_𝑗(𝑡)+^𝑝_𝑟=1𝐶^𝑟𝑑𝑤_𝑟(𝑡),𝑡≥𝑡₀,𝑥𝑡₀=𝑥₀, (2.3) then the 𝑥(𝑡) is expressed by 𝑥(𝑡)≡𝑥𝑡,𝑡₀,𝑥₀1=exp𝐴−2^𝑞_𝑗=1𝐵^𝑗2𝑡−𝑡₀+^𝑞_𝑗=1𝐵^𝑗𝑤_𝑗(𝑡)−𝑤_𝑗𝑡₀𝑥₀+^𝑝_𝑟=1^𝑡_𝑡01exp𝐴−2^𝑞_𝑗=1𝐵^𝑗2(𝑡−𝑠)+^𝑞_𝑗=1𝐵^𝑗𝑤_𝑗(𝑡)−𝑤_𝑗(𝑠)𝐶^𝑟𝑑𝑤_𝑟(𝑠) (2.4) provided that 𝐴𝐵^𝑗=𝐵^𝑗𝐴 and 𝐵^𝑗𝐵^𝑗′=𝐵^𝑗′𝐵^𝑗 for all 𝑗,𝑗′=1,2,…,𝑞.

Let 𝑥(𝑡,𝑇_𝑘,𝑥_𝑘) be the solution to system (2.3) with 𝑡₀∶=𝑇_𝑘, 𝑥₀∶=𝑥_𝑘, 𝐴∶=𝐴_𝑘, 𝐵^𝑗∶=𝐵^𝑗_𝑘, and 𝐶^𝑟∶=𝐶^𝑟_𝑘. Now, we consider the following system of two interconnected stochastic dynamics:𝑑𝑥(𝑡)=𝐴_𝑘−1𝑥(𝑡)𝑑𝑡+^𝑞_𝑘=1𝐵^𝑗_𝑘−1𝑥(𝑡)𝑑𝑤_𝑗(𝑡)+^𝑝_𝑟=1𝐶^𝑟_𝑘−1𝑑𝑤_𝑟(𝑡),𝑇_𝑘−1≤𝑡<𝑇_𝑘,𝑥𝑇_𝑘−1=𝑥_𝑘−1,𝑥_𝑘=𝑧_𝑘𝑥𝑇⁻_𝑘,𝑇_𝑘−1,𝑥_𝑘−1, (2.5) where 𝑧_𝑘, 𝑘=1,2,3,…, are iid positive random variables with 𝑧₀=1, and 𝑥(𝑇⁻_𝑘,𝑇_𝑘−1,𝑥_𝑘−1)=lim_𝑡→𝑇⁻_𝑘𝑥(𝑡,𝑇_𝑘−1,𝑥_𝑘−1). Here, we assume that 𝑁(𝑡), 𝑤(𝑡), 𝑤(𝑡), and 𝑧_𝑘 are independent.

By applying Lemma 2.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. If 𝐴_𝑘𝐵^𝑗_𝑘=𝐵^𝑗_𝑘𝐴_𝑘 and 𝐵^𝑗_𝑘𝐵_𝑘^𝑗′=𝐵_𝑘^𝑗′𝐵^𝑗_𝑖 for all 𝑘∈ℕ∪{0} and 𝑗,𝑗′=1,2,…,𝑞, then the solution to the system (2.5) is given by 𝑥(𝑡)=^𝑁(𝑡)_𝑘=1𝑧_𝑘exp𝐴_𝑁(𝑡)−12^𝑞_j=1𝐵^𝑗_𝑁(𝑡)²𝑡−𝑇_𝑁(𝑡)+^𝑞_𝑗=1𝐵^𝑗_𝑁(𝑡)𝑤_𝑗(𝑡)−𝑤_𝑗𝑇_𝑁(𝑡)×^𝑁(𝑡)_𝑘=1exp𝐴_𝑘−1−12^𝑞_𝑗=1𝐵^𝑗_𝑘−1²𝑇_𝑘−𝑇_𝑘−1+^𝑞_𝑗=1𝐵^𝑗_𝑘−1𝑤_𝑗𝑇_𝑘−𝑤_𝑗𝑇_𝑘−1𝑥₀+^𝑁(𝑡)_𝑖=1^𝑁(𝑡)_𝑘=𝑖𝑧_𝑘exp𝐴_𝑁(𝑡)−12^𝑞_𝑗=1𝐵^𝑗_𝑁(𝑡)²𝑡−𝑇_𝑁(𝑡)+^𝑞_𝑗=1𝐵^𝑗_𝑁(𝑡)𝑤_𝑗(𝑡)−𝑤_𝑗𝑇_𝑁(𝑡)×^𝑁(𝑡)_𝑘=𝑖+1exp𝐴_𝑘−1−12^𝑞_𝑗=1𝐵^𝑗_𝑘−1²𝑇_𝑘−𝑇_𝑘−1+^𝑞_𝑗=1𝐵^𝑗_𝑘−1𝑤_𝑗𝑇_𝑘−𝑤_𝑗𝑇_𝑘−1×^𝑝_𝑟=1^T𝑖_𝑇𝑖−1exp𝐴_𝑖−1−12^𝑞_𝑗=1𝐵^𝑗_𝑖−1²𝑇_𝑖+−𝑠_𝑞𝑗=1𝐵^𝑗_𝑖−1𝑤_𝑗𝑇_𝑖−𝑤_𝑗(𝑠)𝐶^𝑟_𝑖−1𝑑𝑤_𝑟(𝑠)+^𝑝_𝑟=1^𝑡_{𝑇𝑁(𝑡)}exp𝐴_𝑁(𝑡)−12^𝑞_𝑗=1𝐵^𝑗_𝑁(𝑡)²(𝑡−𝑠)+^𝑞_𝑗=1𝐵^𝑗_𝑁(𝑡)𝑤_𝑗(𝑡)−𝑤_𝑗(𝑠)𝐶^𝑟_𝑁(𝑡)𝑑𝑤_𝑟(𝑠). (2.6)

Proof. By applying the result of (2.4) on the subintervals [𝑇_𝑘−1,𝑇_𝑘), for 𝑘=1,…,𝑁(𝑡), and [𝑇_𝑁(𝑡),𝑡), we have the solution to system (2.5) as the following piecewise function: 𝑥(𝑠)=𝑥𝑠,𝑡₀,𝑥₀𝑥𝑡₀=𝑥₀,𝑡₀≤𝑠<𝑇₁,𝑥𝑠,𝑇₁,𝑥₁𝑥𝑇₁=𝑥₁,𝑇₁≤𝑠<𝑇₂,⋯𝑥𝑠,𝑇_𝑖,𝑥_𝑖𝑥𝑇_𝑖=𝑥_𝑘,𝑇_𝑘≤𝑠<𝑇_𝑘+1,⋯𝑥𝑠,𝑇_𝑁(𝑡),𝑥_𝑁(𝑡)𝑥𝑇_𝑁(𝑡)=𝑥_𝑁(𝑡),𝑇_𝑁(𝑡)≤𝑠<𝑡, (2.7) where 𝑥(𝑠,𝑇_𝑘,𝑥_𝑘) is a solution process in (2.4) with 𝑡₀∶=𝑇_𝑘, 𝑥₀∶=𝑥_𝑘, 𝐴∶=𝐴_𝑘, 𝐵^𝑗∶=𝐵^𝑗_𝑘, and 𝐶^𝑟∶=𝐶^𝑟_𝑘, then, we have 𝑥(𝑡)=𝑥𝑡,𝑇_𝑁(𝑡),𝑥_𝑁(𝑡)=exp𝐴_𝑁(𝑡)−12^𝑞_𝑗=1𝐵^𝑗_𝑁(𝑡)²𝑡−𝑇_𝑁(𝑡)+^𝑞_𝑗=1𝐵^𝑗_𝑁(𝑡)𝑤_𝑗(𝑡)−𝑤_𝑗𝑇_𝑁(𝑡)𝑥_𝑁(𝑡)+^𝑝_𝑟=1^𝑡_{𝑇𝑁(𝑡)}exp𝐴_𝑁(𝑡)−12^𝑞_𝑗=1𝐵^𝑗_𝑁(𝑡)²(𝑡−𝑠)+^𝑞_𝑗=1𝐵^𝑗_𝑁(𝑡)𝑤_𝑗(𝑡)−𝑤_𝑗(𝑠)𝐶^𝑟_𝑁(𝑡)𝑑𝑤_𝑟(𝑠)=𝑧_𝑁(𝑡)exp𝐴_𝑁(𝑡)−12^𝑞_𝑗=1𝐵^𝑗_𝑁(𝑡)²𝑡−𝑇_𝑁(𝑡)+^𝑞_𝑗=1𝐵^𝑗_𝑁(𝑡)𝑤_𝑗(𝑡)−𝑤_𝑗𝑇_𝑁(𝑡)×exp𝐴_{𝑁(𝑡)−1}−12^𝑞_𝑗=1𝐵^𝑗_{𝑁(𝑡)−1}²𝑇_𝑁(𝑡)−𝑇_{𝑁(𝑡)−1}+^𝑞_𝑗=1𝐵^𝑗_{𝑁(𝑡)−1}𝑤_𝑗𝑇_𝑁(𝑡)−𝑤_𝑗𝑇_{𝑁(𝑡)−1}𝑥_{𝑁(𝑡)−1}+^𝑝_𝑟=1^{𝑇𝑁(𝑡)}_{𝑇𝑁(𝑡)−1}exp𝐴_{𝑁(𝑡)−1}−12^𝑞_𝑗=1𝐵^𝑗_{𝑁(𝑡)−1}²𝑇_𝑁(𝑡)+−𝑠_𝑞𝑗=1𝐵^𝑗_{𝑁(𝑡)−1}𝑤_𝑗𝑇_𝑁(𝑡)−𝑤_𝑗(𝑠)𝐶^𝑟_{𝑁(𝑡)−1}𝑑𝑤_𝑟(𝑠)+^𝑝_𝑟=1^𝑡_{𝑇𝑁(𝑡)}exp𝐴_𝑁(𝑡)−12^𝑞_𝑗=1𝐵^𝑗_𝑁(𝑡)²(𝑡−𝑠)+^𝑞_𝑗=1𝐵^𝑗_𝑁(𝑡)𝑤_𝑗(𝑡)−𝑤_𝑗(𝑠)𝐶^𝑟_𝑁(𝑡)𝑑𝑤_𝑟(𝑠). (2.8) Next, substitute 𝑥_{𝑁(𝑡)−1}=𝑧_{𝑁(𝑡)−1}𝑥(𝑇⁻_{𝑁(𝑡)−1},𝑇_{𝑁(𝑡)−2},𝑥_{𝑁(𝑡)−2})=𝑧_{𝑁(𝑡)−1}𝑥(𝑇_{𝑁(𝑡)−1},𝑇_{𝑁(𝑡)−2},𝑥_{𝑁(𝑡)−2}). The term 𝑥(𝑇⁻_{𝑁(𝑡)−1},𝑇_{𝑁(𝑡)−2},𝑥_{𝑁(𝑡)−2}) can be replaced by 𝑥(𝑇_{𝑁(𝑡)−1},𝑇_{𝑁(𝑡)−2},𝑥_{𝑁(𝑡)−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 Proposition 2.2. When 𝐶^𝑟_𝑘=0 for all 𝑘 and 𝑟, system (2.5) reduces to the following first-order linear homogeneous system of Itô-Doob type stochastic differential equations with jumps:𝑑𝑥(𝑡)=𝐴_𝑘−1𝑥(𝑡)𝑑𝑡+^𝑞_𝑘=1𝐵^𝑗_𝑘−1𝑥(𝑡)𝑑𝑤_𝑗(𝑡),𝑇_𝑘−1≤𝑡<𝑇_𝑘,𝑥𝑇_𝑘−1=𝑥_𝑘−1,𝑥_𝑘=𝑧_𝑘𝑥𝑇⁻_𝑘,𝑇_𝑘−1,𝑥_𝑘−1. (2.9)

The solution of the above system is given in the following corollary. The result follows from Proposition 2.2 by letting 𝐶^𝑟_𝑘 be zero for all 𝑘 and 𝑟.

Corollary 2.3. If 𝐴_𝑘𝐵^𝑗_𝑘=𝐵^𝑗_𝑘𝐴_𝑘 and 𝐵^𝑗_𝑘𝐵_𝑘^𝑗′=𝐵_𝑘^𝑗′𝐵^𝑗_𝑖 for all 𝑘∈ℕ∪{0} and 𝑗,𝑗′=1,2,…,𝑞, then the solution of system (2.9) is given by 𝑥(𝑡)=^𝑁(𝑡)_𝑘=1𝑧_𝑘exp𝐴_𝑁(𝑡)−12^𝑞_𝑗=1𝐵^𝑗_𝑁(𝑡)²𝑡−𝑇_𝑁(𝑡)+^𝑞_𝑗=1𝐵^𝑗_𝑁(𝑡)𝑤_𝑗(𝑡)−𝑤_𝑗𝑇_𝑁(𝑡)×^𝑁(𝑡)_𝑘=1exp𝐴_𝑘−1−12^𝑞_𝑗=1𝐵^𝑗_𝑘−1²𝑇_𝑘−𝑇_𝑘−1+^𝑞_𝑗=1𝐵^𝑗_𝑘−1𝑤_𝑗𝑇_𝑘−𝑤_𝑗𝑇_𝑘−1𝑥₀. (2.10)

In the case when 𝐵^𝑗_𝑘 are zeros for all 𝑗 and 𝑘, system (2.5) becomes a linear system with additive noise. The continuous dynamics between jumps are now governed by Ornstein-Uhlenbeck equations as follows:𝑑𝑥(𝑡)=𝐴_𝑘−1𝑥(𝑡)𝑑𝑡+^𝑝_𝑟=1𝐶^𝑟_𝑘−1𝑑𝑤_𝑟(𝑡),𝑇_𝑘−1≤𝑡<𝑇_𝑘,𝑥𝑇_𝑘−1=𝑥_𝑘−1,𝑥_𝑘=𝑧_𝑘𝑥𝑇⁻_𝑘,𝑇_𝑘−1,𝑥_𝑘−1. (2.11) Denote 𝐶_𝑘=(𝐶¹_𝑘,𝐶²_𝑘,…,𝐶^𝑝_𝑘) for all 𝑘, then the above system can be rewritten as𝑑𝑥(𝑡)=𝐴_𝑘−1𝑥(𝑡)𝑑𝑡+𝐶_𝑘−1𝑑𝑤(𝑡),𝑇_𝑘−1≤𝑡<𝑇_𝑘,𝑥𝑇_𝑘−1=𝑥_𝑘−1,𝑥_𝑘=𝑧_𝑘𝑥𝑇⁻_𝑘,𝑇_𝑘−1,𝑥_𝑘−1, (2.12) where 𝐶_𝑘's are 𝑛×𝑝 matrices, and 𝑤(𝑡) is a 𝑝-dimensional standard Wiener process.

Corollary 2.4. The solution of system (2.11) is given by 𝑥(𝑡)=^𝑁(𝑡)_𝑘=1𝑧_𝑘exp𝐴_𝑁(𝑡)𝑡−𝑇_𝑁(𝑡)^𝑁(𝑡)_𝑘=1exp𝐴_𝑘−1𝑇_𝑘−𝑇_𝑘−1𝑥₀+^𝑁(𝑡)_𝑖=1^𝑁(𝑡)_𝑘=𝑖𝑧_𝑘exp𝐴_𝑁(𝑡)𝑡−𝑇_𝑁(𝑡)^𝑁(𝑡)_𝑘=𝑖+1exp𝐴_𝑘−1𝑇_𝑘−𝑇_𝑘−1×^𝑝_𝑟=1^𝑇𝑖_𝑇𝑖−1exp𝐴_𝑖−1𝑇_𝑖−𝑠𝐶^𝑟_𝑖−1𝑑𝑤_𝑟(𝑠+_{𝑝𝑟=1𝑡𝑇𝑁(𝑡)})exp𝐴_𝑁(𝑡)(𝑡−𝑠)𝐶^𝑟_𝑁(𝑡)𝑑𝑤_𝑟(𝑠). (2.13)

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, 𝑥(𝑡) takes values in ℝ₊, and in this case, 𝐴_𝑘 and 𝐵_𝑘 are scalars for all 𝑘 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 𝑁(𝑡). 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 process 𝑁(𝑡) with intensity 𝜆(𝑡), the joint density function of the jump times 𝑇₁,𝑇₂,…,𝑇_𝑙 conditioned on 𝑁(𝑡)=𝑙 is given by 𝑓_{𝑇1,𝑇2,…,𝑇𝑙∣𝑁(𝑡)=𝑙}𝑡₁,𝑡₂,…,𝑡_𝑙=𝑙!∏^𝑙_𝑘=1𝜆𝑡_𝑘∫^𝑡₀𝜆(𝑢)𝑑𝑢^𝑙. (3.1)

Next lemma gives the conditional probability density function of 𝑥(𝑡) given the number of jumps and the jump times.

Lemma 3.2. Given that 𝑁(𝑡)=𝑙, and 𝑇₁=𝑡₁,…,𝑇_𝑙=𝑡_𝑙, 𝑥(𝑡) has a probability density function as 𝑓_{𝑥(𝑡)∣𝑁(𝑡)=𝑙,𝑡1,…,𝑡l}(𝑥)=1𝑥^∞_−∞ℎ^∗𝑙(ln𝑥−𝑠)𝜙𝑠;𝜇,𝜎𝑑𝑠,𝑥>0, (3.2) where ℎ^∗𝑙 is the 𝑙th convolution of the common probability density function ℎ of ln𝑧_𝑘, for 𝑘=1,2,…,𝑙, and 𝜙(⋅;𝜇,𝜎) denotes the normal density function with mean 𝜇 and variance 𝜎.

Proof. Given that 𝑁(𝑡)=𝑙, and 𝑇₁=𝑡₁,…, 𝑇_𝑙=𝑡_𝑙, from (2.10), we have ln𝑥(𝑡)=^𝑙_𝑘=1ln𝑧_𝑘+ln𝑥₀+𝐴_𝑙−𝐵²_𝑙2𝑡−𝑡_𝑙+𝐵_𝑙𝑤(𝑡)−𝑤𝑡_𝑙+^𝑙_𝑘=1𝐴_𝑘−1−𝐵²_𝑘−12𝑡_𝑘−𝑡_𝑘−1+𝐵_𝑘−1𝑤𝑡_𝑘−𝑤𝑡_𝑘−1=𝑉+𝑆, (3.3) where we denote 𝑉=∑^𝑙_𝑘=1ln𝑧_𝑘 and 𝑆 as the sum of last three terms in (3.3).
Since ℎ is the common probability density function of ln𝑧_𝑘, 𝑉 as the sum of 𝑙 iid random variables has the probability density function as the 𝑙th convolution ℎ^∗𝑙(𝑣). We further note that 𝑆 is the sum of 𝑙+1 independent normal variables due to the independent increment property of Wiener process. Then 𝑆 is normally distributed with mean and variance as𝜇≡𝜇𝑡,𝑡₁,…,𝑡_𝑙=ln𝑥₀+𝐴_𝑙−𝐵²_𝑙2𝑡−𝑡_𝑙+^𝑙_𝑘=1𝐴_𝑘−1−𝐵²_𝑘−12𝑡_𝑘−𝑡_𝑘−1,𝜎≡𝜎𝑡,𝑡₁,…,𝑡_𝑙=𝐵²_𝑙𝑡−𝑡_𝑙+^𝑙_𝑘=1𝐵²_𝑘−1𝑡_𝑘−𝑡_𝑘−1. (3.4) Since 𝑧_𝑘 and 𝑤(𝑡) are independent, then 𝑉 and 𝑆 are also independent. Using transformation method [26] on ln𝑥(𝑡)=𝑉+𝑆, we can obtain the conditional probability density function of ln𝑥(𝑡) as 𝑓_{ln𝑥(𝑡)∣𝑁(𝑡)=𝑙,𝑡1,…,𝑡𝑙}(𝑥)=^∞_−∞ℎ^∗𝑙(𝑥−𝑠)𝜙𝑠;𝜇,𝜎𝑑𝑠. (3.5) If follows that 𝑓_{𝑥(𝑡)∣𝑁(𝑡)=𝑙,𝑡1,…,𝑡𝑙}(𝑥)=𝑓_{ln𝑥(𝑡)∣𝑁(𝑡)=𝑙,𝑡1,…,𝑡𝑙}()1ln𝑥𝑥=1𝑥^∞_−∞ℎ^∗𝑙(ln𝑥−𝑠)𝜙𝑠;𝜇,𝜎𝑑𝑠,𝑥>0. (3.6)

Having obtained the conditional probability density function of 𝑥(𝑡), 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 process 𝑥(𝑡) to the system (2.5) is given by 𝑓_𝑥(𝑡)(𝑥)=^∞_𝑙=0^𝑡₀⋯^𝑡3₀^𝑡2₀^∞_−∞ℎ^∗𝑙(ln𝑥−𝑠)𝜙𝑠;𝜇,𝜎𝑑𝑠^𝑙_𝑘=1𝜆𝑡_𝑘𝑑𝑡₁𝑑𝑡₂⋯𝑑𝑡_𝑙×1𝑥exp−^𝑡₀𝜆(𝑢)𝑑𝑢. (3.7)

Proof. From (3.1) and (3.2), the joint density function of 𝑥(𝑡), 𝑇₁,…, 𝑇_𝑙 given the condition 𝑁(𝑡)=𝑙 is given by 𝑓_{𝑥(𝑡),𝑇1,…,𝑇𝑙∣𝑁(𝑡)=𝑙}𝑥,𝑡₁,…,𝑡_𝑙=𝑓_{𝑥(𝑡)∣𝑁(𝑡)=𝑙,𝑡1,…,𝑡𝑙}(𝑥)𝑓_{𝑇1,𝑇2,…,𝑇𝑙∣𝑁(𝑡)=𝑙}𝑡₁,𝑡₂,…,𝑡_𝑙=1𝑥^∞_−∞ℎ^∗𝑙(ln𝑥−𝑠)𝜙𝑠;𝜇,𝜎𝑑𝑠𝑙!∏^𝑙_𝑘=1𝜆𝑡_𝑘∫^𝑡₀𝜆(𝑢)𝑑𝑢^𝑙. (3.8) Integrating with respect to 𝑡₁,𝑡₂,…, and 𝑡_𝑙 gives the probability density function of 𝑥(𝑡) given that 𝑁(𝑡)=𝑙 as 𝑓_{𝑥(𝑡)∣𝑁(𝑡)=𝑙}(𝑥)=^𝑡₀⋯^𝑡3₀^𝑡2₀^∞_−∞ℎ^∗𝑙(ln𝑥−𝑠)𝜙𝑠;𝜇,𝜎𝑑𝑠^𝑙_𝑘=1𝜆𝑡_𝑘𝑑𝑡₁𝑑𝑡₂⋯𝑑𝑡_𝑙𝑙!𝑥∫^𝑡₀𝜆(𝑢)𝑑𝑢^𝑙. (3.9) By multiplying the above density by the probability of 𝑙 jumps, namely, 𝑃𝑟(𝑁(𝑡)=𝑙)=exp−^𝑡₀𝜆(𝑢)∫𝑑𝑢^𝑡0𝜆(𝑢)𝑑𝑢^𝑙,𝑙! (3.10) and then taking the summation over 𝑙, the marginal probability density function of 𝑥(𝑡) 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 process 𝑥(𝑡) for the 𝑛-dimensional stochastic system (2.9) under the following assumptions. (i)The drift and diffusion coefficients, 𝐴_𝑘 and 𝐵^𝑗_𝑘 for all 𝑘∈ℕ∪{0} and 𝑗=1,2,…,𝑞, are diagonalizable. (ii)The coefficients in each regime pairwise commute, that is, 𝐴_𝑘𝐵^𝑗_𝑘=𝐵^𝑗_𝑘𝐴_𝑘 and 𝐵^𝑗_𝑘𝐵^𝑙_𝑘=𝐵^𝑙_𝑘𝐵^𝑗_𝑘 for all 𝑘∈ℕ∪{0} and 𝑗,𝑙=1,2,…,𝑞. (iii)𝐴_𝑘∈𝒞 for 𝑘∈ℕ∪{0}, where 𝒞 denotes the set of 𝑛×𝑛 diagonalizable matrices whose eigenvector matrix, 𝑀, is such that M⁻¹𝑥>0 for any 𝑥>0.

We first consider the stochastic system on the interval between jumps. Given that 𝑁(𝑡)=𝑙 and 𝑇₁=𝑡₁, 𝑇₂=𝑡₂,…,𝑇_𝑙=𝑡_𝑙, consider the following SDE on [𝑡_𝑘,𝑡_𝑘+1), for some 𝑘=0,1,…,𝑙:𝑑𝑥(𝑠)=𝐴_𝑘𝑥(𝑠)𝑑𝑡+^𝑞_𝑗=1𝐵^𝑗_𝑘𝑥(𝑠)𝑑𝑤_𝑗(𝑠),𝑡_𝑘≤𝑠<𝑡_𝑘+1,𝑥𝑡_𝑘=𝑥_𝑘. (4.1) 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 𝑥 is said to have a 𝑛-dimensional multivariate normal distribution with mean and covariance matrix, 𝜇_𝑥 and Σ_𝑥, if its probability density function is given by𝑓(𝑥)=(2𝜋)^−𝑛/2Σ_𝑥^−1/2−1exp2𝑥−𝜇_𝑥^𝑇Σ_𝑥⁻¹𝑥−𝜇_𝑥. (4.2) 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]). If 𝑥 has a multivariate normal distribution with mean and covariance matrix, 𝜇_𝑥 and Σ_𝑥, then 𝑦=𝑃𝑥+𝑐 as a linear transformation of 𝑥 follows also multivariate normal distribution with mean and covariance matrix, 𝜇_𝑦=𝑃𝜇_𝑥+𝑐 and Σ_𝑦=𝑃Σ_𝑥𝑃^𝑇.

Proof. Since 𝑦_𝑖=∑^𝑛_𝑗=1𝑃_𝑖𝑗𝑥_𝑗+𝑐_𝑖 and every linear combination of normal random variables is still normal, then 𝑦 follows a multivariate normal distribution. By the linearity of expectation we have 𝜇_𝑦=𝐸𝑦=𝐸𝑃𝜇_𝑥+𝑐=𝑃𝐸[𝑥]+𝑐=𝑃𝜇_𝑥Σ+𝑐,_𝑦=𝐸𝑦−𝜇_𝑦𝑦−𝜇_𝑦^𝑇=𝐸𝑃𝑥−𝜇_𝑥𝑃𝑥−𝜇_𝑥^𝑇=𝐸𝑃𝑥−𝜇_𝑥𝑥−𝜇_𝑥^𝑇𝑃^𝑇=𝑃𝐸𝑥−𝜇_𝑥𝑥−𝜇_𝑥^𝑇𝑃^𝑇=𝑃Σ_𝑥𝑃^𝑇. (4.3)

To find the probability density function of 𝑥(𝑠) 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, 𝐴_𝑘s and 𝐵^𝑗_𝑘s have the same eigenvector matrix, denoted by 𝑀_𝑘. Moreover, 𝐴_𝑘≡𝑀_𝑘⁻¹𝐴_𝑘𝑀_𝑘=diag(𝑎¹_𝑘,…,𝑎^𝑛_𝑘) and 𝐵^𝑗_𝑘≡𝑀_𝑘⁻¹𝐵^𝑗_𝑘𝑀_𝑘=diag(𝑏̃¹_𝑘,𝑗,…,𝑏̃^𝑛_𝑘,𝑗) where {𝑎¹_𝑘,…,𝑎^𝑛_𝑘} and {𝑏̃¹_𝑘,𝑗,…,𝑏̃^𝑛_𝑘,𝑗} are the sets of eigenvalues of 𝐴_𝑘 and 𝐵^𝑗_𝑘, respectively, for all 𝑘,𝑗. Next, we define a linear transformation 𝑦=𝑀_𝑘⁻¹𝑥, then 𝑥=𝑀_𝑘𝑦. Define ln𝑦=(ln𝑦¹,ln𝑦²,…,ln𝑦^𝑛) [28].

The below proposition gives the probability distribution of solution process 𝑥(𝑡) on the interval [𝑡_𝑘,𝑡_𝑘+1) over which the coefficients are constant.

Lemma 4.4. Under assumptions (i)–(iii), the process 𝑥(𝑠) satisfying the SDE (4.1) has a probability density function given by 𝑓_𝑥(𝑠)(𝑥)=(2𝜋)^−𝑛/2Σ_𝑘(𝑠)^−1/2𝑀_𝑘^−1𝑛_𝑖=11𝑀_𝑘⁻¹𝑥^𝑖−1×exp2ln𝑀_𝑘⁻¹𝑥−𝜇_𝑘(𝑠)^𝑇Σ_𝑘(𝑠)⁻¹ln𝑀_𝑘⁻¹𝑥−𝜇_𝑘(𝑠),for𝑠∈t_k,t_k+1, (4.4) where 𝜇_𝑘(𝑠)=ln𝑦¹𝑡_𝑘+𝑎¹_𝑘−12^𝑞_𝑗=1𝑏̃¹_𝑘,𝑗²𝑠−𝑡_𝑘ln𝑦²𝑡_𝑘+𝑎²_𝑘−12^𝑞_𝑗=1𝑏̃²_𝑘,𝑗²𝑠−𝑡_𝑘⋮ln𝑦^𝑛𝑡_𝑘+𝑎^𝑛_𝑘−12^𝑞_𝑗=1𝑏̃^𝑛_𝑘,𝑗²𝑠−𝑡_𝑘,𝐵^∗_𝑘=𝑏̃¹_𝑘,1𝑏̃¹_𝑘,2⋯𝑏̃¹_𝑘,𝑞𝑏̃²_𝑘,1𝑏̃²_𝑘,2⋯𝑏̃²_𝑘,𝑞𝑏̃⋮⋮⋮⋮_𝑛𝑘,1𝑏̃^𝑛_𝑘,2⋯𝑏̃^𝑛_𝑘,𝑞,Σ_𝑘(𝑠)=𝑠−𝑡_𝑘𝐵^∗_𝑘𝐵^∗_𝑘^𝑇. (4.5)

Proof. For 𝑠∈[𝑡_𝑘,𝑡_𝑘+1), we have 𝑦(𝑠)=𝑀_𝑘⁻¹𝑥(𝑠), and 𝑥(𝑠)=𝑀_𝑘𝑦(𝑠). By multiplying 𝑀_𝑘⁻¹ on both sides of the SDE (4.1), we obtain the transformed SDE as follows: 𝑀_𝑘⁻¹𝑑𝑥(𝑠)=𝑀_𝑘⁻¹𝐴_𝑘𝑥(𝑠)𝑑𝑠+^𝑞_𝑗=1𝑀_𝑘⁻¹𝐵^𝑗_𝑘𝑥(𝑠)𝑑𝑤_𝑗(𝑠),⟹𝑑𝑀_𝑘⁻¹𝑥(𝑠)=𝑀_𝑘⁻¹𝐴_𝑘𝑀_𝑘𝑦(𝑠)𝑑𝑠+^𝑞_𝑗=1𝑀_𝑘⁻¹𝐵^𝑗_𝑘𝑀_𝑘𝑦(𝑠)𝑑𝑤_𝑗(𝑠),⟹𝑑𝑦(𝑠)=𝐴_𝑘𝑦(𝑠)𝑑𝑠+^𝑞_𝑗=1𝐵^𝑗_𝑘𝑦(𝑠)𝑑𝑤_𝑗(𝑠), (4.6) where 𝐴_𝑘 and 𝐵^𝑗_𝑘 are diagonal matrices as defined before.
From the application of Lemma 2.1 with 𝐶^𝑟_𝑘=0, the solution process of the transformed system (4.6) is𝑦(𝑠)1=exp𝐴−2^𝑞_𝑗=1𝐵^𝑗_𝑘²𝑠−𝑡_𝑘+^𝑞_𝑗=1𝐵^𝑗_𝑘𝑤_𝑗(𝑠)−𝑤_𝑗𝑡_𝑘𝑦_𝑘, (4.7) for 𝑠∈[𝑡_𝑘,𝑡_𝑘+1).
It follows from assumption (iii) that 𝑦(𝑠)=𝑀_𝑘⁻¹𝑥(𝑠)>0 since 𝑥(𝑠)>0, then, we can rewrite system (4.7) in the following form:ln𝑦(𝑠)=𝜇_𝑘(𝑠)+𝐵^∗_𝑘𝑤(𝑠)−𝑤𝑡_𝑘, (4.8) where 𝜇_𝑘(𝑠) and 𝐵^∗_𝑘 are defined above.
Since 𝑤(𝑡) is a standard Wiener process, then 𝑤(𝑠)−𝑤(𝑡_𝑘) has a multivariate normal distribution with mean zero and covariance matrix (𝑠−𝑡_𝑘)𝐼_𝑛, where 𝐼_𝑛 is the 𝑛×𝑛 identity matrix. Then, by Theorem 4.3, ln𝑦(𝑠) as a linear transformation of 𝑤(𝑠)−𝑤(𝑡_𝑘) is also multivariate normally distributed with mean 𝜇_𝑘(𝑠) and covariance matrix Σ_𝑘(𝑠), whereΣ_𝑘(𝑠)=𝐵^∗_𝑘𝑠−𝑡_𝑘𝐼_𝑛𝐵^∗_𝑘^𝑇=𝑠−𝑡_𝑘𝐵^∗_𝑘𝐵^∗_𝑘^𝑇, (4.9) and the (u,v)th element of Σ_𝑘(𝑠) is (𝑠−𝑡_𝑘)∑^𝑞_𝑗=1𝑏̃^𝑢_𝑘,𝑗𝑏̃^𝑣_𝑘,𝑗, for 𝑢,𝑣=1,2,…,𝑛. The probability density function of ln𝑦(𝑠) is 𝑓_ln𝑦(𝑠)𝑦=(2𝜋)^−𝑛/2Σ_𝑘(𝑠)^−1/2−1exp2𝑦−𝜇_𝑘(𝑠)^𝑇Σ_𝑘(𝑠)⁻¹𝑦−𝜇_𝑘(𝑠). (4.10) We now apply the method of transformation from 𝑦_ln(𝑠) to 𝑦(𝑠), then, 𝑓_𝑦(𝑠)𝑦=𝑓_ln𝑦(𝑠)𝑦_lndet𝐽₁, (4.11) where ln𝑦=(ln𝑦¹,ln𝑦²,…,ln𝑦^𝑛), and 𝐽₁ is the Jacobian matrix. The Jacobian determinant can be computed as det𝐽₁=𝜕ln𝑦¹𝜕𝑦¹𝜕ln𝑦¹𝜕𝑦²⋯𝜕ln𝑦¹𝜕𝑦^𝑛𝜕ln𝑦²𝜕𝑦¹𝜕ln𝑦²𝜕𝑦²⋯𝜕ln𝑦²𝜕𝑦^𝑛⋮⋮⋮⋮𝜕ln𝑦^𝑛𝜕𝑦¹𝜕ln𝑦^𝑛𝜕𝑦²⋯𝜕ln𝑦^𝑛𝜕𝑦^𝑛=1𝑦¹010⋯0𝑦²1⋯0⋮⋮⋱⋮00⋯𝑦^𝑛=^𝑛_𝑖=11𝑦^𝑖. (4.12) Then, from (4.10) and (4.11), we have the probability density function of 𝑦(𝑠) as 𝑓_𝑦(𝑠)𝑦=(2𝜋)^−𝑛/2Σ_𝑘(𝑠)^𝑛−1/2_𝑖=11𝑦^𝑖−1×exp2ln𝑦−𝜇_𝑘(𝑠)^𝑇Σ_𝑘(𝑠)⁻¹ln𝑦−𝜇_𝑘(𝑠). (4.13) Next step is to find the probability density function of 𝑥(𝑠) by using method of transformation from 𝑦(𝑠) to 𝑥(𝑠). Since 𝑥(𝑠)=𝑀_𝑘𝑦(𝑠), or 𝑦(𝑠)=𝑀_𝑘⁻¹𝑥(𝑠), we have 𝑓_𝑥(𝑠)(𝑥)=𝑓_𝑦(𝑠)𝑀_𝑘⁻¹𝑥det𝐽₂, (4.14) where det𝐽₂𝜕=det𝜕𝑥^𝑣𝑛_𝑗=1𝑀_𝑘^−1𝑢𝑗𝑥^𝑗_𝑢,𝑣=det𝑀_𝑘^−1𝑢𝑣_𝑢,𝑣=det𝑀_𝑘⁻¹. (4.15) The result follows from combining (4.13) and (4.14).

Now, we are to develop the conditional probability density function of solution process 𝑥(𝑡), which is a result parallel to the one-dimensional case in Lemma 3.2.

Lemma 4.5. Under assumptions (i)–(iii) and given that 𝑁(𝑡)=𝑙, and 𝑇₁=𝑡₁, 𝑇₂=𝑡₂,…,𝑇_𝑙=𝑡_𝑙, the solution process 𝑥(𝑡) has a conditional probability density function as 𝑓_{𝑥(𝑡)∣𝑁(𝑡)=𝑙,𝑡1,…,𝑡𝑙}(𝑥)=_ℝ𝑛+⋯_ℝ𝑛+𝑓_{𝑥(𝑡),𝑥(𝑡𝑙),…,𝑥(𝑡1)∣𝑁(𝑡)=𝑙,𝑡1,…,𝑡𝑙}𝑥,𝑥_𝑙,…,𝑥₁𝑑𝑥_𝑙⋯𝑑𝑥₁=_ℝ𝑛+⋯_ℝ𝑛+(2𝜋)^−𝑛/2Σ_𝑙(𝑡)^−1/2𝑀_𝑙^−1𝑛_𝑖=11𝑀_𝑙⁻¹𝑥^𝑖−1×exp2ln𝑀_𝑙⁻¹𝑥−𝜇_𝑙(𝑡)^𝑇Σ_𝑙(𝑡)⁻¹ln𝑀_𝑙⁻¹𝑥−𝜇_𝑙(𝑡)×^𝑙−1_𝑘=0^∞₀(2𝜋)^−𝑛/2Σ_𝑘𝑡_𝑘+1^−1/2𝑀_𝑘^−1𝑛_𝑖=11𝑀_𝑘⁻¹𝑥_𝑘+1/𝑧_𝑘+1^𝑖−1×exp2ln𝑀_𝑘⁻¹𝑥_𝑘+1𝑧_𝑘+1−𝜇_𝑘𝑡_𝑘+1^𝑇Σ_𝑘𝑡_𝑘+1⁻¹×ln𝑀_𝑘⁻¹𝑥_𝑘+1𝑧_𝑘+1−𝜇_𝑘𝑡_𝑘+1𝑔𝑧_𝑘+11𝑧^𝑛_𝑘+1𝑑𝑧_𝑘+1𝑑𝑥_𝑙⋯𝑑𝑥₁, (4.16) where 𝑔 is the common probability density function of 𝑧_𝑘, 𝑘=1,2,…,𝑙.

Proof. We will apply the result in Lemma 4.4 piecewisely to the system (2.9) under the conditions 𝑁(𝑡)=𝑙 and 𝑇₁=𝑡₁,𝑇₂=𝑡₂,…,𝑇_𝑙=𝑡_𝑙. First, we note that the joint probability density function of (𝑥(𝑡),𝑥(𝑡_𝑙),…,𝑥(𝑡₁)) can be expressed as 𝑓_{𝑥(𝑡),𝑥(𝑡𝑙),…,𝑥(𝑡1)∣𝑁(𝑡)=𝑙,𝑡1,…,𝑡𝑙}𝑥,𝑥_𝑙,…,𝑥₁=𝑓_{𝑥(𝑡)∣𝑥(𝑡𝑙),𝑁(𝑡)=𝑙,𝑡1,…,𝑡𝑙}(𝑥)𝑓_{𝑥(𝑡𝑙)∣𝑥(𝑡𝑙−1),𝑁(𝑡)=𝑙,𝑡1,…,𝑡𝑙}𝑥_𝑙⋯𝑓_{𝑥(𝑡1)∣𝑁(𝑡)=𝑙,𝑡1,…,𝑡𝑙}𝑥₁, (4.17) then, for 𝑘=0,1,…,𝑙−1, consider that 𝑥(𝑡_𝑘+1)=𝑥(𝑡⁻_𝑘+1)𝑧_𝑘+1 as a product of two random variables where the first one has the probability density function given in (4.4), and 𝑧_𝑘+1 is the random jump factor at time 𝑡_𝑘+1. By the independence of 𝑥(𝑡⁻_𝑘+1) and 𝑧_𝑘+1, we then have 𝑓_{𝑥(𝑡𝑘+1)|𝑥(𝑡𝑘),𝑁(𝑡)=𝑙,𝑡1,…,𝑡𝑙}𝑥_𝑘+1=^∞₀𝑓_{𝑥(𝑡−𝑘+1)|𝑥(𝑡𝑘),𝑁(𝑡)=𝑙,𝑡1,…,𝑡𝑙}𝑥_𝑘+1𝑧_𝑘+1𝑔𝑧_𝑘+11𝑧^𝑛_𝑘+1𝑑𝑧_𝑘+1=^∞₀(2𝜋)^−𝑛/2Σ_𝑘𝑡_𝑘+1^−1/2𝑀_𝑘^−1𝑛_𝑖=11𝑀_𝑘⁻¹𝑥_𝑘+1/𝑧_𝑘+1^𝑖−1×exp2ln𝑀_𝑘⁻¹𝑥_𝑘+1𝑧_𝑘+1−𝜇_𝑘𝑡_𝑘+1^𝑇Σ_𝑘𝑡_𝑘+1⁻¹×ln𝑀_𝑘⁻¹𝑥_𝑘+1𝑧_𝑘+1−𝜇_𝑘𝑡_𝑘+1𝑔𝑧_𝑘+11𝑧^𝑛_𝑘+1𝑑𝑧_𝑘+1, (4.18) then (4.17) can be written as 𝑓_{𝑥(𝑡),𝑥(𝑡𝑙),…,𝑥(𝑡1)∣𝑁(𝑡)=𝑙,𝑡1,…,𝑡𝑙}𝑥,𝑥_𝑙,…,𝑥₁=(2𝜋)^−𝑛/2Σ_𝑙(𝑡)^−1/2𝑀_𝑙^−1𝑛_𝑖=11𝑀_𝑙⁻¹𝑥^𝑖−1×exp2ln𝑀_𝑙⁻¹𝑥−𝜇_𝑙(𝑡)^𝑇Σ_𝑙(𝑡)⁻¹ln𝑀_𝑙⁻¹𝑥−𝜇_𝑙(𝑡)×^𝑙−1_𝑘=0^∞₀(2𝜋)^−𝑛/2Σ_𝑘𝑡_𝑘+1^−1/2𝑀_𝑘^−1𝑛_𝑖=11𝑀_𝑘⁻¹𝑥_𝑘+1/𝑧_𝑘+1^𝑖−1×exp2ln𝑀_𝑘⁻¹𝑥_𝑘+1𝑧_𝑘+1−𝜇_𝑘𝑡_𝑘+1^𝑇Σ_𝑘𝑡_𝑘+1⁻¹×ln𝑀_𝑘⁻¹𝑥_𝑘+1𝑧_𝑘+1−𝜇_𝑘𝑡_𝑘+1𝑔𝑧_𝑘+11𝑧^𝑛_𝑘+1𝑑𝑧_𝑘+1. (4.19) The conditional probability density function of 𝑥(𝑡) given in (4.16) is obtained by integrating (4.19) with respect to 𝑥₁,𝑥₂,…, and 𝑥_𝑙.

Now, we can derive the unconditional probability distribution of the solution process of the 𝑛-dimensional system given in the following proposition.

Proposition 4.6. Under assumptions (i)–(iii), the probability density function of the solution process 𝑥(𝑡) to the n-dimensional system (2.9) is given by 𝑓_𝑥(𝑡)(𝑥)=^∞_𝑙=0^𝑡₀⋯^𝑡3₀^𝑡2₀×_ℝ𝑛+⋯_ℝ𝑛+(2𝜋)^−𝑛/2Σ_𝑙(𝑡)^−1/2𝑀_𝑙^−1𝑛_𝑖=11𝑀_𝑙⁻¹𝑥^𝑖−1×exp2ln𝑀_𝑙⁻¹𝑥−𝜇_𝑙(𝑡)^𝑇Σ_𝑙(𝑡)⁻¹ln𝑀_𝑙⁻¹𝑥−𝜇_𝑙(𝑡)×^𝑙−1_𝑘=0^∞₀(2𝜋)^−𝑛/2Σ_𝑘𝑡_𝑘+1^−1/2𝑀_𝑘^−1𝑛_𝑖=11𝑀_𝑘⁻¹𝑥_𝑘+1/𝑧_𝑘+1^𝑖−1×exp2ln𝑀_𝑘⁻¹𝑥_𝑘+1𝑧_𝑘+1−𝜇_𝑘𝑡_𝑘+1^𝑇Σ_𝑘𝑡_𝑘+1⁻¹×ln𝑀_𝑘⁻¹𝑥_𝑘+1𝑧_𝑘+1−𝜇_𝑘𝑡_𝑘+1×𝑔𝑧_𝑘+11𝑧^𝑛_𝑘+1𝑑𝑧_𝑘+1𝑑𝑥_𝑙⋯𝑑𝑥₁×^𝑙_𝑖=1𝜆𝑡_𝑖𝑑𝑡₁𝑑𝑡₂⋯𝑑𝑡_𝑙exp−^𝑡₀𝜆(𝑢)𝑑𝑢. (4.20)

Proof. The proof follows the argument in Proposition 3.3 by incorporating the random jumps.

Remark 4.7. It is obvious that the result in Proposition 4.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 process 𝑥(𝑡) of the Ornstein-Uhlenbeck equation 𝑑𝑥(𝑡)=𝐴𝑥(𝑡)𝑑𝑡+𝐶𝑑𝑤(𝑡),𝑡≥𝑡₀,𝑥𝑡₀=𝑥₀ (5.1) follows a multivariate normal distribution with mean 𝜇(𝑡) and covariance matrix Σ(𝑡), where 𝜇(𝑡)=exp^{𝐴(𝑡−𝑡0)}𝑥₀,Σ(𝑡)=^𝑡_𝑡0𝑒^{𝐴(𝑡−𝑢)}𝑉𝑒^{𝐴𝑇(𝑡−𝑢)}𝑑𝑢, (5.2) and 𝑉=𝐶𝐶^𝑇.

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 of 𝑥(𝑡).

Lemma 5.2. Under the conditions 𝑁(𝑡)=𝑙, and 𝑇₁=𝑡₁, 𝑇₂=𝑡₂,…,𝑇_𝑙=𝑡_𝑙, the solution process 𝑥(𝑡) for system (2.12) has a conditional probability density function as 𝑓_{𝑥(𝑡)∣𝑁(𝑡)=𝑙,𝑡1,…,𝑡𝑙}(𝑥)=_ℝ𝑛⋯_ℝ𝑛𝑓_{𝑥(𝑡),𝑥(𝑡𝑙),…,𝑥(𝑡1)∣𝑁(𝑡)=𝑙,𝑡1,…,𝑡𝑙}𝑥,𝑥_𝑙,…,𝑥₁𝑑𝑥_𝑙⋯𝑑𝑥₁=_ℝ𝑛⋯_ℝ𝑛(2𝜋)^−𝑛/2Σ_𝑙(𝑡)^−1/2−1exp2𝑥−𝜇_𝑙(𝑡)^𝑇Σ_𝑙(𝑡)⁻¹𝑥−𝜇_𝑙(𝑡×_{𝑙−1𝑘=0∞0})(2𝜋)^−𝑛/2Σ_𝑘𝑡_𝑘+1^−1/2−1×exp2𝑥_𝑘+1𝑧_𝑘+1−𝜇_𝑘𝑡_𝑘+1^𝑇Σ_𝑘𝑡_𝑘+1⁻¹𝑥_𝑘+1𝑧_𝑘+1−𝜇_𝑘𝑡_𝑘+1×𝑔𝑧_𝑘+11𝑧^𝑛_𝑘+1𝑑𝑧_𝑘+1𝑑𝑥_𝑙⋯𝑑𝑥₁, (5.3) where 𝑔 is the common probability density function of 𝑧_𝑘, 𝑘=1,2,…,𝑙.

Proof. As we noted before that the joint probability density function of (𝑥(𝑡),𝑥(𝑡_𝑙),…,𝑥(𝑡₁)), under the conditions 𝑁(𝑡)=𝑙 and 𝑇₁=𝑡₁, 𝑇₂=𝑡₂,…,𝑇_𝑙=t_𝑙, can be expressed as 𝑓_{𝑥(𝑡),𝑥(𝑡𝑙),…,𝑥(𝑡1)∣𝑁(𝑡)=𝑙,𝑡1,…,𝑡𝑙}𝑥,𝑥_𝑙,…,𝑥₁=𝑓_{𝑥(𝑡)∣𝑥(𝑡𝑙),𝑁(𝑡)=𝑙,𝑡1,…,𝑡𝑙}(𝑥)𝑓_{𝑥(𝑡𝑙)∣𝑥(𝑡𝑙−1),𝑁(𝑡)=𝑙,𝑡1,…,𝑡𝑙}𝑥_𝑙⋯𝑓_{𝑥(𝑡1)∣𝑁(𝑡)=𝑙,𝑡1,…,𝑡𝑙}𝑥₁, (5.4) then, by applying the result in Lemma 5.1 on each interval between jumps [𝑡_𝑘,𝑡_𝑘+1), we have, for 𝑠∈[𝑡_𝑘,𝑡_𝑘+1), 𝑓_{𝑥(𝑠)∣𝑥(𝑡𝑘),𝑁(𝑡)=𝑙,𝑡1,…,𝑡𝑙}(𝑥)=𝜙𝑥;𝜇_𝑘(𝑠),Σ_𝑘(𝑠)=(2𝜋)^−𝑛/2Σ_𝑘(𝑠)^−1/2−1exp2𝑥−𝜇_𝑘(𝑠)^𝑇Σ_𝑘(𝑠)⁻¹𝑥−𝜇_𝑘(𝑠), (5.5) where 𝜇_𝑘(𝑠)=exp^{𝐴𝑘(𝑠−𝑡𝑘)}𝑥_𝑘,Σ_𝑘(𝑠)=^𝑠_𝑡𝑘𝑒^{𝐴𝑘(𝑠−𝑢)}𝑉_𝑘𝑒^{𝐴𝑇𝑘(𝑠−𝑢)}𝑑𝑢, (5.6) and 𝑉_𝑘=𝐶_𝑘𝐶^𝑇_𝑘. Then, for 𝑘=0,1,…,𝑙−1, consider that 𝑥(𝑡_𝑘+1)=𝑥(𝑡⁻_𝑘+1)𝑧_𝑘+1 as a product of two random variables where the first one has the probability density function given in (5.5), and 𝑧_𝑘+1 is the random jump factor at time 𝑡_𝑘+1. By the independence of 𝑥(𝑡⁻_𝑘+1) and 𝑧_𝑘+1, then 𝑓_{𝑥(𝑡𝑘+1)∣𝑁(𝑡)=𝑙,𝑡1,…,𝑡𝑙}𝑥_𝑘+1=^∞₀𝑓_{𝑥(𝑡−𝑘+1)∣𝑁(𝑡)=𝑙,𝑡1,…,𝑡𝑙}𝑥_𝑘+1𝑧_𝑘+1𝑔𝑧_𝑘+11𝑧^𝑛_𝑘+1𝑑𝑧_𝑘+1=^∞₀(2𝜋)^−𝑛/2Σ_𝑘𝑡_𝑘+1^−1/2−1×exp2𝑥_𝑘+1𝑧_𝑘+1−𝜇_𝑘𝑡_𝑘+1^𝑇Σ_𝑘𝑡_𝑘+1⁻¹𝑥_𝑘+1𝑧_𝑘+1−𝜇_𝑘𝑡_𝑘+1×𝑔𝑧_𝑘+11𝑧^𝑛_𝑘+1𝑑𝑧_𝑘+1. (5.7) Then the conditional joint probability density function (5.4) can be written as 𝑓_{𝑥(𝑡),𝑥(𝑡𝑙),…,𝑥(𝑡1)∣𝑁(𝑡)=𝑙,𝑡1,…,𝑡𝑙}𝑥,𝑥_𝑙,…,𝑥₁=(2𝜋)^−𝑛/2Σ_𝑙(𝑡)^−1/2−1exp2𝑥−𝜇_𝑙(𝑡)^𝑇Σ_𝑙(𝑡)⁻¹𝑥−𝜇_𝑙(𝑡)×^𝑙−1_𝑘=0^∞₀(2𝜋)^−𝑛/2Σ_𝑘𝑡_𝑘+1^−1/2−1×exp2𝑥_𝑘+1𝑧_𝑘+1−𝜇_𝑘𝑡_𝑘+1^𝑇(Σ_𝑘(𝑡_𝑘+1))⁻¹𝑥_𝑘+1𝑧_𝑘+1−𝜇_𝑘𝑡_𝑘+1×𝑔𝑧_𝑘+11𝑧^𝑛_𝑘+1𝑑𝑧_𝑘+1. (5.8) The conditional probability density function of 𝑥(𝑡) given in (5.3) is obtained by integrating (5.8) with respect to 𝑥₁,𝑥₂,…, and 𝑥_𝑙.

Finally, following the same argument in Proposition 3.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 process 𝑥(𝑡) of system (2.12) is given by 𝑓_𝑥(𝑡)(𝑥)=^∞_𝑙=0^𝑡₀⋯^𝑡3₀^𝑡2_0ℝ𝑛⋯_ℝ𝑛(2𝜋)^−𝑛/2Σ_𝑙(𝑡)^−1/2−1×exp2𝑥−𝜇_𝑙(𝑡)^𝑇Σ_𝑙(𝑡)⁻¹𝑥−𝜇_𝑙(𝑡)×^𝑙−1_𝑘=0^∞₀(2𝜋)^−𝑛/2Σ_𝑘𝑡_𝑘+1^−1/2−1×exp2𝑥_𝑘+1𝑧_𝑘+1−𝜇_𝑘𝑡_𝑘+1^𝑇Σ_𝑘𝑡_𝑘+1⁻¹𝑥_𝑘+1𝑧_𝑘+1−𝜇_𝑘𝑡_𝑘+1×𝑔𝑧_𝑘+11𝑧^𝑛_𝑘+1𝑑𝑧_𝑘+1𝑑𝑥_𝑙⋯𝑑𝑥₁^𝑙_𝑖=1𝜆𝑡_𝑖𝑑𝑡₁𝑑𝑡₂⋯𝑑𝑡_𝑙×exp−^𝑡₀𝜆(𝑢)𝑑𝑢. (5.9)

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 𝑛=1 the Fokker-Planck equation corresponding to system (2.9) even in the absence of discrete time interventions is given by𝜕𝑓𝜕𝜕𝑡=−+𝜕𝜕𝑥𝐴𝑥𝑓²𝜕𝑥²𝐵²𝑥²𝑓=2𝐵²−𝐴𝑓+4𝐵²𝑥𝜕−𝐴𝜕𝑥+𝐵²𝑥²𝜕²𝜕𝑥². (6.1) 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.

Acknowledgments

This research was supported by the Mathematical Sciences Division, U.S. Army Research Office, Grant no. W911NF-07-1-0283. The authors would like to express their gratitude to the referee and the editor for their constructive, scientific thorough review and comments on the paper.