Linear Stochastic Fluid Networks: Rare-Event Simulation and Markov Modulation

TL;DR

This paper develops efficient importance sampling algorithms for rare-event probabilities in Markov-modulated linear stochastic fluid networks, providing theoretical guarantees and analyzing their performance in different regimes.

Contribution

It introduces importance sampling methods with provable efficiency for Markov-modulated fluid networks and extends techniques from light-tailed rare-event simulation.

Findings

Algorithm is asymptotically efficient in slow modulation regime

Number of runs grows polynomially for unmodulated networks

Provides recursion for moments of storage levels

Abstract

We consider a linear stochastic fluid network under Markov modulation, with a focus on the probability that the joint storage level attains a value in a rare set at a given point in time. The main objective is to develop efficient importance sampling algorithms with provable performance guarantees. For linear stochastic fluid networks without modulation, we prove that the number of runs needed (so as to obtain an estimate with a given precision) increases polynomially (whereas the probability under consideration decays essentially exponentially); for networks operating in the slow modulation regime, our algorithm is asymptotically efficient. Our techniques are in the tradition of the rare-event simulation procedures that were developed for the sample-mean of i.i.d. one-dimensional light-tailed random variables, and intensively use the idea of exponential twisting. In passing, we also point out how to set up a recursion to evaluate the (transient and stationary) moments of the joint storage level in Markov-modulated linear stochastic fluid networks.

Figures

Equations

$$X(t)=\sum _{j=1}^{N(t)}{B}_j{\mathrm{e}}^{-r(t-U_j)}\stackrel{\mathrm{d}}{=}\sum _{j=1}^{N(t)}{B}_j{\mathrm{e}}^{-r{U}_j},$$

$$M(\vartheta ):=\mathbb{E}{\mathrm{e}}^{\vartheta X(t)}$$

$$\mathbb{E}X(t)=\frac{\lambda }{r}(1-{\mathrm{e}}^{-rt})\mathbb{E}B=:m(t).$$

$$Y_n(t):=\sum _{i=1}^n{X}_i(t),$$

$$p_n(a)=\mathbb{P}(Y_n(t)⩾na)$$

$$I(a):=\sup_{\vartheta}\big{(}\vartheta a-\log M(\vartheta)\big{)},\>\>\>\>\vartheta^{\star}:=\arg\sup_{\vartheta}\big{(}\vartheta a-\log M(\vartheta)\big{)},$$

$$\underset{n\to \infty }{\lim }\frac{1}{n}\log p_n(a)=-I(a)=-{\vartheta }^{\star }a+\log M({\vartheta }^{\star }).$$

$${\mathbb{E}}_{\mathbb{Q}}{\mathrm{e}}^{\vartheta X(t)}=\frac{\mathbb{E}{\mathrm{e}}^{(\vartheta +{\vartheta }^{\star })X(t)}}{\mathbb{E}{\mathrm{e}}^{{\vartheta }^{\star }X(t)}}=\frac{M(\vartheta +{\vartheta }^{\star })}{M({\vartheta }^{\star })};$$

$${\mathbb{E}}_{\mathbb{Q}}X(t)=\frac{M^{\prime }({\vartheta }^{\star })}{M({\vartheta }^{\star })}=a.$$

$$M(\vartheta +{\vartheta }^{\star })=\sum _{k=0}^{\infty }{\mathrm{e}}^{-\lambda t}\frac{(\lambda t\mathbb{E}\exp ({\vartheta }^{\star }B{\mathrm{e}}^{-rU}){)}^k}{k!}{\left(\frac{\mathbb{E}\exp ((\vartheta +{\vartheta }^{\star })B{\mathrm{e}}^{-rU})}{\mathbb{E}\exp ({\vartheta }^{\star }B{\mathrm{e}}^{-rU})}\right)}^k,$$

$$\sum _{k=0}^{\infty }\exp (-\lambda t\mathbb{E}\exp ({\vartheta }^{\star }B{\mathrm{e}}^{-rU}))\frac{(\lambda t\mathbb{E}\exp ({\vartheta }^{\star }B{\mathrm{e}}^{-rU}){)}^k}{k!}{\left(\frac{\mathbb{E}\exp ((\vartheta +{\vartheta }^{\star })B{\mathrm{e}}^{-rU})}{\mathbb{E}\exp ({\vartheta }^{\star }B{\mathrm{e}}^{-rU})}\right)}^k.$$

$$\lambda t\mathbb{E}\exp ({\vartheta }^{\star }B{\mathrm{e}}^{-rU})=\lambda {\int }_0^t\beta ({\mathrm{e}}^{-ru}{\vartheta }^{\star })\mathrm{d}u,$$

$${\mathbb{E}}_{\mathbb{Q}}\exp ((\vartheta +{\vartheta }^{\star })B{\mathrm{e}}^{-rU})=\frac{\mathbb{E}\exp ((\vartheta +{\vartheta }^{\star })B{\mathrm{e}}^{-rU})}{\mathbb{E}\exp ({\vartheta }^{\star }B{\mathrm{e}}^{-rU})}.$$

$$\mathbb{E}\exp ((\vartheta +{\vartheta }^{\star })B{\mathrm{e}}^{-rU})={\int }_0^t\frac{\beta ({\mathrm{e}}^{-ru}(\vartheta +{\vartheta }^{\star }))}{\beta ({\mathrm{e}}^{-ru}{\vartheta }^{\star })}\frac{1}{t}\beta ({\mathrm{e}}^{-ru}{\vartheta }^{\star })\mathrm{d}u,$$

$${\mathbb{E}}_{\mathbb{Q}}\exp ((\vartheta +{\vartheta }^{\star })B{\mathrm{e}}^{-rU})={\int }_0^t\frac{\beta ({\mathrm{e}}^{-ru}(\vartheta +{\vartheta }^{\star }))}{\beta ({\mathrm{e}}^{-ru}{\vartheta }^{\star })}\frac{\beta ({\mathrm{e}}^{-ru}{\vartheta }^{\star })}{{\int }_0^t\beta ({\mathrm{e}}^{-rv}{\vartheta }^{\star })\mathrm{d}v}\mathrm{d}u.$$

$$f_U^{\mathbb{Q}}(u)=\frac{\beta ({\mathrm{e}}^{-ru}{\vartheta }^{\star })}{{\int }_0^t\beta ({\mathrm{e}}^{-rv}{\vartheta }^{\star })\mathrm{d}v}.$$

$${\int }_0^s\beta ({\mathrm{e}}^{-ru}\vartheta )\mathrm{d}u$$

$$m(t)=a\left(1-\frac{\vartheta }{\mu }\right)\left(1-\frac{\vartheta }{\mu }{\mathrm{e}}^{-rt}\right),$$

$${\vartheta }^{\star }=\frac{\mu {\mathrm{e}}^{rt}}{2}\left((1+{\mathrm{e}}^{-rt})-\sqrt{(1-{\mathrm{e}}^{-rt}{)}^2+4{\mathrm{e}}^{-rt}\frac{m(t)}{a}}\right)$$

$$\frac{\lambda }{r}\log \left(\frac{\mu {\mathrm{e}}^{rt}-{\vartheta }^{\star }}{\mu -{\vartheta }^{\star }}\right)$$

$$F_U^{\mathbb{Q}}(u):=\mathbb{Q}(U⩽u)=\log \left(\frac{\mu {\mathrm{e}}^{ru}-{\vartheta }^{\star }}{\mu -{\vartheta }^{\star }}\right)/\log \left(\frac{\mu {\mathrm{e}}^{rt}-{\vartheta }^{\star }}{\mu -{\vartheta }^{\star }}\right)$$

$$\frac{1}{r}\log \left({\left({\mathrm{e}}^{rt}-\frac{{\vartheta }^{\star }}{\mu }\right)}^H{\left(1-\frac{{\vartheta }^{\star }}{\mu }\right)}^{1-H}+\frac{{\vartheta }^{\star }}{\mu }\right).$$

$$L=\frac{\mathrm{d}\mathbb{P}}{\mathrm{d}\mathbb{Q}}={\mathrm{e}}^{-{\vartheta }^{\star }{Y}_n(t)}{\mathrm{e}}^{n\log M({\vartheta }^{\star })}.$$

$$\frac{1}{p_n(a)}\cdot \frac{T}{\sqrt{N}}\sqrt{\mathbb{V}{\mathrm{ar}}_{\mathbb{Q}}(L^{2}I)},$$

$${\Sigma }_n:=\frac{T^2}{{\epsilon }^2}\frac{\mathbb{V}{\mathrm{ar}}_{\mathbb{Q}}(L^{2}I)}{(p_n(a){)}^2}.$$

$$p_n(a)={\mathbb{E}}_{\mathbb{Q}}(LI)\sim \frac{1}{\sqrt{n}}\frac{1}{{\vartheta }^{\star }\sqrt{2\pi \tau }}{\mathrm{e}}^{-nI(a)},\tau :={\frac{{\mathrm{d}}^2}{\mathrm{d}{\vartheta }^2}\log M(\vartheta )|}_{\vartheta ={\vartheta }^{\star }}.$$

$${\mathbb{E}}_{\mathbb{Q}}(L^{2}I)\sim \frac{1}{\sqrt{n}}\frac{1}{2{\vartheta }^{\star }\sqrt{2\pi \tau }}{\mathrm{e}}^{-2nI(a)};$$

$${\Sigma }_n\sim \alpha \sqrt{n},\alpha =\frac{T^2}{{\epsilon }^2}{\vartheta }^{\star }\cdot \frac{1}{2}\sqrt{2\pi \tau }.$$

$$\tau =\frac{\lambda }{r}\left(\frac{1}{(\mu -{\vartheta }^{\star }{)}^2}-\frac{1}{(\mu {\mathrm{e}}^{rt}-{\vartheta }^{\star }{)}^2}\right)=\mathrm{1.8240.}$$

$$X^{(1)}(t)=\sum _{j=1}^N{B}_j{\mathrm{e}}^{-r_1(t-U_j)}$$

Full text

Affiliations. O. Boxma is with Eurandom and Department of Mathematics and Computer Science, Eindhoven University of Technology, PO Box 513, 5600 MB Eindhoven, the Netherlands.

E. Cahen is with CWI, Science Park 123, 1098 XG Amsterdam, the Netherlands.

D. Koops and M. Mandjes are with Korteweg-de Vries Institute for Mathematics, University of Amsterdam, Science Park 105-107, 1098 XG Amsterdam, the Netherlands.

M. Mandjes is corresponding author. Email: m.r.h.mandjes@uva.nl. Telephone: + 31 20 525 5164. Fax: + 31 20 525 7820.

Acknowledgments. The research of O. Boxma, D. Koops and M. Mandjes was partly funded by the NWO Gravitation Project Networks, Grant Number 024.002.003. The research of O. Boxma was also partly funded by the Belgian Government, via the IAP Bestcom project. The research of E. Cahen was funded by an NWO grant, Grant Number 613.001.352.

1. Introduction

Linear stochastic fluid networks, as introduced in [19], can be informally described as follows. Consider a network consisting of $L$ stations. Jobs, whose sizes are i.i.d. samples from some general $L$-dimensional distribution, arrive at the stations according to a Poisson process. At each of the nodes, in between arrivals the storage level decreases exponentially. Processed traffic is either transferred to the other nodes or leaves the network (according to a given routing matrix). In addition to this basic version of the linear stochastic fluid network, there is also its Markov modulated counterpart [18], in which the arrival rate, the distribution of the job sizes, and the routing matrix depend on the state of an external, autonomously evolving finite-state continuous-time Markov chain (usually referred to as the background process).

Linear stochastic fluid networks can be seen as natural fluid counterparts of corresponding infinite-server queues. As such, they inherit several nice properties of those infinite-server queues. In particular, separate infinitesimally small fluid particles, moving through the network, do not interfere, and are therefore mutually independent. Essentially due to this property, linear stochastic fluid networks allow explicit analysis; in particular, the joint Laplace transform of the storage levels at a given point in time can be expressed in closed form as a function of the arrival rate, the Laplace transform of the job sizes and the routing matrix [19, Thm. 5.1].

When Markov modulation is imposed, the analysis becomes substantially harder. Conditional on the path of the background process, again explicit expressions can be derived, cf. [18, Thm. 1]. Unconditioning, however, cannot be done in a straightforward manner. As a consequence the results found are substantially less explicit than for the non-modulated linear stochastic fluid network. In [18] also a system of ordinary differential equations has been set up that provides the transform of the stationary storage level; in addition, conditions are identified that guarantee the existence of such a stationary distribution.

In this paper we focus on rare events for Markov-modulated linear stochastic fluid networks. More specifically, in a particular scaling regime (parameterized by $n$) we analyze the probability $p_{n}$ that at a given point in time the network storage vector is in a given rare set. By scaling the arrival rate as well as the rare set (which amounts to multiplying them by a scaling parameter $n$), the event of interest becomes increasingly rare. More specifically, under a Cramér-type assumption on the job-size distribution, application of large-deviations theory yields that $p_{n}$ decays (roughly) exponentially. As $p_{n}$ can be characterized only asymptotically, one could consider the option of using simulation to obtain precise estimates. The effectiveness, however, of such an approach is limited due to the rarity of the event under consideration: in order to get a reliable estimate, one needs sufficiently many runs in which the event occurs. This is the reason why one often resorts to simulation using importance sampling (or: change of measure). This is a variance reduction technique in which one replaces the actual probability measure by an alternative measure under which the event under consideration is not rare; correcting the simulation output with appropriate likelihood ratios yields an unbiased estimate.

The crucial issue when setting up an importance sampling procedure concerns the choice of the alternative measure: one would like to select one that provides a substantial variance reduction, or is even (in some sense) optimal. The objective of this paper is to develop a change of measure which performs provably optimally.

Our ultimate goal is to obtain an efficient simulation procedure for Markov-modulated linear stochastic fluid networks. We do so by (i) first considering a single node without modulation, (ii) then multi-node systems, still without modulation, and (iii) finally modulated multi-node systems. There are two reasons for this step-by-step setup:

• $\circ$

For the non-modulated models we have more refined results than for the modulated models. More specifically, for the non-modulated models we have developed estimates for the number of runs $\Sigma_{n}$ required to obtain an estimate with predefined precision (showing that $\Sigma_{n}$ grows polynomially in the rarity parameter $n$), whereas for modulated models we can just prove that $\Sigma_{n}$ grows subexponentially.

• $\circ$

In addition, this approach allows the reader to get gradually familiar with the concepts used in this paper.

The construction and analysis of our importance sampling methodology is based on the ideas developed in [8]; there the focus was on addressing similar issues for a single-node Markov modulated infinite-server system. In line with [8], we consider the regime in which the background process is ‘slow’: while we (linearly) speed up the driving Poisson process, we leave the rates of the Markovian background process unalterned.

A traditional, thoroughly examined, importance sampling problem concerns the sample mean $S_{n}$ of $n$ i.i.d. light-tailed random variables $X_{1},\ldots,X_{n}$; the objective there is to estimate ${\mathbb{P}}(S_{n}\geqslant a)$ for $a>{\mathbb{E}}X_{1}$ and $n$ large. As described in [1, Section VI.2], in this situation importance sampling (i.e., sampling under an alternative measure, and translating the simulation output back by applying appropriate likelihood ratios) works extremely well. To this end, the distribution of the $X_{i}$ s should be exponentially twisted. As it turns out, in our setup, the probability of our interest can be cast in terms of this problem. Compared to the standard setup of sample means of one-dimensional random variables, however, there are a few complications: (i) in our case it is not a priori clear how to sample from the exponentially twisted distributions, (ii) we consider multi-dimensional distributions (i.e., rare-event probabilities that concern the storage levels of all individual buffers in the network), (iii) we impose Markov modulation. We refer to e.g. [14, 20] for earlier work on similar problems.

In passing, we also point out how to set up a recursion to evaluate the (transient and stationary) moments of the joint storage level in Markov-modulated linear stochastic fluid networks (where the results in [18] are restricted to just the first two stationary moments at epochs that the background process jumps).

The single-node model without modulation falls in the class of (one-dimensional) shot-noise models, for which efficient rare-event simulation techniques have been developed over the past, say, two decades. Asmussen and Nielsen [3] and Ganesh et al.  [13] consider the probability that a shot-noise process decreased by a linear drift ever exceeds some given level. Relying on sample-path large deviations results, an asymptotically efficient importance sampling algorithm is developed, under the same scaling as the one we consider in our paper. The major difference with our model (apart from the fact that we deal with considerably more general models, as we focus on networks and allow modulation) is that we focus on a rare-event probability that relates to the position of the process at a fixed point in time; in this setting we succeed in finding accurate estimates of the number of runs needed to get an estimate of given precision.

There is a vast body of literature related to the broader area of rare-event simulation for queueing systems. We refer to the literature overviews [6, 17]; interesting recent papers include [2, 10, 23].

This paper is organized as follows. In Section 2 the focus is on a single-node network, without Markov modulation (addressing complication (i) above), Section 3 addresses the extension to multi-node systems (addressing complication (ii)), and in Section 4 the feature of modulation is added (addressing complication (iii)). In each of these three sections, we propose a change of measure, quantify its performance, and demonstrate its efficiency through simulation experiments. In Section 4.1 we include the explicit expressions for the moments in Markov-modulated linear stochastic fluid networks. A discussion and concluding remarks are found in Section 5.

2. Single resource, no modulation

To introduce the concepts we work with in this paper, we analyze in this section a linear stochastic fluid network consisting of a single node, in which the input is just compound Poisson (so no Markov modulation is imposed). More precisely, in the model considered, jobs arrive according to a Poisson process with rate $\lambda$, bring along i.i.d. amounts of work (represented by the sequence of i.i.d. random variables $(B_{1},B_{2},\ldots)$), and the workload level decays exponentially at a rate $r>0.$ This model belongs to the class of shot-noise processes. As mentioned in the introduction, we gradually extend the model in the next sections.

2.1. Preliminaries

We first present a compact representation for the amount of work in the system at time $t$, which we denote by $X(t)$, through its moment generating function. To this end, let $N(t)$ denote a Poisson random variable with mean $\lambda t$, and $(U_{1},U_{2},\ldots)$ i.i.d. uniformly distributed random variables (on the interval $[0,t]$). Assume in addition that the random objects $(B_{1},B_{2},\ldots)$, $N(t)$, and $(U_{1},U_{2},\ldots)$ are independent. Then it is well-known that the value of our shot-noise process at time $t$ can be expressed as

[TABLE]

where the distributional equality is a consequence of the fact that the distribution of $U$ is symmetric on the interval $[0,t]$. It is easy to compute the moment generating function (mgf) of $X(t)$, by conditioning on the value of $N(t)$:

[TABLE]

where $\beta(\cdot)$ is the mgf corresponding to $B$ (throughout assumed to exist). By differentiating and inserting $\vartheta=0$, it follows immediately that

[TABLE]

Higher moments can be found by repeated differentiation. We note that, as $t$ is held fixed throughout the document, we often write $N$ rather than $N(t).$

2.2. Tail probabilities, change of measure

The next objective is to consider the asymptotics of the random variable $X(t)$ under a particular scaling. In this scaling we let the arrival rate be $n\lambda$ rather than just $\lambda$, for ${n\in{\mathbb{N}}}$. The value of the shot-noise process is now given by

[TABLE]

with the vector $(X_{1}(t),\ldots,X_{n}(t))$ consisting of i.i.d. copies of the random variable $X(t)$ introduced above; here the infinite divisibility of a Compound Poisson distribution is used.

Our goal is to devise techniques to analyze the tail distribution of $Y_{n}(t).$ Standard theory now provides us with the asymptotics of

[TABLE]

for some $a>m(t)$; we are in the classical ‘Cramér setting’ [12, Section 2.2] if it is assumed that $M(\vartheta)$ is finite in a neighborhood around the origin (which requires that the same property is satisfied by $\beta(\cdot)$). Let $I(a)$ and $\vartheta^{\star}\equiv\vartheta^{\star}(a)$, respectively, be defined as

[TABLE]

with $M(\cdot)$ as above. Using ‘Cramér’, we obtain that, under mild conditions,

[TABLE]

More refined asymptotics are available as well; we get back to this issue in Section 2.3.

As these results apply in the regime that $n$ is large, a relevant issue concerns the development of efficient techniques to estimate $p_{n}(a)$ through simulation. An important rare-event simulation technique is importance sampling, relying on the commonly used exponential twisting technique. We now investigate how to construct the exponentially twisted version ${\mathbb{Q}}$ (with twist $\vartheta^{\star}$) of the original probability measure ${\mathbb{P}}.$ The main idea is that under ${\mathbb{Q}}$ the $X_{i}(t)$ have mean $a$, such that under the new measure the event under study is not rare anymore.

More concretely, exponential twisting with parameter $\vartheta^{\star}$ means that under the new measure ${\mathbb{Q}}$, the $X_{i}(t)$ should have the mgf

[TABLE]

under this choice the random variable has the desired mean:

[TABLE]

The question is now: how to sample a random variable that has this mgf? To this end, notice that $M(\vartheta)=\exp(-\lambda t+\lambda t\,{\mathbb{E}}\exp(\vartheta B{\rm e}^{-rU}))$ and

[TABLE]

such that (3) equals

[TABLE]

From this expression we can see how to sample the $X_{i}(t)$ under ${\mathbb{Q}}$, as follows. In the first place we conclude that under ${\mathbb{Q}}$ the number of arrivals becomes Poisson with mean

[TABLE]

rather than $\lambda t$ (which is an increase). Likewise, it entails that under ${\mathbb{Q}}$ the distribution of the $B_{j}{\rm e}^{-rU_{j}}$ should be twisted by $\vartheta^{\star}$, in the sense that these random variables should have under ${\mathbb{Q}}$ the mgf

[TABLE]

We now point out how such a random variable should be sampled. To this end, observe that

[TABLE]

so that

[TABLE]

From this representation two conclusions can be drawn. In the first place, supposing there are $k$ arrivals, then the arrival epochs $U_{1},\ldots,U_{k}$ are i.i.d. under ${\mathbb{Q}}$, with the density given by

[TABLE]

In the second place, given that the $k$-th arrival occurs at time $u$, the density of the corresponding job size $B_{k}$ should be exponentially twisted by ${\rm e}^{-ru}\,\vartheta^{\star}$ (where each of the job sizes is sampled independently of everything else).

Now that we know how to sample from ${\mathbb{Q}}$ it is straightforward to implement the importance sampling. Before we describe its complexity (in terms of the number of runs required to obtain an estimate with given precision), we first provide an example in which we demonstrate how the change of measure can be performed.

Example 1**.**

In this example we consider the case that the $B_{i}$ are exponentially distributed with mean $\mu^{-1}$. Applying the transformation $w:={\rm e}^{-ru}\,\vartheta/\mu,$ it is first seen that

[TABLE]

As $\vartheta^{\star}$ solves the equation $M^{\prime}(\vartheta^{\star})/M(\vartheta^{\star})=a$, we obtain the quadratic equation

[TABLE]

leading to

[TABLE]

(where it is readily checked that $\vartheta^{\star}\in(0,\mu)$).

Now we compute what the alternative measure ${\mathbb{Q}}$ amounts to. In the first place, the number of arrivals should become Poisson with parameter

[TABLE]

(which is larger than $\lambda t$). In addition, we can check that

[TABLE]

(rather than $u/t$). The function $F_{U}^{\mathbb{Q}}(u)$ has the value [math] for $u=0$ and the value $1$ for $u=t$, and is concave. This concavity reflects that the arrival epochs of the shots tend to be closer to [math] under ${\mathbb{Q}}$ than under ${\mathbb{P}}$. This is because we identified each of the $U_{i}$ with $t$ minus the actual corresponding arrival epoch in (1); along the most likely path of $Y_{n}(t)$ itself the shots will be typically closer to $t$ under ${\mathbb{Q}}$. Observe that one can sample $U$ under ${\mathbb{Q}}$ using the classical inverse distribution function method [1, Section II.2a]: with $H$ denoting a uniform number on $[0,1)$, we obtain such a sample by

[TABLE]

Also, conditional on a $U_{i}$ having attained the value $u$, the jobs $B_{i}$ should be sampled from an exponential distribution with mean $(\mu-{\rm e}^{-ru}\,\vartheta^{\star})^{-1}$.

Remark 1**.**

In the model we study in this section, the input of the linear stochastic fluid network is a compound Poisson process. As pointed out in [19] the class of inputs can be extended to the more general class of increasing Lévy processes in a straightforward manner. **

2.3. Efficiency properties of importance sampling procedure

In this subsection we analyze the performance of the procedure introduced in the previous section. The focus is on a characterization of the number of runs needed to obtain an estimate with a given precision (at a given confidence level).

In every run $Y_{n}(t)$ is sampled under ${\mathbb{Q}}$, as pointed out above. As ${\mathbb{Q}}$ is an implementation of an exponential twist (with twist $\vartheta^{\star}$), the likelihood ratio (of sampling $Y_{n}(t)$ under ${\mathbb{P}}$ relative to ${\mathbb{Q}}$) is given by

[TABLE]

In addition, define $I$ as the indicator function of the event $\{Y_{n}(t)\geqslant na\}$. Clearly, ${\mathbb{E}}_{\mathbb{Q}}(LI)=p_{n}(a)$. We keep generating samples $LI$ (under ${\mathbb{Q}}$), and estimate $p_{n}(a)$ by the corresponding sample mean, until the ratio of the half-width of the confidence interval (with critical value $T$) and the estimator drops below some predefined $\varepsilon$ (say, 10%). Under ${\mathbb{P}}$ the number of runs needed is effectively inversely proportional to $p_{n}(a)$, hence exponentially increasing in $n$. We now focus on quantifying the reduction of the number of runs when using the importance sampling procedure we described above, i.e., the one based on the measure ${\mathbb{Q}}$.

Using a Normal approximation, it is a standard reasoning that when performing $N$ runs the ratio of the half-width of the confidence interval and the estimator is approximately

[TABLE]

and hence the number of runs needed is roughly

[TABLE]

We now analyze how $\Sigma_{n}$ behaves as a function of the ‘rarity parameter’ $n$. Due to the Bahadur-Rao result [4], with $f_{n}\sim g_{n}$ denoting $f_{n}/g_{n}\to 1$ as $n\to\infty$,

[TABLE]

Using the same proof technique as in [4], it can be shown that

[TABLE]

see Appendix A for the underlying computation. It also follows that ${\mathbb{E}}_{\mathbb{Q}}(L^{2}I)\sim{\mathbb{V}}{\rm ar}_{\mathbb{Q}}(L^{2}I)$.

We can use these asymptotics, to conclude that under ${\mathbb{Q}}$ the number of runs required grows slowly in $n$. More specifically, $\Sigma_{n}$ is essentially proportional to $\sqrt{n}$ for $n$ large. This leads to the following result; cf. [5, Section 2] for related findings in a more general context.

Proposition 1**.**

As $n\to\infty$,

[TABLE]

2.4. Simulation experiments

In this subsection we present numerical results for the single-node model without Markov modulation. We focus on the case of exponential jobs, as in Example 1. We simulate until the estimate has reached the precision $\varepsilon=0.1$, with confidence level $0.95$ (such that the critical value is $T=1.96$). The parameters chosen are: $t=1$, $r=1$, $\lambda=1$, and $\mu=1$. We set $a=1$ (which is larger than $m(t)=1-e^{-1}$). As it turns out, $\vartheta^{\star}=0.2918$ and

[TABLE]

The top-left panel of Fig. 1 confirms the exponential decay of the probability of interest, as a function of $n$. In the top-right panel we verify that the number of runs indeed grows proportionally to $\sqrt{n}$; the value of $\alpha$, as defined in (7), is $198.7$, which is depicted by the horizontal line. The bottom-left panel shows the density of the arrival epochs, which confirms that the arrival epochs tend to be closer to [math] under ${\mathbb{Q}}$ than under ${\mathbb{P}}$; recall that under ${\mathbb{P}}$ these epochs are uniformly distributed on $[0,t]$. Recall that we reversed time in (1): for the actual shot-noise system that we are considering, it means that in order to reach the desired level at time $t$, the arrival epochs tend to be closer to $t$ under ${\mathbb{Q}}$ than under ${\mathbb{P}}$. The bottom-right panel presents the rate of the exponential job sizes as a function of $u$. Using (4), the arrival rate under ${\mathbb{Q}}$ turns out to be $1.2315$.

3. Multi-node systems, no modulation

In this section we consider multi-node stochastic fluid linear stochastic fluid networks, of the type analyzed in the work by Kella and Whitt [19]. It is instructive to first consider the simplest multi-node system: a tandem network without external input in the downstream node and no traffic leaving after having been served by the upstream node (and rate $r_{\ell}$ for node $\ell$, $\ell=1,2$); later we extend the ideas developed to general linear stochastic fluid networks.

3.1. Preliminaries

As mentioned above, we first consider the two-node tandem. The content of the first node is, as before,

[TABLE]

(with $N$ having a Poisson distribution with mean $\lambda t$), but it can be argued that the content of the second node satisfies a similar representation. More specifically, using the machinery developed in [19], it turns out that

[TABLE]

As before, perform the scaling by $n$, meaning that the arrival rate $\lambda$ is inflated by a factor $n$. It leads to the random vectors $(X^{(1)}_{1}(t),\ldots,X^{(1)}_{n}(t))$ and $(X^{(2)}_{1}(t),\ldots,X^{(2)}_{n}(t))$. With these vectors we can define $Y_{n}^{(1)}(t)$ and $Y_{n}^{(2)}(t)$, analogously to how this was done in the single-node case; these two random quantities represent the contents of the upstream resource and the downstream resource, respectively.

The state of this tandem system can be uniquely characterized in terms of its (bivariate) moment generating function. The technique to derive an explicit expression is by relying on the above distributional equality (8). Again, the key step is to condition on the number of shots that have arrived in the interval $[0,t]$: with ${\boldsymbol{\vartheta}}=(\vartheta_{1},\vartheta_{2})$,

[TABLE]

The above computation is for the two-node tandem system, but the underlying procedure can be extended to the case of networks with more than 2 nodes, and external input in each of the nodes. To this end, we consider the following network consisting of $L$ nodes. Jobs are generated according to a Poisson process. At an arrival epoch, an amount is added to the content of each of the resources $\ell\in\{1,\ldots,L\}$, where the amount added to resource $\ell$ is distributed as the (non-negative) random variable $B^{(\ell)}$; $\beta({\boldsymbol{\vartheta}})$, with ${\boldsymbol{\vartheta}}\in{\mathbb{R}}^{L}$, is the joint mgf of $B^{(1)}$ up to $B^{(L)}$ (note that the components are not assumed independent). In addition, let the traffic level at node $\ell$ decay exponentially with rate $r_{\ell}$ (i.e., the value of the output rate is linear in the current level, with proportionality constant $r_{\ell}$). A deterministic fraction $p_{\ell\ell^{\prime}}\geqslant 0$ ($\ell\not=\ell^{\prime}$) is then fed into node $\ell^{\prime}$, whereas a fraction $p_{\ell\ell}\geqslant 0$ leaves the network (with $\sum_{\ell^{\prime}=1}^{L}p_{\ell\ell^{\prime}}=1$). We denote $r_{\ell\ell^{\prime}}:=r_{\ell}p_{\ell\ell^{\prime}}.$ As an aside we mention that this general model covers models in which some arrivals (of the Poisson process with parameter $\lambda$) actually lead to arrivals at only a subset of the $L$ queues (since the job sizes $B^{(1)},\ldots,B^{(L)}$ are allowed to equal [math]).

We now point out how the joint buffer content process can be analyzed. Again our objective is to evaluate the moment generating function. Define the matrix $R$ as follows: its $(\ell,\ell)$-th entry is $r_{\ell\ell}+\sum_{\ell^{\prime}\not=\ell}r_{\ell\ell^{\prime}}$, whereas its $(\ell,\ell^{\prime})$-th entry (with $\ell\not=\ell^{\prime}$) is $-r_{\ell\ell^{\prime}}$. We have, according to Kella and Whitt [19], with $N$ again Poisson with mean $\lambda t$, the following distributional equality: for any $\ell\in\{1,\ldots,L\}$,

[TABLE]

It means we can compute the joint mgf of $X^{(1)}(t)$ up to $X^{(L)}(t)$ as follows, cf. [19, Thm. 5.1]:

[TABLE]

which is the matrix/vector-counterpart of the expression (2) that we found in the single-node case; for the two-node case the special form (9) applies.

3.2. Tail probabilities, change of measure

In this subsection we introduce the change of measure that we use in our importance sampling approach. Many of the concepts are analogous to concepts used for the single-node case in Section 2.

Define (in self-evident notation)

[TABLE]

Due to the multivariate version of Cramér’s theorem, with $A:=[a_{1},\infty)\times\cdots\times[a_{L},\infty)$,

[TABLE]

More refined asymptotics than the logarithmic asymptotics of (10) are available as well, but these are not yet relevant in the context of the present subsection; we return to these ‘exact asymptotics’ in Section 3.3.

We assume that the set $A$ is ‘rare’, in the sense that

[TABLE]

Let us now construct the importance sampling measure. Let ${\boldsymbol{\vartheta}}^{\star}$ be the optimizing ${\boldsymbol{\vartheta}}$ in the decay rate of $p_{n}({\boldsymbol{a}}).$ Mimicking the reasoning we used in the single-node case, the number of arrivals becomes Poisson with mean

[TABLE]

The arrival epochs should be drawn using the density

[TABLE]

Given an arrival at time $u$, $(B^{(1)},\ldots,B^{(L)})$ should be exponentially twisted by

[TABLE]

3.3. Efficiency properties of importance sampling procedure

We now consider the efficiency properties of the change of measure proposed in the previous subsection. To this end, we first argue that the vector ${\boldsymbol{\vartheta}}$ generally has some (at least one) strictly positive entries, whereas the other entries equal 0; i.e., there are no negative entries. To this end, we first denote by ${\boldsymbol{b}}^{\star}$ the ‘most likely point’ in $A$:

[TABLE]

so that ${\boldsymbol{\vartheta}}^{\star}={\boldsymbol{\vartheta}}({\boldsymbol{b}}).$ It is a standard result from convex optimization that

[TABLE]

Suppose now that $\vartheta_{i}({\boldsymbol{b}}^{\star})<0$. Increasing the $i$-th component of the ${\boldsymbol{b}}^{\star}$ (while leaving all other components unchanged) would lead to a vector that is still in $A$, but that by virtue of (11) corresponds to a lower value of the objective function $I(\cdot)$, thus yielding that ${\boldsymbol{b}}$ was not the optimizer; we have thus found a contradiction. Similarly, when $\vartheta_{i}({\boldsymbol{b}}^{\star})=0$ we have that $b^{\star}_{i}>a_{i}$ (as otherwise a reduction of the objective function value would be possible, which contradicts ${\boldsymbol{b}}^{\star}$ being minimizer).

Now define $\Theta$ as the subset of $i\in\{1,\ldots,L\}$ such that $\vartheta_{i}>0$, and let $D\in\{1,\ldots,L\}$ the number of elements of $\Theta.$ We now argue that the number of runs needed to obtain an estimate of predefined precision scales as $n^{D/2}.$ Relying on the results from [11] (in particular their Thm. 3.4), it follows that $p_{n}({\boldsymbol{a}})$ behaves (for $n$ large) proportionally to $n^{-D/2}\exp(-nI({\boldsymbol{b}}))$; using the same machinery, ${\mathbb{E}}_{\mathbb{Q}}(L^{2}I)$ behaves proportionally to $n^{-D/2}\exp(-2nI({\boldsymbol{b}}))$. Mimicking the line of reasoning of Section 2.3, we conclude that the number of runs needed is essentially proportional to $n^{D/2}$. The formal statement is as follows.

Proposition 2**.**

As $n\to\infty$,

[TABLE]

where $\tau$ is the determinant of the Hessian of $\log M({\boldsymbol{\vartheta}})$ in ${\boldsymbol{\vartheta}}^{\star}.$

We further illustrate the ideas and intuition behind the qualitative result described in the above proposition by considering the case $L=2.$ It is noted that three cases may arise: (i) $\Theta=\{1,2\}$, (ii) $\Theta=\{1\}$, (iii) $\Theta=\{2\}$; as case (iii) can be dealt with in the same way as case (ii), we concentrate on the cases (i) and (ii) only. In case (i), where $D=2$, the necessary condition [11, Eqn. (3.4)] is fulfilled as ${\boldsymbol{\vartheta}}>0$ componentwise. As in addition the conditions A–C of [11] are in place, it is concluded that [11, Thm. 3.4] can be applied, leading to ${\boldsymbol{b}}^{\star}={\boldsymbol{a}},$ and

[TABLE]

where $\tau$ is the determinant of the Hessian of $\log M({\boldsymbol{\vartheta}})$ in ${\boldsymbol{\vartheta}}^{\star}.$ Along the same lines, it can be shown that

[TABLE]

It now follows that $\Sigma_{n}$ is roughly linear in $n$: with $\varepsilon$ and $T$ as introduced in Section 2.3,

[TABLE]

In case (ii), we do not have that ${\boldsymbol{\vartheta}}>0$ componentwise, and hence [11, Thm. 3.4] does not apply; in the above terminology, $D=1<2=L.$ Observe that in this case the exponential decay rate of the event $\{Y^{(1)}_{n}(t)\geqslant na_{1},Y^{(2)}_{n}(t)<na_{2}\}$ strictly majorizes that of $\{Y^{(1)}_{n}(t)\geqslant na_{1}\}$ (informally: the former event is substantially less likely than the latter). It thus follows that $b_{1}^{\star}=a_{1}$ and $b_{2}^{\star}>a_{2}$, and

[TABLE]

and in addition

[TABLE]

As a consequence in this regime $\Sigma_{n}$ grows essentially proportional to $\sqrt{n}$ for $n$ large:

[TABLE]

In case (iii) $\Sigma_{n}$ behaves proportionally to $\sqrt{n}$ as well.

3.4. Simulation experiments

We conclude this section by providing a few numerical illustrations. In the first set we focus on the downstream queue only (i.e., we analyze $p_{n}(0,a_{2})$), whereas in the second set we consider the joint exceedance probability $p_{n}({\boldsymbol{a}}).$ The precision and confidence have been chosen as in Example 1. Throughout we take $t=1$, $r_{1}=2$, $r_{2}=1$, $\lambda=1$, and $\mu=1$.

In the first set of experiments we take $a_{1}=0$ and $a_{2}=1$. Elementary numerical analysis yields that $\vartheta^{\star}=0.8104$ and $\tau=1.4774$, leading to $\alpha$, as defined in (13), equalling 474.3. The four panels of Fig. 2 should be interpreted as their counterparts in Fig. 1. Interestingly, the bottom-left panel indicates that in the tandem system it does not pay off to let jobs arrive right before $t$ (as they first have to go through the first resource to end up in the second resource), as reflected by the shape of the density of the arrival epochs under ${\mathbb{Q}}$; to this end, recall that we reversed time in (8), so that a low density at $u=0$ in the graph corresponds to a high density at $u=t$ in the actual system. The arrival rate under ${\mathbb{Q}}$ is $1.5103$.

In the second set of experiments we take $a_{1}=1.2$ and $a_{2}=1.1$; all other parameters are the same as in the first set. As mentioned above, we now consider the joint exceedance probability. As it turns out, $\vartheta_{1}^{\star}=0.1367$ and $\vartheta^{\star}_{2}=0.2225$. The top-right panel of Fig. 3 shows that for this specific parameter setting $\Sigma_{n}/n$ converges to the limiting constant rather slowly. Concerning the bottom-left panel, note that in Section 2 we saw that to make sure the first queue gets large it helps to have arrivals at the end of the interval, whereas above we observed that to make the second queue large arrivals should occur relatively early. We now focus on the event that both queues are large, and consequently the arrival distribution becomes relatively uniform again, as shown in the bottom-left panel. The arrival rate under ${\mathbb{Q}}$ is $2.3478$.

4. Multi-node systems under Markov modulation

In this section consider the networks analyzed in the previous section, but now in a random environment. More specifically, the type of random environment we focus on here is known as Markov modulation: the system dynamics are affected by the state of an external finite-state irreducible Markov process $J(\cdot)$ with generator matrix $Q=(q_{jj^{\prime}})_{j,j^{\prime}=1}^{d}$. When this Markov process (usually referred to as the background process) is in state $j$, arrivals occur according to a Poisson process with rate $\lambda_{j}$, the mgf of the job size is $\beta_{j}({\boldsymbol{\vartheta}})$, and the routing matrix is $R_{j}$. Analogously to the definitions used in the case without Markov modulation, this routing matrix’ $(i,i)$-th entry is

[TABLE]

which can be interpreted as the rate at which fluid leaves server $i$ when the background process is in $j$. Likewise, for $i\not=i^{\prime}$,

[TABLE]

which is the rate at which fluid flows from server $i$ to $i^{\prime}$ when the background process is in $j$.

Below we assume that $J(0)=j_{0}$ for a fixed state $j_{0}\in\{1,\ldots,d\}$; it is seen that all results generalize to an arbitrary initial distribution in a straightforward manner.

The structure of the section is in line with that of the previous two sections: we subsequently describe general results for the model under consideration (extending earlier results from [18]), propose an importance sampling measure, establish efficiency properties of the corresponding estimator, and present a number of numerical experiments.

4.1. Exact expressions for moments

Before focusing on simulation-based techniques, this subsection (which can be read independently of the rest of the section) shows that various moment-related quantities can be computed in closed form.

Multi-node systems under Markov modulation have been studied in detail by Kella and Stadje [18]. We start this subsection by providing a compact derivation of a pde characterizing the system’s transient behavior, which was not included in that paper. To this end, we define, for $j\in\{1,\ldots,d\}$,

[TABLE]

with $1_{j}(t)$ the indicator function of the event that $J(t)=j$. Using the standard ‘Markov machinery’, $\Xi_{j}({\boldsymbol{\vartheta}},t+\Delta t)$ equals (up to $o(\Delta t)$ terms) the sum of a contribution

[TABLE]

due to the scenario that an arrival occurs between $t$ and $t+\Delta t$, a contribution

[TABLE]

due to the scenario that the modulating Markov process jumps between $t$ and $t+\Delta t$, and a contribution

[TABLE]

with $q_{j}:=-q_{jj}$; regarding the last term, observe that when the background process is in state $j$, and no new job arrives between $t$ and $t+\Delta t$,

[TABLE]

We thus find that

[TABLE]

This immediately leads to (by subsequently subtracting $\Xi_{j}({\boldsymbol{\vartheta}},t)$ from both sides, dividing by $\Delta t$, and letting $\Delta t\downarrow 0$)

[TABLE]

Let us now compactly summarize the relation (14), in vector/matrix notation. This notation will prove practical when computing higher moments; in other (but related) contexts, similar procedures have been proposed in e.g. [15, 22]. Let ${\mathscr{M}}^{n_{1}\times n_{2}}$ be the set of ${\mathbb{R}}$-valued matrices of dimension $n_{1}\times n_{2}$ (for generic $n_{1},n_{2}\in{\mathbb{N}}$). In addition, $I_{n}$ is the identity matrix of dimension $n\in{\mathbb{N}}$. We introduce the following three matrices in ${\mathscr{M}}^{d\times d}$, ${\mathscr{M}}^{d\times d}$, and ${\mathscr{M}}^{Ld\times Ld}$, respectively:

[TABLE]

We use the conventional notation $\otimes$ for the Kronecker product. Recall that the Kronecker product is bilinear, associative and distributive with respect to addition; these properties we will use in the sequel without mentioning. It also satisfies the mixed product property $(A\otimes B)(C\otimes D)=(AC)\otimes(BD)$. Furthermore, note that $I_{n_{1}}\otimes I_{n_{2}}=I_{n_{1}n_{2}}$.

We now consider some differentiation rules for matrix-valued functions which will allow us to iteratively evaluate moments. In the first place we define the operator $\nabla_{\boldsymbol{\vartheta}}$ for ${\boldsymbol{\vartheta}}\in{\mathbb{R}}^{L}$; to keep notation compact, we often suppress the subscript ${\boldsymbol{\vartheta}}$, and write just $\nabla$. Let $f\equiv f({\boldsymbol{\vartheta}})$ be a mapping of ${\mathbb{R}}^{L}$ to $\mathscr{M}^{n_{1}\times n_{2}}$. Then $\nabla f\equiv\nabla f({\boldsymbol{\vartheta}})\in{\mathscr{M}}^{n_{1}L\times n_{2}}$ is defined by

[TABLE]

In the above definition $\nabla f_{ij}\equiv\nabla f_{ij}({\boldsymbol{\vartheta}})$ is to be understood as the usual gradient; the symbol $\partial_{i}$ is used to denote the partial derivative with respect to the $i$-th variable, in the sense of

[TABLE]

Furthermore, we define inductively $\nabla^{k}f\equiv\nabla^{k}f({\boldsymbol{\vartheta}}):=\nabla(\nabla^{k-1}f)$, $k\in{\mathbb{N}}$, with $\nabla^{0}f:=f$. It is checked that $\nabla^{k}f({\boldsymbol{\vartheta}})$ is a mapping of ${\mathbb{R}}^{L}$ to $\mathscr{M}^{L^{k}n_{1}\times n_{2}}$.

In the sequel we use a couple of differentiation rules, that we have listed below. Let $A(\cdot)$ be a matrix-valued function from ${\mathbb{R}}^{L}$ to ${\mathscr{M}}^{n_{1}\times n_{2}}$, and $B(\cdot)$ a matrix-valued function from ${\mathbb{R}}^{L}$ to ${\mathscr{M}}^{n_{2}\times n_{3}}$, and let $I_{q}$ be a $q\times q$ identity matrix (for some $q\in\mathbb{N}$). Then,

• –

Product rule:

[TABLE]

being an element of ${\mathscr{M}}^{Ln_{1}\times n_{3}}$.

• –

Differentiation of Kronecker product (1):

[TABLE]

• –

Differentiation of Kronecker product (2):

[TABLE]

where $K_{m,n}$ is the commutation matrix defined by

[TABLE]

and $H_{ij}\in\mathscr{M}^{m\times n}$ denotes a matrix with a $1$ at its $(i,j)$-th position and zeros elsewhere, cf. [21].

The first rule can be checked componentwise and the second rule is trivial. The third rule follows from the first and second rule in combination with the fact that the Kronecker product commutes after a correction with the commutation matrices. Moreover, we use the property $K_{m,n}^{-1}=K_{n,m}$. An overview of the properties of commutation matrices can be found in [21].

In the introduced terminology, it follows that (14) can be written as

[TABLE]

We now point out how (transient and stationary) moments can be evaluated; note that [18] focuses on the first two stationary moments at epochs that the background process jumps. We throughout use the notation ${\boldsymbol{z}}_{i}(t)$ for the $i$-th derivative of ${\boldsymbol{\Xi}}({\boldsymbol{\vartheta}},t)$ in $({\boldsymbol{0}},t)$, for $t\geqslant 0$:

[TABLE]

for $i\in{\mathbb{N}}$. Note that, with $\pi_{j}(t)=(\exp(Qt))_{j_{0},j}$,

[TABLE]

$\circ$   We start by characterizing the first moments. Applying the operator $\nabla\equiv\nabla_{\boldsymbol{\vartheta}}$ to the differential equation (15) yields

[TABLE]

using standard properties of the Kronecker product in combination with

[TABLE]

where $\boldsymbol{e}_{i}$ denotes the $L$-dimensional column vector in which component $i$ equals $1$ and all other components are [math]. Then, inserting ${\boldsymbol{\vartheta}}=\boldsymbol{0}$ into (16) yields the system of (non-homogeneous) linear differential equations

[TABLE]

In the stationary case, we obtain

[TABLE]

with $\boldsymbol{\pi}:=\lim_{t\to\infty}\boldsymbol{\pi}(t)$ (being the unique non-negative solution of ${\boldsymbol{\pi}}^{\rm T}Q={\boldsymbol{0}}^{\rm T}$ such that the entries of ${\boldsymbol{\pi}}$ sum to $1$).

$\circ$   We now move to second moments. Applying the $\nabla_{\boldsymbol{\vartheta}}$-operator to (16),

[TABLE]

in which the factor $\nabla_{\boldsymbol{\vartheta}}(\Lambda B(\boldsymbol{\vartheta})\otimes I_{L})$ can be expressed more explicitly as

[TABLE]

and the factor $\nabla_{\boldsymbol{\vartheta}}(((I_{d}\otimes\boldsymbol{\vartheta}^{\rm T})R^{\rm T})\otimes I_{L})$ simplifies to $(K_{d,L}\otimes I_{L})K_{L,dL}(R^{\rm T}\otimes I_{L})$. Inserting $\boldsymbol{\vartheta}=0$ yields the system of linear differential equations

[TABLE]

where $\boldsymbol{z}_{1}(t)$ solves Eqn. (17). As before, the stationary quantities can be easily derived (by equating $\boldsymbol{z}_{2}^{\prime}(t)$ to [math]). One has to keep in mind, however, that some of the mixed partial derivatives occur multiple times in $\boldsymbol{z}_{k}$, for $k\in\{2,3,\ldots\}$, and therefore the solution will only be unique after removing the corresponding redundant rows. Alternatively, the system can be completed by including equations which state that these mixed partial derivatives are equal.

$\circ$   It now follows that higher moments can be found recursively, using the three differentiation rules that we stated above.

Remark 2**.**

Various variants of our model can be dealt with similarly. In this remark we consider the slightly adapted model in which shots only occur simultaneously with a jump in the modulating Markov chain. Then (up to $o(\Delta t)$ terms) $\Xi_{j}(\boldsymbol{\vartheta},t+\Delta t)$ is the sum of a contribution

[TABLE]

due to the scenario that there is a jump in the modulating chain in the interval $[t,t+\Delta t]$ (which also induces a shot), and a contribution of

[TABLE]

with $q_{j}:=-q_{jj}$, in the scenario that there is no jump. Performing the same steps as above, we obtain

[TABLE]

which has a similar structure as (14). It follows that the moments can be found as before. With $\widetilde{Q}:={\rm diag}\{q_{1},\ldots,q_{d}\}$, it turns out that the transient means are given by

[TABLE]

In particular, the stationary first moment equals

[TABLE]

Consider the following numerical example for the computation of the expected values and variances, in which the technique described above is illustrated.

Example 2**.**

In this example, we choose the parameters in such a way that we see non-monotonic behavior. Our example corresponds to a case in which the system does not start empty, which is dealt with by imposing suitable starting conditions. We consider a two-dimensional ($L=2$) queueing system, with a two-dimensional state space of the Markov modulating process ($d=2$). We pick $q_{12}=q_{21}=1$, $\lambda_{1}=\lambda_{2}=1$, ${\mathbb{E}}\,B_{1}={\mathbb{E}}\,B_{2}={\mathbb{E}}\,B_{1}^{2}={\mathbb{E}}\,B_{2}^{2}=1$, $J(0)=1$, $(X_{1}(0),X_{2}(0))=(3,3)$, and the rate matrices

[TABLE]

Due to the symmetry in the choice of the parameters, one can expect that for both states of the background process expected value tends (as $t$ grows large) to the same steady-state value; the same is anticipated for the stationary variance. This is confirmed by Fig. 4. For $t$ small, the two queues behave differently due to $J(0)=1$, which implies that queue $1$ drains faster. Note that ${\mathbb{E}}\,X_{2}(t)$ even increases for $t$ small, due to the fact that the flow from node 1 to 2 equals the flow from 2 to 1, constituting a net flow of zero, so that the additional contribution of external output to node 2 leads to a net increase of ${\mathbb{E}}\,X_{2}(t)$. The transient correlation is plotted in Fig. 5. At time $t=0$ the queues are perfectly correlated, since the starting state is deterministic. Then the correlation decreases due to the asymmetric flow rates until around $t=1$, which is when the Markov chain $J$ is expected to switch, after which the correlation monotonously tends to the steady state.**

4.2. Tail probabilities, change of measure

We now characterize the decay rate of the rare-event probability under study, and we propose a change of measure to efficiently estimate it. In the notation we have been using so far, we again focus on

[TABLE]

where ${\boldsymbol{Y}}_{n}(t)=(Y_{n}^{(1)}(t),\ldots,Y_{n}^{(L)}(t)).$ It is stressed that, following [8], we consider the regime in which the background process is ‘slow’. In concrete terms, this means that we linearly speed up the driving Poisson process (i.e., we replace the arrival rates $\lambda_{j}$ by $n\lambda_{j}$), but leave the rates of the Markovian background process unaltered.

First we find an alternative characterization of the state of the system at time $t$. Let ${\mathscr{F}}_{t}$ denote the set of all functions from $[0,t]$ onto the states $\{1,\ldots,d\}$. Consider a path $f\in{\mathscr{F}}_{t}$. Let $f$ have $K{(f)}$ jumps between [math] and $t$, whose epochs we denote by $t_{1}{(f)}$ up to $t_{K{(f)}}{(f)}$ (and in addition $t_{0}{(f)}:=0$ and $t_{K{(f)}+1}{(f)}:=t$). Let

[TABLE]

(i.e., the state of $f$ immediately after the $i$-th jump). We also introduce

[TABLE]

Suppose now that the Markov process $J(\cdot)$ follows the path $f\in{\mathscr{F}}_{t}$. Then the contribution to the mgf of ${\boldsymbol{X}}(t)$ due to shots that arrived between $t_{i}(f)$ and $t_{i+1}(f)$ is, mimicking the arguments that we used in Section 3.2 for non-modulated networks,

[TABLE]

As a consequence, the mgf of $X(t)$ given the path $f$ is

[TABLE]

First conditioning on the path of $J(\cdot)\in{\mathscr{F}}_{t}$ between [math] and $t$ and then unconditioning, it then immediately follows that the mgf of ${\boldsymbol{X}}(t)$ is given by

[TABLE]

Then, precisely as is shown in [8] for a related stochastic system, the decay rate can be characterized as follows:

[TABLE]

The argumentation to show this is analogous to the one in [8, Thm. 1], and can be summarized as follows. In the first place, let $f^{\star}$ be the optimizing path in (19). Then, as $J(\cdot)$ does not depend on $n$, we can choose a ‘ball’ ${\mathscr{B}}_{t}({f^{\star}})$ around $f^{\star}$ such that the decay rate of the probability of $J(\cdot)$ being in that ball is [math]. The lower bound follows by only taking into account the contribution due to paths in ${\mathscr{B}}_{t}({f^{\star}})$. The upper bound follows by showing that the contribution of all $f\in{\mathscr{F}}_{t}\setminus{\mathscr{B}}_{t}({f^{\star}})$ is negligible.

Informally, the path $f^{\star}$ has the interpretation of the most likely path of $J(\cdot)$ given that the rare event under consideration happens. To make sure that the event under consideration is rare, we assume that for all $f\in{\mathscr{F}}_{t}$

[TABLE]

The change of measure we propose is the following. In every run we first sample the path $J(s)$ for $s\in[0,t]$ under the original measure ${\mathbb{P}}$ (i.e., with $J(0)=j_{0}$, and then using the generator matrix $Q$). We call the resulting path $f\in{\mathscr{F}}_{t}$. For this path, define ${\boldsymbol{\vartheta}}^{\star}_{f}\geqslant{\bf 0}$ as the optimizing $\boldsymbol{\vartheta}$ in the definition of ${\mathbb{I}}(f)$ in (19); ${\boldsymbol{b}}^{\star}_{f}\in A$ is the optimizing ${\boldsymbol{b}}.$

Conditional on the path $f$ of the background process, under the new measure ${\mathbb{Q}}$ the number of external arrivals between $t_{i}(f)$ and $t_{i+1}(f)$ is Poisson with parameter

[TABLE]

where $P_{i}(u,f):=D_{i}(u,f)D_{i+1}(f)\cdots D_{K(f)}(f).$ The arrival epochs between $t_{i}(f)$ and $t_{i+1}(f)$ should be drawn using the density

[TABLE]

Given an arrival at time $u$ between $t_{i}(f)$ and $t_{i+1}(f)$, the job sizes $(B^{(1)},\ldots,B^{(L)})$ should be sampled from a distribution with mgf $\beta_{j_{i}(f)}({\boldsymbol{\vartheta}})$, but then exponentially twisted by

[TABLE]

Remark 3**.**

As mentioned above, the background process is sampled under the original measure, whereas an alternative measure is used for the number of arrivals, the arrival epochs, and the job sizes. The intuition behind this, is that the rare event under consideration is caused by two effects:

• $\circ$

In the first place, samples of the background process $J$ should be close to $f^{\star}$. Under ${\mathbb{P}}$ a reasonable fraction ends up close to $f^{\star}$ — more precisely, the event of $J$ being close to $f^{\star}$ does not become increasingly rare when $n$ grows. As a consequence, no change of measure is needed here.

• $\circ$

In the second place, given the path of the background process, the $Y^{(\ell)}_{n}(t)$ should exceed the values $na_{\ell}$, for $\ell=1,\ldots,L$. This event does become exponentially rare as $n$ grows, so importance sampling is to be applied here.

4.3. Efficiency properties of importance sampling procedure

We now analyze the speed up realized by the change of measure introduced in the previous subsection. Unlike our results for the non-modulated systems, now we cannot find the precise rate of growth of $\Sigma_{n}.$ What is possible though, is proving asymptotic efficiency (also sometimes referred to as logarithmic efficiency), in the sense that we can show that

[TABLE]

(where the second equality is a consequence of (19)). This equality is proven as follows. As by Jensen’s inequality ${\mathbb{E}}_{\mathbb{Q}}(L^{2}I)\geqslant({\mathbb{E}}_{\mathbb{Q}}(LI))^{2}=(p_{n}({\boldsymbol{a}}))^{2},$ we are left to prove the upper bound:

[TABLE]

If the path of $J(\cdot)$ equals $f\in{\mathscr{F}}_{t}$, it follows by an elementary computation that we have constructed the measure ${\mathbb{Q}}$ such that

[TABLE]

The fact that ${\boldsymbol{\vartheta}}^{\star}_{f}$ is componentwise non-negative, in combination with the fact that ${\boldsymbol{Y}}_{n}(t)\geqslant{\boldsymbol{a}}$ when $I=1$, entails that

[TABLE]

noting that ${\boldsymbol{a}}$ and ${\boldsymbol{b}}^{\star}_{f}$ may only differ if the corresponding entry of ${\boldsymbol{\vartheta}}^{\star}_{f}$ equals [math] (that is, $\langle{\boldsymbol{a}}-{\boldsymbol{b}}_{f}^{\star},{\boldsymbol{\vartheta}}_{f}^{\star}\rangle=0$). The upper bound thus follows: with $f^{\star}$ the minimizing path in (19), recalling that $J(\cdot)$ is sampled under ${\mathbb{P}}$,

[TABLE]

We have established the following result.

Proposition 3**.**

As $n\to\infty$, the proposed importance sampling procedure is asymptotically efficient. This means that the number of runs needed grows subexponentially:

[TABLE]

Remark 4**.**

In the scaling considered, for both the logarithmic asymptotics of $p_{n}({\boldsymbol{a}})$ and our importance sampling algorithm, the precise transition rates $q_{ij}$ do not matter; the only crucial element is that the background process is irreducible. Observe that the probability of interest does depend on the rates $q_{ij}$, and so do the exact asymptotics. We refer to [7] for the exact asymptotics of a related infinite-server model; it is noted that the derivation of such precise asymptotics is typically highly involved.**

The above reasoning indicates that the proposed procedure remains valid under more general conditions: the ideas carry over to any situation in which the rates are piecewise constant along the most likely path. **

4.4. Simulation experiments

We performed experiments featuring a single-node system under Markov modulation. In our example the job sizes stem from an exponential distribution. When the background process is in state $i$, the arrival rate is $\lambda_{i}$, the job-size distribution is exponential with parameter $\mu_{i}$, and the rate at which the storage level decays is $r_{i}$, for $i\in\{1,\ldots,d\}.$

The change of measure is then implemented as follows. As pointed out in Section 4.2, per run a path $f$ of the background process is sampled under the original measure ${\mathbb{P}}.$ Suppose along this path there are $K$ transitions (remarking that, for compactness, we leave out the argument $f$ here), say at times $t_{1}$ up to $t_{K}$; with $t_{0}=0$ and $t_{K+1}=t$, the state between $t_{i}$ and $t_{i+1}$ is denoted by $j_{i}$, for $i=0,\ldots,K.$ Per run a specific change of measure is to be computed, parametrized by the $t_{i}$ and $j_{i}$, as follows.

We define

[TABLE]

the product in this expression should be interpreted as $1$ if $i+1>K$. It is readily checked that

[TABLE]

Let $\vartheta^{\star}$ be the maximizing argument of $\vartheta a-\log M(\vartheta)$.

We can now provide the alternative measure ${\mathbb{Q}}$ for this path of the background process. The number of arrivals between $t_{i}$ and $t_{i+1}$ (for $i=0,\ldots,K$) becomes Poisson with parameter

[TABLE]

(where it is noted that this expression is larger than $\lambda_{j_{i}}(t_{i+1}-t_{i})$, which was the parameter under ${\mathbb{P}}$). The density of each of the arrivals between $t_{i}$ and $t_{i+1}$ becomes

[TABLE]

(rather than a uniform distribution, as was the case under ${\mathbb{P}}$); sampling from this distribution is easy, since the inverse distribution function can be determined in closed form. Given an arrival that takes place at time $u$ between $t_{i}$ and $t_{i+1}$, the job size is exponential with parameter $\mu_{j_{i}}-P_{i}(u)\,\vartheta^{\star}$ (rather than exponential with parameter $\mu_{j_{i}}$).

We now describe two examples in which the dimension of the background process is $d=2$, $q_{12}=q_{21}=2$, and $t=1$. In the first example we fix $a=3$, ${\boldsymbol{\lambda}}=(2,1)$, ${\boldsymbol{\mu}}=(\frac{1}{2},1),$ and ${\boldsymbol{r}}=(5,1)$, in the second example $a=0.8$, ${\boldsymbol{\lambda}}=(0.9,1)$, ${\boldsymbol{\mu}}=(0.9^{-1},1),$ and ${\boldsymbol{r}}=(0.3,0.6)$. As before, we simulate until the precision of the estimate has reached $\varepsilon=0.1$. The top two panels in Figs. 6–7 should be read as those in Figs. 1–3; the bottom two panels correspond to the density of the arrival epochs and the rate of the exponential job sizes, respectively, for $f$ the ‘empirical maximizer’ of ${\mathbb{I}}_{f}(a)$ (i.e., the maximizer of ${\mathbb{I}}_{f}(a)$ over all paths $f$ of the background process that were sampled in the simulation experiment).

In the first example the thus obtained ‘optimal path’ subsequently visits states 1, 2, and 1, where the corresponding jump times are $t_{1}^{\star}=0.654$ and $t_{2}^{\star}=0.739$, and the decay rate is $0.573$. The mean numbers of arrivals in the three parts of the optimal path are $1.392$, $0.090$ and $0.963$ respectively, whereas for Monte Carlo sampling these are $1.308$, $0.085$ and $0.522$ respectively.

In the second example the optimal path subsequently visits states 2 and 1, where the corresponding jump time is $t_{1}^{\star}=0.790$. In this case the decay rate has the value $0.000806$. The mean numbers of arrivals in the two parts of the optimal path are $0.812$ and $0.195$ respectively, which are slightly higher than the corresponding values under Monte Carlo sampling ($0.790$ and $0.189$ respectively). Observe that in this example the difference between the two measures is relative small, also reflected by the small value of the decay rate; the event under consideration technically qualifies as ‘rare’ in that $p_{n}(0.8)\to 0$ as $n\to\infty$, but has a relatively high likelihood (e.g. as compared to the first example). As a consequence of the fact that both measures almost coincide, the two densities in the bottom-left panel can hardly be distinguished.

We observe that the top panels confirm that in both examples (i) $p_{n}(a)$ decays roughly exponentially in $n$, (ii) the number of runs needed grows roughly linearly in $n$ (in the first example slightly sublinearly).

5. Discussion and concluding remarks

In this paper we have considered the probability of attaining a value in a rare set $A$ at a fixed point in time $t$: with $A=[a_{1},\infty)\times\cdots\times[a_{L},\infty)$,

[TABLE]

A relevant related quantity is the probability of having reached the set $A$ before $t$:

[TABLE]

observe that this probability increases to $1$ as $t\to\infty$. Alternatively, one could study the probability that all $a_{\ell}$ (for $\ell=1,\ldots,L$) are exceeded before $t$, but not necessarily at the same time:

[TABLE]

Powerful novel sample-path large deviations results by Budhiraja and Nyquist [9], which deal with a general class of multi-dimensional shot-noise processes, may facilitate the development of efficient importance sampling algorithms for non-modulated linear stochastic fluid networks. The results in [9] do not cover Markov modulation, though.

In the current setup of Section 4 the speed of the background process is kept fixed, i.e., not scaled by $n$. For modulated diffusions a sample-path large deviation principle has been recently established in [16] for the case that the background process is sped up by a factor $n$ (which amounts to multiplying the generator matrix $Q$ by $n$); the rate function decouples into (i) a part concerning the rare-event behavior of the background process and (ii) a part concerning the rare-event behavior of the diffusion (conditional on the path of the background process). With a similar result for the Markov-modulated linear stochastic fluid networks that we have studied in this paper, one could potentially set up an efficient importance sampling procedure for the probabilities (20) and (21) under this scaling.

Appendix A

We here point out how (6) can be established; the line of reasoning is precisely the same as in the derivation of (5) in [12, Thm. 3.7.4]. First write

[TABLE]

which, with $Z_{n}:=(Y_{n}(t)-na)/\sqrt{n}$, equals

[TABLE]

Observe that ${\mathbb{E}}_{\mathbb{Q}}\,Y_{n}=na$, due to the very choice of ${\mathbb{Q}}.$ This entails that $Z_{n}$ converges in distribution to a centered Normal random variable; as can be verified, the corresponding variance is $\tau$ (where $\tau$ is defined in (5)). Using the Berry-Esseen-based justification presented in [12, page 111], we conclude that, as $n\to\infty$,

[TABLE]

Completing the square, the right-hand side of the previous display equals, with ${\mathscr{N}}(\mbox{\sc m},\mbox{\sc v})$ a normal random variable with mean m and variance v,

[TABLE]

Now we use the standard equivalence (as $x\to\infty$)

[TABLE]

to obtain

[TABLE]

Combining the above, we derive the claim:

[TABLE]