{'fraction_non_alphanumeric': 0.0533704583511819, 'fraction_numerical': 0.014252891467736282, 'mean_word_length': 3.98534265191225, 'pattern_counts': {'":': 0, '<': 40, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 24, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hypothesis testing for a L\'evy-driven storage system by Poisson sampling

This paper focuses hypothesis testing for the input of a L\'evy-driven storage system by sampling of the storage level. As the likelihood is not explicit we propose two tests that rely on transformation of the data. The first approach uses i.i.d. `quasi-busy-periods' between observations of zero workload. The distribution of the duration of quasi-busy-periods is determined. The second method is a conditional likelihood ratio test based on the Bernoulli events of observing a zero or positive workload, conditional on the previous workload. Performance analysis is presented for both tests along with speed-of-convergence results, that are of independent interest.


Introduction
In statistical hypothesis testing, one wishes to distinguish between a null hypothesis H 0 and an alternative hypothesis H 1 by observing a series of random variables. In the common setup the hypotheses directly relate to the observations: for instance, observing a series of random variables, the hypotheses could correspond to these observations stemming from two specific distributions. In the statistical literature a vast body of results has been established that facilitate such tests. In many applications, however, the hypotheses relate to the observations in a more involved manner. A prominent example of such a situation lies in the domain of storage systems and queues: one has periodic observations of the storage level, but the hypotheses are in terms of the system's input process. For instance, by periodically observing the workload, one would like to distinguish between two values of the arrival rate. While some work on hypothesis testing for this context has been done, a general framework is still lacking, and, as a consequence, various open questions remain. In this paper we consider the workload process of a rather broad set of storage systems. We focus on the situation of a resource that is fed by an increasing Lévy process (often referred to as a 'subordinator') which is uniquely characterized by its Laplace exponent ϕ(·), and that is emptied at a deterministic rate. This class of workload models covers the intensively studied storage system with compound Poisson input (often referred to as the M/G/1 queue), but it allows the driving Lévy process to be any subordinator (for instance a Gamma process or an inverse Gaussian process). The main objective of this paper is to develop methods for distinguishing between two characteristic exponent functions, say ϕ 0 (·) and ϕ 1 (·), based on observations of the corresponding workload process, rather than on observations of the Lévy input processes themselves. It is assumed that the workload is observed at Poisson instants. A major complication is that, although in our Lévy-input context the increments of the cumulative input process are independent, subsequent workload observations are not, so that the likelihood cannot be evaluated and consequently conventional tests cannot be applied.
Contributions. As mentioned, this paper develops tests for distinguishing between Laplace exponents ϕ 0 (·) and ϕ 1 (·) based on workload observations. We present two approaches that succeed in resolving the complications identified above. A key feature of both approaches is that the sequential test has power one; if the null hypothesis is wrong, then this will be detected with probability one as long as there is no restriction on the number of observations. However, a type-I error is still possible, i.e., a false rejection of the null hypothesis. To assess the performance of the test, one wishes to compute (or approximate) the type-I error probability.
• In the first approach we consider the so-called 'quasi busy period', which is defined as the number of Poisson observations until the workload hits 0 again. It is an inherent feature of the model that subsequent quasi busy periods are i.i.d. As a consequence, a traditional likelihood ratio test can be used. To evaluate the likelihood, the distribution of the duration of quasi busy periods should be determined; we point out how this can be done. In addition, we present results that assess the performance of the test. • The second approach works with a conditional likelihood ratio test. It is based on the fact that we can explicitly compute the probability of the events of observing a zero or positive workload, conditional on the value of the previous observation. In this setup the increments of the log-likelihood are not i.i.d., so that a rather delicate analysis is needed to analyze the performance of the test.
To quantify the performance of the test pertaining to the second approach, we require a number of results describing the convergence of a Lévy-driven storage system to its stationary version.
To this end, we derive a series of speed-of-convergence results, that are also of independent interest.
Background and related literature. In the situation we are considering the systems's input is a non-decreasing Lévy process minus a deterministic drift, while the workload is sampled according at Poisson epochs. This setup has been considered before in [12], where a method was developed for consistent and asymptotically normal semi-parametric estimation of the Laplace exponent, based on workload observations at Poisson epochs. The present paper can be seen as the hypothesis testing counterpart of [12]. The two approaches we pursue are attempts to remain as closely as possible to conventional likelihood ratio tests, both in terms of the evaluation of the likelihood and the assessment of the test's performance. However, as pointed out above, due to workload's intricate dynamics, various complications need to be overcome. For more general background on hypothesis testing we refer to the textbooks [13,17]. Our methods naturally extend to a stability test for a queue with unknown input as well as changepoint detection.
As mentioned above, in our storage system setting a general theory for hypothesis testing is lacking. A review of the classical literature on this topic is given in [4,Section 5]. Without attempting to provide an exhaustive overview, we mention a few specific references. If Markov chains can be embedded in the queueing process (which is the case in e.g. systems of the type M/G/1 and GI/M/s), and if the corresponding jump process is observed, the likelihood can be evaluated in closed-form [4, Section 5.e]; cf. the sequential test proposed in [3]. A large sample asymptotic test for the traffic intensity of a G/G/s system in which interarrival and service times are observed, is presented in [14]; relying on a delta-method type argument the authors construct a normal approximation for the error probabilities of the test. A test based on the distribution of the number of arrivals during a service period is proposed in [11].
Organization. This paper is organized as follows. Section 2 describes the setup considered in this paper, presents some preliminaries on Lévy-driven storage systems and hypothesis testing, and formally states our objectives. Then in Section 3 we detail the quasi-busy-period based approach, including the analysis of its performance. Then, in Section 4 we provide a series of results on convergence of the storage-level process to its stationary version, which are applied in Section 5 but which are relevant in their own right as well. Then Section 5 presents the conditional likelihood ratio approach and its performance analysis. Numerical experiments are described in Section 6, while Section 7 provides a brief discussion and some concluding remarks.

Model, preliminaries, objectives
In this section we provide a model description (and introduce the notation that will be used throughout the paper), present some preliminaries on Lévy-driven storage systems, and state our objectives.
2.1. Model. We consider a storage system fed by a non-decreasing Lévy process (J(t)) t 0 ≡ J(·). The output of the system is a unit-rate linear drift. The system's net input process (X(t)) t 0 ≡ X(·) is therefore given through X(t) = J(t) − t, which is a spectrally-positive Lévy process characterized by its Laplace exponent where ν is a Lévy jump measure such that ν(−∞, 0) = 0.
In the sequel we let (V (t)) t 0 ≡ V (·) denote the corresponding workload process. It can be represented as the net input process reflected at zero, in that V (t) = X(t) + max{V (0), L(t)}, where L(t) := − inf 0 s t X(s). Under the stability condition ϕ ′ (0) = −EX(1) > 0 the workload has a stationary distribution V := V (∞) with an LST given by the generalized Pollaczek-Khintchine formula [6, p. 27]: otherwise the workload process is unstable, meaning that V (∞) = ∞ almost surely. The first two moments of the stationary distribution are given by where Let ρ := EJ(1) = (0,∞) xν(dx) denote the expected input per unit of time. Then the expected net input per unit of time is Therefore, EX(1) < 0 is equivalent to ρ < 1, which is the typical form of the stability condition in queueing theory.

2.2.
Poisson sampling. In this paper the workload process will be sampled at Poisson epochs. Our analysis strongly relies on the availability of explicit expressions of the workload after an exponentially distributed time as a function of the initial workload level. The workload process starts, at time 0, at some known level V (0). It is observed according to an independent Poisson process with rate ξ > 0. Let T 1 , T 2 , . . . be i.i.d. exponentially distributed with parameter ξ > 0. We denote by S k := T 1 + . . . + T k the epoch of the k-th observation, i.e., S k has an Erlang distribution with scale parameter k and and shape parameter ξ. From now on we use the compact notation V k := V (S k ) to denote the workload process at the k-th Poisson epoch, and in addition V 0 := V (0). If the input process is spectrally positive, then the distribution of the workload at sample i ∈ {1, 2, . . .} conditional on the workload at sample i − 1 is characterized through where ψ(x) := ϕ −1 (x); see e.g. [8] and [6,Ch. IV]. If, furthermore, the input process of the queue is a subordinator, then the workload process attains the value zero with positive probability. This probability can be computed by taking α → ∞ in (3): This probability will be the main building block of the conditional likelihood ratio test that will be introduced in Section 5. In Section 3 we will present new results on the distribution of the number of samples taken between consecutive observations of zero workload. Note that when ξ ↓ 0, if ρ < 1 then ψ(ξ) → 0 whereas otherwise ψ(ξ) → ψ(0) > 0. If ρ < 1 then, for all v 0, This relation has an intuitive backing: when decreasing the sampling rate the events of finding an idle server become 'increasingly independent', and therefore the probability approaches the steady-state idle-server probability 1 − ρ.

Hypothesis testing.
Our goal is to provide a framework for testing hypotheses related to the distribution of the net input of the queue: with some abuse of notation, An important special case concerns tests that correspond to the traffic intensity: and in particular detecting stability if ϕ ′ 0 (0) > 0 and ϕ ′ 1 (0) 0. Another objective concerns providing a procedure for changepoint detection. In this context the system may start with Laplace exponent ϕ 0 (·), but a change may occur at a random time after which the Laplace exponent becomes ϕ 1 (·). The goal is to identify if (and when) this change occurs. In many applications one is interested in a dynamic test as observations are collected, as opposed to a static test with a given sample of size n. Then the test is defined by a test statistic g n (V 0 , . . . , V n ), in combination with two disjoint decision sets, G 0 and G 1 . Let N k := inf{n : g n (V 0 , . . . , V n ) ∈ G k } for k ∈ {0, 1}. The stopping rule is: collect samples until the stopping time N := min{N 0 , N 1 }, and reject (accept, respectively) the null hypothesis if N 1 < N 0 (N 0 < N 1 , respectively). As usual, there are two types of possible errors: the type-I error is quantified as α := P H 0 (N 1 < N 0 ), whereas the type-II error is 1 − π := P H 1 (N 1 > N 0 ), where π is typically referred to as the power of the test. In some cases the sampling is stopped only if the null hypothesis is rejected, for example in the context of change-point detection. If P H 1 (N < ∞) = 1, then the test is called a power-one sequential test. The significance level of the test is P H 0 (N < ∞). The most common approach for hypothesis testing is the likelihood-ratio-test (LRT): reject H 0 if L n x for some x > 0, where This test is proven to have certain optimal properties, such as having the highest power for any value of α. In our setting, however, serious complications arise, as we do not have closed forms expressions for the likelihoods P H i (V 1 , . . . , V n ), i = 0, 1. This is a consequence of the fact that the observations V i are in general not identically distributed (as we do not necessarily start at time 0 with a stationary workload), and, more importantly, have a rather elaborate dependence structure. Therefore, if we would like to apply the above LRT, then a first idea would be to rely on numerical or approximate techniques to evaluate the likelihoods P H i (V 1 , . . . , V n ). However, the remedies that we propose, and which will be described in detail later, are of a different nature: we transform the observations V 1 , . . . , V n into a vector of which we can compute the likelihood, effectively bringing us back into a classical LRT framework.
Remark 1. A (naïve) benchmark test on the mean rate generated by the driving Lévy process is the following. Perform a simple threshold test on the average workload: x for a suitably chosen x. In case one is interested in detecting stability, then this is a power-one test as in the unstable setting the average workload will grow unbounded. However, if both hypotheses correspond to stable queues, then there is a positive probability of a type-II error (but this can be controlled by choosing the test's parameters appropriately). Large-sample asymptotics for the test statistic are readily available, such as those in e.g. [7]: after centering and normalizing by √ n it converges to a normal random variable. This CLT could be used for an approximative likelihood-ratio test corresponding to i.i.d. normal random variables, but this would ignore the dependence between the observations. The aim of this paper is to develop new tests that combine the advantages of the likelihood-ratio test with tractable methods for quantifying its performance. We get back to the naïve test described above in the numerical analysis of Section 6. ♦

Approach I: quasi busy periods
As mentioned in the previous section, the approach we take is to transform the vector of workload observations (performed at Poisson instances) into an alternative vector for which LRT-type tests can be performed. In this section we focus on using the quasi busy period (QBP), being the number of observations until the workload hits zero again. Evidently, these QBPs constitute a sequence of i.i.d. random variables. In Section 3.1 we derive distributional properties of the QBP, which enable us to evaluate the likelihood. Section 3.2 then describes how the LRT should be set up.
3.1. Distribution of quasi busy period. In this section we concentrate on the random quantity, conditional on V 0 = 0, R thus records how many Poisson epochs it takes until the buffer is empty again. The objective of this subsection is to devise a procedure that facilitates the computation of the distribution of R, through the probabilities We do so by first computing the probabilities p k (ξ) := P(V k = 0), after which we express the r k (ξ) in terms of the p k (ξ). Let S k be an Erlang random variable with scale parameter n and shape parameter ξ, i.e., the sum of n i.i.d. exponentially distributed random variables with mean ξ −1 . We first point out how to compute, with f S k (·) denoting the density of S k , To this end, observe that In other words, we obtain the recursion, for k = 2, 3, . . . , The initialization of the recursion follows from p 1 (ξ) = ξ/ψ(ξ). The next term is This recursion can be further expanded, so as to obtain the following result. where In addition, ̺ ℓ (ξ) can be found recursively from Proof. We prove (6) inductively. The validity of the expression for p k (ξ) is obvious for k = 1. Now suppose the claim holds for some k ∈ N. Then, where the third step follows by recognizing a telescopic series. This proves the first claim. The second claim follows trivially by ℓ times differentiating ξ/ψ(ξ) (where we use the binomial expansion for higher derivatives of products of functions). Regarding the third claim, observe that, for any ℓ ∈ N, This immediately yields the stated recursion.
Remark 2. A sanity check of the above formula is that it should yield, as a consequence of the celebrated PASTA property ('Poisson arrivals see time averages'; see [2, Section VII.6.1]), that lim t→∞ P(V (t) = 0) = lim k→∞ P(V k = 0). This relation indeed holds, as follows from the leftmost expression equals lim k→∞ P(V k = 0), whereas the rightmost expression, as a direct consequence of the generalized Pollaczek-Khinchine formula equals lim t→∞ P(V (t) = 0) (see [6,Thm. 3.2]). The first equality in (7) is due to our explicit expression for p 1 (ξ). ♦ Remark 3. Remarkably, inspecting the proof of Proposition 2 reveals that this approach provides us with a general devise to translate a transform at an exponential epoch (with mean ξ −1 ) into its counterpart at an Erlang epoch (with parameters k and ξ). Indeed, defining α k (ξ | s) := E e −sĀ k , with (A t ) t∈R some stochastic process andĀ k := A S k , we obtain that This idea can e.g. be used to find the transform of the workload process in a Lévy-driven queue, say with spectrally-positive input, at an Erlang epoch. Once more assuming V 0 = 0, in this case we should take [6, Thm. 4.1] In this case the evaluation of the derivatives α (ℓ) 1 (ξ | s) is more involved than in the setting described in Proposition 2; it can be checked that α 2 (ξ | s) agrees with the formula in the last display of [6, p. 54] (choosing x = 0 there), that was derived in an entirely different manner. The setup considered in [15] is related; there the focus is on the case that the inter-event times T k have different means (where the case with equal means can in principle be dealt with applying a limiting argument). ♦ In our setting, we typically have ϕ(·) at our disposal (and all its derivatives), and in addition we can numerically evaluate ψ(ξ) (e.g. by bisection), but (as we lack a closed-form expression for ψ(ξ)) we cannot easily evaluate the derivatives of ψ(ξ). (An exception is the case that (X t ) t∈R is a compound Poisson process with exponential jumps; then the inverse ψ(·) allows an explicit expression.) To remedy this, the classical Faà di Bruno formula is helpful. Regarding the first derivative, we know that ϕ(ψ(ξ)) = ξ, so that by differentiation we obtain ϕ ′ (ψ(ξ))ψ ′ (ξ) = 1, and hence .
The above relation is a genuine recursion, due to the fact that m 1 + · · · + m k < k for all m ∈ M • k , entailing that when evaluating ψ (k) (ξ) all quantities appearing in the right-hand side of (8) are known. The next step is to compute the r k (ξ) from the p k (ξ), which can be done recursively. We partition the event of having a positive workload at observations 1 up k − 1 (i.e., V 1 > 0, . . . , V k−1 > 0) and a workload 0 at the k-th observation (i.e., V k = 0), as follows. The main idea is that the event under consideration can be written as the difference between (A) the event We thus obtain that r k (ξ) can be evaluated recursively through the following relation, providing us with the distribution of the QBP.
where the empty sum is defined as 0. Here the probabilities p j (ξ) directly follow from Proposition 2 and Lemma 3.

3.2.
LRT for quasi busy periods. In this subsection we point out how the results that were presented in Section 3.1 can be used to develop a test for the setting (5), relying on a sample of QBPs (R 1 , . . . , R n ). The log-likelihood of a sample (R 1 , . . . , R n ) is now, in self-evident notation, , which can be computed relying on the expressions derived in Section 3.1. We can now construct standard LRT. In a two-sided test there are two thresholds, say x 0 and x 1 such that The decision rule is based on N := inf{n 1 : ℓ n / ∈ [x 0 , x 1 ]}: reject the null hypothesis if ℓ N > x 1 and accept the null hypothesis if ℓ N < x 0 . If x 0 = −∞, then this is a power-one test and the type-I error probability is given by We proceed by reflecting on the pros and cons of this test. The main disadvantage of this method is that a substantial amount of information is lost when transforming the workload observations into QBPs: in fact it is only used whether a observation is zero or positive (i.e., its precise value is ignored). There are, however, two important attractive properties: (1) The test allows distinguishing between any pair of two Lévy subordinators with ϕ 0 (·) = ϕ 1 (·) and an arbitrary sampling rate ξ. This is because the distribution function r k depends on all derivatives of the inverse function ψ(·). As we will see later, the test that is presented in Section 5 does not have this property. (2) As mentioned, the sample (R 1 , . . . , R n ) is i.i.d. As a consequence, standard methods for approximating the error α x are readily available; see, e.g., the textbooks [13,17].
From a practical point of view, as the number of operations required to evaluate r is of the order i!, there may be computational issues. This means that, particularly for larger k, the computations may become time consuming. This problem can be overcome by a truncation: for some appropriately chosen K, all QBP durations of at least K are lumped together, and have probability We conclude this section by providing a quantification of the type-1 error. Consider the random walk ℓ n = n i=0 P i . For any given threshold x we define the first passage time through N x := inf{n : ℓ n x}. We are interested in the hitting probability α(x) := P H 0 (N x < ∞), and the expected first passage time ; an explicit expression for C can be found in e.g. [9]. In addition, by standard arguments,

Rate of convergence to steady state
In the previous section we have set up a test based on quasi busy periods. In Section 5 we propose an alternative approach, which we call a conditional likelihood ratio test (CLRT). As a preparation to the performance analysis of this CLRT, in the present section we provide results describing the speed of convergence of a Lévy-driven storage system. More specifically, they (i) facilitate the computation of the asymptotic variance of the likelihood ratio, and (ii) enable the construction of a functional limit theorem for the likelihood-ratio that can then be used in order to approximate the test's error probability by a Brownian motion hitting probability. Importantly, however, these speed-of-convergence results are, to the best of our knowledge, new, and relevant in their own right.
The convergence rate conditions we deal with here are weaker than those required for geometric ergodicity, in the sense of [10]. Informally speaking, we do not require that the observed workload process converges exponentially fast to the stationary distribution, but rather that it does so at a rate such that the sum of absolute deviations of certain functions of the workload (relative to their expectations according to the stationary distribution, that is) converges.
In what follows we make use of both the continuous time workload V (t) = X(t) + L(t), where L(t) = − inf 0 s t X(s), and the corresponding discretely observed workload process V n = V (S n ), where S n is Erlang distributed with scale parameter n and shape parameter ξ > 0. Recall that V denotes the steady-state workload, where the stability condition ϕ ′ (0) > 0 is assumed throughout. The first lemma states a PASTA-type equivalence result for comparing the accumulated deviation from the stationary expectation in continuous time to the corresponding deviation at Poisson epochs.
Proof. We can write the left-hand side of (9) as which equals, with T ξ denoting an exponentially distributed random variable with mean ξ −1 , Now recall that a geometrically distributed (with success parameter q) number of exponentially distributed random variables (with mean ξ −1 ) is exponentially distributed (with mean (qξ) −1 ). This means that we can rewrite the expression in the previous display as which equals the right-hand side of (9).
The following theorem states a number of results that describe the workload (observed at Poisson epochs) convergence to stationarity, conditionally on the initial workload V (0).

Theorem 7.
If P(V (0) < ∞) = 1, then the following equations hold almost surely: (i) assuming that X(1) has a finite third moment (i.e., |ϕ k (0)| < ∞), and defining for any v 0 (ii) defining for any v, α 0 (iii) defining for any v, α 0 we have, for any α 0, Proof. We first evaluate the series (13) for a given V (0) = v by considering the continuous-time analog, where (3) has been used. The limit can computed by applying L'Hôpital's rule twice and yields that d v (α equals k 2 (v, α), as defined above. By Lemma 5 we thus obtain As P(V (0) < ∞) = 1 we can condition on V (0) to obtain (13): To establish (12) and (14) we follow similar arguments, also using Lemma 6. For any V (0) = v, let Observe that Evaluating the above derivatives and taking the conditional expectation with respect to V (0) immediately yields (12) and (14).
The following technical lemma will play a crucial role in Section 5.
By the definition of ϕ(α) we have that lim α→∞ ϕ(α)/α = 1, and thus, Therefore, for any v 0. We conclude that, for any v 0, |k 2 (v, α)| is a bounded and continuous function, with respect to α ∈ (0, ∞], and it therefore admits a maximal value at some α ⋆ < ∞ or approaches an upper bound as α → ∞. Finally, E[k ⋆ 2 (V )] < ∞, as a consequence of the fact that (2) in combination with the finite second moment assumption implies that EV < ∞.
For the case of the initial workload V (0) being 0 we refine the results of Theorem 7 to absolute convergence of the series. This result will be useful in establishing the asymptotic variance of the likelihood ratio in Section 5.2.
Lemma 9. If X(1) has a finite second moment (i.e, ϕ (2) (0) < ∞), then for any α > 0, (15) If X(1) has a finite third moment (i.e, |ϕ and, for any α > 0, Recall that V is a random variable that is distributed as the stationary workload (whose transform is given by Eqn. (1)), so that V ⋆ (t) is the workload at time t starting with a stationary workload at time 0 (implying that V ⋆ (t) = d V for all t 0). Note that V (t) is its counterpart, but starting with an empty workload at time 0. Because of V (0) = 0, we have that V (t) V ⋆ (t) for all t 0 for every sample path of X(·), and in particular V (S n ) V ⋆ (S n ) for any Erlang distributed sampling time S n . Hence, Similarly we also have that e −αV (t) e −αV ⋆ (t) for any α 0 and all t 0 for every sample path of X(·), hence and applying (13) for V (0) = 0 yields (18). If g(·) is a Lipschitz continuous function, there exists a constant K ∈ (0, ∞) such that, for all t 0, By Jensen's inequality, the Lipschitz assumption, and Eqn. (9) in Lemma 5 , Therefore, the inequality in (17) follows by verifying that x → xe −αx is a function with Lipschitz constant K = 1.

Remark 4. In [16, Thm. 2] it was shown that the M/G/1 virtual waiting time is geometrically ergodic if and only if 1 − G(x)
ce −µx for some c, µ > 0. In Theorem 7 above we assumed weaker conditions: finiteness of specific moments rather than light-tailed jumps. Note, however, that geometric ergodicity is stronger than the property that we found, because we do not specify a rate of convergence, but just the finiteness of the integral. For the Brownian approximation of the error probability in the hypothesis testing application presented in Section 5.3 we will see that a geometric convergence rate is a sufficient, but not necessary, condition. ♦ Corollary 10. If X(1) has a finite second moment (i.e, ϕ (2) (0) < ∞), then for any α > 0, Proof. Recall that ϕ(α) αϕ ′ (0). Observing that the claim follows.

Approach II: conditional likelihood ratio test
As in our first approach, in our second approach we use a transform of the observation that allows us to evaluate the likelihood. The underlying idea is that we construct an LRT based on Bernoulli variables of observing either a zero workload or a positive workload, conditional on the value of the previous workload observation. To this end, we introduce Y i := 1(V i = 0) to denote a sequence of idle-period indicators, for i = 1, . . . , n. Conditional on V := (V 0 , . . . , V n ), Y = (Y 1 , . . . , Y n ) is distributed as a sequence of independent, but not identically distributed, Bernoulli random variables. More specifically, the corresponding likelihood reads, using (4), The likelihood function (19) depends on the input distribution only through the constant θ := ψ(ξ) ∈ (ξ, ∞). We aim at developing a test for the simple hypothesis testing problem From now on we assume that the hypothesis H 0 and H 1 correspond to the Laplace exponents ϕ 0 (·) and ϕ 1 (·) that can be distinguished by the parameters θ 0 = ψ 0 (ξ) and θ 1 = ψ 1 (ξ), given the sampling rate ξ. For a sample of workload observations, V = (V 0 , . . . , V n ), the conditional likelihood ratio test (CLRT) is based on the statistic can be evaluated using (4).
Before proceeding to the analysis of this test, a few remarks are in place. As can be seen from the above expressions, the test can only distinguish between distributions if ψ 0 (ξ) = ψ 1 (ξ) for the current sampling rate ξ. This is a minor problem though: one could always choose a sampling rate ξ for which ψ 0 (ξ) and ψ 1 (ξ) do not coincide. In the three bullets below we consider this issue in greater detail.
• For an M/M/1 queue with known service rate µ, we can test for hypotheses on the arrival rates: λ 0 = λ 1 . There is the advantage that the inverse of the exponent function is known in closed-form: This function is clearly monotone in λ (for any given ξ). A hypothesis on λ is therefore equivalent to an hypothesis on ψ(ξ). Conclude that for any ξ the values of ψ 0 (ξ) and ψ 1 (ξ) do not coincide (as long as λ 0 = λ 1 ). Note that the same argument holds if we replace the roles of arrival and service rates, i.e., fix λ and test for µ 0 = µ 1 . The monotonicity of ψ(ξ) with respect to µ can be verified by straightforward algebra. • For an M/M/1 queue with unknown arrival and service rates λ and µ, there may be values of ξ for which we cannot always test hypotheses on the input intensity ρ = λ/µ. For example, λ 1 = 0.8 and µ 1 = 1 yield ρ 1 = 0.8, and λ 0 = 0.235 and µ 0 = 2 yield ρ 0 = 0.51 < ρ 1 . However, ψ 1 (2) = ψ 2 (2) = 2.576. We conclude that in this case even if the traffic intensity is very different, then the conditional likelihood cannot distinguish between the two instances by testing for the parameter ψ(ξ). The obvious remedy is, as mentioned, to pick another sampling rate ξ. • Similar considerations play a role for more general Lévy subordinators. Consider for instance a storage system in which the driving Lévy process is a Gamma process with parameters (β, γ). A Gamma process has increments that are, per time unit, Γ(β, γ) distributed. Similar to the M/M/1 case, for any sampling rate ξ a test for either β or γ can be constructed, but when setting up a test for the traffic intensity ρ = β/γ one has to check whether the chosen ξ is such that ψ 0 (ξ) = ψ 1 (ξ); if not, then ξ has to be adapted.

5.1.
Power-one sequential test. In the sequential test based on the conditional likelihood ratio, the null hypothesis is rejected when the test statistic attains a high value: for some threshold x, reject H 0 if L n x. Note that this test either rejects the null hypothesis or does not terminate. The latter can only occur if H 0 holds, as L n will reach any threshold x < ∞ with probability one under H 1 . Formally, with N x = inf{n 1 : L n x}, we have P H 1 (N x < ∞) = 1, making it a power-one sequential test. On the other hand, there is a positive probability that the threshold x will be reached under H 0 , i.e., a false rejection of the null hypothesis. We denote the probability of this type-I error by our goal is to compute or approximate this quantity. Applying (19) yields In what follows we work with the log-likelihood ratio (LLR), defined as ℓ n := log L n . The LLR ℓ n is more convenient than L n because of its additive structure. Concretely, we can write ℓ n = n i=1 Z i for random variables Z i . It should be kept in mind, however, that ℓ n is not a classical random walk, as the increments are neither identically distributed (unless one starts off in stationarity) nor independent. From now on we consider the test with the stopping rule N x = inf{n 1 : ℓ n x}. There are two key performance measures for this test. In the first place one commonly considers the significance level The second key metric is the expected number of samples until rejection of the null hypothesis From now we use the short notations P k and E k for probabilities and expectations under the null hypothesis (k = 0) and alternative hypothesis (k = 1).

5.2.
Large sample asymptotics. In the sequel we assume that the workload process starts in stationarity. This makes the increments Z i identically distributed, but they are obviously not independent. As a consequence, standard methods from sequential analysis are not directly applicable. Nevertheless, as we will show below, we can provide approximations for the test's performance measures. In this subsection we establish the asymptotic behavior of the log-likelihood ratio ℓ n as the sample size n grows large. Specifically, we establish a strong law of large numbers (SLLN) for the sample mean and a functional central limit theorem (FCLT) for the corresponding centered empirical process, assuming that the workload process starts in stationarity. These results will be used in the sequel to approximate the performance of the CLRT: we approximate α(x) and τ x by the hitting probability and expected hitting time, respectively, of a Brownian motion with an appropriately chosen drift and variance. Taking the logarithm of (20) we have ℓ n = n i=1 Z i , where If the workload process is stable under H k for k ∈ {0, 1}, i.e., E k X(1) = −ϕ ′ k (0) < 0, then by the PASTA property the stationary distribution V is also the limit of V n as n → ∞ with respect to P k , so that V n = d V ). We denote a stationary increment of the LLR by Z; because we start off in stationarity, Z n = d Z.
In the following lemma we establish a SLLN for the stationary first and second moment of the sequence (Z n ) n∈N . These will be used later to establish a FCLT for the LLR process. A complication lies in the aforementioned fact that the Z n are not independent, so that standard LLN and FCLT techniques cannot be applied directly. To remedy this, we resort in the rest of this subsection to utilizing the special structure of the workload sampled at Poisson times, in combination with a martingale FCLT for weakly dependent stationary random variables. Define so that, by (21), . and Proof. We start by proving claim (i). Applying PASTA, we have that 1 n ℓ n = 1 n n i=1 Z i → as E k Z. The stationary sample average can be computed as follows. Recalling that V (0) = d V , the mean of the first increment Z 1 (and hence also the mean of all other increments) equals where in the last equality (4) has been used. Claim (ii) follows in the same manner.
We next turn our attention to the asymptotic distribution of √ n (ℓ n /n − m k ) = (ℓ n − nm k ) / √ n as n → ∞, for k = 0, 1. Let z i := E k Z i and M n : As V 0 (and Z 1 ) is stationary, we have z i = m k for all i, so that we can focus on deriving the limiting distribution of M n / √ n as n → ∞. Let BM(d, σ 2 ) be a Brownian motion with drift d and variance coefficient σ 2 , i.e., at time t having a normal distribution with mean dt and variance σ 2 t. The main result of this subsection is the following FCLT.
The proof of this theorem relies on the methodology presented in [5,. In particular, we decompose M n into a martingale difference process and an additional term that vanishes when scaled by 1/ √ n, and then apply the FCLT for stationary ergodic martingale difference processes (see [5,Thm. 18.3]). The decomposition is similar to the one used in the proof of [5,Thm. 19.1], but we utilize the specific structure of our problem and do not impose the stronger assumptions that are required there. Informally, this essentially boils down to showing that the dependence between the increments in the LLR diminishes sufficiently fast. To this end, we apply the results of Section 4 that describe the convergence rate of the LLR to stationarity. This will then be used to verify the FCLT conditions of the martingale approximation, and in particular to show that σ 2 k < ∞. We proceed by sketching the structure of the remainder of this subsection, geared towards proving Theorem 12.
• Lemma 13, which builds on Theorem 7 and Lemma 8, establishes that the LLR converges fast enough to satisfy a sufficient condition for the FCLT of Eqn. (25). • Then, in Proposition 14, we rely on Lemma 9 and Corollary 10 to show that σ 2 k := lim n→∞ n −1 Var k M n is finite, and moreover the absolute convergence of the series of covariance terms (Cov k (Z 1 , Z i+1 )) i=0,1,... . • Then we are in a position to prove Theorem 12. We use the martingale decomposition and apply Lemma 13 to verify that M n can be written as a sum of a martingale difference processM n with stationary increments that have a finite second moment (thus satisfying the conditions of [5,Thm. 18.3]) and an error term R n that is almost surely finite. Then Proposition 14 is used to show that the asymptotic variance ofM n / √ n equals σ 2 k . The proofs of Lemma 13 and Proposition 14 are quite lengthy and are therefore relegated to the appendix. Lemma 13. Let k ∈ {0, 1}. If ϕ ′ k (0) > 0, and the initial workload is stationary (i.e., V 0 = d V ), then the following claims hold: (a) With respect to P k we have that almost surely, where The following proposition establishes that the variance of the LLR is finite, and in addition that the series of covariance terms (Cov k (Z 1 , Z i+1 )) i=0,1,... converges absolutely. The second part of this statement is important not just for establishing the FCLT approximation but also from a computational perspective. It implies that the covariance series, and subsequently σ 2 k , can be evaluated efficiently by means of truncation as the remainder of the series vanishes in absolute terms.
k (0)| < ∞, and the initial workload is stationary (i.e., V 0 = d V ), then as n → ∞, where c ki := Cov k (Z 1 , Z i+1 ) is given by Moreover, ∞ i=1 |c ki | < ∞ and therefore σ 2 k < ∞. Proof of Theorem 12. In Lemma 14 it is established that this is, however, a necessary but not sufficient condition for the FCLT (25) to hold. We will show that M n can be decomposed into two terms; one that is almost surely finite and the other has the same asymptotic distribution as a martingale difference sum that satisfies the conditions of the martingale FCLT [5,Thm. 18.3].
The decomposition relies on iterating expectations, as follows. Define F n := σ{V 0 , . . . , V n } and consider Then observe that the first term, i.e., D i,1 : Similarly, the second term can be decomposed into where the first term in the right-hand side of the previous display is now a F i−2 -martingale difference. Continuing along these lines, we readily obtain, for j ∈ {1, . . . , i}, that is a F i−j -martingale difference. By performing i iterations, we thus find that which implies that Let us consider the first term in the right-hand side of (31). By changing the order of summation and some relabelling, we obtain Note that for now we have assumed that the infinite series are finite almost surely, and we will later verify that this indeed the case. The stationarity of V 0 implies that for any i > j, Therefore, for every j ∈ N, with (D j ) j∈N being a sequence of stationary random variables. Now Lemma 13(a) implies that, almost surely, Under the assumption of |ϕ By applying similar arguments E k [z(V 0 )z(V 1 )] < ∞ as well (in particular, note that one can condition on V 0 and use (3) and (4) to obtain a decomposition of E k [z(V 0 )z(V 1 )] as linear and quadratic terms of V 0 ). Therefore, E k [D 2 j ] < ∞. Applying the telescopic structure in (30) yields Upon combining the above, we find from (31) that we can write Observe thatM n is a martingale difference sum. Also, R n is almost surely finite due to Lemma 13(b), so that R n / √ n → 0 almost surely as n → ∞. Combining the above, we conclude notice that the right-hand side involves a series of the stationary martingale differences (D j ) j∈N . As E k [D 2 j ] < ∞, the conditions of the martingale FCLT [5,Thm. 18.3] are satisfied and we conclude that (25) holds. Furthermore, E k [D 2 j ] = lim n→∞ n −1 Var kMn . Using the Cauchy-Schwarz inequality together with n −1 Var k M n → σ 2 k (see Lemma 14) and n −1 Var k R n → 0, we conclude that as n → ∞, This concludes the proof.

Brownian approximation.
A standard way to approximate α(x) is by relying on a Brownian approximation, in our case facilitated by Theorem 12. Applying this theorem, with (B(t)) t 0 a standard Brownian motion, realizing that m 0 < 0, where the last expression is a standard result for the maximum of a Brownian motion with negative drift. Similarly, we can approximate the expected rejection time if the null hypothesis is wrong by

5.4.
Change of measure and approximation. Apart from a Brownian approximation, another standard technique to approximately evaluate α(x) is by using exponential-change-ofmeasure techniques. Again, the key difficulty is the intricate dependence between increments of the log-likelihood process. As before, ℓ n = n i=1 Z i , where the increments Z i are given by (21), and N x = inf{n 1 : n i=1 Z i x}. Let the cumulant generating function be given by Then M n (β) := exp(βℓ n − nκ n (β)) is a martingale with respect to F n := σ{V 0 , . . . , V n } such that EM n (β) = 1 for all n ∈ N. Let γ n denote the sequence of solutions to κ n (γ n ) = 0, and consider the change of measure Lemma 15. If ϕ ′ 0 (0) > 0 and the initial workload is stationary (i.e., V 0 = d V ), then Proof. Again, applying PASTA and the SLLN, in combination with V 0 = d V , yields the stated. First observe that which equals which, by (1), equals (33).
We now point out how α(x) can be approximated in the regime that the sampling rate is very slow, i.e., ξ → 0. Then the Z i are (almost) independent and identically distributed (as the random variable Z). Hence, M n (β) roughly equals exp(βℓ n − nκ 1 (β)), where κ 1 (β) is given by (33). The Lundberg coefficient γ > 0 is the solution to κ 1 (γ) = 0, which is unique because κ ′ 1 (0) = E 0 Z 1 = m 0 < 0. Consider the change of measure Then applying [2, Thm. III.5.1] we have that, with C defined as lim x→∞ E γ e −γ(Nx−x) , a more explicit characterization of C is given in [9]. We now investigate to what extent the correlation between the Z i can be taken into account. Ideally, we would like to relate the stationary coefficient γ to the sequence of coefficients γ n that satisfies κ n (γ n ) = 0. Analyzing the sequence γ n is challenging because it involves the joint distribution of (Z 1 , . . . , Z n ). The following procedure can be followed. Let h n (β) := E 0 e β n i=1 Z i , and observe that κ n (β) = 0 is equivalent to h n (β) = 1. We have that h 2 (β) = E 0 e βZ 1 +βZ 2 = Cov 0 e βZ 1 , e βZ 2 + E 0 e βZ 1 E 0 e βZ 2 , and as Z 1 and Z 2 both have the same marginal (and stationary) distribution we can write Eqn. (34) can be generalized to The conclusion is that the sequence γ n can be numerically approximated by simulating the solutions to (35). Such a simulated sequence of solutions is illustrated in Figure 1. In Section 6 we explore the performance of the error approximations based on γ 1 and γ 20 .

Simulation analysis
In this section we present numerical results on tests corresponding to an M/M/1 queue, distinguishing between two values of the arrival rate λ (λ 0 = 0.6 and λ 1 = 0.8), for a given value of the service rate (µ = 10). 'QPBT' refers to the test based on the quasi busy periods (Approach I), and 'CLRT' to the conditional likelihood ratio test (Approach II). From our experiments, the following observations can be made: • The QBPT and CLRT tests presented in Sections 3 and 5 agree in over 80% of the sample-path realizations. That is to say, in the vast majority of all cases, if one test rejects, then so does the other. The thresholds for both tests were chosen so that they have approximately the same significance level α(x). • Both tests greatly outperform the naïve mean test n −1 n i=1 V i x that was discussed in Remark 1. This naïve test has a very lower power for the same significance levels.
• All approximation techniques presented in this paper perform reasonably well. See for example Figure 2 that compares the type-I error probability of the CLRT, based on simulation with three approximations: the Brownian approximation of Section 5.3, the change-of-measure approximation of Section 5.4 using γ 1 , and the change-of-measure approximation of Section 5.4 using γ 20 . • In Figure 3 the type-I error probabilities are plotted for increasing sampling rates ξ. The thresholds are chosen so that the approximated error probability equals 5%. For the QBPT, x = −γ −1 log(0.05), where γ corresponds to the change-of-measure approximation presented in Section 3. For the CLRT, x = −γ 20 −1 log(0.05), where γ 20 corresponds to the change-of-measure approximation presented in Section 5.4. We observe that the approximation improves as the sampling rate increases.
• In Figure 4 the type-II error probabilities β are plotted for increasing sampling rates ξ, using the same parameters as in Figure 3. The tests were truncated at n = 10 3 . For a low sampling rate, the truncation has no effect; we observe that β ≈ 1, as is the case for n → ∞. However, for lower sampling rates, many more observations are required, and in addition the test does not have power one. This is particularly true for the QBPT, were the power drops to almost 70%, while for the CLRT it remains above 90%. • Figure 5 plots the corresponding expected number of samples until rejection (conditional on rejection before truncation). The QBPT test generally rejects faster, but as we saw in Figure 4 this comes with a lower power, meaning that if the test was not truncated these expectations will be much higher.
x α(x) Figure 2. Probability α(CLRT) of type-I errors as a function of the threshold x for the CLRT, together with three approximations. The (non-approximative) probabilities were computed using simulation with n = 10 3 samples. The sampling rate is ξ = 3. Figure 5. Expected number of samples until rejection of the null-hypothesis (under H 1 ), using both approaches (τ (CLRT) and τ (QBPT)). The (nonapproximative) probabilities were computed using simulation with n = 10 3 samples.

Discussion and concluding remarks
This paper has focused on hypothesis testing for Lévy-driven storage systems. Tests are developed that are capable of distinguishing a Laplace exponent ϕ 0 (·) from an alternative Laplace exponent ϕ 1 (·). In the setup considered the driving Lévy process is a subordinator minus a deterministic drift. In both approaches proposed, we exploit the property that for such Lévy processes the workload level is zero with positive probability. An open challenge is to set up tests for more general Lévy input. We remark that reflected Brownian motion is not covered by the current framework, but it can be dealt with due to the fact that for this specific model the likelihood can be evaluated in closed form.
In both approaches information is 'lost': in the first approach it is only used whether a workload sample is positive or zero (thus leading to i.i.d. quasi busy periods), in the second approach it is only used whether a workload sample is 0 conditional on the previous observation. One thus wonders whether one could develop a test in which 'more information is used'. Due to the lack of an explicit expression for the joint density of V 1 , . . . , V n , inevitably, either some information provided by the sample will be lost, or the likelihood of the sample is to be evaluated in an approximate manner. Poisson sampling has evident computational advantages, as pointed out in our paper, but in some practical settings one may prefer equidistant sampling. The techniques presented in [15] could potentially be used to approximate the corresponding density, but computational challenges are anticipated. Alternatively, the likelihood can be evaluated by numerical inversion of (3); this concerns an inversion to obtain a density from a transform as well as an inversion to translate the exponentially distributed time into its deterministic counterpart. In a related setting, such (double) inversion techniques have been explored in [1].
Using the upper bounds established in part (b) we have that for some positive and finite constants A 1 , A 2 , A 3 . Hence, by the definitions of k 2 (v, α) and k 3 (v, α) in Theorem 7 and k ⋆ 2 (V 0 ) in Lemma 8, z(V ) 2 is upper bounded by a linear combination of terms of the form V , V 2 , e −αV , V e −αV and V 2 e −αV (with α ∈ (0, 2 max{θ 0 , θ 1 , α ⋆ }]). Therefore, appealing to (2), if ϕ ′ (0) > 0 and |ϕ Proof of Proposition 14. If Z 0 is stationary then Z i = d Z 0 for any i ∈ N. The first claim then follows from the standard identity the decomposition of c ki in (29) is obtained by plugging in (21), and using the assumption that V 0 is stationary. It is left to prove that ∞ i=1 |c ki | < ∞, which we show by using that V n converges fast enough to its stationary distribution, relying on Lemma 9. Below we denote the density of the workload level at sample n by f n (·); likewise, the density of the conditional workload at sample n given an initial workload V 0 = v 0 is denoted by by f n (· | v 0 ). We treat the four individual terms in the right-hand side of (29) separately. For ease, we leave out the subscript k throughout the proof.
-First term. Recall that Y n = 1(V n = 0). Thus, As V i+1 and V 0 are independent conditional on V 1 for any i ∈ N, we have that As a consequence, the first term of (29) can be written as ] . Next, applying (4) yields, for i ∈ N, Applying (17) in Lemma 9 we conclude that if |ϕ (3) k (0)| < ∞, then ∞ i=1 |a i | < ∞. -Second term. Applying similar arguments, the second term term of (29) is given by where Then it is a matter of straightforwardly combining ideas from the second and third term to prove that ∞ i=1 |f i | < ∞. Now that we have bounds on all four terms, by applying the triangle inequality once more we conclude that there exists a constant 0 < κ < ∞ such that and we thus conclude that σ 2 k < ∞.