Sensitivities of some performance measures of quasi-birth-and-death processes

Abstract.

Quasi-birth-and-death processes (QBDs) form a class of Markovian models with many applications in which the level-variable X(t) records the number of individuals in some system and a phase-variable φ(t) records some additional information about the system. We study the sensitivity of various performance measures of QBDs with respect to the perturbation of a generator matrix. Thereby, we extend the sensitivity analysis of level-dependent QBDs in Gomez-Corral and Lopez-Garcia[4] and develop results with a focus on applications in modelling healthcare systems.

Keywords:

• Matrix analytic methods
• performance measures
• quasi-birth-and-death process
• sensitivities

1. Introduction

1.1. Motivating example: the hospital

The functioning of a hospital can be successfully modeled via a quasi-birth-and-death process. Such a model of a hospital captures random arrivals and departures, the random length of stay of an individual patient, and possibly different types of patients (i.e., complex and non-complex). Based on the occupation rates, the hospital must adjust their decisions to deliver health services to the community. Undoubtedly, hospital managers face the hard problem of improving hospital performance. The challenge also comes from the fact that “hospital performance” is a concept that is not clearly defined and based on not only (or even primarily) financial criteria. Hospital performance evaluation takes into account many factors, such as mortality, safety of care, readmission rate, patient experience, and service times. These measures depend on the number of patients in the hospital, where both too low and too high numbers come with challenges for a hospital. On the one hand, there may not be enough demand to ensure a wide array of specialists, and on the other hand, there is a risk of congestion. Both situations result in worse outcomes for patients, greater expense for the hospital, and a heavier load on hospital staff.

Hospitals run at various levels of urgency depending on current occupancy levels. Applying different policies, like scheduling on-call physicians, planning/cancelling elective surgeries, etc., allows for adjusting to the urgency level and impacting the arrival and departure rates. It is crucial to know how sensible various performance metrics are with respect to these hospital decisions in order to help stabilize the hospital at the optimal level.

1.2. Motivating example: two bacterial strains

We recall an example from Lipsitch et al.[10] and Gómez-Corral and López-García[4] about two bacterial strains: antibiotic-sensitive (strain 1) and antibiotic-resistant (strain 2). These bacterial strains spread among patients in a hospital such that the infection by one bacterial strain provides immunity against the other. Patients in the hospital are routinely provided antibiotics 1 and 2, regardless of whether these patients are infected or not by bacteria; antibiotic 1 is only effective against the AS bacterial strain, whereas antibiotic 2 is effective against both strains of bacteria. There might be effects of spontaneous clearance of bacteria and contributions of antibiotics 1 and 2 to this recovery. In such a situation, we might be interested in the outbreak for the type-k strain and global outbreak (length). Similarly as in Section 1.1, the transition matrix, and hence performance measures of interest, can be impacted by changes in administering antibiotics.

1.3. Aim and organization of the paper

Quasi-birth-and-death processes (QBDs) are the fundamental class of Markovian models in the theory of matrix-analytic methods, with a level-variable X(t) and a phase-variable φ(t) forming a two-dimensional state space. In many applications of the QBDs, the phase variable φ(t) is used to model information about the underlying environment that drives the evolution of some systems. A QBD is a model that lends itself to representing healthcare systems in a natural, intuitive manner (see Figure 1 and the above motivating example), and so the application potential of the QBDs in this area is immense, as demonstrated by Heydar et al.[6].

Figure 1. Evolution of a QBD $\{(X(t),\phi (t)):t\ge 0\}$ modeling the total number of patients X(t) and some information φ(t) about the system e.g. (here) the number of patients of a particular class.

The aim of this article is to study the sensitivities of various performance measures of QBDs with respect to a parameter dependent generator matrix. Thereby, we extend the sensitivity analysis of level-dependent QBDs in Gomez-Corral and Lopez-Garcia[4] to a wide range of metrics of interest, and focus on applications in healthcare systems. This work also lays the foundation or computing and analyzing such sensitivities given a hospital data set[15]. We note that the ideas presented in this article about sensitivities of continuous time QBDs, can be successfully applied in discrete time case.

This article is organized as follows. In Section 2, we describe the model in use, namely level-dependent QBD. In Section 3, we present performance metrics of interest and the first results concerning these metrics. In Section 4, we derive sensitivities and illustrate these quantities in several examples.

2. Level-dependent QBD model

Suppose that $\{(X(t),\phi (t)):t\ge 0\}$ is a continuous-time Markov chain (CTMC) with a two-dimensional state (X(t), φ(t)) consisting of the level variable X(t) and the phase variable, φ(t), taking values in an irreducible state space given by $$\begin{array}{l}\mathcal{S}=\{(n,i):n=0,1,2,\dots ,N;i=0,1,\dots ,m_n\},\end{array}$$ and with transition rates recorded in the generator matrix $\mathbf{Q}=[q_{(n,i)(n^{\prime },i^{\prime })}{]}_{(n,i),(n^{\prime },i^{\prime })\in \mathcal{S}}$ made of block matrices ${\mathbf{Q}}^{[n,n^{\prime }]}=[q_{(n,i)(n^{\prime },i^{\prime })}{]}_{i=0,1,\dots ,m_n,i^{\prime }=0,1,\dots ,m_{n^{\prime }}}$ such that $$\begin{array}{ll}\mathbf{Q}& =[{\mathbf{Q}}^{[n,n^{\prime }]}]\\ & =\left[\begin{array}{cccccccc}{\mathbf{Q}}^{[0,0]}& {\mathbf{Q}}^{[0,1]}& \mathbf{0}& \cdots & \cdots & \cdots & \cdots & \mathbf{0}\\ {\mathbf{Q}}^{[1,0]}& {\mathbf{Q}}^{[1,1]}& {\mathbf{Q}}^{[1,2]}& \mathbf{0}& \cdots & \cdots & \cdots & \mathbf{0}\\ \mathbf{0}& {\mathbf{Q}}^{[2,1]}& {\mathbf{Q}}^{[2,2]}& {\mathbf{Q}}^{[2,3]}& \cdots & \cdots & \cdots & \mathbf{0}\\ ⋮& ⋮& ⋮& ⋮& \cdots & \cdots & \cdots & ⋮\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \cdots & {\mathbf{Q}}^{[N-1,N-2]}& {\mathbf{Q}}^{[N-1,N-1]}& {\mathbf{Q}}^{[N-1,N]}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \cdots & \mathbf{0}& {\mathbf{Q}}^{[N,N-1]}& {\mathbf{Q}}^{[N,N]}\end{array}\right]\end{array}$$ so that only transitions to the neighboring levels are possible. We refer to such a process as a continuous-time level-dependent quasi-birth-and-death process (LD-QBD), bounded from above by level N. The level variable X(t) may be used to record the number of individuals in some system at time t, while φ(t) may be used to record some additional information about the system at time t.

Such a LD-QBD is a continuous-time Markov chain (CTMC), and so standard expressions from the theory of CTMCs apply. However, as the state space 𝒮 may be very large, we apply methods from the theory of matrix-analytic methods instead, resulting in efficient computational methods. We also consider the case when the LD-QBD does not have an upper boundary N.

We note that a range of various transient and stationary performance measures of such defined LD-QBD can be readily derived using the existing methods in the literature of matrix-analytic methods, see Grant[5], Ramaswami[13], and Joyner and Fralix[8]. Since these performance measures depend on the parameters of the model, we are interested in the sensitivity analysis of these measures.

Suppose that the generator 𝐐(𝜽) of a LD-QBD depends on some parameters recorded in a row vector $\mathit{\theta }=[{\theta }_i{]}_{i=1,\dots ,k}$. Various stationary (long-run) and transient (time-dependent) quantities in the analysis of the LD-QBDs, recorded as matrices 𝐀(𝜽) = [A_ij(𝜽)], can then be expressed using expressions involving 𝐐(𝜽), and so also depend on 𝜽. We derive the sensitivity analysis of relevant quantities, where given a vector of parameters $\mathit{\theta }=[{\theta }_i{]}_{i=1,\dots ,k}$ and matrix 𝐀(𝜽), we write $$\begin{array}{l}\frac{\partial \mathbf{A}(\mathit{\theta })}{\partial \mathit{\theta }}=\left[\begin{array}{c}\frac{\partial \mathbf{A}(\mathit{\theta })}{\partial {\theta }_1}\dots \frac{\partial \mathbf{A}(\mathit{\theta })}{\partial {\theta }_k}\end{array}\right].\end{array}$$

3. Performance measures

In this section, we consider the following key quantities in the analysis of the LD-QBDs:

1. the stationary distribution 𝝅 in Subsection 3.1;

2. first hitting times to lower levels in Subsection 3.2;

3. first hitting times to higher levels in Subsection 3.3;

4. sojourn times in specified sets in Subsection 3.4;

5. the distribution of the process observed at time t in Subsection 3.5.

We follow the approach summarized in Grant[5], which is built on the results in Ramaswami[13], Joyner and Fralix[8], and Phung-Duc et al.[12]. We also state expressions for the relevant Laplace-Stieltjes transforms (LSTs) of various quantities. These can be inverted using numerical inversion techniques by Abate and Whitt[1], Den Iseger[3], or Horváth et al.[7], to compute the corresponding quantities. We refer a reader to [2, 6, 11] for computations of stationary distribution, sojourn times, and applications.

3.1. Stationary distribution

For all $n=0,1,\dots ,N$, $i=0,1,\dots ,m_n$, define the limiting probabilities $$\begin{array}{l}{\pi }_{(n,i)}=\underset{t\to \mathrm{\infty }}{lim}\mathbb{P}(X(t)=n,\phi (t)=i)\end{array}$$ recording the long-run proportions of time spent in states (n, i), and collecting these in a row vector $\mathit{\pi }=[{\mathit{\pi }}_n{]}_{n=0,1,\dots ,N}$, where ${\mathit{\pi }}_n=[{\pi }_{(n,i)}{]}_{i=0,1,\dots ,m_n}$.

To evaluate 𝝅, we follow the approach in Grant[5], and write the expressions for 𝝅_n in terms of 𝝅_N (rather than in terms of 𝝅₀, since potential close-to-zero values in 𝝅₀ may lead to computational errors). We consider matrices $$\begin{array}{l}{\hat{\mathbf{R}}}^{(n)}=[{\hat{R}}_{ij}^{(n)}{]}_{i=0,1,\dots ,m_{n+1},j=0,1,\dots ,m_n}\end{array}$$ recording the expected times ${\hat{R}}_{ij}^{(n)}$ spent in states (n, j), given the process starts in state $(n+1,i)$, before returning to level n+1. Then, for $n=0,1,\dots ,N-1$, we know that $$\begin{array}{ll}{\hat{\mathbf{R}}}^{(0)}(s)& =-{\mathbf{Q}}^{[1,0]}({\mathbf{Q}}^{[0,0]}-s\mathbf{I}{)}^{-1},\\ {\hat{\mathbf{R}}}^{(n)}(s)& =-{\mathbf{Q}}^{[n+1,n]}({\hat{\mathbf{R}}}^{(n-1)}(s){\mathbf{Q}}^{[n-1,n]}+{\mathbf{Q}}^{[n,n]}-s\mathbf{I}{)}^{-1}\end{array}$$ with ${\hat{\mathbf{R}}}^{(n)}={\hat{\mathbf{R}}}^{(n)}(0)$. Let (1) $$\begin{array}{ll}{\mathit{\pi }}_n& ={\mathit{\pi }}_{n+1}{\hat{\mathbf{R}}}^{(n)}={\mathit{\pi }}_N\prod _{k=N-1}^n{\hat{\mathbf{R}}}^{(k)}.\end{array}$$(1)

Finally, 𝝅_N is the solution of the set of equations, (2) $$\begin{array}{l}\phantom{.}\{\begin{array}{l}{\mathit{\pi }}_N({\hat{\mathbf{R}}}^{(N-1)}{\mathbf{Q}}^{[N-1,N]}+{\mathbf{Q}}^{[N,N]})=\mathbf{0}\\ {\mathit{\pi }}_N(\mathbf{1}+\sum _{n=0}^{N-1}\prod _{k=N-1}^n{\hat{\mathbf{R}}}^{(k)}\mathbf{1})=1.\end{array}\end{array}$$(2)

Remark 1.

Similar methods may be applied to derive 𝝅 for a LD-QBD in which the level variable X(t)≥0 has no upper boundary N. First, following Phung-Duc et al.[12], find a sufficiently large truncation level L such that $||{\hat{\mathbf{R}}}_L^{(n)}-{\hat{\mathbf{R}}}_{L-1}^{(n)}||<ϵ$ for a required criterion ϵ > 0, where ${\hat{\mathbf{R}}}_L^{(n)}$ and ${\hat{\mathbf{R}}}_{L-1}^{(n)}$ denote matrix ${\hat{\mathbf{R}}}^{(n)}$ computed for a LD-QBD with an upper boundary N = L and $N=L-1$, respectively. Next, apply the approximation ${\hat{\mathbf{R}}}^{(n)}\approx {\hat{\mathbf{R}}}_L^{(n)}$. This technique may be applied to the remaining quantities.

3.2. First hitting times to lower levels

We analyze here the distribution of time for the process to first reach level n−k and does so in state $(n-k,j)$ given start from state (n, i). Let ${\theta }_n=inf\{t>0:X(t)=n\}$ and, for any $n=1,\dots ,N$, $k=1,\dots ,n$, let $$\begin{array}{l}{\mathbf{g}}^{n,n-k}(t)=[g_{ij}^{n,n-k}(t){]}_{i=1,\dots ,m_n;j=1,\dots ,m_{n-k}}\end{array}$$ be the matrix such that the entry (3) $$\begin{array}{ll}{g}_{ij}^{n,n-k}(t)& =\frac{\partial }{\partial t}\mathbb{P}({\theta }_{n-k}\le t,\phi ({\theta }_{n-k})=j|X(0)=n,\phi (0)=i)\end{array}$$(3) records the probability density that the process first visits level n−k at time t and does so in phase j, given that it starts from level n in phase i. We also define $$\begin{array}{l}{\mathbf{G}}^{n,n-k}(s)=[G_{ij}^{n,n-k}(s){]}_{i=1,\dots ,m_n;j=1,\dots ,m_{n-k}}\end{array}$$ where (4) $$\begin{array}{ll}{G}_{ij}^{n,n-k}(s)& ={\int }_0^{\mathrm{\infty }}{e}^{-st}{g}_{ij}^{n,n-k}(t)dt\end{array}$$(4) is the corresponding Laplace–Stieltjes transform (LST) matrix. Further, let ${\mathbf{G}}^{n,n-k}={\mathbf{G}}^{n,n-k}(0)$ and ${\mathbb{E}}^{n,n-k}=-\frac{\partial }{\partial s}{\mathbf{G}}^{n,n-k}(s){|}_{s=0}$ be the probability and expectation matrix, respectively.

To evaluate ${\mathbf{G}}^{n,n-k}(s)$ and ${\mathbf{G}}^{n,n-k}$, we apply the following recursion summarized in Grant[5], (5) $$\begin{array}{ll}{\mathbf{G}}^{N,N-1}(s)& =-({\mathbf{Q}}^{[N,N]}-s\mathbf{I}{)}^{-1}{\mathbf{Q}}^{[N,N-1]},\end{array}$$(5) (6) $$\begin{array}{ll}{\mathbf{G}}^{n,n-1}(s)& =-({\mathbf{Q}}^{[n,n]}-s\mathbf{I}+{\mathbf{Q}}^{[n,n+1]}{\mathbf{G}}^{n+1,n}(s){)}^{-1}{\mathbf{Q}}^{[n,n-1]}\\ & \text{for }n=N-1,\dots ,1\end{array}$$(6) and (7) $$\begin{array}{ll}{\mathbf{G}}^{n,n-k}(s)& ={\mathbf{G}}^{n,n-1}(s){\mathbf{G}}^{n-1,n-2}(s)\times \cdots \times {\mathbf{G}}^{n-k+1,n-k}(s).\end{array}$$(7)

By the above recursion, it follows that $$\begin{array}{ll}{\mathbb{E}}^{n,n-k}& ={\mathbb{E}}^{n,n-1}{\mathbf{G}}^{n-1,n-2}\times \cdots \times {\mathbf{G}}^{n-k+1,n-k}+\cdots \\ & +{\mathbf{G}}^{n,n-1}{\mathbf{G}}^{n-1,n-2}\times \cdots \times {\mathbf{G}}^{n-k+2,n-k+1}{\mathbb{E}}^{n-k+1,n-k}\\ & =\sum _{m=n}^{n-k+1}\left(\prod _{\ell =n}^{m+1}{\mathbf{G}}^{\ell ,\ell -1}\right){\mathbb{E}}^{m,m-1}\left(\prod _{\ell =m-1}^{n-k+1}{\mathbf{G}}^{\ell ,\ell -1}\right).\end{array}$$

Then, by rearranging (5(5) $$\begin{array}{ll}{\mathbf{G}}^{N,N-1}(s)& =-({\mathbf{Q}}^{[N,N]}-s\mathbf{I}{)}^{-1}{\mathbf{Q}}^{[N,N-1]},\end{array}$$(5) ) and then taking $\partial /\partial s$ at s = 0, we get $$\begin{array}{ll}& ({\mathbf{Q}}^{[N,N]}-s\mathbf{I}){\mathbf{G}}^{N,N-1}(s)=-{\mathbf{Q}}^{[N,N-1]}\\ & \Rightarrow -{\mathbf{G}}^{N,N-1}+{\mathbf{Q}}^{[N,N]}(-{\mathbb{E}}^{N,N-1})=\mathbf{0}\\ & \Rightarrow {\mathbf{Q}}^{[N,N]}{\mathbb{E}}^{N,N-1}=-{\mathbf{G}}^{N,N-1}\end{array}$$ which implies that (8) $$\begin{array}{ll}{\mathbb{E}}^{N,N-1}& =-({\mathbf{Q}}^{[N,N]}{)}^{-1}{\mathbf{G}}^{N,N-1}.\end{array}$$(8)

Furthermore, by (6(6) $$\begin{array}{ll}{\mathbf{G}}^{n,n-1}(s)& =-({\mathbf{Q}}^{[n,n]}-s\mathbf{I}+{\mathbf{Q}}^{[n,n+1]}{\mathbf{G}}^{n+1,n}(s){)}^{-1}{\mathbf{Q}}^{[n,n-1]}\\ & \text{for }n=N-1,\dots ,1\end{array}$$(6) ), for $n=N-1,\dots ,1$ $$\begin{array}{ll}& ({\mathbf{Q}}^{[n,n]}-s\mathbf{I}+{\mathbf{Q}}^{[n,n+1]}{\mathbf{G}}^{n+1,n}(s)){\mathbf{G}}^{n,n-1}(s)=-{\mathbf{Q}}^{[n,n-1]}\\ & \Rightarrow (-\mathbf{I}-{\mathbf{Q}}^{[n,n+1]}{\mathbb{E}}^{n+1,n}){\mathbf{G}}^{n,n-1}+({\mathbf{Q}}^{[n,n]}\\ & +{\mathbf{Q}}^{[n,n+1]}{\mathbf{G}}^{n+1,n})(-{\mathbb{E}}^{n,n-1})=\mathbf{0}\\ & \Rightarrow ({\mathbf{Q}}^{[n,n]}+{\mathbf{Q}}^{[n,n+1]}{\mathbf{G}}^{n+1,n}){\mathbb{E}}^{n,n-1}=(-\mathbf{I}-{\mathbf{Q}}^{[n,n+1]}{\mathbb{E}}^{n+1,n}){\mathbf{G}}^{n,n-1}\end{array}$$ and we conclude that (9) $$\begin{array}{ll}{\mathbb{E}}^{n,n-1}& =-({\mathbf{Q}}^{[n,n]}+{\mathbf{Q}}^{[n,n+1]}{\mathbf{G}}^{n+1,n}{)}^{-1}(\mathbf{I}+{\mathbf{Q}}^{[n,n+1]}{\mathbb{E}}^{n+1,n}){\mathbf{G}}^{n,n-1}.\end{array}$$(9)

3.3. First hitting times to higher levels

We analyze here the distribution of time for the process to first reach level n+k and does so in state $(n+k,j)$ given start from state (n₀, i). For any $n=0,1,\dots ,N-1$, $k=1,\dots ,N-n$, let $$\begin{array}{l}{\mathbf{h}}^{n,n+k}(t)=[h_{ij}^{n,n+k}(t){]}_{i=1,\dots ,m_n;j=1,\dots ,m_{n+k}}\end{array}$$ be the matrix such that (10) $$\begin{array}{ll}{h}_{ij}^{n,n+k}(t)& =\frac{\partial }{\partial t}\mathbb{P}({\theta }_{n+k}\le t,\phi ({\theta }_{n+k})=j|X(0)=n,\phi (0)=i)\end{array}$$(10) records the probability density that the process first visits level n+k at time t and does so in phase j, given that it starts from level n in phase i. Let $$\begin{array}{l}{\mathbf{H}}^{n,n+k}(s)=[H_{ij}^{n,n+k}(s){]}_{i=1,\dots ,m_n;j=1,\dots ,m_{n+k}}\end{array}$$ given by (11) $$\begin{array}{ll}{H}_{ij}^{n,n+k}(s)& ={\int }_0^{\mathrm{\infty }}{e}^{-st}{h}_{ij}^{n,n+k}(t)dt,\end{array}$$(11) be the corresponding LST matrix, and ${\mathbf{H}}^{n,n+k}={\mathbf{H}}^{n,n+k}(0)$ and ${\mathbb{E}}^{n,n+k}=-\frac{\partial }{\partial s}{\mathbf{H}}^{n,n+k}(s){|}_{s=0}$ be the probability and expectation matrix, respectively.

To evaluate ${\mathbf{H}}^{n,n+k}(s)$ and ${\mathbf{H}}^{n,n+k}$, we apply the following recursion summarized in Grant[5], (12) $$\begin{array}{ll}{\mathbf{H}}^{0,1}(s)& =-({\mathbf{Q}}^{[0,0]}-s\mathbf{I}{)}^{-1}{\mathbf{Q}}^{[0,1]},\end{array}$$(12) (13) $$\begin{array}{ll}{\mathbf{H}}^{n,n+1}(s)& =-({\mathbf{Q}}^{[n,n]}-s\mathbf{I}+{\mathbf{Q}}^{[n,n-1]}{\mathbf{H}}^{n-1,n}(s){)}^{-1}{\mathbf{Q}}^{[n,n+1]}\\ & \text{for }n=N-1,\dots ,1\end{array}$$(13) and then let (14) $$\begin{array}{ll}{\mathbf{H}}^{n,n+k}(s)& ={\mathbf{H}}^{n,n+1}(s){\mathbf{H}}^{n+1,n+2}(s)\times \cdots \times {\mathbf{H}}^{n+k-1,n+k}(s).\end{array}$$(14)

By the above recursion, it follows that (15) $$\begin{array}{ll}{\mathbb{E}}^{n,n+k}& ={\mathbb{E}}^{n,n+1}{\mathbf{H}}^{n+1,n+2}\times \cdots \times {\mathbf{H}}^{n+k-1,n+k}+\cdots \\ & +{\mathbf{H}}^{n,n+1}{\mathbf{H}}^{n+1,n+2}\times \cdots \times {\mathbf{H}}^{n+k-2,n+k-1}{\mathbb{E}}^{n+k-1,n+k}\\ & =\sum _{m=n}^{n+k-1}\left(\prod _{\ell =n}^{m-1}{\mathbf{H}}^{\ell ,\ell +1}\right){\mathbb{E}}^{m,m+1}\left(\prod _{\ell =m+1}^{n+k-1}{\mathbf{H}}^{\ell ,\ell +1}\right),\end{array}$$(15) where, by (12(12) $$\begin{array}{ll}{\mathbf{H}}^{0,1}(s)& =-({\mathbf{Q}}^{[0,0]}-s\mathbf{I}{)}^{-1}{\mathbf{Q}}^{[0,1]},\end{array}$$(12) ), analogly to the corresponding computations for (8(8) $$\begin{array}{ll}{\mathbb{E}}^{N,N-1}& =-({\mathbf{Q}}^{[N,N]}{)}^{-1}{\mathbf{G}}^{N,N-1}.\end{array}$$(8) ), (16) $$\begin{array}{ll}{\mathbb{E}}^{0,1}& =-({\mathbf{Q}}^{[0,0]}{)}^{-1}{\mathbf{H}}^{0,1}\end{array}$$(16) and by (13(13) $$\begin{array}{ll}{\mathbf{H}}^{n,n+1}(s)& =-({\mathbf{Q}}^{[n,n]}-s\mathbf{I}+{\mathbf{Q}}^{[n,n-1]}{\mathbf{H}}^{n-1,n}(s){)}^{-1}{\mathbf{Q}}^{[n,n+1]}\\ & \text{for }n=N-1,\dots ,1\end{array}$$(13) ), for $n=1,\dots ,N-1$, (17) $$\begin{array}{ll}{\mathbb{E}}^{n,n+1}& =-({\mathbf{Q}}^{[n,n]}+{\mathbf{Q}}^{[n,n-1]}{\mathbf{H}}^{n-1,n}{)}^{-1}(\mathbf{I}+{\mathbf{Q}}^{[n,n-1]}{\mathbb{E}}^{n-1,n}){\mathbf{G}}^{n,n+1}.\end{array}$$(17)

3.4. Sojourn times in specified sets

Suppose that the process starts from state (n, i) at time 0. Let $\mathcal{A}\subset \{0,1,\dots ,N\}$ be some set of desirable or undesirable levels, such as $\mathcal{A}=\{0,1,\dots ,A\}$ or $\mathcal{A}=\{B,B+1,\dots ,N\}$ where A and B are some thresholds.

Let ${\mathbf{G}}_{\mathcal{A}}^{n,n-k}(s)=[G_{\mathcal{A};ij}^{n,n-k}(s){]}_{i=1,\dots ,m_n;j=1,\dots ,m_{n-k}}$ be the LST matrix of the distribution of the total time spent in levels in the set 𝒜 during a sample path corresponding to $G_{ij}^{n,n-k}$. Let ${\mathbf{G}}_{\mathcal{A}}^{n,n-k}={\mathbf{G}}_{\mathcal{A}}^{n,n-k}(0)$ and ${\mathbb{E}}_{\mathcal{A}}^{n,n-k}=-\frac{\partial }{\partial s}{\mathbf{G}}_{\mathcal{A}}^{n,n-k}(s){|}_{s=0}$ be the probability and expectation matrix, respectively.

Similarly, let ${\mathbf{H}}_{\mathcal{A}}^{n,n+k}(s)=[H_{\mathcal{A};ij}^{n,n+k}(s){]}_{i=1,\dots ,m_n;j=1,\dots ,m_{n+k}}$ be the LST matrix of the distribution of the total time spent at levels in the set 𝒜 during a sample path corresponding to $H_{ij}^{n,n+k}$. Let ${\mathbf{H}}_{\mathcal{A}}^{n,n+k}={\mathbf{H}}_{\mathcal{A}}^{n,n+k}(0)$ and ${\mathbb{E}}_{\mathcal{A}}^{n,n+k}=-\frac{\partial }{\partial s}{\mathbf{H}}_{\mathcal{A}}^{n,n+k}(s){|}_{s=0}$ be the probability and expectation matrix, respectively.

By the decomposition of a sample path analogous to Section 3.2, we have, (18) $$\begin{array}{ll}{\mathbf{G}}_{\mathcal{A}}^{n,n-k}(s)& ={\mathbf{G}}_{\mathcal{A}}^{n,n-1}(s){\mathbf{G}}_{\mathcal{A}}^{n-1,n-2}(s)\times \cdots \times {\mathbf{G}}_{\mathcal{A}}^{n-k+1,n-k}(s),\end{array}$$(18) where (19) $$\begin{array}{ll}{\mathbf{G}}_{\mathcal{A}}^{N,N-1}(s)& =-({\mathbf{Q}}^{[N,N]}-s\mathbf{I}\times I(N\in \mathcal{A}){)}^{-1}{\mathbf{Q}}^{[N,N-1]},\end{array}$$(19) and for $n=N-1,\dots ,1$, (20) $$\begin{array}{ll}{\mathbf{G}}_{\mathcal{A}}^{n,n-1}(s)& =-({\mathbf{Q}}^{[n,n]}-s\mathbf{I}\times I(n\in \mathcal{A})+{\mathbf{Q}}^{[n,n+1]}{\mathbf{G}}_{\mathcal{A}}^{n+1,n}(s){)}^{-1}{\mathbf{Q}}^{[n,n-1]},\end{array}$$(20) where I(⋅) is an indicator function.

Similarly, by the decomposition of a sample path analogous to Section 3.2, we have, (21) $$\begin{array}{ll}{\mathbf{H}}_{\mathcal{A}}^{n,n+k}(s)& ={\mathbf{H}}_{\mathcal{A}}^{n,n+1}(s){\mathbf{H}}_{\mathcal{A}}^{n+1,n+2}(s)\times \cdots \times {\mathbf{H}}_{\mathcal{A}}^{n+k-1,n+k}(s).\end{array}$$(21) where (22) $$\begin{array}{ll}{\mathbf{H}}_{\mathcal{A}}^{0,1}(s)& =-({\mathbf{Q}}^{[0,0]}-s\mathbf{I}\times I(0\in \mathcal{A}){)}^{-1}{\mathbf{Q}}^{[0,1]},\end{array}$$(22) and for $n=N-1,\dots ,1$, (23) $$\begin{array}{ll}{\mathbf{H}}_{\mathcal{A}}^{n,n+1}(s)& =-({\mathbf{Q}}^{[n,n]}-s\mathbf{I}\times I(n\in \mathcal{A})+{\mathbf{Q}}^{[n,n-1]}{\mathbf{H}}^{n-1,n}(s){)}^{-1}{\mathbf{Q}}^{[n,n+1]}.\end{array}$$(23)

3.5. Distribution at time t

Suppose that the QBD starts at some level $n_0=0,1,2\dots ,N$ in some phase $i_0=0,1,\dots ,m_{n_0}$ according to the initial distribution of phases $\mathit{\alpha }=[{\alpha }_j{]}_{j\in 0,1,\dots ,m_{n_0}}$ such that ${\alpha }_j=\mathbb{P}(\phi (0)=j)$.

Define vector $\mathbf{f}(t)=[{\mathbf{f}}_n(t){]}_{n=0,1,\dots ,N}$ such that (24) $$\begin{array}{ll}[{\mathbf{f}}_n(t){]}_j& ={\mathbb{P}}_{\mathit{\alpha }}(X_t=n,\phi (t)=j)\end{array}$$(24) is the probability that at time t the process is on level n and in phase j, given 𝜶; and the corresponding Laplace Transform (LT) vector $\tilde{\mathbf{f}}(s)=[{\tilde{\mathbf{f}}}_n(s){]}_{n=0,1,\dots ,N}$ such that (25) $$\begin{array}{ll}[{\tilde{\mathbf{f}}}_n(s){]}_j& ={\int }_0^{\mathrm{\infty }}{e}^{-st}{\mathbb{P}}_{\mathit{\alpha }}(X_t=n,\phi (t)=j)dt.\end{array}$$(25)

The standard Kolmogorov differential equations (DEs) of the process are, (26) $$\begin{array}{ll}\frac{\partial {\mathbf{f}}_0(t)}{\partial t}& ={\mathbf{f}}_0(t){\mathbf{Q}}^{[0,0]}+{\mathbf{f}}_1(t){\mathbf{Q}}^{[1,0]},\\ \frac{\partial {\mathbf{f}}_n(t)}{\partial t}& ={\mathbf{f}}_n(t){\mathbf{Q}}^{[n,n]}+{\mathbf{f}}_{n-1}(t){\mathbf{Q}}^{[n-1,n]}+{\mathbf{f}}_{n+1}(t){\mathbf{Q}}^{[n+1,n]}\\ & \text{for }n=1,\dots ,N-1\\ \frac{\partial {\mathbf{f}}_N(t)}{\partial t}& ={\mathbf{f}}_N(t){\mathbf{Q}}^{[N,N]}+{\mathbf{f}}_{N-1}(t){\mathbf{Q}}^{[N-1,N]}\end{array}$$(26) with the initial condition ${\mathbf{f}}_{n_0}(0)=\mathit{\alpha }$, and 𝐟_n(0) = 𝟎 for n≠n₀. To evaluate $\tilde{\mathbf{f}}(s)$, we apply the following recursion summarized in Grant[5], (27) $$\begin{array}{ll}{\tilde{\mathbf{f}}}_{n_0}(s)& =-\mathit{\alpha }{(({\mathbf{Q}}^{[n_0,n_0]}-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n_0-1)}(s){\mathbf{Q}}^{[n_0-1,n_0]}+{\tilde{\mathbf{R}}}^{(n_0+1)}(s){\mathbf{Q}}^{[n_0+1,n_0]})}^{-1}\end{array}$$(27) and then for n > n₀, (28) $$\begin{array}{ll}{\tilde{\mathbf{f}}}_n(s)& =-\mathit{\alpha }{\mathbf{H}}^{n_0,n}(s){(({\mathbf{Q}}^{[n,n]}-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n-1)}(s){\mathbf{Q}}^{[n-1,n]}+{\tilde{\mathbf{R}}}^{(n+1)}(s){\mathbf{Q}}^{[n+1,n]})}^{-1}\end{array}$$(28) and for n < n₀, $$\begin{array}{ll}{\tilde{\mathbf{f}}}_n(s)& =-\mathit{\alpha }{\mathbf{G}}^{n_0,n}(s){(({\mathbf{Q}}^{[n,n]}-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n-1)}(s){\mathbf{Q}}^{[n-1,n]}+{\tilde{\mathbf{R}}}^{(n+1)}(s){\mathbf{Q}}^{[n+1,n]})}^{-1}\end{array}$$ with $$\begin{array}{ll}{\tilde{\mathbf{R}}}^{(n)}(s)& ={\tilde{R}}_n\circ {\tilde{R}}_{n+1}\circ \dots \circ {\tilde{R}}_N(\mathbf{0}),\\ {\tilde{R}}_n(\mathbf{X})& =-{\mathbf{Q}}^{[n-1,n]}({\mathbf{Q}}^{[n,n]}-s\mathbf{I}+\mathbf{X}{\mathbf{Q}}^{[n+1,n]}{)}^{-1}\end{array}$$ and $$\begin{array}{ll}{\hat{\mathbf{R}}}^{(n)}(s)& ={\hat{R}}_n\circ {\hat{R}}_{n-1}\circ \dots \circ {\hat{R}}_0(\mathbf{0}),\\ {\hat{R}}_n(\mathbf{X})& =-{\mathbf{Q}}^{[n+1,n]}(\mathbf{X}{\mathbf{Q}}^{[n-1,n]}+{\mathbf{Q}}^{[n,n]}-s\mathbf{I}{)}^{-1}.\end{array}$$

4. Sensitivity analysis

4.1. Examples

We consider simple examples to illustrate the theory. We state the expressions for the derivatives $\frac{\partial }{\partial \mathit{\theta }}\mathbf{Q}(\mathit{\theta })$ since the derivatives $\frac{\partial }{\partial \mathit{\theta }}$ of quantities of interest can be expressed in terms of these.

Example 1.

As a simple example, suppose that X(t) records the total number of customers in the system at time t. Let $\phi (t)\in \{1,\dots ,k\}$ be the phase of the environment that drives the evolution of the system, so that λ_φ(t) > 0 is the arrival rate to the system (provided X(t) < N), and μ_φ(t) > 0 is the service rate (provided X(t) > 0). Assume that $\{\phi (t):t\ge 0\}$ is a continuous-time Markov chain with generator 𝐓 = [T_ij]. The system can be modeled as a LD-QBD with generator $\mathbf{Q}(\mathit{\theta })=[q(\mathit{\theta }{)}_{(n,i)(m,j)}]$ that depends on the vector of parameters $\mathit{\theta }=[{\lambda }_1,\dots ,{\lambda }_k,{\mu }_1,\dots ,{\mu }_k]$, such that $$\begin{array}{ll}q(\mathit{\theta }{)}_{(n,i)(m,j)}& =\{\begin{array}{ll}{T}_{ij}& \text{when }j\ne i,m=n;\\ {\lambda }_i& \text{when }j=i,m=n+1,n<N;\\ n{\mu }_i& \text{when }j=i,m=n-1,n>0;\\ T_{ii}-{\lambda }_i-n{\mu }_i& \text{when }j=i,m=n,0<n<N;\\ T_{ii}-{\lambda }_i& \text{when }j=i,m=n=0;\\ T_{ii}-n{\mu }_i& \text{when }j=i,m=n=N;\\ 0& \text{otherwise.}\end{array}\end{array}$$

Then $\frac{\partial }{\partial \mathit{\theta }}\mathbf{Q}(\mathit{\theta })$ is given by $$\begin{array}{ll}\frac{\partial }{\partial {\lambda }_i}q(\mathit{\theta }{)}_{(n,i)(m,j)}& =\{\begin{array}{ll}1& \text{when }j=i,m=n+1,n<N;\\ -1& \text{when }j=i,m=n,0\le n<N;\\ 0& \text{otherwise;}\end{array}\end{array}$$ and $$\begin{array}{ll}\frac{\partial }{\partial {\mu }_i}q(\mathit{\theta }{)}_{(n,i)(m,j)}& =\{\begin{array}{ll}n& \text{when }j=i,m=n-1,n>0;\\ -n& \text{when }j=i,m=n,0<n\le N;\\ 0& \text{otherwise.}\end{array}\end{array}$$

Example 2.

Suppose that customers of type $k\in \{1,2\}$ arrive at the system with capacity N at the total rate λ_k > 0, and are served at rate μ_k > 0 per customer. The system can be modeled as a LD-QBD $\{(X(t),\phi (t)):t\ge 0\}$ where X(t) records the total number of customers and φ(t)≤X(t) records the number of customers of type 1 in the system at time t. Then the generator $\mathbf{Q}(\mathit{\theta })=[q(\mathit{\theta }{)}_{(n,i)(m,j)}]$ that depends on the vector of parameters $\mathit{\theta }=[{\lambda }_1,{\lambda }_2,{\mu }_1,{\mu }_2]$, is such that $$\begin{array}{ll}q(\mathit{\theta }{)}_{(n,i)(m,j)}& =\{\begin{array}{ll}{\lambda }_2& \text{when }j=i,m=n+1,n<N;\\ {\lambda }_1& \text{when }j=i+1,m=n+1,n<N;\\ (n-i){\mu }_2& \text{when }j=i,m=n-1,n>0;\\ i{\mu }_1& \text{when }j=i-1,m=n-1,n>0;\\ -({\lambda }_1+{\lambda }_2+(n-i){\mu }_2+i{\mu }_1)& \text{when }j=i,m=n,0<n<N;\\ -({\lambda }_1+{\lambda }_2)& \text{when }j=i,m=n=0;\\ -((n-i){\mu }_2+i{\mu }_1)& \text{when }j=i,m=n=N;\\ 0& \text{otherwise.}\end{array}\end{array}$$

Then $\frac{\partial }{\partial \mathit{\theta }}\mathbf{Q}(\mathit{\theta })$ is given by $$\begin{array}{ll}\frac{\partial }{\partial {\lambda }_1}q(\mathit{\theta }{)}_{(n,i)(m,j)}& =\{\begin{array}{ll}1& \text{when }j=i+1,m=n+1,n<N;\\ -1& \text{when }j=i,m=n,0\le n<N;\\ 0& \text{otherwise;}\end{array}\end{array}$$ and $$\begin{array}{ll}\frac{\partial }{\partial {\lambda }_2}q(\mathit{\theta }{)}_{(n,i)(m,j)}& =\{\begin{array}{ll}1& \text{when }j=i,m=n+1,n<N;\\ -1& \text{when }j=i,m=n,0\le n<N;\\ 0& \text{otherwise;}\end{array}\end{array}$$ and $$\begin{array}{ll}\frac{\partial }{\partial {\mu }_1}q(\mathit{\theta }{)}_{(n,i)(m,j)}& =\{\begin{array}{ll}i& \text{when }j=i-1,m=n-1,n>0;\\ -i& \text{when }j=i,m=n,0<n\le N;\\ 0& \text{otherwise.}\end{array}\end{array}$$ and $$\begin{array}{ll}\frac{\partial }{\partial {\mu }_2}q(\mathit{\theta }{)}_{(n,i)(m,j)}& =\{\begin{array}{ll}(n-i)& \text{when }j=i,m=n-1,n>0;\\ -(n-i)& \text{when }j=i,m=n,0<n\le N;\\ 0& \text{otherwise.}\end{array}\end{array}$$

Example 3.

In this example, we consider a model in which the generator is perturbed. Suppose that the generator of a LD-QBD is a function of $\mathit{ϵ}=[{ϵ}_i{]}_{i=1,\dots ,k}>\mathbf{0}$ such that (29) $$\begin{array}{ll}\mathbf{Q}(\mathit{ϵ})& =\mathbf{Q}+\sum _{i=1}^k{ϵ}_i\times {\tilde{\mathbf{Q}}}_i\end{array}$$(29) is a generator for sufficiently small $||\mathit{ϵ}||>\mathbf{0}$. Then $$\begin{array}{ll}\frac{\partial \mathbf{Q}(\mathit{ϵ})}{\partial \mathit{ϵ}}& =\left[\begin{array}{ccc}\frac{\partial \mathbf{Q}(\mathit{ϵ})}{\partial {ϵ}_1}& \dots & \frac{\partial \mathbf{Q}(\mathit{ϵ})}{\partial {ϵ}_k}\end{array}\right]=\left[\begin{array}{ccc}{\tilde{\mathbf{Q}}}_1& \dots & {\tilde{\mathbf{Q}}}_k\end{array}\right],\end{array}$$ and the derivatives $\frac{\partial }{\partial \mathit{ϵ}}$ of various quantities are computed at ϵ = 𝟎.

4.2. Stationary distribution

In computing sensitivities, we will use these key expressions for the matrix derivative calculus, namely the derivative of a matrix product in (30(30) $$\begin{array}{ll}\frac{\partial }{\partial \mathit{\theta }}(\mathbf{A}(\mathit{\theta })\times \mathbf{B}(\mathit{\theta }))& =\frac{\partial \mathbf{A}(\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes \mathbf{B}(\mathit{\theta }))+\mathbf{A}(\mathit{\theta })\times \frac{\partial \mathbf{B}(\mathit{\theta })}{\partial \mathit{\theta }}\end{array}$$(30) ) and the derivative of the matrix inverse in (31(31) $$\begin{array}{ll}\frac{\partial (\mathbf{A}(\mathit{\theta }){)}^{-1}}{\partial \mathit{\theta }}& =-(\mathbf{A}(\mathit{\theta }){)}^{-1}\times \frac{\partial \mathbf{A}(\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes (\mathbf{A}(\mathit{\theta }){)}^{-1})\end{array}$$(31) ) (see e.g.Kowal[9]): (30) $$\begin{array}{ll}\frac{\partial }{\partial \mathit{\theta }}(\mathbf{A}(\mathit{\theta })\times \mathbf{B}(\mathit{\theta }))& =\frac{\partial \mathbf{A}(\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes \mathbf{B}(\mathit{\theta }))+\mathbf{A}(\mathit{\theta })\times \frac{\partial \mathbf{B}(\mathit{\theta })}{\partial \mathit{\theta }}\end{array}$$(30) (31) $$\begin{array}{ll}\frac{\partial (\mathbf{A}(\mathit{\theta }){)}^{-1}}{\partial \mathit{\theta }}& =-(\mathbf{A}(\mathit{\theta }){)}^{-1}\times \frac{\partial \mathbf{A}(\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes (\mathbf{A}(\mathit{\theta }){)}^{-1})\end{array}$$(31) with 𝐈 denoting an identity matrix of an appropriate size.

To derive the stationary distribution sensitivity, we apply (30(30) $$\begin{array}{ll}\frac{\partial }{\partial \mathit{\theta }}(\mathbf{A}(\mathit{\theta })\times \mathbf{B}(\mathit{\theta }))& =\frac{\partial \mathbf{A}(\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes \mathbf{B}(\mathit{\theta }))+\mathbf{A}(\mathit{\theta })\times \frac{\partial \mathbf{B}(\mathit{\theta })}{\partial \mathit{\theta }}\end{array}$$(30) ) in taking $\frac{\partial }{\partial \mathit{\theta }}$ in (1(1) $$\begin{array}{ll}{\mathit{\pi }}_n& ={\mathit{\pi }}_{n+1}{\hat{\mathbf{R}}}^{(n)}={\mathit{\pi }}_N\prod _{k=N-1}^n{\hat{\mathbf{R}}}^{(k)}.\end{array}$$(1) ). That leads to the following recursion: $$\begin{array}{ll}\frac{\partial {\mathit{\pi }}_n(\mathit{\theta })}{\partial \mathit{\theta }}& =\frac{\partial {\mathit{\pi }}_{n+1}(\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes {\hat{\mathbf{R}}}^{(n)}(\mathit{\theta }))+{\mathit{\pi }}_{n+1}(\mathit{\theta })\times \frac{\partial {\hat{\mathbf{R}}}^{(n)}(\mathit{\theta })}{\partial \mathit{\theta }}.\end{array}$$

Then, by taking partial derivative $\frac{\partial }{\partial \mathit{\theta }}$ in the second equation of (2(2) $$\begin{array}{l}\phantom{.}\{\begin{array}{l}{\mathit{\pi }}_N({\hat{\mathbf{R}}}^{(N-1)}{\mathbf{Q}}^{[N-1,N]}+{\mathbf{Q}}^{[N,N]})=\mathbf{0}\\ {\mathit{\pi }}_N(\mathbf{1}+\sum _{n=0}^{N-1}\prod _{k=N-1}^n{\hat{\mathbf{R}}}^{(k)}\mathbf{1})=1.\end{array}\end{array}$$(2) ), we obtain $$\begin{array}{ll}\frac{\partial {\mathit{\pi }}_N(\mathit{\theta })}{\partial \mathit{\theta }}& =-{\mathit{\pi }}_N(\mathit{\theta })\times \frac{\partial ({\hat{\mathbf{R}}}^{(N-1)}(\mathit{\theta }){\mathbf{Q}}^{[N-1,N]}(\mathit{\theta })+{\mathbf{Q}}^{[N,N]}(\mathit{\theta }))}{\partial \mathit{\theta }}\\ & \times ({\mathbf{I}}_k\otimes {({\hat{\mathbf{R}}}^{(N-1)}(\mathit{\theta }){\mathbf{Q}}^{[N-1,N]}(\mathit{\theta })+{\mathbf{Q}}^{[N,N]}(\mathit{\theta }))}^{-1}).\end{array}$$

In the above expressions, we evaluate the derivatives of ${\hat{\mathbf{R}}}^{(n)}(s)$ which are of the form $$\begin{array}{ll}\frac{\partial }{\partial \mathit{\theta }}{\hat{\mathbf{R}}}^{(0)}(s)& =\frac{\partial }{\partial \mathit{\theta }}(-{\mathbf{Q}}^{[1,0]}({\mathbf{Q}}^{[0,0]}-s\mathbf{I}{)}^{-1}),\end{array}$$ and $$\begin{array}{ll}\frac{\partial }{\partial \mathit{\theta }}{\hat{\mathbf{R}}}^{(n)}(s)& =\frac{\partial }{\partial \mathit{\theta }}(-{\mathbf{Q}}^{[n+1,n]}({\hat{\mathbf{R}}}^{(n-1)}(s){\mathbf{Q}}^{[n-1,n]}+{\mathbf{Q}}^{[n,n]}-s\mathbf{I}{)}^{-1}).\end{array}$$

4.3. First hitting times to lower levels

Using the results of Section 3.2 and the expressions in (5(5) $$\begin{array}{ll}{\mathbf{G}}^{N,N-1}(s)& =-({\mathbf{Q}}^{[N,N]}-s\mathbf{I}{)}^{-1}{\mathbf{Q}}^{[N,N-1]},\end{array}$$(5) )–(7(7) $$\begin{array}{ll}{\mathbf{G}}^{n,n-k}(s)& ={\mathbf{G}}^{n,n-1}(s){\mathbf{G}}^{n-1,n-2}(s)\times \cdots \times {\mathbf{G}}^{n-k+1,n-k}(s).\end{array}$$(7) ), to evaluate $\frac{\partial }{\partial \mathit{\theta }}{\mathbf{G}}^{n,n-k}(\mathit{\theta })$, we obtain the following recursion. We have: $$\begin{array}{ll}\frac{\partial }{\partial \mathit{\theta }}{\mathbf{G}}^{N,N-1}(\mathit{\theta })& =-\frac{\partial }{\partial \mathit{\theta }}(({\mathbf{Q}}^{[N,N]}(\mathit{\theta }){)}^{-1}{\mathbf{Q}}^{[N,N-1]}(\mathit{\theta }))\\ & =-\frac{\partial ({\mathbf{Q}}^{[N,N]}(\mathit{\theta }){)}^{-1}}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[N,N-1]}(\mathit{\theta }))\\ & -({\mathbf{Q}}^{[N,N]}(\mathit{\theta }){)}^{-1}\times \frac{\partial {\mathbf{Q}}^{[N,N-1]}(\mathit{\theta })}{\partial \mathit{\theta }}\\ & =({\mathbf{Q}}^{[N,N]}(\mathit{\theta }){)}^{-1}\times \frac{\partial {\mathbf{Q}}^{[N,N]}(\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes ({\mathbf{Q}}^{[N,N]}(\mathit{\theta }){)}^{-1})\\ & \times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[N,N-1]}(\mathit{\theta }))-({\mathbf{Q}}^{[N,N]}(\mathit{\theta }){)}^{-1}\times \frac{\partial {\mathbf{Q}}^{[N,N-1]}(\mathit{\theta })}{\partial \mathit{\theta }},\end{array}$$ and for $n=N-1,\dots ,1$, $$\begin{array}{ll}\frac{\partial }{\partial \mathit{\theta }}{\mathbf{G}}^{n,n-1}(\mathit{\theta })& =-\frac{\partial }{\partial \mathit{\theta }}(({\mathbf{Q}}^{[n,n]}(\mathit{\theta })+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta }){\mathbf{G}}^{n+1,n}(\mathit{\theta }){)}^{-1}{\mathbf{Q}}^{[n,n-1]}(\mathit{\theta }))\\ & =-\frac{\partial ({\mathbf{Q}}^{[n,n]}(\mathit{\theta })+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta }){\mathbf{G}}^{n+1,n}(\mathit{\theta }){)}^{-1}}{\partial \mathit{\theta }}\\ & \times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[n,n-1]}(\mathit{\theta }))-({\mathbf{Q}}^{[n,n]}(\mathit{\theta })+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta }){\mathbf{G}}^{n+1,n}(\mathit{\theta }){)}^{-1}\\ & \times \frac{\partial {\mathbf{Q}}^{[n,n-1]}(\mathit{\theta })}{\partial \mathit{\theta }}\end{array}$$ where (32) $$\begin{array}{ll}& \frac{\partial ({\mathbf{Q}}^{[n,n]}(\mathit{\theta })+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta }){\mathbf{G}}^{n+1,n}(\mathit{\theta }){)}^{-1}}{\partial \mathit{\theta }}\\ & =-({\mathbf{Q}}^{[n,n]}(\mathit{\theta })+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta }){\mathbf{G}}^{n+1,n}(\mathit{\theta }){)}^{-1}\\ & \times \frac{\partial ({\mathbf{Q}}^{[n,n]}(\mathit{\theta })+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta }){\mathbf{G}}^{n+1,n}(\mathit{\theta }))}{\partial \mathit{\theta }}\\ & \times ({\mathbf{I}}_k\otimes ({\mathbf{Q}}^{[n,n]}(\mathit{\theta })+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta }){\mathbf{G}}^{n+1,n}(\mathit{\theta }){)}^{-1})\\ & =-({\mathbf{Q}}^{[n,n]}(\mathit{\theta })+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta }){\mathbf{G}}^{n+1,n}(\mathit{\theta }){)}^{-1}\\ & \times (\frac{\partial ({\mathbf{Q}}^{[n,n]}(\mathit{\theta })}{\partial \mathit{\theta }}+\frac{\partial {\mathbf{Q}}^{[n,n+1]}(\mathit{\theta })}{\partial \mathit{\theta }}\\ & \times ({\mathbf{I}}_k\otimes {\mathbf{G}}^{n+1,n}(\mathit{\theta }))+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta })\times \frac{\partial {\mathbf{G}}^{n+1,n}(\mathit{\theta })}{\partial \mathit{\theta }})\\ & \times ({\mathbf{I}}_k\otimes ({\mathbf{Q}}^{[n,n]}(\mathit{\theta })+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta }){\mathbf{G}}^{n+1,n}(\mathit{\theta }){)}^{-1})\end{array}$$(32) and for k≥2 we obtain the recursion, $$\begin{array}{ll}\frac{\partial }{\partial \mathit{\theta }}{\mathbf{G}}^{n,n-k}(\mathit{\theta })& =\frac{\partial }{\partial \mathit{\theta }}({\mathbf{G}}^{n,n-1}(\mathit{\theta }){\mathbf{G}}^{n-1,n-2}(\mathit{\theta })\times \cdots \times {\mathbf{G}}^{n-k+1,n-k}(\mathit{\theta }))\\ & =\frac{\partial }{\partial \mathit{\theta }}({\mathbf{G}}^{n,n-k+1}(\mathit{\theta }){\mathbf{G}}^{n-k+1,n-k}(\mathit{\theta }))\\ & =\frac{\partial {\mathbf{G}}^{n,n-k+1}(\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes {\mathbf{G}}^{n-k+1,n-k}(\mathit{\theta }))\\ & +{\mathbf{G}}^{n,n-k+1}(\mathit{\theta })\times \frac{\partial {\mathbf{G}}^{n-k+1,n-k}(\mathit{\theta })}{\partial \mathit{\theta }}.\end{array}$$

Further, $$\begin{array}{ll}\frac{\partial {\mathbb{E}}^{N,N-1}(\mathit{\theta })}{\partial \mathit{\theta }}& =-\frac{\partial ({\mathbf{Q}}^{[N,N]}(\mathit{\theta }){)}^{-1}{\mathbf{G}}^{N,N-1}(\mathit{\theta })}{\partial \mathit{\theta }}\\ & =({\mathbf{Q}}^{[N,N]}(\mathit{\theta }){)}^{-1}\times \frac{\partial {\mathbf{Q}}^{[N,N]}(\mathit{\theta })}{\partial \mathit{\theta }}\\ & \times ({\mathbf{I}}_k\otimes ({\mathbf{Q}}^{[N,N]}(\mathit{\theta }){)}^{-1})\times ({\mathbf{I}}_k\otimes {\mathbf{G}}^{[N,N-1]}(\mathit{\theta }))\\ & -({\mathbf{Q}}^{[N,N]}(\mathit{\theta }){)}^{-1}\times \frac{\partial {\mathbf{G}}^{[N,N-1]}(\mathit{\theta })}{\partial \mathit{\theta }},\end{array}$$ and for $n=N-1,\dots ,1$, $$\begin{array}{ll}\frac{\partial {\mathbb{E}}^{n,n-1}(\mathit{\theta })}{\partial \mathit{\theta }}& =-\frac{\begin{array}{c}\partial ({\mathbf{Q}}^{[n,n]}(\mathit{\theta })+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta }){\mathbf{G}}^{n+1,n}(\mathit{\theta }){)}^{-1}\\ (\mathbf{I}+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta }){\mathbb{E}}^{n+1,n}(\mathit{\theta })){\mathbf{G}}^{n,n-1}(\mathit{\theta })\end{array}}{\partial \mathit{\theta }}\\ & =-\frac{\partial ({\mathbf{Q}}^{[n,n]}(\mathit{\theta })+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta }){\mathbf{G}}^{n+1,n}(\mathit{\theta }){)}^{-1}}{\partial \mathit{\theta }}\\ & \times ({\mathbf{I}}_k\otimes (\mathbf{I}-{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta }){\mathbb{E}}^{n+1,n}(\mathit{\theta })){\mathbf{G}}^{n,n-1}(\mathit{\theta }))\\ & +({\mathbf{Q}}^{[n,n]}(\mathit{\theta })+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta }){\mathbf{G}}^{n+1,n}(\mathit{\theta }){)}^{-1}\\ & \times \frac{\partial (\mathbf{I}+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta }){\mathbb{E}}^{n+1,n}(\mathit{\theta })){\mathbf{G}}^{n,n-1}(\mathit{\theta })}{\partial \mathit{\theta }},\end{array}$$ where the first derivative in the above is given in (32(32) $$\begin{array}{ll}& \frac{\partial ({\mathbf{Q}}^{[n,n]}(\mathit{\theta })+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta }){\mathbf{G}}^{n+1,n}(\mathit{\theta }){)}^{-1}}{\partial \mathit{\theta }}\\ & =-({\mathbf{Q}}^{[n,n]}(\mathit{\theta })+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta }){\mathbf{G}}^{n+1,n}(\mathit{\theta }){)}^{-1}\\ & \times \frac{\partial ({\mathbf{Q}}^{[n,n]}(\mathit{\theta })+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta }){\mathbf{G}}^{n+1,n}(\mathit{\theta }))}{\partial \mathit{\theta }}\\ & \times ({\mathbf{I}}_k\otimes ({\mathbf{Q}}^{[n,n]}(\mathit{\theta })+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta }){\mathbf{G}}^{n+1,n}(\mathit{\theta }){)}^{-1})\\ & =-({\mathbf{Q}}^{[n,n]}(\mathit{\theta })+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta }){\mathbf{G}}^{n+1,n}(\mathit{\theta }){)}^{-1}\\ & \times (\frac{\partial ({\mathbf{Q}}^{[n,n]}(\mathit{\theta })}{\partial \mathit{\theta }}+\frac{\partial {\mathbf{Q}}^{[n,n+1]}(\mathit{\theta })}{\partial \mathit{\theta }}\\ & \times ({\mathbf{I}}_k\otimes {\mathbf{G}}^{n+1,n}(\mathit{\theta }))+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta })\times \frac{\partial {\mathbf{G}}^{n+1,n}(\mathit{\theta })}{\partial \mathit{\theta }})\\ & \times ({\mathbf{I}}_k\otimes ({\mathbf{Q}}^{[n,n]}(\mathit{\theta })+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta }){\mathbf{G}}^{n+1,n}(\mathit{\theta }){)}^{-1})\end{array}$$(32) ), the second derivative is given by $$\begin{array}{ll}& \frac{\partial (\mathbf{I}+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta }){\mathbb{E}}^{n+1,n}(\mathit{\theta })){\mathbf{G}}^{n,n-1}(\mathit{\theta })}{\partial \mathit{\theta }}\\ & =\frac{\partial ({\mathbf{Q}}^{[n,n+1]}(\mathit{\theta }){\mathbb{E}}^{n+1,n}(\mathit{\theta }))}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes {\mathbf{G}}^{n,n-1}(\mathit{\theta }))\\ & +(\mathbf{I}+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta }){\mathbb{E}}^{n+1,n}(\mathit{\theta }))\times \frac{\partial {\mathbf{G}}^{n,n-1}(\mathit{\theta })}{\partial \mathit{\theta }}\\ & =(\frac{\partial {\mathbf{Q}}^{[n,n+1]}(\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes {\mathbb{E}}^{n+1,n}(\mathit{\theta }))+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta })\times \frac{\partial {\mathbb{E}}^{n+1,n}(\mathit{\theta })}{\partial \mathit{\theta }})\\ & \times ({\mathbf{I}}_k\otimes {\mathbf{G}}^{n,n-1}(\mathit{\theta }))\\ & +(\mathbf{I}+{\mathbf{Q}}^{[n,n+1]}(\mathit{\theta }){\mathbb{E}}^{n+1,n}(\mathit{\theta }))\times \frac{\partial {\mathbf{G}}^{n,n-1}(\mathit{\theta })}{\partial \mathit{\theta }},\end{array}$$ and for k≥2 we obtain the recursion, $$\begin{array}{ll}& \frac{\partial }{\partial \mathit{\theta }}{\mathbb{E}}^{n,n-k}(\mathit{\theta })\\ & =\frac{\partial }{\partial \mathit{\theta }}({\mathbb{E}}^{n,n-1}(\mathit{\theta }){\mathbf{G}}^{n-1,n-2}(\mathit{\theta })\times \cdots \times {\mathbf{G}}^{n-k+1,n-k}(\mathit{\theta })+\cdots \\ & +{\mathbf{G}}^{n,n-1}(\mathit{\theta }){\mathbf{G}}^{n-1,n-2}(\mathit{\theta })\times \cdots \times {\mathbf{G}}^{n-k+2,n-k+1}(\mathit{\theta }){\mathbb{E}}^{n-k+1,n-k}(\mathit{\theta }))\\ & =\frac{\partial }{\partial \mathit{\theta }}({\mathbb{E}}^{n,n-k+1}(\mathit{\theta }){\mathbb{G}}^{n-k+1,n-k}(\mathit{\theta })+{\mathbb{G}}^{n,n-k+1}(\mathit{\theta }){\mathbb{E}}^{n-k+1,n-k}(\mathit{\theta }))\\ & =\frac{\partial {\mathbb{E}}^{n,n-k+1}(\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes {\mathbb{G}}^{n-k+1,n-k}(\mathit{\theta }))+{\mathbb{E}}^{n,n-k+1}(\mathit{\theta })\\ & \times \frac{\partial {\mathbb{G}}^{n-k+1,n-k}(\mathit{\theta })}{\partial \mathit{\theta }}+\frac{\partial {\mathbb{G}}^{n,n-k+1}(\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes {\mathbb{E}}^{n-k+1,n-k}(\mathit{\theta }))\\ & +{\mathbb{G}}^{n,n-k+1}(\mathit{\theta })\times \frac{\partial {\mathbb{E}}^{n-k+1,n-k}(\mathit{\theta })}{\partial \mathit{\theta }}.\end{array}$$

Remark 2.

The computations of sensitivities in this section can be easily adapted to compute the following:

1. First hitting times to higher levels (see Subsection 3.3); we obtain similar expressions.

2. Sojourn times in specified sets (see Subsection 3.4); again, the analysis is similar.

4.4. Distribution at time t

In order to compute the distribution at time t, we shall apply the Laplace transform. We note that: (33) $$\begin{array}{ll}{\int }_0^{\mathrm{\infty }}{e}^{-st}\frac{\partial }{\partial \mathit{\theta }}{\mathbf{f}}_n(t;\mathit{\theta })dt& =\frac{\partial }{\partial \mathit{\theta }}{\int }_0^{\mathrm{\infty }}{e}^{-st}{\mathbf{f}}_n(t;\mathit{\theta })dt=\frac{\partial }{\partial \mathit{\theta }}{\tilde{\mathbf{f}}}_n(s;\mathit{\theta }).\end{array}$$(33)

The right-hand side of the above expression can be computed using the results from the earlier sections, and then inverted to obtain the quantity of interest $\frac{\partial }{\partial \mathit{\theta }}{\mathbf{f}}_n(t;\mathit{\theta })$.

As example, for n = n₀, by (27(27) $$\begin{array}{ll}{\tilde{\mathbf{f}}}_{n_0}(s)& =-\mathit{\alpha }{(({\mathbf{Q}}^{[n_0,n_0]}-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n_0-1)}(s){\mathbf{Q}}^{[n_0-1,n_0]}+{\tilde{\mathbf{R}}}^{(n_0+1)}(s){\mathbf{Q}}^{[n_0+1,n_0]})}^{-1}\end{array}$$(27) ) we have: $$\begin{array}{ll}& {\tilde{\mathbf{f}}}_{n_0}(s;\mathit{\theta })(({\mathbf{Q}}^{[n_0,n_0]}(\mathit{\theta })-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n_0-1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n_0-1,n_0]}(\mathit{\theta })\\ & +{\tilde{\mathbf{R}}}^{(n_0+1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n_0+1,n_0]}(\mathit{\theta }))=-\mathit{\alpha },\end{array}$$ and hence, the following sum is equal to zero, $$\begin{array}{ll}& \frac{\partial {\tilde{\mathbf{f}}}_{n_0}(s;\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes (({\mathbf{Q}}^{[n_0,n_0]}(\mathit{\theta })-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n_0-1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n_0-1,n_0]}(\mathit{\theta })\\ & +{\tilde{\mathbf{R}}}^{(n_0+1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n_0+1,n_0]}(\mathit{\theta })))\\ & +{\tilde{\mathbf{f}}}_{n_0}(s;\mathit{\theta })\times \frac{\begin{array}{c}\partial (({\mathbf{Q}}^{[n_0,n_0]}(\mathit{\theta })-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n_0-1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n_0-1,n_0]}(\mathit{\theta })\\ +{\tilde{\mathbf{R}}}^{(n_0+1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n_0+1,n_0]}(\mathit{\theta }).)\end{array}}{\partial \mathit{\theta }}\end{array}$$

This gives, $$\begin{array}{ll}\frac{\partial {\tilde{\mathbf{f}}}_{n_0}(s;\mathit{\theta })}{\partial \mathit{\theta }}& =-{\tilde{\mathbf{f}}}_{n_0}(s;\mathit{\theta })\\ & \times \frac{\partial }{\partial \mathit{\theta }}(({\mathbf{Q}}^{[n_0,n_0]}(\mathit{\theta })-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n_0-1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n_0-1,n_0]}(\mathit{\theta })\\ & +{\tilde{\mathbf{R}}}^{(n_0+1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n_0+1,n_0]}(\mathit{\theta }))\\ & \times ({\mathbf{I}}_k\otimes (({\mathbf{Q}}^{[n_0,n_0]}(\mathit{\theta })-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n_0-1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n_0-1,n_0]}(\mathit{\theta })\\ & +{{\tilde{\mathbf{R}}}^{(n_0+1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n_0+1,n_0]}(\mathit{\theta }))}^{-1})\end{array}$$ where the right-hand side can be readily evaluated.

Alternatively, directly from (27(27) $$\begin{array}{ll}{\tilde{\mathbf{f}}}_{n_0}(s)& =-\mathit{\alpha }{(({\mathbf{Q}}^{[n_0,n_0]}-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n_0-1)}(s){\mathbf{Q}}^{[n_0-1,n_0]}+{\tilde{\mathbf{R}}}^{(n_0+1)}(s){\mathbf{Q}}^{[n_0+1,n_0]})}^{-1}\end{array}$$(27) ), $$\begin{array}{ll}\frac{\partial }{\partial \mathit{\theta }}{\tilde{\mathbf{f}}}_{n_0}(s)& =-\mathit{\alpha }\frac{\partial }{\partial \mathit{\theta }}(({\mathbf{Q}}^{[n_0,n_0]}(\mathit{\theta })-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n_0-1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n_0-1,n_0]}(\mathit{\theta })\\ & +{{\tilde{\mathbf{R}}}^{(n_0+1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n_0+1,n_0]}(\mathit{\theta }))}^{-1}\\ & =\mathit{\alpha }(({\mathbf{Q}}^{[n_0,n_0]}(\mathit{\theta })-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n_0-1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n_0-1,n_0]}(\mathit{\theta })\\ & +{{\tilde{\mathbf{R}}}^{(n_0+1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n_0+1,n_0]}(\mathit{\theta }))}^{-1}\\ & \times \frac{\partial }{\partial \mathit{\theta }}(({\mathbf{Q}}^{[n_0,n_0]}(\mathit{\theta })-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n_0-1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n_0-1,n_0]}(\mathit{\theta })\\ & +{\tilde{\mathbf{R}}}^{(n_0+1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n_0+1,n_0]}(\mathit{\theta }))\\ & \times ({\mathbf{I}}_k\otimes (({\mathbf{Q}}^{[n_0,n_0]}(\mathit{\theta })-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n_0-1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n_0-1,n_0]}(\mathit{\theta })\\ & +{{\tilde{\mathbf{R}}}^{(n_0+1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n_0+1,n_0]}(\mathit{\theta }))}^{-1})\end{array}$$ and so we get the same form, since the right-hand side of second line is just $-{\tilde{\mathbf{f}}}_{n_0}(s;\mathit{\theta })$.

For n > n₀, by (28(28) $$\begin{array}{ll}{\tilde{\mathbf{f}}}_n(s)& =-\mathit{\alpha }{\mathbf{H}}^{n_0,n}(s){(({\mathbf{Q}}^{[n,n]}-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n-1)}(s){\mathbf{Q}}^{[n-1,n]}+{\tilde{\mathbf{R}}}^{(n+1)}(s){\mathbf{Q}}^{[n+1,n]})}^{-1}\end{array}$$(28) ), the expression $-\mathit{\alpha }{\mathbf{H}}^{n_0,n}(s;\mathit{\theta })$ is equal to $$\begin{array}{l}{\tilde{\mathbf{f}}}_n(s;\mathit{\theta })(({\mathbf{Q}}^{[n,n]}(\mathit{\theta })-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n-1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n-1,n]}(\mathit{\theta })+{\tilde{\mathbf{R}}}^{(n+1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n+1,n]}(\mathit{\theta }))\end{array}$$ and therefore, the expression $-\mathit{\alpha }\frac{\partial {\mathbf{H}}^{n_0,n}(s;\mathit{\theta })}{\partial \mathit{\theta }}$ is equal to the sum: $$\begin{array}{ll}& \frac{\partial {\tilde{\mathbf{f}}}_n(s;\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes (({\mathbf{Q}}^{[n,n]}(\mathit{\theta })-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n-1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n-1,n]}(\mathit{\theta })\\ & +{\tilde{\mathbf{R}}}^{(n+1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n+1,n]}(\mathit{\theta })))\\ & +{\tilde{\mathbf{f}}}_n(s;\mathit{\theta })\times \frac{\partial }{\partial \mathit{\theta }}(({\mathbf{Q}}^{[n,n]}(\mathit{\theta })-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n-1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n-1,n]}(\mathit{\theta })\\ & +{\tilde{\mathbf{R}}}^{(n+1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n+1,n]}(\mathit{\theta })).\end{array}$$

By the above, we obtain that the sensitivity $\frac{\partial {\tilde{\mathbf{f}}}_n(s;\mathit{\theta })}{\partial \mathit{\theta }}$ can be expressed as $$\begin{array}{ll}& (-\mathit{\alpha }\frac{\partial {\mathbf{H}}^{n_0,n}(s;\mathit{\theta })}{\partial \mathit{\theta }}-{\tilde{\mathbf{f}}}_n(s;\mathit{\theta })\\ & \times \frac{\partial }{\partial \mathit{\theta }}(({\mathbf{Q}}^{[n,n]}(\mathit{\theta })-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n-1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n-1,n]}(\mathit{\theta })\\ & +{\tilde{\mathbf{R}}}^{(n+1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n+1,n]}(\mathit{\theta }))\\ & \times ({\mathbf{I}}_k\otimes (({\mathbf{Q}}^{[n,n]}(\mathit{\theta })-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n-1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n-1,n]}(\mathit{\theta })\\ & +{{\tilde{\mathbf{R}}}^{(n+1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n+1,n]}(\mathit{\theta }))}^{-1}).\end{array}$$

Alternatively, directly from (28(28) $$\begin{array}{ll}{\tilde{\mathbf{f}}}_n(s)& =-\mathit{\alpha }{\mathbf{H}}^{n_0,n}(s){(({\mathbf{Q}}^{[n,n]}-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n-1)}(s){\mathbf{Q}}^{[n-1,n]}+{\tilde{\mathbf{R}}}^{(n+1)}(s){\mathbf{Q}}^{[n+1,n]})}^{-1}\end{array}$$(28) ), $\frac{\partial }{\partial \mathit{\theta }}{\tilde{\mathbf{f}}}_n(s;\mathit{\theta })$ is equal to: $$\begin{array}{ll}& -\mathit{\alpha }\frac{\partial }{\partial \mathit{\theta }}{\mathbf{H}}^{n_0,n}(s;\mathit{\theta })\times ({\mathbf{I}}_k\otimes (({\mathbf{Q}}^{[n,n]}(\mathit{\theta })-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n-1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n-1,n]}(\mathit{\theta })\\ & +{{\tilde{\mathbf{R}}}^{(n+1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n+1,n]}(\mathit{\theta }))}^{-1})\\ & -\mathit{\alpha }{\mathbf{H}}^{n_0,n}(s;\mathit{\theta })\times \frac{\partial }{\partial \mathit{\theta }}(({\mathbf{Q}}^{[n,n]}(\mathit{\theta })-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n-1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n-1,n]}(\mathit{\theta })\\ & +{{\tilde{\mathbf{R}}}^{(n+1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n+1,n]}(\mathit{\theta }))}^{-1}\end{array}$$ which gives, $$\begin{array}{ll}& \frac{\partial }{\partial \mathit{\theta }}{\tilde{\mathbf{f}}}_n(s;\mathit{\theta })\\ & =-\mathit{\alpha }\frac{\partial }{\partial \mathit{\theta }}{\mathbf{H}}^{n_0,n}(s;\mathit{\theta })\\ & \times ({\mathbf{I}}_k\otimes {(({\mathbf{Q}}^{[n,n]}(\mathit{\theta })-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n-1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n-1,n]}(\mathit{\theta })+{\tilde{\mathbf{R}}}^{(n+1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n+1,n]}(\mathit{\theta }))}^{-1})\\ & +\mathit{\alpha }{\mathbf{H}}^{n_0,n}(s;\mathit{\theta })\\ & \times {(({\mathbf{Q}}^{[n,n]}(\mathit{\theta })-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n-1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n-1,n]}(\mathit{\theta })+{\tilde{\mathbf{R}}}^{(n+1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n+1,n]}(\mathit{\theta }))}^{-1}\\ & \times \frac{\partial }{\partial \mathit{\theta }}(({\mathbf{Q}}^{[n,n]}(\mathit{\theta })-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n-1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n-1,n]}(\mathit{\theta })+{\tilde{\mathbf{R}}}^{(n+1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n+1,n]}(\mathit{\theta }))\\ & \times ({\mathbf{I}}_k\otimes {(({\mathbf{Q}}^{[n,n]}(\mathit{\theta })-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n-1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n-1,n]}(\mathit{\theta })+{\tilde{\mathbf{R}}}^{(n+1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n+1,n]}(\mathit{\theta }))}^{-1})\\ & =-\mathit{\alpha }\frac{\partial }{\partial \mathit{\theta }}{\mathbf{H}}^{n_0,n}(s;\mathit{\theta })\\ & \times ({\mathbf{I}}_k\otimes {(({\mathbf{Q}}^{[n,n]}(\mathit{\theta })-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n-1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n-1,n]}(\mathit{\theta })+{\tilde{\mathbf{R}}}^{(n+1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n+1,n]}(\mathit{\theta }))}^{-1})\\ & -{\tilde{\mathbf{f}}}_n(s;\mathit{\theta })\\ & \times \frac{\partial }{\partial \mathit{\theta }}(({\mathbf{Q}}^{[n,n]}(\mathit{\theta })-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n-1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n-1,n]}(\mathit{\theta })+{\tilde{\mathbf{R}}}^{(n+1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n+1,n]}(\mathit{\theta }))\\ & \times ({\mathbf{I}}_k\otimes {(({\mathbf{Q}}^{[n,n]}(\mathit{\theta })-s\mathbf{I})+{\hat{\mathbf{R}}}^{(n-1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n-1,n]}(\mathit{\theta })+{\tilde{\mathbf{R}}}^{(n+1)}(s;\mathit{\theta }){\mathbf{Q}}^{[n+1,n]}(\mathit{\theta }))}^{-1})\end{array}$$ and hence, again, we derive the same result.

4.5. Distribution at time t via the DEs (26)

We consider the DEs, as another approach to the analysis in Subsection 4.4.

4.5.1. General case

By taking the derivatives $\frac{\partial }{\partial \mathit{\theta }}$ in the DEs in (26(26) $$\begin{array}{ll}\frac{\partial {\mathbf{f}}_0(t)}{\partial t}& ={\mathbf{f}}_0(t){\mathbf{Q}}^{[0,0]}+{\mathbf{f}}_1(t){\mathbf{Q}}^{[1,0]},\\ \frac{\partial {\mathbf{f}}_n(t)}{\partial t}& ={\mathbf{f}}_n(t){\mathbf{Q}}^{[n,n]}+{\mathbf{f}}_{n-1}(t){\mathbf{Q}}^{[n-1,n]}+{\mathbf{f}}_{n+1}(t){\mathbf{Q}}^{[n+1,n]}\\ & \text{for }n=1,\dots ,N-1\\ \frac{\partial {\mathbf{f}}_N(t)}{\partial t}& ={\mathbf{f}}_N(t){\mathbf{Q}}^{[N,N]}+{\mathbf{f}}_{N-1}(t){\mathbf{Q}}^{[N-1,N]}\end{array}$$(26) ), we have $$\begin{array}{ll}\frac{\partial }{\partial \mathit{\theta }}\frac{\partial {\mathbf{f}}_0(t;\mathit{\theta })}{\partial t}& =\frac{\partial }{\partial \mathit{\theta }}{\mathbf{f}}_0(t;\mathit{\theta }){\mathbf{Q}}^{[0,0]}(\mathit{\theta })+\frac{\partial }{\partial \mathit{\theta }}{\mathbf{f}}_1(t;\mathit{\theta }){\mathbf{Q}}^{[1,0]}(\mathit{\theta })\\ \frac{\partial }{\partial \mathit{\theta }}\frac{\partial {\mathbf{f}}_n(t;\mathit{\theta })}{\partial t}& =\frac{\partial }{\partial \mathit{\theta }}{\mathbf{f}}_n(t;\mathit{\theta }){\mathbf{Q}}^{[n,n]}(\mathit{\theta })+\frac{\partial }{\partial \mathit{\theta }}{\mathbf{f}}_{n-1}(t;\mathit{\theta }){\mathbf{Q}}^{[n-1,n]}(\mathit{\theta })\\ & +\frac{\partial }{\partial \mathit{\theta }}{\mathbf{f}}_{n+1}(t;\mathit{\theta }){\mathbf{Q}}^{[n+1,n]}(\mathit{\theta })\text{for }n=1,\dots ,N-1,\\ \frac{\partial }{\partial \mathit{\theta }}\frac{\partial {\mathbf{f}}_N(t;\mathit{\theta })}{\partial t}& =\frac{\partial }{\partial \mathit{\theta }}{\mathbf{f}}_N(t;\mathit{\theta }){\mathbf{Q}}^{[N,N]}(\mathit{\theta })+\frac{\partial }{\partial \mathit{\theta }}{\mathbf{f}}_{N-1}(t;\mathit{\theta }){\mathbf{Q}}^{[N-1,N]}(\mathit{\theta })\end{array}$$ with the initial condition ${\mathbf{f}}_{n_0}(0;\mathit{\theta })=\mathit{\alpha }$, and ${\mathbf{f}}_n(0;\mathit{\theta })=\mathbf{0}$ for n≠n₀.

Next, we apply expressions from the matrix calculus (30(30) $$\begin{array}{ll}\frac{\partial }{\partial \mathit{\theta }}(\mathbf{A}(\mathit{\theta })\times \mathbf{B}(\mathit{\theta }))& =\frac{\partial \mathbf{A}(\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes \mathbf{B}(\mathit{\theta }))+\mathbf{A}(\mathit{\theta })\times \frac{\partial \mathbf{B}(\mathit{\theta })}{\partial \mathit{\theta }}\end{array}$$(30) )–(31(31) $$\begin{array}{ll}\frac{\partial (\mathbf{A}(\mathit{\theta }){)}^{-1}}{\partial \mathit{\theta }}& =-(\mathbf{A}(\mathit{\theta }){)}^{-1}\times \frac{\partial \mathbf{A}(\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes (\mathbf{A}(\mathit{\theta }){)}^{-1})\end{array}$$(31) ) to the right-hand side, to obtain (34) $$\begin{array}{ll}\frac{\partial }{\partial \mathit{\theta }}\frac{\partial {\mathbf{f}}_0(t;\mathit{\theta })}{\partial t}& =\frac{\partial {\mathbf{f}}_0(t;\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[0,0]}(\mathit{\theta }))+{\mathbf{f}}_0(t;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[0,0]}(\mathit{\theta })}{\partial \mathit{\theta }}\\ & +\frac{\partial {\mathbf{f}}_1(t;\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[1,0]}(\mathit{\theta }))+{\mathbf{f}}_1(t;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[1,0]}(\mathit{\theta })}{\partial \mathit{\theta }};\\ \frac{\partial }{\partial \mathit{\theta }}\frac{\partial {\mathbf{f}}_n(t;\mathit{\theta })}{\partial t}& =\frac{\partial {\mathbf{f}}_n(t;\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[n,n]}(\mathit{\theta }))+{\mathbf{f}}_n(t;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[n,n]}(\mathit{\theta })}{\partial \mathit{\theta }}\\ & +\frac{\partial {\mathbf{f}}_{n-1}(t;\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[n-1,n]}(\mathit{\theta }))\\ & +{\mathbf{f}}_{n-1}(t;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[n-1,n]}(\mathit{\theta })}{\partial \mathit{\theta }}\\ & +\frac{\partial {\mathbf{f}}_{n+1}(t;\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[n+1,n]}(\mathit{\theta }))\\ & +{\mathbf{f}}_{n+1}(t;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[n+1,n]}(\mathit{\theta })}{\partial \mathit{\theta }}\text{ for }n=1,\dots ,N-1;\\ \frac{\partial }{\partial \mathit{\theta }}\frac{\partial {\mathbf{f}}_N(t;\mathit{\theta })}{\partial t}& =\frac{\partial {\mathbf{f}}_N(t;\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[N,N]}(\mathit{\theta }))+{\mathbf{f}}_N(t;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[N,N]}(\mathit{\theta })}{\partial \mathit{\theta }}\\ & +\frac{\partial {\mathbf{f}}_{N-1}(t;\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[N-1,N]}(\mathit{\theta }))\\ & +{\mathbf{f}}_{N-1}(t;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[N-1,N]}(\mathit{\theta })}{\partial \mathit{\theta }}.\end{array}$$(34)

Further, we multiply both sides in the above by e^−st and then integrate, to evaluate the Laplace transforms: $$\begin{array}{ll}{\tilde{\mathbf{f}}}_n(s;\partial \mathit{\theta })& ={\int }_{t=0}^{\mathrm{\infty }}{e}^{-st}\frac{\partial {\mathbf{f}}_n(t;\mathit{\theta })}{\partial \mathit{\theta }}dt=\frac{\partial }{\partial \mathit{\theta }}{\int }_{t=0}^{\mathrm{\infty }}{e}^{-st}{\mathbf{f}}_n(t;\mathit{\theta })dt\\ {\tilde{\mathbf{f}}}_n(s;\mathit{\theta })& ={\int }_{t=0}^{\mathrm{\infty }}{e}^{-st}{\mathbf{f}}_n(t;\mathit{\theta })dt\end{array}$$

In the following, we shall apply the following identity for the Laplace transform ℒ(F(t)): $$\begin{array}{ll}\mathcal{L}(F^{\prime }(t))& =s\mathcal{L}(F(t))-F(0).\end{array}$$

Our key aim here is to evaluate $\frac{\partial {\mathbf{f}}_n(t;\mathit{\theta })}{\partial \mathit{\theta }}$ by inverting ${\tilde{\mathbf{f}}}_n(s;\partial \mathit{\theta })$. This can be achieved once we get a formula for ${\tilde{\mathbf{f}}}_n(s;\partial \mathit{\theta })$. Therefore, we focus on the formula for ${\tilde{\mathbf{f}}}_n(s;\partial \mathit{\theta })$ next. We have: $$\begin{array}{ll}{\tilde{\mathbf{f}}}_0(s;\partial \mathit{\theta })& =(-s\frac{\partial {\tilde{\mathbf{f}}}_0(s)}{\partial \mathit{\theta }}+{\tilde{\mathbf{f}}}_0(s;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[0,0]}(\mathit{\theta })}{\partial \mathit{\theta }}\\ & +{\tilde{\mathbf{f}}}_1(s;\partial \mathit{\theta })\times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[1,0]}(\mathit{\theta }))+{\tilde{\mathbf{f}}}_1(s;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[1,0]}(\mathit{\theta })}{\partial \mathit{\theta }})\\ & \times {({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[0,0]}(\mathit{\theta }))}^{-1}.\end{array}$$

Indeed, by (34(34) $$\begin{array}{ll}\frac{\partial }{\partial \mathit{\theta }}\frac{\partial {\mathbf{f}}_0(t;\mathit{\theta })}{\partial t}& =\frac{\partial {\mathbf{f}}_0(t;\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[0,0]}(\mathit{\theta }))+{\mathbf{f}}_0(t;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[0,0]}(\mathit{\theta })}{\partial \mathit{\theta }}\\ & +\frac{\partial {\mathbf{f}}_1(t;\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[1,0]}(\mathit{\theta }))+{\mathbf{f}}_1(t;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[1,0]}(\mathit{\theta })}{\partial \mathit{\theta }};\\ \frac{\partial }{\partial \mathit{\theta }}\frac{\partial {\mathbf{f}}_n(t;\mathit{\theta })}{\partial t}& =\frac{\partial {\mathbf{f}}_n(t;\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[n,n]}(\mathit{\theta }))+{\mathbf{f}}_n(t;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[n,n]}(\mathit{\theta })}{\partial \mathit{\theta }}\\ & +\frac{\partial {\mathbf{f}}_{n-1}(t;\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[n-1,n]}(\mathit{\theta }))\\ & +{\mathbf{f}}_{n-1}(t;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[n-1,n]}(\mathit{\theta })}{\partial \mathit{\theta }}\\ & +\frac{\partial {\mathbf{f}}_{n+1}(t;\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[n+1,n]}(\mathit{\theta }))\\ & +{\mathbf{f}}_{n+1}(t;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[n+1,n]}(\mathit{\theta })}{\partial \mathit{\theta }}\text{ for }n=1,\dots ,N-1;\\ \frac{\partial }{\partial \mathit{\theta }}\frac{\partial {\mathbf{f}}_N(t;\mathit{\theta })}{\partial t}& =\frac{\partial {\mathbf{f}}_N(t;\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[N,N]}(\mathit{\theta }))+{\mathbf{f}}_N(t;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[N,N]}(\mathit{\theta })}{\partial \mathit{\theta }}\\ & +\frac{\partial {\mathbf{f}}_{N-1}(t;\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[N-1,N]}(\mathit{\theta }))\\ & +{\mathbf{f}}_{N-1}(t;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[N-1,N]}(\mathit{\theta })}{\partial \mathit{\theta }}.\end{array}$$(34) ), $$\begin{array}{ll}{\int }_0^{\mathrm{\infty }}{e}^{-st}\frac{\partial }{\partial \mathit{\theta }}\frac{\partial {\mathbf{f}}_0(t;\mathit{\theta })}{\partial t}dt& =\frac{\partial }{\partial \mathit{\theta }}{\int }_0^{\mathrm{\infty }}{e}^{-st}\frac{\partial {\mathbf{f}}_0(t;\mathit{\theta })}{\partial t}dt\\ & =\frac{\partial }{\partial \mathit{\theta }}(s{\tilde{\mathbf{f}}}_0(s)-{\tilde{\mathbf{f}}}_0)\\ & =s\frac{\partial {\tilde{\mathbf{f}}}_0(s)}{\partial \mathit{\theta }}\\ & ={\int }_0^{\mathrm{\infty }}{e}^{-st}\frac{\partial {\mathbf{f}}_0(t;\mathit{\theta })}{\partial \mathit{\theta }}dt\times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[0,0]}(\mathit{\theta }))\\ & +{\int }_0^{\mathrm{\infty }}{e}^{-st}{\mathbf{f}}_0(t;\mathit{\theta })dt\times \frac{\partial {\mathbf{Q}}^{[0,0]}(\mathit{\theta })}{\partial \mathit{\theta }}\\ & +{\int }_0^{\mathrm{\infty }}{e}^{-st}\frac{\partial {\mathbf{f}}_1(t;\mathit{\theta })}{\partial \mathit{\theta }}dt\times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[1,0]}(\mathit{\theta }))\\ & +{\int }_0^{\mathrm{\infty }}{e}^{-st}{\mathbf{f}}_1(t;\mathit{\theta })dt\times \frac{\partial {\mathbf{Q}}^{[1,0]}(\mathit{\theta })}{\partial \mathit{\theta }},\end{array}$$ and so $$\begin{array}{ll}s\frac{\partial {\tilde{\mathbf{f}}}_0(s)}{\partial \mathit{\theta }}& ={\tilde{\mathbf{f}}}_0(s;\partial \mathit{\theta })\times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[0,0]}(\mathit{\theta }))+{\tilde{\mathbf{f}}}_0(s;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[0,0]}(\mathit{\theta })}{\partial \mathit{\theta }}\\ & +{\tilde{\mathbf{f}}}_1(s;\partial \mathit{\theta })\times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[1,0]}(\mathit{\theta }))+{\tilde{\mathbf{f}}}_1(s;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[1,0]}(\mathit{\theta })}{\partial \mathit{\theta }}.\end{array}$$

This gives $$\begin{array}{ll}{\tilde{\mathbf{f}}}_0(s;\partial \mathit{\theta })& =-s\frac{\partial {\tilde{\mathbf{f}}}_0(s)}{\partial \mathit{\theta }}+{\tilde{\mathbf{f}}}_0(s;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[0,0]}(\mathit{\theta })}{\partial \mathit{\theta }}\\ & +{\tilde{\mathbf{f}}}_1(s;\partial \mathit{\theta })\times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[1,0]}(\mathit{\theta }))+{\tilde{\mathbf{f}}}_1(s;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[1,0]}(\mathit{\theta })}{\partial \mathit{\theta }},\end{array}$$ and thus $$\begin{array}{ll}{\tilde{\mathbf{f}}}_0(s;\partial \mathit{\theta })& =(-s\frac{\partial {\tilde{\mathbf{f}}}_0(s)}{\partial \mathit{\theta }}+{\tilde{\mathbf{f}}}_0(s;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[0,0]}(\mathit{\theta })}{\partial \mathit{\theta }}\\ & +{\tilde{\mathbf{f}}}_1(s;\partial \mathit{\theta })\times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[1,0]}(\mathit{\theta }))+{\tilde{\mathbf{f}}}_1(s;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[1,0]}(\mathit{\theta })}{\partial \mathit{\theta }})\\ & \times {({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[0,0]}(\mathit{\theta }))}^{-1}\end{array}$$ as required.

By a similar argument, $$\begin{array}{ll}{\tilde{\mathbf{f}}}_n(s;\partial \mathit{\theta })& =(-s\frac{\partial {\tilde{\mathbf{f}}}_n(s)}{\partial \mathit{\theta }}+{\tilde{\mathbf{f}}}_n(s;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[n,n]}(\mathit{\theta })}{\partial \mathit{\theta }}\\ & +{\tilde{\mathbf{f}}}_{n-1}(s;\partial \mathit{\theta })\times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[n-1,n]}(\mathit{\theta }))+{\tilde{\mathbf{f}}}_{n-1}(s;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[n-1,n]}(\mathit{\theta })}{\partial \mathit{\theta }}\\ & +{\tilde{\mathbf{f}}}_{n+1}(s;\partial \mathit{\theta })\times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[n+1,n]}(\mathit{\theta }))+{\tilde{\mathbf{f}}}_{n+1}(s;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[n+1,n]}(\mathit{\theta })}{\partial \mathit{\theta }})\\ & \times {(\mathbf{I}-{\mathbf{I}}_k\otimes {\mathbf{Q}}^{[n,n]}(\mathit{\theta }))}^{-1}\end{array}$$ for $n=1,\dots ,N-1$, and ${\tilde{\mathbf{f}}}_N(s;\partial \mathit{\theta })$ is equal to $$\begin{array}{ll}& (-s\frac{\partial {\tilde{\mathbf{f}}}_N(s)}{\partial \mathit{\theta }}+{\tilde{\mathbf{f}}}_N(s;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[N,N]}(\mathit{\theta })}{\partial \mathit{\theta }}\\ & +{\tilde{\mathbf{f}}}_{N-1}(s;\partial \mathit{\theta })\times ({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[N-1,N]}(\mathit{\theta }))+{\tilde{\mathbf{f}}}_{N-1}(s;\mathit{\theta })\times \frac{\partial {\mathbf{Q}}^{[N-1,N]}(\mathit{\theta })}{\partial \mathit{\theta }})\\ & \times {({\mathbf{I}}_k\otimes {\mathbf{Q}}^{[N,N]}(\mathit{\theta }))}^{-1}.\end{array}$$

Thus, we have derived a recursion.

Inverting the above Laplace transforms will give us the required $$\begin{array}{l}\frac{\partial {\mathbf{f}}_n(t,\mathit{\theta })}{\partial \mathit{\theta }}\text{for }n=0,1,\dots ,N.\end{array}$$

4.5.2. A special form of the DEs

Let us consider a simpler version of the DEs: $$\begin{array}{ll}\frac{\partial \mathbf{f}(t;\mathit{\theta })}{\partial t}& =\mathbf{f}(t;\mathit{\theta })\mathbf{Q}(\mathit{\theta }).\end{array}$$

In such a case, we obtain: $$\begin{array}{ll}\frac{\partial }{\partial \mathit{\theta }}\frac{\partial \mathbf{f}(t;(\mathit{\theta }))}{\partial t}& =\frac{\partial }{\partial \mathit{\theta }}\mathbf{f}(t;\mathit{\theta })\mathbf{Q}(\mathit{\theta })\\ & =\frac{\partial \mathbf{f}(t;\mathit{\theta })}{\partial \mathit{\theta }}\times ({\mathbf{I}}_k\otimes \mathbf{Q}(\mathit{\theta }))+\mathbf{f}(t;\mathit{\theta })\times \frac{\partial \mathbf{Q}(\mathit{\theta })}{\partial \mathit{\theta }}\end{array}$$ and so the Laplace transform is expressed by $$\begin{array}{ll}{\int }_0^{\mathrm{\infty }}{e}^{-st}\frac{\partial }{\partial \mathit{\theta }}\frac{\partial \mathbf{f}(t;\mathit{\theta })}{\partial t}dt& =\frac{\partial }{\partial \mathit{\theta }}{\int }_0^{\mathrm{\infty }}{e}^{-st}\frac{\partial \mathbf{f}(t;\mathit{\theta })}{\partial t}dt\\ & =\frac{\partial }{\partial \mathit{\theta }}(s\times {\int }_0^{\mathrm{\infty }}{e}^{-st}\mathbf{f}(t;\mathit{\theta })dt-\mathbf{f}(0;\mathit{\theta }))\\ & =s\times \frac{\partial }{\partial \mathit{\theta }}({\int }_0^{\mathrm{\infty }}{e}^{-st}\mathbf{f}(t;\mathit{\theta })dt)\\ & ={\int }_0^{\mathrm{\infty }}{e}^{-st}\frac{\partial \mathbf{f}(t;\mathit{\theta })}{\partial \mathit{\theta }}dt\times ({\mathbf{I}}_k\otimes \mathbf{Q}(\mathit{\theta }))+{\int }_0^{\mathrm{\infty }}{e}^{-st}\mathbf{f}(t;\mathit{\theta })dt\\ & \times \frac{\partial \mathbf{Q}(\mathit{\theta })}{\partial \mathit{\theta }}.\end{array}$$

Therefore $$\begin{array}{ll}& {\int }_0^{\mathrm{\infty }}{e}^{-st}\frac{\partial \mathbf{f}(t;\mathit{\theta })}{\partial \mathit{\theta }}dt\\ & =(-s\times \frac{\partial }{\partial \mathit{\theta }}({\int }_0^{\mathrm{\infty }}{e}^{-st}\mathbf{f}(t;\mathit{\theta })dt)+{\int }_0^{\mathrm{\infty }}{e}^{-st}\mathbf{f}(t;\mathit{\theta })dt\times \frac{\partial \mathbf{Q}(\mathit{\theta })}{\partial \mathit{\theta }})\\ & \times {({\mathbf{I}}_k\otimes \mathbf{Q}(\mathit{\theta }))}^{-1}.\end{array}$$

5. Numerical example

To illustrate the application of the theory, consider the LD-QBD $\{(X(t),\phi (t)):t\ge 0\}$ from Example 2 modeling a system with capacity N = 100 and two types of customers arriving at rates λ₁ = 5 and λ₂ = 8, and departing at rates ${\mu }_1=1/10$ and ${\mu }_2=1/5$, respectively. We evaluate the stationary distribution 𝝅 of this process using expressions (1(1) $$\begin{array}{ll}{\mathit{\pi }}_n& ={\mathit{\pi }}_{n+1}{\hat{\mathbf{R}}}^{(n)}={\mathit{\pi }}_N\prod _{k=N-1}^n{\hat{\mathbf{R}}}^{(k)}.\end{array}$$(1) )–(2(2) $$\begin{array}{l}\phantom{.}\{\begin{array}{l}{\mathit{\pi }}_N({\hat{\mathbf{R}}}^{(N-1)}{\mathbf{Q}}^{[N-1,N]}+{\mathbf{Q}}^{[N,N]})=\mathbf{0}\\ {\mathit{\pi }}_N(\mathbf{1}+\sum _{n=0}^{N-1}\prod _{k=N-1}^n{\hat{\mathbf{R}}}^{(k)}\mathbf{1})=1.\end{array}\end{array}$$(2) ) and present the output in Figure 2.

Figure 2. Stationary distribution of the LD-QBD $\{(X(t),\phi (t)):t\ge 0\}$ with N = 100, λ₁ = 5, λ₂ = 8, ${\mu }_1=1/10$ and ${\mu }_2=1/5$in Example 2.

We note from Figure 2 that the stationary distribution is supported on the set n≥60 and $i\in [30,70]$ with median for n is high 80 and median for i in high 40.

Next, we evaluate the derivative $\frac{\partial }{\partial \mathit{\theta }}\mathit{\pi }$ of the stationary distribution of this process using expressions in Section 4.2 and present the output in Figure 3. We note that sensitivities of π_{n, i} with respect to λ₁ and μ₂ are non-negative and with respect to λ₂ and μ₁ are non-positive, which is consistent with intuition. More precisely, recall that φ(t) denotes number of patients of type 1 in the system. Therefore, if λ₁ increases the number of patients of type 1 increases and the proportion of type 1 patients π_{n, i} must increase. Similarly, if the service rate of type 2 patients increases, the number of patients of type 2 decreases and the proportion of type 1 patients π_{n, i} must increase.

Figure 3. The derivatives of the stationary distribution of the LD-QBD $\{(X(t),\phi (t)):t\ge 0\}$ with N = 100, λ₁ = 5, λ₂ = 8, ${\mu }_1=1/10$ and ${\mu }_2=1/5$ in Example 2. Color of the line decodes specific value of n: blue 65, orange 74, yellow 85 and purple 95.

Conversely, if λ₂ increases the number of patients of type 2 increases and the proportion of type 1 patients π_{n, i} must decrease. Similarly, if service rate of type 1 patients increases,the number of patients of type 1 decreases and the proportion of type 1 patients π_{n, i} must decrease too.

We also remark that the amplitude of change in π_{n, i} depends both on the level of n and the level of i. Absolute value of the respective derivatives for fixed n are uni-modal in i (as we read from the Figure 3). Looking at the color code in Figure 3, we also note that the impact of all parameters on the stationary distribution is stronger around median values of 𝝅.

Further, we focus on first hitting times to lower levels, and using results from Section 3.2 and Section 3.4, we evaluate $$\begin{array}{l}\sum _j[{\mathbf{g}}^{n,n-k}(t){]}_{ij},\sum _j[\frac{\partial }{\partial \mathit{\theta }}{\mathbf{g}}^{n,n-k}(t){]}_{ij},\sum _j[{\mathbf{g}}_{\mathcal{A}}^{n,n-k}(t){]}_{ij}\text{and}\sum _j[\frac{\partial }{\partial \mathit{\theta }}{\mathbf{g}}_{\mathcal{A}}^{n,n-k}(t){]}_{ij}\end{array}$$ given the initial state $(n,i)=(100,100)$ and the destination level $n-k=90$, with $\mathcal{A}=[95,100]$. The output is presented in Figures 4–6.

Figure 4. Scalar quantities $\sum _j[{\mathbf{g}}^{n,n-k}(t){]}_{ij}$ and $\sum _j[{\mathbf{g}}_{\mathcal{A}}^{n,n-k}(t){]}_{ij}$, given $(n,i)=(100,100)$, $n-k=90$ and $\mathcal{A}=[95,100]$, of the LD-QBD $\{(X(t),\phi (t)):t\ge 0\}$ with N = 100, λ₁ = 5, λ₂ = 8, ${\mu }_1=1/10$ and ${\mu }_2=1/5$ in Example 2.

Figure 5. Scalar quantities $\sum _j[\frac{\partial }{\partial \mathit{\theta }}{\mathbf{g}}^{n,n-k}(t){]}_{ij}$, given $(n,i)=(100,100)$ and $n-k=90$ of the LD-QBD $\{(X(t),\phi (t)):t\ge 0\}$ with N = 100, λ₁ = 5, λ₂ = 8, ${\mu }_1=1/10$ and ${\mu }_2=1/5$ in Example 2.

Figure 6. Scalar quantities $\sum _j[\frac{\partial }{\partial \mathit{\theta }}{\mathbf{g}}_{\mathcal{A}}^{n,n-k}(t){]}_{ij}$, given $(n,i)=(100,100)$, $n-k=90$ and $\mathcal{A}=[95,100]$ of the LD-QBD $\{(X(t),\phi (t)):t\ge 0\}$ with N = 100, λ₁ = 5, λ₂ = 8, ${\mu }_1=1/10$ and ${\mu }_2=1/5$ in Example 2.

In Figure 4, we present the density of the distribution of time before first arriving at level $n-k=90$ and any phase given initial state (100, 100) (full system, only patients of type 1) on the left graph; and the density of distribution of time spent at levels [95, 100] before first arriving at level $n-k=90$ and any phase given initial state (100, 100) on the right graph. Both graphs have the same shape representing that the time spent in the intervals [91, 100] and [95, 100] increases until some threshold time (at t = 5 and just before 5, respectively), from which it decreases and becomes extremely unlikely (values 30 and 25, respectively).

We note that to obtain the derivatives $$\begin{array}{l}\frac{\partial }{\partial \mathit{\theta }}{\mathbf{g}}^{n,n-k}(t)\text{and }\frac{\partial }{\partial \mathit{\theta }}{\mathbf{g}}_{\mathcal{A}}^{n,n-k}(t),\end{array}$$ we have applied the algorithm in Den Iseger[3], its convenient translation to a code in Toutain et al.[14], and an observation analogous to (33(33) $$\begin{array}{ll}{\int }_0^{\mathrm{\infty }}{e}^{-st}\frac{\partial }{\partial \mathit{\theta }}{\mathbf{f}}_n(t;\mathit{\theta })dt& =\frac{\partial }{\partial \mathit{\theta }}{\int }_0^{\mathrm{\infty }}{e}^{-st}{\mathbf{f}}_n(t;\mathit{\theta })dt=\frac{\partial }{\partial \mathit{\theta }}{\tilde{\mathbf{f}}}_n(s;\mathit{\theta }).\end{array}$$(33) ), that is, $$\begin{array}{ll}{\int }_0^{\mathrm{\infty }}{e}^{-st}\frac{\partial }{\partial \mathit{\theta }}{\mathbf{g}}^{n,n-k}(t;\mathit{\theta })dt& =\frac{\partial }{\partial \mathit{\theta }}{\mathbf{G}}^{n,n-k}(s;\mathit{\theta })\text{and}\\ {\int }_0^{\mathrm{\infty }}{e}^{-st}\frac{\partial }{\partial \mathit{\theta }}{\mathbf{g}}_{\mathcal{A}}^{n,n-k}(t;\mathit{\theta })dt& =\frac{\partial }{\partial \mathit{\theta }}{\mathbf{G}}_{\mathcal{A}}^{n,n-k}(s;\mathit{\theta }).\end{array}$$

The recursions for $$\begin{array}{l}\frac{\partial }{\partial \mathit{\theta }}{\mathbf{G}}^{n,n-k}(s;\mathit{\theta })\text{and}\frac{\partial }{\partial \mathit{\theta }}{\mathbf{G}}_{\mathcal{A}}^{n,n-k}(s;\mathit{\theta })\end{array}$$ are analogous to the recursion for $\frac{\partial }{\partial \mathit{\theta }}{\mathbf{G}}^{n,n-k}(\mathit{\theta })$ in Section 4.3, with the suitable adaptation in the formula being the use of $({\mathbf{Q}}^{[n,n]}(\mathit{\theta })-s\mathbf{I})$ and $({\mathbf{Q}}^{[n,n]}(\mathit{\theta })-s\mathbf{I}\times I(n\in \mathcal{A}))$, respectively, in place of 𝐐^{[n, n]}(𝜽), for all n.

We note that derivatives presented in Figure 5 (A) and (B), and Figure 6 (A) and (B) exhibit the same behavior. When λ₁ or λ₂ increase, there are more patients in the system, and it takes more time to first reduce the number of patients from level 100 to level 90. Hence, both graphs in Figure 4 move to the right when λ_i for $i=1,2$ increases. That means that for small times, density decreases, for moderate times density increases and for large times, density does not change from 0.

Of course, we have observed opposite behavior in Figure 5 (C) and (D), and Figure 6 (C) and (D). It is due to the fact that the increase of μ₁ or μ₂ results in fewer patients in the system, and it takes less time to first reduce the number of patients from level 100 to level 90. Hence, both graphs in Figure 4 move to the left when μ_i for $i=1,2$ increases, meaning that for small times density increases, for moderate times density decreases, and for large times density does not change from 0.

6. Conclusions

In this study, we developed novel results for the sensitivity analysis of level-dependent quasi-birth-and-death processes (LD-QBDs), which have potential applications in many complex systems, including healthcare systems, as discussed throughout the article.

We derive recursive expressions for several key performance metrics of LD-QBDs. We demonstrate the application of the theory through numerical examples where we highlight relevant intuitions and the physical interpretations of these metrics.

This analysis forms the foundation for future applications in healthcare systems and other important processes within queueing theory. The ideas developed here can be applied to continuous and discrete time frameworks.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

Anna Aksamit is supported by the Australian Research Council Early Career Researcher Award DE200100896. Małgorzata O’Reilly is supported by the Australian Research Council Discovery Project DP180100352. Zbigniew Palmowski was partially supported by the National Science Centre (Poland) under the grant 2021/41/B/HS4/00599. Małgorzata O’Reilly is grateful for the support and hospitality of the Sydney Mathematical Research Institute (SMRI) in June 2023.