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

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, themajor 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.


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 piecewisedeterministic Markov process, where transitions between discrete modes are triggered by random events and deterministic conditions for hitting the boundary, while the continuousstate 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

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 x t be a real n-dimensional process.A k and B j k are n×n matrices for any k ∈ N∪{0} and j 1, 2, . . ., q.Let C r k be n-dimensional vectors for any k ∈ N ∪ {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: C r N t dw r t , t ≥ t 0 , x t 0 x 0 > 0, 2.1 where w t w 1 t , . . ., w q 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 x 1 , x 2 , . . ., x n > 0 as x i > 0 for all i 1, 2, . . ., n.When C r k 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 Here, A k and B j k 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, B j , and C r 's are fixed over time, that is, the solution process on a subinterval between jumps.Under the condition that the matrices A, B 1 , B 2 , . . ., B q 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 x t ≡ x t, t 0 , x 0 be the solution of the following initial value problem (IVP): then the x t is expressed by provided that AB j B j A and B j B j B j B j for all j, j 1, 2, . . ., q.
Let x t, T k , x k be the solution to system 2.3 with t 0 : T k , x 0 : x k , A : A k , B j : B j k , and C r : C r k .Now, we consider the following system of two interconnected stochastic dynamics: where z k , k 1, 2, 3, . .., are iid positive random variables with z 0 1, and Here, we assume that N t , w t , w t , and z k 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
i for all k ∈ N ∪ {0} and j, j 1, 2, . . ., q, then the solution to the system 2.5 is given by

2.6
Proof.By applying the result of 2.4 on the subintervals T k−1 , T k , for k 1, . . ., N t , and T N t , t , we have the solution to system 2.5 as the following piecewise function: where x s, T k , x k is a solution process in 2.4 with t 0 : T k , x 0 : x k , A : A k , B j : B j k , and C r : C r k , then, we have

2.8
Next, substitute 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 C r k 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
The solution of the above system is given in the following corollary.The result follows from Proposition 2.2 by letting C r k be zero for all k and r.
i for all k ∈ N ∪ {0} and j, j 1, 2, . . ., q, then the solution of system 2.9 is given by

2.10
In the case when B j k 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
Denote p k for all k, then the above system can be rewritten as

2.12
where C k '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

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 R , and in this case, A k and B k 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 Proof.Given that N t l, and T 1 t 1 , . .., T l t l , from 2.10 , we have

3.3
where we denote V l k 1 ln z k and S as the sum of last three terms in 3.3 .Since h is the common probability density function of ln z k , 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
Since z k and w t are independent, then V and S are also independent.Using transformation method 26 on ln x t V S, we can obtain the conditional probability density function of ln x t as If follows that

3.6
Having obtained the conditional probability density function of x 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 process x t to the system 2.5 is given by

3.7
Proof.From 3.1 and 3.2 , the joint density function of x t , T 1 , . .., T l given the condition N t l is given by

3.8
Integrating with respect to t 1 , t 2 , . .., and t l gives the probability density function of x t given that N t l as 3.9 By multiplying the above density by the probability of l jumps, namely, and then taking the summation over l, the marginal probability density function of x t given in 3.7 is established.

Probability Distribution of Multivariate Linear Homogeneous Models
In this section, we will derive the probability density function of the solution process x t for the n-dimensional stochastic system 2.9 under the following assumptions.
i The drift and diffusion coefficients, A k and B j k for all k ∈ N ∪ {0} and j 1, 2, . . ., q, are diagonalizable.
ii The coefficients in each regime pairwise commute, that is, k for all k ∈ N ∪ {0} and j, l 1, 2, . . ., q. iii A k ∈ C for k ∈ N ∪ {0}, where C denotes the set of n × n diagonalizable matrices whose eigenvector matrix, M, is such that M −1 x > 0 for any x > 0.
We first consider the stochastic system on the interval between jumps.Given that N t l and T 1 t 1 , T 2 t 2 , . . ., T l t l , consider the following SDE on t k , t k 1 , for some k 0, 1, . . ., l: 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 The following lemma gives the useful fact that the linear transformation of a multivariate normal random vector is again multivariate normally distributed.

4.3
To find the probability density function of x 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, A k s and B j k s have the same eigenvector matrix, denoted by M k .Moreover, where { a 1 k , . . ., a n k } and { b 1 k,j , . . ., b n k,j } are the sets of eigenvalues of A k and B j k , respectively, for all k, j.Next, we define a linear transformation y M −1 k x, then x M k y.Define ln y ln y 1 , ln y 2 , . . ., ln y n 28 .The below proposition gives the probability distribution of solution process x t on the interval t k , t k 1 over which the coefficients are constant.

Lemma 4.4. Under assumptions (i)-(iii), the process x s satisfying the SDE 4.1 has a probability density function given by
where

4.5
Proof.For s ∈ t k , t k 1 , we have y s M −1 k x s , and x s M k y s .By multiplying M −1 k on both sides of the SDE 4.1 , we obtain the transformed SDE as follows: ⇒ dy s A k y s ds q j 1 B j k y s dw j s ,

4.6
where A k and B j k are diagonal matrices as defined before.From the application of Lemma 2.1 with C r k 0, the solution process of the transformed system 4.6 is It follows from assumption iii that y s M −1 k x s > 0 since x s > 0, then, we can rewrite system 4.7 in the following form: where μ k s and B * k are defined above.
Since w t is a standard Wiener process, then w s − w t k has a multivariate normal distribution with mean zero and covariance matrix s − t k I n , where I n is the n × n identity matrix.Then, by Theorem 4.3, ln y s as a linear transformation of w s − w t k is also multivariate normally distributed with mean μ k s and covariance matrix Σ k s , where The probability density function of ln y s is We now apply the method of transformation from y ln s to y s , then, Then, from 4.10 and 4.11 , we have the probability density function of y s as

4.13
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 M k y s , or y s M −1 k x s , we have The result follows from combining 4.13 and 4.14 .Now, we are to develop the conditional probability density function of solution process x t , which is a result parallel to the one-dimensional case in Lemma 3.2.

Lemma 4.5. Under assumptions (i)-(iii) and given that N t
l, and T 1 t 1 , T 2 t 2 , . . ., T l t l , the solution process x t has a conditional probability density function as where g is the common probability density function of z k , k 1, 2, . . ., l.
Proof.We will apply the result in Lemma 4.4 piecewisely to the system 2.9 under the conditions N t l and T 1 t 1 , T 2 t 2 , . . ., T l t l .First, we note that the joint probability density function of x t , x t l , . . ., x t 1 can be expressed as 4.17 then, for k 0, 1, . . ., l − 1, consider that x t k 1 x t − k 1 z k 1 as a product of two random variables where the first one has the probability density function given in 4.4 , and z k 1 is the random jump factor at time t k 1 .By the independence of x t − k 1 and z k 1 , we then have 4.18 then 4.17 can be written as f x t ,x t l ,...,x t 1 |N t l,t 1 ,...,t l x, x l , . . ., x

4.19
The conditional probability density function of x t given in 4.16 is obtained by integrating 4.19 with respect to x 1 , x 2 , . .., and x l .Now, we can derive the unconditional probability distribution of the solution process of the n-dimensional system given in the following proposition.Proposition 4.6.Under assumptions (i)-(iii), the probability density function of the solution process x t to the n-dimensional system 2.9 is given by

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.

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 .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 x t .Lemma 5.2.Under the conditions N t l, and T 1 t 1 , T 2 t 2 , . . ., T l t l , the solution process x t for system 2.12 has a conditional probability density function as where g is the common probability density function of z k , k 1, 2, . . ., l.
Proof.As we noted before that the joint probability density function of x t , x t l , . . ., x t 1 , under the conditions N t l and T 1 t 1 , T 2 t 2 , . . ., T l t l , can be expressed as where Then, for k 0, 1, . . ., l − 1, consider that x t k 1 x t − k 1 z k 1 as a product of two random variables where the first one has the probability density function given in 5.5 , and z k 1 is the random jump factor at time t k 1 .By the independence of x t − k 1 and z k 1 , then

5.8
The conditional probability density function of x t given in 5.3 is obtained by integrating 5.8 with respect to x 1 , x 2 , . .., and x l .

Journal of Probability and Statistics
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 x t of system 2.12 is given by 5.9

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
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.
Cox and Lewis 24 or Crowder et al. 25 .For a non-homogeneous Poisson process N t with intensity λ t , the joint density function of the jump times T 1 , T 2 , . . ., T l conditioned on N t l is given by f T 1 ,T 2 ,...,T l |N t l t 1 , t 2 , ...l .3.1Next lemma gives the conditional probability density function of x t given the number of jumps and the jump times.Lemma 3.2.Given that N t l, and T 1 t 1 , ..., T l t l , x t has a probability density function asf x t |N t l,t 1 ,...,t l x 1 x ∞ −∞ h * l ln x − s φ s; μ, σ ds, x > 0, 3.2where h * l is the lth convolution of the common probability density function h of ln z k , for k 1, 2, . . ., l, and φ •; μ, σ denotes the normal density function with mean μ and variance σ.
Lemma 5.1 see 18, 29 .The solution process x t of the Ornstein-Uhlenbeck equation