Estimating Customer Impatience in a Service System With Unobserved Balking

1. Introduction

In many service systems, potential customers are provided with real-time information on the delay they will experience (Ibrahim 2018, Dong et al. 2019). This information, typically in terms of anticipated delay or current workload, is used by the potential customers to decide to either join the system or balk. In general, one could view balking as an implicit form of admission control; the load on the service system does not explode as a consequence of the fact that during busy periods some of the potential customers decide to leave, thus mitigating the overall congestion. This means that the system is in effect self-regulating even if the potential service demand is unknown to the administrator. To make sound decisions in designing the service system, however, it is crucial to have knowledge of the volume of the potential demand and the mechanism based on which potential customers decide to join or to leave. With this knowledge, the service provider can, for example, determine the optimal admission price or decide whether it is economically viable to increase the service capacity.

This paper is motivated by the fact that it is challenging to extract this knowledge from data when balking customers cannot be observed by the administrator. Such complication arises naturally in many service systems. The most obvious case is a physical service facility, where a visible queue (e.g., shops, parking lots, museums), sometimes even outside of the facility, will deter some of the potential customers. This has become even more prevalent recently, with social distancing measures limiting the number of people that can be inside a public space at a given time. This form of data censorship occurs also in systems that provide real-time delay data via electronic means such as dedicated apps for the service. Some examples are expected travel times in transportation services (e.g., Google Maps and Uber) and expected delivery times in food delivery services (e.g., UberEats in the United States or Deliveroo in Europe). Note that expected delay is often listed for the possible services (e.g., estimated arrival time of a taxi or delivery times for a list of Italian restaurants) even if a user does not indicate the specific service he or she is interested in, and therefore, the system cannot know whether a user balked from the service because of the delay information.

In this paper, we model such service systems with unobservable balking as an M/G/s + H queue, where the last symbol stands for the cumulative distribution function of the customer patience. More specifically, we consider a first-come, first-served (FCFS) queue with s ($s⩾1$) servers and potential customers arriving according to a Poisson process with rate $\lambda >0$. Potential customers have independent and identically distributed (i.i.d.) patience levels that are distributed as the nonnegative random variable $Y\sim H$ and bring i.i.d. service requirements that are distributed as the nonnegative random variable $B\sim G$. Customer i joins the system if Y_i is greater or equal to the virtual waiting time at the moment of his or her arrival and otherwise balks without being observed. Note that although customers are homogeneous in a statistical sense, they have individual (thus heterogeneous) patience levels, so this model can clearly represent the types of customers that will effectively join the system.

We set up an estimation problem by parametrizing the patience distribution as $H_{\theta }$ ($\theta \in \mathrm{\Theta }$), where Θ denotes a compact set satisfying some regularity conditions. The goal of this work is to estimate the pair $(\lambda ,\theta )$ based on observation of the effective queueing process, which is constructed from records of arrival and departure times not including balking customers. As some of the customers balk, the corresponding effective interarrival times are not the usual i.i.d. exponentially distributed random variables. The effective arrival process is in fact a nonhomogeneous Poisson process whose rate depends on the virtual waiting time process. In particular, the effective arrival rate to the system when the virtual waiting time $v⩾0$ is $\lambda (1-H(v))$. Note that the marginal effective arrival process is not a nonhomogeneous Poisson process with respect to time. Therefore, the dependence structure of the observations of interarrival times is directly determined by that of the virtual waiting time process. The workload-dependent arrival dynamics are utilized to derive a maximum-likelihood estimator (MLE) based on the effective interarrival observations. Specifically, the density of an interarrival time A conditional on initial virtual waiting time $v⩾0$ is given by

The Markovian dynamics of the queue further imply that the likelihood of the interarrival times conditional on the waiting times and job sizes is given by a product of densities with the above form for any sample of observations.

We further study the asymptotic properties of the proposed estimator. Importantly, the asymptotic performance of the estimator is not given by standard results due to the intricate underlying dependence structure. Indeed, so as to prove consistency and asymptotic normality, a subtle reasoning is required to make sure that the situation at hand fulfills specific regularity conditions. It is further shown that for homogeneous (deterministic) patience levels, that is, $\text{P}(Y=\theta )=1$ for some $\theta >0$, the asymptotic distribution of the errors is not normal but rather exponential (with rate n and not $\sqrt{n}$). The main contribution of this paper is the development of procedures for consistent estimation of the total arrival rate (that is, corresponding to customers joining the system and balking customers) and the patience-level distribution, observing the queue-length process only. Hence, somewhat counterintuitively, we can estimate load that has never been observed.

The main purpose of this paper is to lay foundations for statistical inference for congested service systems with unobserved balking. Although the queueing model considered in this paper is of course not a perfect replica of any specific real system, it does provide a methodological statistical approach for service systems where unobserved balking due to congestion is a key feature. Furthermore, the framework presented here can be extended to capture additional features such as noisy waiting time information, as will be discussed in the concluding remarks.

It is also important to stress that the framework we develop in this paper is essentially different from the extensively studied problem of estimating patience parameters from observed abandonment data (see Brown et al. 2005, Mandelbaum and Zeltyn 2013). If balking customers are observed, then the data are comprised of waiting times and corresponding indicator variables stating whether a customer joined or balked. Specifically, when it is observed that a customer i joined (resp. balked) given the waiting time W_i, the observation is of the form $Y⩾W_i$ (resp. $Y<W_i$). These data can be directly applied to construct an estimator for θ, for example, by applying an MLE method. If abandonment times are directly observed, then even more information is available, namely the realization of Y_i for every abandoning customer, and $Y⩾W_i$ for customers who obtained service. Note that in the above scenarios estimation of the total arrival rate is straightforward because all arrivals are observed. However, when balking is unobserved, the only information available is that of interarrival and waiting times, and therefore, indirect estimation methods are called for. Of course, there may be applications where some of the balking or reneging customers are observed, whereas others are not. For such systems, it may be reasonable to adopt a hybrid approach that uses classical patience estimation methods together with those presented in this paper. For the sake of brevity, we assume throughout that balking happens immediately and is not observed by the system operator so that no data regarding the balking customers is available for estimation purposes.

1.1. Contributions

We proceed by detailing our paper’s contributions. The focus is on developing estimators with provable performance guarantees. Notably, setting up our estimator is straightforward in the sense that it relies on a closed-form expression for the likelihood and does not require any queueing-theoretic analysis. As a result, even for models for which stationary performance analysis is involved (or intractable), parameter inference can be done in a relatively direct manner. Our contributions relate to (i) our estimator and its asymptotic properties, (ii) results and experiments for specific patience-level distributions, and (iii) extensions and ramifications.

• Framing the system described above as a queueing model with impatient customers, we can evaluate the corresponding likelihood function (in terms of the observed quantities). This means that, in a parametric context, we can estimate the unknown parameters — that is, λ and the parameter(s) θ corresponding to the patience level Y — relying on a maximum-likelihood procedure. The estimation method has several attractive features, the most prominent one being that it does not require simultaneous estimation of the service requirement distribution or even making any distributional assumption about it. In addition, our methodology uses the fundamental queueing dynamics and does not, as is common in the literature, rely on any fluid and/or diffusion approximations.

• In the case of a continuous patience-level distribution, we prove, under appropriate regularity conditions, strong consistency and asymptotic normality of the resulting estimation errors (scaled by $\sqrt{n}$). The proof of the asymptotic normality relies on an application of a suitable version of the martingale central limit theorem.

• In the case that the patience level Y is constant (i.e., $Y=\theta$ for some θ), more refined results can be obtained. Most notably, we show that the estimation error (scaled by n) converges to an exponential distribution with a parameter that depends on the stationary workload density at θ and the loss probability of customers.

• A number of other special cases are dealt with, covering exponentially and (generalized) hyperexponentially distributed patience levels; we present closed-form expressions for the asymptotic variance of the estimation error. The class of generalized hyperexponentially distributions is highly relevant for various reasons, most importantly because these can be used to approximate any positive continuous distribution arbitrarily well (Botta and Harris 1986). In addition, in the call center literature (Roubos and Jouini 2013), it has been observed that such distributions yield a good fit for observed customer patience in terms of the time until abandonment. We provide several numerical examples that illustrate the robustness of our methodology. Indeed, our experiments show that even if the patience-level distribution is misspecified (i.e., is not generalized as hyperexponential itself), the distribution function of the estimated generalized hyperexponential distribution is still remarkably close to the true one.

• We also discuss several ramifications and model extensions that can be handled using our framework. These include extending the framework to a system with noisy delay information, estimating the proportion of impatient customers in the population, and estimating specific utility parameters such as service value and waiting time sensitivity.

1.2. Related literature

There is a substantial body of literature on queues with impatient customers. Without attempting to provide an exhaustive overview, we include a brief account of the various research lines. Naor’s influential paper (Naor 1969) can be seen as a pioneering contribution to the field of “behavioral queues.” It presented a powerful stylized model for a queueing resource in which potential customers decide to join or to balk based on delay information. Since then, performance and economic analysis of queueing systems with balking customers has been studied extensively. In particular, the effect of providing workload information to strategic customers was investigated in, for example, Guo and Zipkin (2007, 2008). Many other model variations are reviewed in Hassin and Haviv (2003) and Hassin (2016). More background on applied probability and queueing systems in a general context is provided in, for example, Assmussen (2003) and Dębicki and Mandjes (2015).

For various specific types of queueing systems with impatient customers, explicit performance results are known; see, for example, Boxma et al. (2010, 2011), Liu and Kulkarni (2006, 2008), and Mandelbaum and Shimkin (2000, 2004). Importantly, however, the estimation techniques presented in our paper do not require knowledge of, for example, the queue’s stationary distribution. It is also worth stressing that despite the fact that queues with impatient customers are well studied, to the best of our knowledge, hardly any workload-based estimation techniques have been developed. A notable exception is Gorst-Rasmussen and Hansen (2009), featuring an asymptotic analysis of statistical inference of waiting times truncated by abandonments; our setting, however, is crucially different from the one in Gorst-Rasmussen and Hansen (2009), because in the latter, abandoning customers and their sojourn times are observed.

There are various papers on estimation problems in the setting where all customers join but leave when getting impatient; examples are Courcoubetis et al. (1995), Akşin et al. (2013), Mandelbaum and Zeltyn (2013), and Yefenof et al. (2022). Importantly, this setting is crucially different from ours, in which the balking customers are not observed. In this branch of the literature, the challenge lies in the fact that for customers who are served, we have just a lower bound on their impatience time, making this a “right-censored” estimation problem.

Another related research direction concerns the so-called queue inference engine, a setting in which one does not observe the arrival process of the customers but rather the ordered service entry and service completion times (so-called “transactional data”); see, for example, Larson (1990) and Daily and Servi (1998). One comes across this situation when considering, for example, automatic teller machines, where the queue of customers waiting for the machine (and in particular their arrival times) is not observable, but the transactional data are recorded. In Daley and Servi (1998), a multiserver queue is considered, but with the specific feature that balking is assumed to occur when all servers are busy upon arrival, with a fixed probability independent of the queue length, making the underlying queueing dynamics fundamentally different from ours.

Strictly speaking, Negahban (2019, 2022) did not deal with transactional data; it is assumed that the interarrival times corresponding to the nonbalking customers are observed (so that its scope is somewhat similar to ours). At the methodological level, however, Negahban (2019, 2022) crucially differed from our approach; the proposed estimation procedure relies on discrete-event simulation rather than a rigorously backed maximum-likelihood procedure.

Our work fits in the larger branch of the literature on inference for queues. This domain focuses on estimation problems where the queue is observed (for instance periodically), with the objective to estimate its input parameters. We refer to Asanjarani et al. (2021) for a broad recent overview of this area. Examples are the procedures for the M/G/1 queue, as presented in Hansen and Pitts (2006) and den Boer and Mandjes (2017). We also refer to the Poisson-probing based approach (Ravner et al. 2019) for Lévy-driven queues (Dębicki and Mandjes 2015) and the general framework presented in Duffie and Glynn (2004), as well as the corresponding hypothesis-testing problem (Mandjes and Ravner 2020). Some approaches rely on exploiting knowledge of the queue’s tail probabilities (Courcoubetis et al. 1995, Glynn and Zeevi 2000, Duffy and Meyn 2011). A related line of research concerns the detection of stability or instability of queueing networks; see the Monte Carlo-based methods of Leahu et al. (2017) and Mandjes et al. (2017).

The main conclusion from the above is that, despite the mature literature on queues with impatient customers and related estimation problems, there is a lack of methodologies that can deal with the situation in which balking customers are not observed.

1.3. Paper Organization

This paper is organized as follows. In Section 2 we describe the model featuring in this paper and introduce the notation used. In addition, we present a number of known results on queues with customer impatience that are relied on later. Then, in Section 3, we construct and analyze our maximum-likelihood estimator. For continuous patience-level distributions, we prove strong consistency and asymptotic normality of the estimation error. The case of deterministic patience is treated separately, leading to a similar conclusion regarding strong consistency, but with the notable difference that in this case the limiting distribution of the estimation error is exponential. Section 4 provides detailed analysis of the MLE for exponential, hyperexponential, and generalized hyperexponential patience distributions, the latter cases being relevant because they can be used to approximate the distributions of arbitrary nonnegative random variables. This includes explicit analysis of the likelihood function and its corresponding asymptotic variance as well as simulation experiments assessing the performance of the estimators. Section 5 illustrates how the method can be utilized for several applications. In Section 6, we provide concluding remarks and directions for followup research, including a brief discussion of a nonparametric estimation procedure.

2. Model and Preliminaries

In this section, we provide a detailed model description, introduce the notation used throughout this paper, state the objective of the paper, and present preliminaries.

2.1. Model Description

We consider a service system with s servers ($s\in \mathbb{N}$). Potential customers arrive at the service system according to a Poisson process with rate λ. Each customer has a service requirement with cumulative distribution function (cdf) $G(·)$ and a patience level with cdf $H(·)$. We define B and Y as the generic random variables corresponding to service requirements and impatience levels, respectively. The interarrival times ${(T_i)}_{i=1,2\dots }$, the service requirements ${(B_i)}_{i=1,2,\dots }$, and patience levels ${(Y_i)}_{i=1,2,\dots }$ are independent sequences of i.i.d. distributed random variables. Customer i joins the queue if Y_i is smaller than or equal to the waiting time he or she will experience. More specifically, let V(t) (for $t⩾0$) denote the virtual waiting time at time t, which is defined as the time it would take before at least one server becomes idle if no customers have joined the system after time t:

where Q(t) denotes the number of customers in the system (including those in service) at time t and $N_{\text{D}}(t)$ (for $t⩾0$) denotes the total number of customers finishing their service in a time interval $(0,t]$. For the single-server case (i.e., if s = 1), the virtual waiting time equals to the workload in system, that is, the sum of remaining service requirements of customers. Note that the virtual waiting time in the multiserver case ($s=2,3,\dots$) does not coincide with the workload in system, unlike in the single-server case:

• (i) When $i\in \{0,1,\dots ,s-1\}$ servers are busy, the virtual waiting time equals zero, whereas the workload equals some positive value and decreases with slope i as time passes without effective arrivals.

• (ii) When all s servers are busy, the virtual waiting time and the workload decrease with slope 1 and s, respectively, as time passes without effective arrivals.

Figure 1 illustrates an example of the realization of the virtual waiting time process for s = 2. The actual waiting time of a customer arriving at time t equals $V(t-)$ so that a customer balks if the virtual waiting time just before the arrival instant exceeds the customer’s patience threshold. The resulting system can be considered as an M/G/s queue with impatient customers.

Because some customers will balk, there will be a difference between the interarrival times of the customers (covering both the balking and nonbalking customers) on one hand, and the effective interarrival times (covering just the nonbalking customers) on the other hand. Throughout this paper, the sequence of effective interarrival times is denoted by ${(A_i)}_{i=1,2,\dots }$. Importantly, observe that the random variables A_i are generally neither identically distributed nor independent. We can reconstruct them from the sequences ${(T_i)}_{i=1,2,\dots }$, ${(B_i)}_{i=1,2,\dots }$, and ${(Y_i)}_{i=1,2,\dots }$ in the following manner. For the following construction, we let a nonbalking customer arrive at time t = 0, as will be further motivated in Remark 1. Denoting the virtual waiting time process at time t by V(t), we have that $A_1={\sum }_{k=1}^{j_1}{T}_k$, where

Similarly, $A_i={\sum }_{k=j_{i-1}+1}^{j_i}{T}_k$, where

Observe that j_i has the interpretation of being the index of the i-th nonbalking customer. Let ${\tilde{A}}_i≔{\sum }_{k=1}^i{A}_k$ denote the i-th effective arrival time. We denote the upward-jump size in the virtual waiting time caused by the i-th joining customer by $X_i≔V({\tilde{A}}_i)-V({\tilde{A}}_i-)$. In the single-server case (s = 1), we have $X_i=B_{j_i}$. In the multiserver case ($s⩾2$), on the other hand, X_i takes a complicated form, depending on both the queue-length and the residual service times of customers in service seen on the arrival instant. For a single-server system (s = 1), Figure 2 illustrates the workload process during a single busy period for the case of deterministic patience, that is, $Y=\theta$ for some $\theta >0$. The building blocks of our framework, viz. the interarrival and service requirements and their effective counterparts, are highlighted in the figure. It is stressed that the distribution of the effective interarrival times can significantly deviate from the (exponential) distribution of the interarrival times. Figure 3 illustrates this effect for two different patience-level distributions. In addition, as mentioned, the effective interarrival times are not independent, as opposed to the interarrival times.

We focus on the case that the patience-level distribution is determined by a parameter θ that is a vector in ${\mathbb{R}}^p$ for some (known) $p\in \mathbb{N}$, entailing the estimation problem is parametric. In addition, the arrival rate λ is to be estimated, which is evidently of a parametric nature as well. Alternatively, θ could be a function from some more general space, rendering the estimation problem nonparametric; in Section 6, we briefly discuss an approach that can be used in this setting. Throughout this paper, we assume that we observe the full queue-length process. This in particular means that we observe the arrival and departure epochs of the nonbalking customers (and that we do not observe the balking customers at all). The objective is to use this information to somehow learn the true arrival rate and the parameters of the patience-level distribution. More concretely, we wish to devise statistical procedures for estimating the arrival rate λ and the distribution of the patience Y (both corresponding to nonobservable quantities) with provable performance properties.

2.2. Stationary Distribution

We assume that the mean service time $\text{E}[B]$ is finite and that, with $\rho ≔\lambda \text{E}[B]/s$ denoting the traffic intensity,

Under these assumptions, we have that the queue is stable and that a stationary queue-length $Q≔Q(\infty )$ and a stationary virtual waiting time $V≔V(\infty )$ exist (see Baccelli and Hebuterne 1981 and Baccelli et al. 1984). Note that ${\lim }_{x\to \infty }H(x)<1$ holds only if there are customers with infinite patience, that is, $\text{P}(Y=\infty )>0$. If $\text{P}(Y<\infty )=1$, on the other hand, then we have ${\lim }_{x\to \infty }H(x)=1$, so that the system is stable irrespective of the value of the traffic intensity ρ. As we will see, the MLE method to be developed relies only on the transient dynamics of the system so that it can be applied even if the stability condition (2) is not satisfied. On the other hand, we will utilize the system stability to discuss the asymptotic properties of the proposed estimator.

For the M/G/s queue with impatient customers, exact expressions for the distributions of the stationary queue-length Q and virtual waiting time V are not known in the literature. However, we can derive the following relations they satisfy in a similar way to the single-server case (Kovalenko 1961). Let $q_n≔\text{P}(Q=n)$ ($n=0,1,\dots$) denote the probability mass function of the stationary queue length. The stationary virtual waiting time distribution has probability mass ${\pi }_0≔\text{P}(V=0)={\sum }_{n=0}^{s-1}{q}_n$ at zero, and it is absolutely continuous on $(0,\infty )$. Let $X_{|y}$ denote a generic random variable for the stationary upward-jump size X conditioned that it takes a positive value and that the immediately preceding virtual waiting time equals y:

With the level-crossing argument, we can verify that the density function v(x) of the stationary virtual waiting time satisfies the following integral equation (cf. Boxma et al. 2011):

where ${\overline{J}}_{|y}(x)≔\text{P}(X_{|y}>x)$ denotes the complementary cdf of the conditional upward-jump size $X_{|y}$.

Owing to the PASTA property, the virtual waiting time seen by an arriving customer has the same distribution as the stationary virtual waiting time V. The stationary waiting time W of a nonbalking customer is thus given by a conditional random variable $[V|V⩽Y]$. We can thus identify $\text{P}(W=0)$ and the probability density function (pdf) $w(·)$ of W in terms of the stationary virtual waiting time distribution $v(·)$:

where $P_{\ell }$ denotes the loss probability given by

For s = 1, (3) simplifies to

where $\overline{G}(x)≔1-G(x)$ denotes the complementary cdf of service requirements. A solution of the integral Equation (6) is, as can be found in Baccelli et al. (1984) and Inoue and Takine (2015), given in terms of the pdf $g_{\text{e}}(x)≔\overline{G}(x)/\text{E}[X]$ of the equilibrium distribution (i.e., the residual lifetime distribution) corresponding to the service requirements:

In some special cases, the expression (7) for the virtual waiting time density $v(·)$ can be further simplified, as shown in the examples below.

Example 1

(Constant Patience Levels). We consider the case that the patience levels take a constant value θ₀ (for some ${\theta }_0>0$). Let $g_{\text{e}}^{(n)}(·)$ (for $x⩾0,n=1,2,\dots$) denote the n-fold convolution of the pdf $g_{\text{e}}(·)$ of the equilibrium distribution of the service times and $G_{\text{e}}^{(n)}(x)$ the corresponding cdf. Also, define $g_{\text{e},{\theta }_0}^{(n)}(·)$ by

$$g_{\text{e},{\theta }_0}^{(n)}(x)≔\{\begin{array}{lll}{g}_{\text{e}}^{(n)}(x)/G_{\text{e}}^{(n)}({\theta }_0),\hfill & \hfill & x⩽{\theta }_0,\hfill \\ 0,\hfill & \hfill & x>{\theta }_0.\hfill \end{array}$$

We then have, by de Kok and Tijms (1985) and Inoue and Takine (2015), with as usual “$\mathrm{\star }$” denoting the convolution operator,

$${\pi }_0={\left(1+\rho +\rho {\displaystyle \sum _{n=1}^{\infty }{\rho }^n}{G}_{\text{e}}^{(n)}({\theta }_0)\right)}^{-1},$$(8)

$$v(x)=\{\begin{array}{ll}{\pi }_0{\displaystyle \sum _{n=1}^{\infty }{\rho }^n}{g}_{\text{e}}^{(n)}(x),& 0⩽x⩽{\theta }_0,\\ {\pi }_0\rho g_{\text{e}}(x)+{\pi }_0{\displaystyle \sum _{n=2}^{\infty }{\rho }^n}{G}_{\text{e}}^{(n-1)}({\theta }_0)·[g_{\text{e},{\theta }_0}^{(n-1)}\mathrm{\star }{g}_{\text{e}}](x),& x>{\theta }_0,\end{array}$$(9)

Example 2

(Exponential Service Times). If the service times follow an exponential distribution with mean $1/\mathit{\mu }$, we have, by Stanford (1979) and Inoue and Takine (2015), that

$$\begin{array}{l}{\pi }_0={\left(1+\lambda {\displaystyle {\int }_0^{\infty }\exp }\left[-\mathit{\mu }x+\lambda {\displaystyle {\int }_0^x\text{P}}(Y⩾y)\text{d}y\right]\text{d}x\right)}^{-1},\\ v(x)={\pi }_0\lambda \exp \left[-\mathit{\mu }x+\lambda {\displaystyle {\int }_0^x\text{P}}(Y⩾y)\text{d}y\right],x⩾0.\end{array}$$

Observe that if the patience distribution is also exponential with rate θ, then the expressions can be simplified by substituting ${\displaystyle {\int }_0^x\text{P}}(Y⩾y)\text{d}y=\exp (-\theta x)$.

3. Parametric Estimation Procedure

In the setting considered, we focus on the first nonbalking n + 1 customers. We set the time origin t = 0 to the time instant that a nonbalking customer joins the system. We record the effective arrival times ${\tilde{A}}_1,{\tilde{A}}_2,\dots ,{\tilde{A}}_n$ and departure times. From this information, we can fully reconstruct the virtual waiting time process using (1); in particular, we obtain the sequence of virtual waiting times $V(0),V({\tilde{A}}_1),\dots ,V({\tilde{A}}_n)$ observed by the nonbalking customers and the sizes of upward jumps $X_0,X_1,\dots ,X_n$ caused by them. In this section, the goal is to estimate in a parametric context the arrival rate and patience-level distribution from the data.

The crucial observation is that we can construct a likelihood function of the patience-level distribution using an observed sample of effective interarrival times $A_1,A_2,\dots ,A_n$, waiting times $W_0=V(0),W_1=V({\tilde{A}}_1),W_2=V({\tilde{A}}_2),\dots ,W_n=V({\tilde{A}}_n)$, and upward-jump sizes $X_0,X_1,\dots ,X_n$, which can be reconstructed from the observed sequences mentioned in the previous paragraph. We work under the natural assumption that the observed system is in stationarity, so that W₀ is the stationary waiting time of a nonbalking customer and X₀ is a stationary upward-jump size.

Remark 1.

To make the sequence of waiting times $W_1,W_2,\dots ,W_n$ stationary, we should let the time origin coincide with an arrival instant of a nonbalking customer. The validity of this statement can be argued as follows.

We first point out that the sequence of waiting times is in general not a Markov process; only for s = 1 it is. We obtain a Markov process when considering the vector of residual service times at effective arrival instants. The waiting times and upward jump sizes are in fact measurable functions of this process. Therefore, in what follows, we let the underlying Markov process of residuals be stationary.

First, it is observed that assuming the virtual waiting time at time t = 0 being stationary does not imply that the waiting time W₁ of the first customer follows the stationary distribution W. To see this, consider (for simplicity) an ordinary M/G/s queue (that is, without customer impatience). If the virtual waiting time at time t = 0 follows the stationary virtual waiting time distribution V, then the first customer arriving after time t = 0 finds an idle server with probability

$$\text{P}(W_1=0)=\text{P}(V=0)+{\displaystyle {\int }_0^{\infty }{e}^{-\lambda x}}v(x)\text{d}x>\text{P}(V=0).$$

Because the stationary waiting time W has the same distribution as the stationary virtual waiting time V in the ordinary M/G/s queue (due to PASTA), we find the inequality

$$\text{P}(W_1=0)>\text{P}(V=0)=\text{P}(W=0).$$

This shows that the waiting time W₁ of the first-arriving customer is biased in that it is not distributed as W.

Intuitively, the bias is a consequence of the fact that a customer arriving after a long interarrival time is more likely chosen as the first arriving customer in our experiment (cf. the well-known inspection paradox in renewal theory), and such a customer tends to observe the system less congested than time average. It is easily seen that letting time t = 0 correspond to the arrival instant of a nonbalking customer makes the sequence of waiting times $W_1,W_2,\dots ,W_n$ stationary. □

Let $\overline{H}(x)≔1-H(x)$ denote the complementary cdf of patience levels. Also, we denote by $\tilde{H}(x)$ the probability of a customer joining when the virtual waiting time equals $x⩾0$:

Note that $\tilde{H}(x)=\overline{H}(x)=\text{P}(Y>x)$ for continuous patience distributions.

We next describe the construction of the likelihood function. We first recall that the effective interarrival times $A_1,A_2,\dots ,A_n$ are not exponentially distributed because of the fact that between two effective arrival instants there may have been arriving customers who observed a virtual waiting time level exceeding their patience level (as depicted in Figure 3). Despite this, we can still characterize the effective interarrival time distribution in terms of the observed quantities. To see this, suppose that, for some $t⩾0$, there have been no effective arrivals in the interval $({\tilde{A}}_{i-1},{\tilde{A}}_{i-1}+t]$. An arrival in $({\tilde{A}}_{i-1}+t,{\tilde{A}}_{i-1}+t+\mathrm{\Delta }t]$ then occurs with probability $\lambda \mathrm{\Delta }t+o(\mathrm{\Delta }t)$, as $\mathrm{\Delta }t\downarrow 0$. The corresponding customer joins the system (i.e., becomes effective) with probability $\tilde{H}(V({\tilde{A}}_{i-1})-t)$ because (i) the virtual waiting time seen on the arrival equals $\max \{0,V({\tilde{A}}_{i-1})-t\}$ and (ii) $\tilde{H}(x)=1$ for x < 0. Therefore, the occurrence of the next (i.e., the i-th) effective arrival follows a time-inhomogeneous Poisson process with time-dependent intensity $\lambda \tilde{H}(V({\tilde{A}}_{i-1})-t)$. Noting that $V({\tilde{A}}_{i-1})=W_{i-1}+X_{i-1}$ (for $i=1,2,\dots$), we thus conclude that

This representation of the distribution of the effective interarrival times facilitates the evaluation of the likelihood. The (conditional) likelihood function for $\mathit{A}≔(A_1,A_2,\dots ,A_n)$ given $\mathit{W}≔(W_1,W_2,\dots ,W_n)$ and $\mathit{X}≔(X_1,X_2,\dots ,X_n)$ is given by the product of the conditional densities

This yields the conditional likelihood function

where we used a Lindley-type recursion, generalized to the multiserver queue:

In the rest of this section, we assume that the patience-level distribution is characterized by a finite-dimensional vector of parameters $\theta \in \mathrm{\Theta }\subseteq {\mathbb{R}}^p$. For $\theta \in \mathrm{\Theta }$, let $H_{\theta }(·),{\overline{H}}_{\theta }(·)$, and ${\tilde{H}}_{\theta }(·)$ denote $H(·),\overline{H}(·)$, and $\tilde{H}(·)$ given parameter θ, respectively. Let any $H_{\theta }(·)$, with $\theta \in \mathrm{\Theta }$, be identifiable in the conventional Kullback-Leibler sense. The maximum-likelihood estimator (MLE) of $(\lambda ,\theta )$ is then given by

where

In the sequel, we provide an asymptotic analysis describing the performance of this MLE. We first focus, in Sections 3.1 and 3.2, on the parameters pertaining to the patience only; that is, the estimation does not cover the arrival rate λ. The object of study is, in self-evident notation, for $\lambda >0$ given,

Extending this to a procedure to also include estimation of the arrival rate is relatively straightforward; we get back to it in Section 3.3.

Remark 2.

Observe that L is a conditional likelihood and not a full likelihood. This is due to the fact that the upward jump-sizes X_i have an elaborate distribution both marginally and jointly with the waiting times. The exception is the single-server case where $X_i\sim B$ and is independent of the previous waiting and arrival times. Hence, for s = 1, the maximization of $L_n(H;\mathit{A},\mathit{W}|\mathit{X})$ with respect to them is equivalent to that of the unconditional likelihood function for $(\mathit{A},\mathit{W},\mathit{X})$. It is noted that the procedure to estimate θ does not require the simultaneous estimation of the service requirement or upward-jump size distributions. $\diamond$

In more detail, the remainder of this section is organized as follows. First, we establish (in Section 3.1) conditions for strong consistency and asymptotic normality of our estimator of the patience parameters, with the errors scaled by $\sqrt{n}$, focusing on continuous patience-level distributions. We then consider (in Section 3.2) the case of deterministic patience and establish consistency (independently) for this case as well. It is noted that in the latter case the asymptotic errors are not normally but rather exponentially distributed, with the errors scaled by n. This is due to the fact that for deterministic patience the MLE is obtained on the boundary of the sample data, much like in the well-known case of estimating the parameter θ of a uniform distribution on $[0,\theta ]$. As mentioned, Section 3.3 discusses the estimation of the arrival rate.

3.1. Continuous Patience-Level Distribution

This subsection covers the asymptotic performance of the MLE for the case that the patience-level distribution is continuous and parametric. Let θ₀ denote the true parameter. Throughout the following analysis, the underlying probability measure is ${\text{P}}_{{\theta }_0}$, that is, the probability measure corresponding to the true patience-level distribution. For the asymptotic results, we make the following assumptions, which will be discussed in Section 3.1.1.

Assumption.

The following assumptions are imposed:

• (A1) The parameter space $\mathrm{\Theta }\subset {\mathbb{R}}^p$ is a compact set such that the true parameter lies in the interior of the set, that is, ${\theta }_0\in {\mathrm{\Theta }}^{\text{o}}$.

• (A2) The observation period commences at t = 0, which corresponds to the stationary arrival instant of a nonbalking customer, and thus the sequence of waiting times $W_1,W_2,\dots ,W_n$ is stationary.

• (A3) Let ${\overline{H}}_{\inf }(x)≔{\inf }_{\theta \in \mathrm{\Theta }}{\overline{H}}_{\theta }(x)$. Depending on whether the right endpoint

  $$h_{\sup }≔\inf \{x⩾0:{\overline{H}}_{\inf }(x)=0\}$$
  takes a finite value or not, we assume one of the following properties: (i) If $h_{\sup }=\infty$, then there exists a positive nondecreasing function $f(·)$ such that ${\lim }_{x\to \infty }f(x)\in (0,\infty ]$ and for some constants $c_1,c_2\in [0,\infty )$,
  $$\underset{x\to \infty }{\lim }\frac{e^{-f(x)}}{{\overline{H}}_{\inf }(x)}=c_1,\underset{x\to \infty }{\lim }f(x){\overline{H}}_{{\theta }_0}(x)=c_2.$$
  (ii) If $h_{\sup }<\infty$, then there exists a positive nondecreasing function $f(·)$ such that ${\lim }_{x\to h_{\sup }-}f(x)\in (0,\infty ]$ and for some constants $c_1,c_2\in [0,\infty )$,
  $$\underset{x\to h_{\sup }-}{\lim }\frac{e^{-f(x)}}{{\overline{H}}_{\inf }(x)}=c_1,\underset{x\to h_{\sup }-}{\lim }f(x){\overline{H}}_{{\theta }_0}(x)=c_2.$$

• (A4) The gradient vector and Hessian matrix of $H_{\theta }$ are continuous with respect to θ. With ${\mathbf{\Psi }}_1(\theta )$ denoting the Hessian matrix of the log-likelihood corresponding to a single waiting time $($in stationarity$),\text{\hspace{0.17em}E}{\mathbf{\Psi }}_1(\theta )$ has finite elements in all coordinates for any $\theta \in \mathrm{\Theta }$.

• (A5) The collection of functions $\{H_{\theta }(x):\theta \in \mathrm{\Theta }\}$ is equicontinuous; for any $x⩾0$ and $\mathrm{ϵ}>0$ there exists a $\delta >0$ such that $|H_{\theta }(x)-H_{\theta }(y)|<\mathrm{ϵ}$ for any y such that $|x-y|<\delta$ uniformly over all $\theta \in \mathrm{\Theta }$.

We now state the main results of this subsection. In the sequel, $\text{N}(\mathit{\mu },\mathbf{\Sigma })$ denotes a normally distributed random variable with mean vector μ and covariance matrix $\mathbf{\Sigma }.$

Theorem 1.

If Assumptions (A1)–(A3) and (A5) are satisfied, then, as $n\to \infty$,

$${\widehat{\theta }}_n{\to }_{\text{as}}{\theta }_0.$$(15)

Theorem 2.

If Assumptions (A1)–(A5) hold, then, because $n\to \infty$,

$$\sqrt{n}({\widehat{\theta }}_n-{\theta }_0){\to }_{\text{d}}\text{N}(0,{\mathit{I}}^{-1}({\theta }_0)),$$(16)

where $\mathit{I}({\theta }_0)≔-\text{E}{\mathbf{\Psi }}_1({\theta }_0)$, and ${\mathit{I}}^{-1}({\theta }_0)$ denotes the inverse of $\mathit{I}({\theta }_0)$.

The proofs of the above theorems are provided in Section 3.1.2 after the detailed discussion of (A1)–(A5) that we give in Section 3.1.1.

3.1.1. Discussion of Assumptions.

We next provide more background on the assumptions imposed. In addition, we discuss their possible relaxation.

• (1) Assumption (A1) is natural, because it requires that the parameter space is big enough so that it contains the true parameter. In practice, the parameter space can be adjusted on the fly, for example, if for large n we obtain boundary solutions for the first-order conditions.

• (2) Assumption (A2) facilitates the use of known results on the convergence of stationary dependent sequences. Note that, as long as a stationary virtual waiting time distribution exists (for which we have given the condition as Equation (2) in Section 2), the same results should hold without making this assumption because of the regenerative nature of the process. However, without this assumption, the conditions for both consistency and asymptotic normality are harder to verify. In particular, the crucial step for consistency is the uniform convergence of the log-likelihood, established in Lemma 1 below. More elaborate conditions for uniform convergence are detailed in, for example, Ranga Rao (1962), Andrews (1987), and Pötscher and Prucha (1989).

• (3) If there are patient customers who do not balk regardless of the waiting time, that is, ${\lim }_{x\to \infty }{\overline{H}}_{{\theta }_0}(x)>0$, then Assumption (A3) requires that the parameter space Θ is chosen so that the existence of patient customers (${\lim }_{x\to \infty }{\overline{H}}_{\theta }(x)>0$) is assumed for all $\theta \in \mathrm{\Theta }$, which is a reasonable assumption in modeling a service system where both patient and impatient customers exist. In the case of ${\lim }_{x\to \infty }{\overline{H}}_{{\theta }_0}(x)=0$, on the other hand, Assumption (A3) requires that the decay rate of ${\overline{H}}_{\theta }(x)$ does not vary too strongly among the distributions $H_{\theta }$, with θ in the parameter space Θ. In practice, this assumption is seldom violated. For example, if the true patience-level distribution decays exponentially, that is, ${\lim }_{x\to \infty }{e}^{\nu x}{\overline{H}}_{{\theta }_0}(x)=c$ for some $\nu >0$ and c > 0, then (A3) is satisfied if the infimum tail function ${\overline{H}}_{\inf }(x)$ does not decay faster than the doubly exponential function $e^{-e^{\nu x}}$ (which evidently decays exceptionally fast). As another example, supposing that the true patience-level distribution obeys a power law, that is, ${\lim }_{x\to \infty }{x}^k{\overline{H}}_{{\theta }_0}(x)=c$ for some k > 0 and c > 0, then (A3) is satisfied if ${\overline{H}}_{\inf }(x)$ does not decay faster than $e^{-x^k}$.

• (4) Assumption (A4) enables the construction of a standard martingale CLT for the asymptotic distribution of the estimation error of the MLE. It is not a necessary condition, but in cases where the assumption is not satisfied, one is typically required to apply ad hoc analysis to derive an asymptotic distribution. In Section 3.2, we show that for a deterministic patience level that does not satisfy (A4), the asymptotic distribution of the error is exponential and not normal. Assumption (A2) implies that the expectation of the Hessian matrix of the log-likelihood ${\mathbf{\Psi }}_1(\theta )$ is the covariance matrix of the gradient, and thus it is always positive definite. As a consequence, one just needs to verify that the coordinates are finite for (A4) to hold.

• (5) The equicontinuity assumption (A5) enables concise analysis and can be verified for many continuous patience-level distributions. For example, if there exists a bounded density for every $H_{\theta }(·)$, then it is Lipschitz continuous, and a uniform bound is given by the supremum of the constants in the compact set Θ. The assumption does not hold if the support of the distribution depends on θ. If $H_{\theta }(x)$ has some discontinuities (with respect to x), one could pursue replacing (A5) by an appropriate upper semicontinuity assumption; see theorem 16b in Ferguson (1996) and Pötscher and Prucha (1989). We will present such an example in Section 5.1. Our asymptotic results may hold by replacing (A5) with (A4) and verifying additional measurability conditions on ${\inf }_{\eta \in {\mathcal{B}}_{\theta }}{\overline{H}}_{\eta }(x)$ for a neighborhood ${ℬ}_{\theta }$ of any $\theta \in \mathrm{\Theta }$ (see Andrews 1987, corr. 2).

3.1.2. Proofs.

From (14), first note that we can express the log-likelihood ${\ell }_n(\theta ;\mathit{A},\mathit{W}|\mathit{X})≔\log L_n(\theta ;\mathit{A},\mathit{W}|\mathit{X})$ by

Observe that we replaced ${\tilde{H}}_{\theta }(·)$ with ${\overline{H}}_{\theta }(·)$ because $H_{\theta }(·)$ is continuous.

With $W≔{\lim }_{i\to \infty }{W}_i$ existing, because $n\to \infty$, the continuous mapping theorem yields

where A₁ is the effective interarrival time when the initial virtual waiting time is $W_0+X_0$ and $W_0{=}_{\text{d}}W$, that is, the stationary waiting time. From now on, we use the compact notations ${\ell }_n(\theta )≔{\ell }_n(\theta ;\mathit{A},\mathit{W}|\mathit{X})$ and $\ell (\theta )≔\text{E}[{\ell }_1(\theta ;A_1,W_1,W_0|X_0)]$. The density of the effective interarrival times (11) is uniquely determined by the function $H_{\theta }(·)$, and thus so is the likelihood (14). This is because the functionis uniquely determined by θ as a consequence of the fact that for every $\theta \in \mathrm{\Theta }$ the measure ${\nu }_{\theta }(·)$ corresponds to a different distribution (that is, in the almost-everywhere sense). As a consequence, $\{H_{\theta }(·):\theta \in \mathrm{\Theta }\}$ is a parametric collection of distribution functions such that there is no pair ${\theta }_1,{\theta }_2\in \mathrm{\Theta }$ for which $H_{{\theta }_1}(·)=H_{{\theta }_2}(·)$ almost everywhere. This entails that the model is identifiable in the Kullback-Leibler sense (see, e.g., Ferguson 1996, chapter 17), and hence,

The key step in the proof of Theorem 1 is establishing a uniform version of (18); then, strong consistency follows by the methodology of Wald (Ferguson 1996, chapters 16 and 17). The virtual waiting time observations are not independent, but by (A2) they are stationary so that we can apply a uniform law for stationary sequences that is commonly used in the econometrics literature. We rely specifically on a theorem taken from Ranga Rao (1962) and its extensions, in particular the ones developed in Andrews (1987) and Pötscher and Prucha (1989). The uniform convergence is established in Lemma 1 (proven in the Appendix).

Lemma 1.

If Assumptions (A1)–(A3) and (A5) hold, then, because $n\to \infty$,

$$\underset{\theta \in \mathrm{\Theta }}{\sup }\left|\frac{1}{n}{\ell }_n(\theta )-\ell (\theta )\right|{\to }_{\text{as}}\hspace{0.5em}0.$$(20)

Proof of Theorem 1.

Observe that the strong consistency follows from the identifiability property (19) and from Lemma 1. Note in particular that the proof of Ferguson (1996) (theorem 17) does not rely on i.i.d. observations once uniform convergence has been established, and so the same steps can be applied here. In particular, a sufficient condition (relying on uniform convergence for strong consistency when observations are dependent) can also be found in Heijmans and Magnus (1986); for every $\theta \ne {\theta }_0$, there exists a neighborhood ${ℬ}_{\delta }(\theta )$ such that

$$\underset{n\to \infty }{\lim }\underset{\gamma \in {ℬ}_{\delta }(\theta )}{\sup }\left(\frac{1}{n}{\ell }_n(\gamma )-\frac{1}{n}{\ell }_n({\theta }_0)\right)<0,$$(21)

almost surely. This follows directly from (19) and (20). □

We continue by proving asymptotic normality that was stated in Theorem 2. This amounts to showing that the estimation error $\sqrt{n}({\widehat{\theta }}_n-{\theta }_0)$ is asymptotically normal because $n\to \infty$, assuming that (A1)–(A5) are satisfied. To this end, we apply the well-known delta method and an appropriate version of the martingale CLT (see, e.g., Heyde 1997, theorem 12.6).

Proof of Theorem 2.

For any given v, we let $\nabla {\overline{H}}_{\theta }(v)\in {\mathbb{R}}^p$ and ${\nabla }^2{\overline{H}}_{\theta }(v)\in {\mathbb{R}}^{p\times p}$ denote the gradient and Hessian, respectively, of ${\overline{H}}_{\theta }$. These are both continuous with respect to any coordinate θ_k as a consequence of (A4). Thus, both the gradient and the Hessian of the log-likelihood are continuous functions. These are given by $\nabla {\ell }_n(\theta )≔{\dot{\mathbf{\ell }}}_n(\theta )$ and ${\nabla }^2{\ell }_n(\theta )≔{\sum }_{i=1}^n{\mathbf{\Psi }}_i(\theta )$, respectively, where, for $k=1,\dots ,p$,

$${({\dot{\ell }}_n(\theta ))}_k={\displaystyle \sum _{i=1}^n\left(\frac{{\nabla }_k{\overline{H}}_{\theta }(W_i)}{{\overline{H}}_{\theta }(W_i)}-\lambda {\displaystyle {\int }_0^{A_i}{\nabla }_k}{\overline{H}}_{\theta }(W_{i-1}+X_{i-1}-u)\text{d}u\right)}{\mathbf{1}}_{\{{\overline{H}}_{\theta }(W_i)>0\}},$$(22)

and, for $k,l=1,\dots ,p$,

$$\begin{array}{ll}{\left({\mathbf{\Psi }}_i(\theta )\right)}_{kl}\hfill & =(\frac{{\overline{H}}_{\theta }(W_i){\nabla }_{kl}^2{\overline{H}}_{\theta }(W_i)-{\nabla }_k{\overline{H}}_{\theta }(W_i){\nabla }_l{\overline{H}}_{\theta }(W_i)}{{\overline{H}}_{\theta }{(W_i)}^2}-\lambda {\displaystyle {\int }_0^{A_i}{\nabla }_{kl}^2}{\overline{H}}_{\theta }(W_{i-1}+X_{i-1}-u)\text{d}u){\mathbf{1}}_{\{{\overline{H}}_{\theta }(W_i)>0\}}.\hfill \end{array}$$(23)

Following the lines of the standard delta method (van der Vaart 1998, chapter 3), we consider the expansion of ${\dot{\mathbf{\ell }}}_n(·)$ at the MLE ${\widehat{\theta }}_n$ around the true parameter:

$$\frac{1}{n}{\dot{\mathbf{\ell }}}_n({\widehat{\theta }}_n)=\frac{1}{n}{\dot{\mathbf{\ell }}}_n({\theta }_0)+({\widehat{\theta }}_n-{\theta }_0){\displaystyle {\int }_0^1\frac{1}{n}}{\displaystyle \sum _{i=1}^n{\mathbf{\Psi }}_i}(\alpha {\widehat{\theta }}_n+(1-\alpha ){\theta }_0)\text{d}\alpha .$$(24)

If ${\theta }_0\in \mathrm{\Theta }$, then because of the strong consistency that we found in Theorem 1, we have, because $n\to \infty$, that

$${\widehat{\theta }}_n{\to }_{\text{as}}{\theta }_0.$$

Furthermore, because ${\theta }_0\in \mathrm{\Theta }$ by assumption (A1), the smoothness assumption (A4) implies that, as n grows large, ${\ell }_n(\theta )$ converges to a concave function, and the MLE is given by a sequence of roots ${\widehat{\theta }}_n$ satisfying ${\dot{\mathbf{\ell }}}_n({\widehat{\theta }}_n)=0$; that is, from some N on there is no boundary solution for any $n⩾N$ with probability one. Therefore, (24) yields, with

$$B_n≔-{\displaystyle {\int }_0^1\frac{1}{n}}{\displaystyle \sum _{i=1}^n{\mathbf{\Psi }}_i}(\alpha {\widehat{\theta }}_n+(1-\alpha ){\theta }_0)\text{d}\alpha ,$$

that

$$\underset{n\to \infty }{\lim }\sqrt{n}({\widehat{\theta }}_n-{\theta }_0)B_n\stackrel{\text{d}}{=}\underset{n\to \infty }{\lim }\frac{1}{\sqrt{n}}{\dot{\mathbf{\ell }}}_n({\theta }_0),$$(25)

that is, both sides of (25) converge to the same distribution (if the limits exist). In Lemma 2, we apply the stationarity of W_i in (A2) to show that $B_n{\to }_{\text{as}}\mathit{I}({\theta }_0)$, where $\mathit{I}({\theta }_0)=-\text{E}{\mathbf{\Psi }}_1({\theta }_0)$. Furthermore, ${\dot{\mathbf{\ell }}}_n({\theta }_0)/\sqrt{n}$ is shown to satisfy a martingale CLT with asymptotic variance $\mathit{I}({\theta }_0)$. Because $\text{E}\nabla {\ell }_1({\theta }_0)=0$, Assumption (A2) implies that the asymptotic variance equals the stationary covariance of the gradient

$$\mathit{I}({\theta }_0)=\text{E}{\nabla }^2({\ell }_1({\theta }_0))=\text{E}[\nabla ({\ell }_1({\theta }_0))\nabla {({\ell }_1({\theta }_0))}^{\mathrm{\top }}],$$(26)

and hence, $\mathit{I}({\theta }_0)$ is a positive definite matrix (see Ferguson 1996, chapter 18, for more details). Assumption (A4) further demands that the elements of $\mathit{I}({\theta }_0)$ are all finite, and then combining the above and applying Slutsky’s theorem to (25) yields Theorem 2. □

We are thus left with showing Lemma 2 below; its proof is given in the Appendix.

Lemma 2.

If (A1)–(A5) hold and $\mathit{I}({\theta }_0)=-\text{E}{\mathbf{\Psi }}_1({\theta }_0)$ is a positive definite matrix with finite elements, then (a) as $n\to \infty ,B_n{\to }_{\text{as}}\mathit{I}({\theta }_0)$, and (b) as $n\to \infty$,

$$\frac{1}{\sqrt{n}}{\dot{\ell }}_n({\theta }_0){\to }_{\text{d}}\text{N}(0,\mathit{I}({\theta }_0)).$$(27)

3.2. Constant Patience-Level MLE

We next consider the case where all customers have the same patience level θ₀, that is, $H(y)={\mathbf{1}}_{\{y⩾{\theta }_0\}}$; observe that in this case the stability condition (2) is always satisfied, but the continuity assumptions of Section 3.1 do not apply. As we will see, in this case the properties of the MLE are markedly different from those identified in Section 3.1.

In this setting, the likelihood function (14) reduces to, with the event ${\text{ℰ}}_i$ denoting $\{W_{i-1}+X_{i-1}⩽\theta \}$ and ${\text{ℰ}}_i^{\text{c}}$ its complement,

Observe that $L(\theta ;\mathit{A},\mathit{W}|\mathit{X})$ (i) equals zero for $\theta \in [0,{\max }_{i=1,2,\dots ,n}{W}_i)$, (ii) has an upward discontinuity at $\theta ={\max }_{i=1,2,\dots ,n}{W}_i$ and decreases for

and (iii) takes a constant value for $\theta ⩾{\max }_{i=1,2,\dots ,n}\{W_i+X_i\}$. We thus find the intuitively appealing property that the MLE is given bythat is, the maximum virtual waiting time at jump times. Note that this estimator resembles the MLE of the parameter θ when the observations are uniformly distributed on $[0,\theta ]$.

Although the estimator ${\widehat{\theta }}_n$ is clearly biased (which follows from ${\text{E}}_{{\theta }_0}[{\widehat{\theta }}_n]<{\theta }_0$ for all n), we show that it converges almost surely to θ₀ as $n\to \infty$. In addition, we prove that the estimation error scaled at rate n converges to an exponential random variable.

Theorem 3.

Because $n\to \infty$,

$${\widehat{\theta }}_n{\to }_{\text{as}}{\theta }_0,$$(29)

and

$$\text{P}(n({\theta }_0-{\widehat{\theta }}_n)⩾x)\to {\text{e}}^{-v({\theta }_0)x/(1-P_{\ell })},$$(30)

where $P_{\ell }$ and $v({\theta }_0)$ denote the stationary loss probability and the virtual waiting time density at level θ₀. In addition, the asymptotic variance of the estimation error agrees with the variance of the limiting exponential distribution and is given by

$$\underset{n\to \infty }{\lim }{n}^2\text{Var}\left[{\theta }_0-{\widehat{\theta }}_n\right]=\underset{n\to \infty }{\lim }\left\{\text{E}\left[{\{n({\theta }_0-{\widehat{\theta }}_n)\}}^2\right]-{\left(\text{E}\left[n({\theta }_0-{\widehat{\theta }}_n)\right]\right)}^2\right\}$$(31)

$$={\left(\frac{1-P_{\ell }}{v({\theta }_0)}\right)}^2.$$(32)

The proof of Theorem 3, which is given in the Appendix, constructs lower and upper bounds for the MLE and establishes that they both converge to θ₀. The same bounds are also used to characterize the asymptotic distribution of the estimation error.

3.3. Estimating the Arrival Rate

Where in the preceding subsections we focused on estimating the patience-level distribution (in a parametric context) for a given value of the arrival rate λ, we now discuss the estimation of λ. Let λ₀ denote the true value of λ. There are several ways to estimate it, the most basic one relying on the idle period observations. Denote the duration of the k-th idle period by I_k, where in the multiserver setting (s > 1) an idle period refers to time intervals such that at least one server is idle. Denote by E_k the total number of effective arrivals during idle period $k=1,2,\dots$ (meaning, for instance, E_k = 1 for all $k=1,2,\dots$ when s = 1). Let C_n denote the number of idle periods observed up until (and including) the n-th effective arrival. Then we propose the estimator

The rationale behind this estimator is that there is no balking during idle periods because the virtual waiting time is zero. Hence, the arrival process during these idle times is homogeneous Poisson with rate λ. Therefore, the estimator ${\widehat{\lambda }}_n$ is a standard MLE of the rate parameter of an exponentially distributed random variable and satisfies all desired asymptotic properties.

Of course, the above procedure does not exploit a substantial amount of potentially useful data that is collected during busy periods. For the case of a deterministic patience level θ₀, the above estimator is easily improved upon. Recall that ${\widehat{\theta }}_n$ was defined as $\max \{W_1,W_2,\dots ,W_n\}$, with the immediate consequence that ${\widehat{\theta }}_n⩽{\theta }_0$. Supposing we observe an arrival such that immediately after this arrival the virtual waiting time level is still below the current value of the estimator, then the next arrival is an effective arrival and hence, occurs after an exponentially distributed random variable with rate λ. In this way, we generate more observations, thus allowing us to estimate λ with a better precision. Observe that the new observations and the idle times form an i.i.d. sequence.

For the case of a continuous patience-level distribution, one may use (33) or, alternatively, set up a joint MLE for λ and θ. The latter has clear advantages, in particular for small samples or heavily loaded systems. Because the log-likelihood function (17) is a smooth and concave function with respect to λ, the MLE of the arrival rate for any estimator ${\breve{\theta }}_n$ of θ is

The asymptotic results of Section 3.1 can, therefore, be extended in a straightforward manner so as to cover the joint estimation of $(\lambda ,\theta )$. We demonstrate this in Section 4.1, where we detail a procedure for jointly estimating the arrival rate and the patience parameter for the case of an exponentially distributed patience level.

4. Exponential and Generalized Hyperexponential Patience

In this section, we discuss a robust and practical approach for estimating continuous patience distributions. In our approach, this is achieved by fitting the MLE of a generalized hyperexponential (GHE) distribution. This approach is attractive because the class of GHE distributions is known to be dense in the space of nonnegative continuous distributions (see, for instance, Botta and Harris 1986), which, in practical terms, means that any nonnegative continuous distribution can be approximated arbitrarily closely by a GHE distribution. In our simulation-based experiments, highly accurate estimates are obtained, even when the baseline patience distribution itself is not GHE.

We first provide a detailed analysis of the joint MLE for the arrival rate and the single parameter of an exponential patience distribution to then move to the cases of hyperexponential and generalized hyperexponential patience distributions. We also present a heuristic search method that fits an approximate GHE distribution for any continuous patience distribution.

4.1. Exponentially Distributed Patience

Suppose that the arrival rate λ₀ is unknown and that the patience-level distribution is exponential, that is, $H_{{\theta }_0}(x)=1-{\text{e}}^{-{\theta }_{0}x}$ for an unknown θ₀ (and $x⩾0$). Because of (17), it can be verified, by distinguishing between the cases $W_0+X_0⩾A_1$ and $W_0+X_0<A_1$, that the log-likelihood for a single observation is

here the Lindley recursion (13) has been used.

Our objective is to analyze the MLE of both λ and θ for a sample of size n, which we denote by $({\widehat{\lambda }}_n,{\widehat{\theta }}_n)$. We assume that (A1) and (A2) hold, that is, compact parameter space and stationary W₀. In addition, (A3) holds because $H_{\theta }$ is exponential for all $\theta \in \mathrm{\Theta }$. The log-likelihood ${\ell }_1(\lambda ,\theta )$ is Lipschitz continuous (with respect to the observations), and so (A5) holds. We thus have that Theorem 1 holds and the MLE for $(\lambda ,\theta )$ is strongly consistent.

Clearly, this is a smooth function. We proceed by computing the gradient and Hessian with respect to $(\theta ,\lambda )$. It takes some elementary calculus to verify that gradient is given by

For any given value of θ, ${\ell }_1(\lambda ,\theta )$ is a concave function in λ, and therefore, the MLE of (34) is given by

The MLE ${\widehat{\theta }}_n$ now follows by maximizing $\frac{1}{n}{\sum }_{i=1}^n{\ell }_1(\theta ,{\widehat{\lambda }}_n(\theta ))$ for $\theta \in \mathrm{\Theta }$. The optimizing θ is obtained either on the boundary of the parameter space Θ or by solving the first-order condition $\nabla ({\ell }_1{(\theta ,{\widehat{\lambda }}_n(\theta ))}_{\theta }=0$. Because $n\to \infty$, we are guaranteed to find an interior solution as long as ${\theta }_0\in {\mathrm{\Theta }}^{\text{o}}$.

Taking second derivatives, we find that the entries of the Hessian matrix are given by

Because $w^k{\text{e}}^{-\theta w}$ is a bounded function on $w\in [0,\infty )$ for any $\theta >0$ and $k\in \mathbb{N}$, we conclude that $-\text{E}{\nabla }^2({\ell }_1({\lambda }_0,{\theta }_0))<\infty$, and therefore, (A4) applies, which implies that Theorem 2 holds. In particular, because $n\to \infty$,

and the asymptotic covariance is given by the inverse of the Hessian,

We performed numerical experiments to assess the performance of the estimation procedure. Table 1 presents approximated confidence intervals for the maximum-likelihood estimators $({\widehat{\lambda }}_n,{\widehat{\theta }}_n)$; in addition, we also evaluated the estimator of the arrival rate based on idle periods, denoted by ${\tilde{\lambda }}_n$, as discussed in Section 3.3. The confidence intervals are evaluated for three congestion levels, namely $\rho =\lambda \text{\hspace{0.17em}E}X\in \{0.5,1,2\}$. In all experiments, the sample size (of observed waiting times) was $n=1000$. Evidently, however, the heavier the system load, the fewer the number of idle periods; on average, $C_n=C_{1000}$ equals 665, 457, and 267 for the three congestion levels.

Table 1. Confidence Intervals for the MLE s $({\widehat{\theta }}_n,{\widehat{\lambda }}_n)$ as Well as the Arrival Rate Estimated from Idle Periods ${\tilde{\lambda }}_n$ for Different Confidence Levels

The numerical output is summarized in Table 1. A first observation is that the MLE for ${\widehat{\theta }}_n$ is more accurate as the system load increases. This is because the patience threshold is reached more often, and therefore, the data are more informative. Furthermore, for a lightly loaded system, the MLE ${\widehat{\lambda }}_n$ and idle-period based estimator ${\tilde{\lambda }}_n$ yield similar (accurate) results, as opposed to the high-load regime, in which the two estimators behave quite differently:

• As a consequence of the fact that idle periods are observed considerably less frequently for a high load, ${\tilde{\lambda }}_n$ becomes substantially less accurate.

• The accuracy of ${\widehat{\lambda }}_n$, however, is only slightly reduced in the high-load regime. This is potentially because of the better accuracy of the estimation of θ, being jointly estimated with λ in the MLE procedure.

In the most heavily loaded example (ρ = 2, that is), the 95% confidence interval of the MLE is very similar to the 80% confidence interval of the idle period-based estimator. We conclude that even though the observations of the waiting times and the effective interarrival observations are highly dependent, the MLE ${\widehat{\lambda }}_n$ provides substantially better confidence intervals than the idle period-based estimator ${\tilde{\lambda }}_n$.

Table 2 presents empirical confidence intervals for the maximum-likelihood estimators of $(\lambda ,\theta )$ for a multiserver system with s = 5. The confidence intervals are wider than in the single-server case, even though the sample size was taken twice as big, being indicative of the fact that the variance of the estimation error in a system with multiple servers is higher. As in the single-server case, lower system load increases the variance of the estimation error for the patience parameter θ. This effect is even stronger in the multiserver setting because of the fact that balking occurs only when all five servers are working, and such a state is not frequently observed. For higher load, the estimation of the patience parameter is much more accurate, as was the case for the single-server system. For the arrival rate, we observe high accuracy for all load levels.

Table 2. Empirical Confidence Intervals for the MLEs $({\widehat{\theta }}_n,{\widehat{\lambda }}_n)$ for Different Confidence Levels

4.2. (Generalized) Hyperexponential Patience

This subsection focuses on the case of the patience-level distribution being generalized hyperexponential (GHE), meaning that the corresponding cdf can be written as a mixture of exponential cdfs with weights that sum to 1 but that are not necessarily positive; this is in contrast with the standard hyperexponential (HE) distribution, where the weights are assumed to be positive. This case is particularly relevant because, as argued in Roubos and Jouini (2013), this distribution has a good empirical fit-to-patience data. Another motivation for considering this distribution lies in the known fact that the cdf of any continuous positive random variable can be approximated arbitrarily accurately (in terms of a suitably defined metric) by a GHE cdf (Botta and Harris 1986).

For convenience, we now assume the arrival rate λ to be known, but it is noted that it can be estimated in a similar manner as in the exponential case discussed in Section 4.1. Suppose the degree of the GHE distribution is $p\in \mathbb{N}$:

where $\theta =({\alpha }_1,\dots ,{\alpha }_p,{\beta }_1,\dots ,{\beta }_p)\in {\mathbb{R}}^{2p},{\sum }_{k=1}^p{\alpha }_k=1$ (where, importantly, we do not assume the α_k to be positive) and ${\beta }_k>0$ for $k=1,\dots ,p$. Without loss of generality, we assume that ${\alpha }_1>{\alpha }_2>\cdots >{\alpha }_p$. Denote the coordinates of the true parameter θ₀ by $({\alpha }_{0k},{\beta }_{0k})$ for $k=1,\dots ,p$. It means that we are to identify the 2p-dimensional parameter vector

Similarly to the previous example of the exponential distribution, Assumptions (A1)–(A3) and (A5) are satisfied when assuming stationary waiting times and a compact parameter space Θ. Therefore, the MLE is strongly consistent by Theorem 1. One should be cautious when fitting a GHE distribution because of further conditions to be imposed on the parameters to make sure that ${\overline{H}}_{\theta }(·)$ is a proper cdf. Even though convergence to the true parameters of the distribution is eventually guaranteed, for finite samples the estimated parameters may not yield a proper distribution. For an in-depth discussion on these conditions and the identifiability of GHE distributions, we refer to Botta and Harris (1986) (Section 3).

Using the representation (17), the log-likelihood pertaining to a single observation is calculated in a similar manner as in the exponential case, yielding

As follows with some elementary algebra, the gradient is then given by, for $j=1,\dots ,p$,

The Hessian can be derived by computing the matrix of second derivatives. For the evaluation of the asymptotic covariance, a convenient alternative is to apply (26), that is,

so as to avoid the symbolic evaluation of the matrix of second derivatives. The entries of $\mathit{I}({\theta }_0)$ are finite because these are combinations of terms of the following types:

• Products of polynomials (of degree at most 2) and exponentials. For example, we come across a term that is, up to a multiplicative constant, $W_0^2{\text{e}}^{-{\beta }_j{W}_0}$. Observe that the mapping $x\mapsto x^k{e}^{-x}$ is bounded for $x⩾0.$

• Ratios of the form, for example, up to a multiplicative constant

  $$\left|\frac{W_1{\text{e}}^{-{\beta }_j{W}_1}}{{\sum }_{k=1}^p{\alpha }_k{\text{e}}^{-{\beta }_k{W}_1}}\right|\mathrm{⩽}{W}_1.$$

In addition, it is verified that if the patience-level distribution is light-tailed (which is the case for the GHE distribution), then the stationary waiting time W is also light-tailed and has finite moments. Indeed, note that (3) implies

so that we have from (4) that the density of W is bounded as $w(x)⩽\lambda \overline{H}(x)$; that is, if the patience level distribution is light-tailed, so is W. As a consequence, from the above and Theorem 2, the estimation errors (that is, scaled by $\sqrt{n}$) converge to a multivariate normal distribution.

Although there are theoretical guarantees for the asymptotic performance of the MLE, computation of the MLE is not straightforward, even for small parameter spaces such as p = 3. It requires maximizing ${\ell }_n(\theta )/n$, that is, solving a nonlinear and nonconcave program in $p\times (p-1)$ variables, with both equality and inequality constraints. This is computationally highly challenging, and standard optimization methods may lead to local maxima. In particular, observe that $\nabla {({\ell }_1(\theta ))}_{{\beta }_j}$ is very small for large values of β_j, which implies that ${\ell }_n(\theta )/n$ displays very “flat” behavior for large values of these β_j. Our experiments revealed that a direct implementation using standard optimization packages often led to points that were even not local maxima. As a consequence, we decided to write dedicated code to compute the MLE.

In light of the inherent complexity of maximizing the likelihood function, in our numerical experiments we have applied the following nested two-step heuristic optimization method.

• (a) For each vector of $\mathit{\alpha }=({\alpha }_1,\dots ,{\alpha }_p)$, the objective function ${\ell }_n(\theta )/n$ was maximized with respect to $\mathit{\beta }=({\beta }_1,\dots ,{\beta }_p)$ using a conventional coordinate descent algorithm. This step was typically fast and accurate; as for given $\mathit{\alpha }$, the objective function behaves nicely. Now we have reduced the problem to an optimization over $\mathit{\alpha }.$

• (b) Then, a standard L-BFGS quasi-Newton method (see e.g., Byrd et al. 1995) is applied to ${\ell }_n(\theta )/n$, so as to search for the optimal vector $\mathit{\alpha }$, with $\mathit{\beta }$ being parameterized by $\mathit{\alpha }$. The optimization is carried out with the following constraints that ensure that the parameters yield a proper distribution (Botta and Harris 1986):

  $$\left\{\mathit{\alpha }:{\alpha }_1>{\alpha }_2>\cdots >{\alpha }_p,{\displaystyle \sum _{i=1}^p{\alpha }_i}=1,{\displaystyle \sum _{i=1}^p{\alpha }_i}{\beta }_i>0\right\}.$$

Note that there is no firm guarantee that this heuristic method converges to the optimal parameters. To overcome this, the search for the optimizing vector $\mathit{\alpha }$ was conducted multiple times for different initial values. The resulting method turned out to be time-consuming, especially for a relatively large sample size n and/or a relatively large dimension p of the parameter space. However, it typically returns considerably more robust results than off-the-shelf optimization routines (that in addition tended to converge very slowly).

Table 3 presents the marginal confidence intervals obtained by the normal approximation for an example of a HE cdf $H_{\theta }(·)$ with p = 2. Observe that the variance of the estimation error is large, even for a substantial sample size. The accuracy is much higher for the low rate of ${\beta }_1=0.25$ than the higher rate of ${\beta }_2=1$. This may be explained by the fact that the likelihood function is almost flat for high values of β₂.

Table 3. Confidence Intervals for the MLEs $({\widehat{\alpha }}_{1n},{\widehat{\beta }}_{1n},{\widehat{\beta }}_{2n})$ Based on a Normal Approximation for a Sample Size of $n=100000$ for Different Confidence Levels

The results presented in Table 3 show that the estimates are concentrated around the true values. The variance of the estimation error, however, is high, even more so given the large number of observations used. It should be realized, though, that HE distributions are not always easily distinguishable, in the sense that seemingly different HE distributions may lead to a very similar cdf. Clearly, from a practical standpoint, the main question is not whether the correct parameters are recovered but rather whether the cdf corresponding to the estimated parameters provides a good fit with the true cdf. In this respect, the main conclusion of our experiments is that our MLE procedure performs remarkably well. For example, in Figure 4, we illustrate this by presenting the true cdf and the estimated cdf based on four random samples (using $n=10000$ observations). We have used the same parameters as in Table 3, with ρ = 1. We observe that, although the fitted parameters greatly differ from the correct ones, the fit of the (complementary) cdf is highly accurate.

In the remainder of this subsection, we further study the performance of the HE MLE. Our experiments lead to the conclusion that HE MLE provides a good fit even in cases where the true distribution is not HE (i.e., in misspecified scenarios). The examples focus mainly on fitting the target distribution by a (conventional) HE distribution, rather than a GHE distribution, by computing the MLE. The motivation behind this choice is that, when working with a GHE distribution, there is the additional complication that the estimated parameters may not yield a proper cdf. More specifically, the space of valid parameters is hard to characterize, let alone to be coded in terms of constraints of an optimization problem. Nevertheless, we would still like to exploit the extra versatility that the GHE class offers so as to improve the fit (when compared with the HE MLE). This is especially relevant for scenarios in which the model is misspecified, bearing in mind that, as mentioned earlier, in principle any distribution can be arbitrarily accurately approximated by a GHE distribution.

With the above considerations in mind, we implemented the following (seemingly naïve) heuristic model selection method. We generate various random weight vectors $\mathit{\alpha }$ (equipped with a random size $p\in \mathbb{N}$). Then, for each of them, we optimize the likelihood function over $\mathit{\beta }$ using only step (a) above. The best model is then chosen by comparing the various combinations of $\mathit{\alpha }$ and p, relying on the Akaike Information Criterion (AIC). The AIC encompasses both the log-likelihood and, in order to avoid overfitting, a penalty for the number of parameters (Burnham and Anderson 2002). Because step (a) can be performed highly efficiently, a main advantage of this heuristic is that it works very fast, even for bigger values of p, thus providing us with a technique to fit general continuous distributions. Recall from Botta and Harris (1986) that if we order the components of the weight vector (i.e., ${\alpha }_1>{\alpha }_2>\cdots >{\alpha }_p$), then a sufficient condition for $(\mathit{\alpha },\mathit{\beta })$ to yield a cdf is that ${\alpha }_1>0$ and ${\sum }_{i=1}^p{\alpha }_i{\beta }_i>0$; we use this principle to select feasible solutions produced by the above heuristic.

Figure 5, (a) and (c), illustrates the estimated survival function of using the HE MLE with varying weights for three examples of patience-level distributions. In the first example, the patience-level distribution is indeed HE with p = 4, whereas the second and third examples are misspecified (corresponding to a lognormal and gamma patience-level distribution, respectively). For each distribution, the MLE was computed using $p=1,$ 2, and 4 (where p = 3 has been left out because it is barely distinguishable from p = 4).

• As was the case in the setting of Figure 4, the experiments corresponding to the HE patience-level distribution show that, even for p = 4 and as many as $n=10000$ waiting-time observations, the MLE does not accurately capture the true parameters. Nevertheless, the fit in terms of the cdf, as displayed in Figure 5(a), is remarkably good. We in addition performed the GHE fitting heuristic, which also provides a highly accurate fit and is considerably faster to compute than the MLE for p = 4.

• The case of lognormal patience, as illustrated by Figure 5(b), shows that the fit is quite good for all estimators, except for p = 1 (i.e., an exponential distribution).

• Figure 5(c) presents the fitted distributions for gamma patience. The fit of the HE MLE is decent, but importantly, in this case the GHE heuristic performs considerably better than the HE MLE (for p = 1, 2, 4). The GHE heuristic selects a model with p = 10 weights.

The above observations are replicated for a multiserver system with s = 10. Figure 6, (a) and (b), plots the true and estimated patience distribution for lognormal and gamma distributions, respectively, using the misspecified hyperexponential MLE. As in the single-server case, the GHE approximation yields an accurate fit.

5. Applications

In this section, we discuss a series of applications in which our estimation methodology can be directly used.

5.1. Unknown Proportion of Impatient Customers

In this subsection, we consider a system where an unknown proportion $\theta \in (0,1)$ of the arriving customers has a known deterministic patience threshold w > 0. The other customers are patient and always join the queue. Suppose one wishes to estimate the fraction θ from data.

To this end, first observe that

Assumption (A4) is not satisfied because $H_{\theta }(x)$ is not continuous (in x) at x = w. However, with respect to θ, we do have that $H_{\theta }(x)$ is continuous. In addition, ${\ell }_n(\theta )$ is concave, and therefore, the MLE is given by the first-order condition ${\ell }_n^{\prime }({\widehat{\theta }}_n)=0$. In this case, we can apply results from Heijmans and Magnus (1986) to establish strong consistency. In particular, the smoothness of ${\ell }_n^{\prime }(\theta )$ enables direct verification of the sufficient condition (21). Asymptotic normality then follows by verifying the conditions of Theorem 2 directly.

5.2. Noisy Delay Messages

Suppose that the patience threshold is a constant θ, but the customers do not observe their exact waiting time but rather some noisy estimate $W_e=e(W)$ for some random function e such that $\text{E}[W_e|W]=W$. We assume that customers join based on this “perturbed delay” W_e, that is, if $W_e⩽\theta$. Suppose we are in the context that the parameters underlying the noise distribution are known, but the threshold θ is not.

The probability of joining at virtual waiting time level W can be computed given the specific noise distribution. A few examples are as follows:

• Additive perturbations. In this case, $W_e=W+\mathrm{ϵ}$. One could, for instance, consider normally distributed perturbations: $\mathrm{ϵ}\sim \text{N}(0,{\sigma }^2)$, independent of W, with $\sigma >0$. Let customers facing $W_e<\theta$ join the system. We thus have, with $\mathrm{\Phi }(·)$ the cdf of the standard normal distribution,

  $${\overline{H}}_{\theta }(v)=\text{P}(W_e⩽\theta |W=v)=\text{P}(\mathrm{ϵ}⩽\theta -v)=\mathrm{\Phi }\left(\frac{\theta -v}{\sigma }\right).$$

  The asymptotic properties of Section 3.1 hold in this case because $\mathrm{\Phi }(·)$ satisfies the regularity conditions.

• Multiplicative perturbations. Now, W_e = WG, with G nonnegative unit-mean and independent of W. In this case,

  $${\overline{H}}_{\theta }(v)=\text{P}(W_e⩽\theta |W=v)=\text{P}\left(G⩽\frac{\theta }{v}\right).$$

Now, if the random variable G is such that the regularity conditions of Section 3.1 are met, then it follows that the asymptotic properties hold in this case as well (with an asymptotic variance that can be expressed in terms of G).

5.3. Admission Pricing

Our estimation procedure can be exploited in the context of various problems rooted in operations research. Evidently, when having estimates of the arrival rate and the patience-level distribution, one could consider the option of increasing the service rate so as to potentially raise profits. Thus, the operational decision to be made is whether the increased revenues outweigh the cost of speeding up the service rate. In this subsection, we consider another operations research-related problem.

Suppose that the queue has an admission price of p and that customers are homogenous with a utility function as featured in the Naor (1969) model. More concretely, let there be constants r and c such that a customer will join the queue only if the virtual waiting time w upon the customer’s arrival satisfies

or, alternatively, w is smaller than the threshold value $(r-p)/c.$ Now observe that for any fixed price p the threshold $\theta (p)=(r-p)/c$ can be estimated using the MLE procedure presented in Section 3.2. If one of the cost function parameters r and c is known, then the other parameter can be estimated directly. If both are unknown, then their estimation can be performed by an exploration procedure; set two prices, say, p₁ and p₂, and observe the system for each price. Supposing each of these two experiments is done with n clients, we let ${\widehat{\theta }}_n(p_i)$ be the MLE for price i = 1, 2. Then the estimators for the cost function parameters r and c are given by, respectively,

Theorem 3, in combination with the continuous mapping theorem, implies that, because $n\to \infty$, both ${\widehat{r}}_n{\to }_{\text{as}}r$ and ${\widehat{c}}_n{\to }_{\text{as}}c$. Furthermore, for large n, confidence intervals for the stationary average revenue per unit of time can be approximated using (30). In particular, observing that the loss probability depends on the price, in self-evident notation,

because $n\to \infty .$ This opens the possibility of approximating, for n large, the (joint) distribution of $n\left({\widehat{r}}_n-r_0,{\widehat{c}}_n-c_0\right))$. Note that the loss probability $P_{\ell }$ needs to be estimated as well. For the single-server case, this can be done by estimating the idle probability upon arrival for a given price and using (5).

6. Concluding Remarks

This paper has considered a service system in which clients potentially balk based on the virtual waiting time level they face at arrival. The main objective concerned the development of a framework for estimating the arrival rate and patience-level distribution. Our approach resolves the complication that in our setup only nonbalking clients are observed. Distinguishing between the case of a continuous patience-level distribution and constant patience, we developed MLE estimators and quantified their asymptotic properties. Through a sequence of examples and ramifications, we have illustrated the performance and broad applicability of our findings.

An important next step could concern the extension to a nonstationary arrival process. For example, one may assume that the potential arrival process is a nonhomogeneous Poisson process with an arrival rate function that depends on time, that is, $\{\lambda (t)\text{}:t⩾0\}$. This means that when the virtual waiting time is v, the effective arrival rate at time t is given by $\lambda (t)(1-H(v))$, and the likelihood function in Section 2 can be updated accordingly. If the arrival rate function is known, then our estimation procedure for θ essentially carries over, including its performance guarantees. If the arrival rate function is unknown and parametric assumptions about it are made, then the joint estimation of the arrival rate and patience parameters is possible. Furthermore, in the practically relevant case that the time-dependent arrival rate is cyclic (i.e, $\lambda (t)=\lambda (t+s)$ for all $t⩾0$ and some cycle length $s⩾0$; think of daily or weekly patterns), then the queueing process will still have regenerative dynamics that can be exploited for asymptotic analysis. It is thus anticipated that the framework laid out in this paper can serve as a basis for developing estimation techniques for a wide class of nonstationary models.

Although we have developed estimation procedures assuming the independence between the customers’ service requirements and patience levels, one could consider an extension that allows dependence. Such an extension could, for example, model the situation in which customers with larger service requirements can be assumed to have more patience. Multiclass queueing models can be used to capture such dependence in service requirements and patience levels (Sakuma and Takine 2017). More specifically, suppose that there are two customer classes 1 and 2, where class k (k = 1, 2) customers arrive with rate λ_k and have service requirements (resp. patience levels), i.i.d. according to a class-dependent cdf $G_k(·)$ (resp. $H_k(·)$). It is natural to assume that the customer class is unobservable by the estimator, thus yielding a sequence of dependent service requirements and patience levels. In this setting, we have to estimate the service-requirement cdf $G_k(·)$ jointly with the patience-level cdf $H_k(·)$ (and the arrival rate λ_k) to compute the likelihood function corresponding to an observed sequence $(\mathit{A},\mathit{W},\mathit{X})$.

Another possible direction for future work concerns the situation in which each customer decides to balk based on his or her sojourn time, that is, the virtual waiting time seen by this customer just after (instead of before) his or her arrival. If the customer precisely knows his or her service requirement, this problem is easily reduced to the one considered in this paper. The other obvious option is that the customer balks if the virtual waiting time just before his or her arrival increased by a “guess” of his or her service requirement exceeds the customer’s patience threshold.

An alternative to the approach followed in the paper would be to pursue nonparametric estimation, cf. the results in Hansen and Pitts (2006) for the conventional case without balking. In the single-server case, it may be of help that we can write, with W(x) and D(x) denoting the cdf s of the waiting times and sojourn times of nonbalking customers, respectively,

if w(x) > 0 and 0 else; the validity of (38) follows from combining (4) with a level-crossing identity. Whereas it is clear how to estimate W(x) and D(x) using the empirical distribution, estimation of w(x) is more challenging; kernel-based techniques may turn out useful in this context.

A final issue to consider is that in many systems customers do not observe exact waiting times but rather queue lengths. Specifically, this can be modeled as an M/G/s + H system with customers that join or balk after having observed the number of customers in the system. In this case, the patience level Y is a discrete random variable indicating at what queue length a customer is willing to join. Therefore, the interarrival process at any given queue length $q⩾0$ is a Poisson process with rate ${\lambda }_q=\lambda {\overline{H}}_{\theta }(q)$. As was done in this paper, and MLE for the arrival rates (and the corresponding parameter θ) can be derived from a sample of interarrival times and the respective queue lengths.