Call center service process analysis: Bayesian parametric and semi-parametric mixture modeling

European Journal of Operational Research 234 (2014) 709–719

Contents lists available at ScienceDirect

European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor

Stochastics and Statistics

Call center service process analysis: Bayesian parametric and semi-parametric mixture modeling Tevfik Aktekin ⇑ Department of Decision Sciences, University of New Hampshire, Durham, NH 03824, USA

a r t i c l e

i n f o

Article history: Received 10 December 2012 Accepted 27 October 2013 Available online 6 November 2013 Keywords: Call center Bayesian finite mixture Dirichlet processes Service modeling Queuing Staffing

a b s t r a c t The main focus of the call center research has been on models that assume all input distributions are known in queuing theory which gives birth to staffing and the estimation of operating characteristics. Studies investigating uncertainty of the input distributions and its implications on call center management are scarce. This study attempts to fill this gap by analyzing the call center service distribution behavior by using Bayesian parametric and semi-parametric mixture models that are capable of exhibiting non-standard behavior such as multi-modality, skewness and excess kurtosis motivated by real call center data. The study is motivated by the observation that different customer profiles might require different agent skill sets which can create additional sources of uncertainty in the behavior of service distributions. In estimating model parameters, Markov chain Monte Carlo methods such as the Gibbs sampler and the reversible jump algorithms are presented and the implications of using such models on system performance and staffing are discussed. Ó 2013 Elsevier B.V. All rights reserved.

1. Introduction In the last decade call center operations became an inevitable part of many global businesses. Thanks to the technological advancement in the communications field, call centers grew to be the milieu through which firms interact with their potential and existing customers. The complexities involved in managing call center operations gave birth to managerial challenges worth investigating. A call center operation is typically characterized by three fundamental processes of arrival, service and abandonment. Each of these processes offer challenges from both modeling and practical points of view. One of the major concerns for call center customers is fast and efficient service which can effect the perceived quality of a call center system via the fraction of customers that abandon. Call center customers are usually served by service agents by a phone or an online-chat feature if their needs are not satisfied by automatic voice response units. The challenge from a managerial point of view is to determine the optimal number of agents without compromising the perceived quality of the operation. The main focus of the call center related research has been on models that assume all input distributions are known in queuing theory. Studies which address statistical inference issues are scarce with the exception of Brown et al. (2005) and Aldor-Noiman, Feigin, and Mandelbaum (2009). As pointed out by Mandelbaum ⇑ Tel.: +1 603 862 3387. E-mail address: tevfi[email protected] 0377-2217/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.ejor.2013.10.064

and Zeltyn (2007, chap.), modeling the uncertainty of the system processes is of concern to call center practitioners with possible implications in the estimation of operating characteristics, staffing and service levels. In this paper, we focus specifically on the modeling and the inference of the service process which was not previously studied in detail in the literature. Our study is motivated by the observation that different customer profiles might require different agent skill sets which can create additional sources of uncertainty in the behavior of service distributions. For instance, the behavior exhibited by agents that are capable of answering quick stock exchange related questions would be quite different than those that are capable of handling day to day business type customers. Being able to model such heterogeneous behavior of the service process will be crucial as part of a queuing system that is capable of capturing general service behavior. One such queuing system is denoted by M t =G=st þ G where the arrival times are modeled by a non-homogeneous Poisson process with the arrival rate function kt , the service times and the abandonment times are distributed according to a general density and there are st agents for each t time interval as discussed in Whitt (2007). Furthermore, being able to quickly update service uncertainty as new service information becomes available for adequately implementing dynamic staffing policies (i.e. updating st ) is another area where some of our proposed models can provide some guidance on. The workload of a queuing system is commonly used as a measure in call center research and requires the estimation of the service process in addition to the arrival process. There are

710

T. Aktekin / European Journal of Operational Research 234 (2014) 709–719

many studies that focus solely on the forecasting of the arrival process such as Weinberg, D.Brown, and Stroud (2007), Shen and Huang (2008) and Aktekin and Soyer (2011). Thus, we focus our attention in the modeling of the service process which in turn coupled with the arrival models can be used to estimate the offered workloadpffiffiffi (R) to determine the optimal number of agents via s ¼ R þ b R as discussed by Whitt (2007). One of the commonly made assumptions in call center queuing models is that the service process follows an exponential distribution with fixed rate. However, studies such as Brown et al. (2005) and Aldor-Noiman et al. (2009) argue that the service process roughly follows a lognormal distribution. In this paper, we investigate the potential behavior of the service process using mixture models that are capable of capturing non-standard behavior such as multi-modality, skewness and excess kurtosis motivated by real call center data for three distinct customer profiles; service provided to stock exchange customers (NE), service provided to new/potential customers (NW) that are interested in acquiring services and the service provided for regular day to day business customers (PS). Another motivation for using the mixture models is the potential unobserved heterogeneity of the call center customers and their servers due to the lack of available characteristics that makes modeling the service process using traditional methods unrealistic. We propose three classes of inference strategies in modeling the service process. The first two are parametric approaches using mixtures with known and unknown number of components. The third one is a semi-parametric approach for mixtures of Dirichlet processes. We consider their Bayesian posterior inference using Markov chain Monte Carlo (MCMC) methods such as the Gibbs sampler and the reversible jump algorithms. In addition, we formally test and investigate the existence of lognormal service density behavior against mixture behavior. To the best of our knowledge, such Bayesian mixture models have not been considered in the analysis of call center service processes and represent the novelty of our approach in the call center literature. The stream of operational research in call centers mainly focuses on optimal staffing (capacity sizing) and queuing theoretic implications. From a staffing perspective, Aksin and Harker (2003) study the capacity sizing problem using heuristics for different types of calls where resource sharing is allowed. Garnett, Mandelbaum, and Reiman (2002) investigate the effects of ignoring the abandonment behavior of costumers on staffing and Aguir, Aksin, Karaesmen, and Dallery (2008) study the effects of ignoring the possibility of retrials on staffing. Batta, Berman, and Wang (2007) consider the problem of balancing staffing and switching costs in a call center with several customer types with time-dependent service. Avramidis, Chan, Gendreau, L’Ecuyer, and Pisacane (2010) propose simulation-based algorithms for multi-skill call center staffing. From a queuing perspective, Hoffman and Harris (1986) introduce a model to estimate the caller retrial rate for a call center. So and Tang (1996) develop a model that analyzes the benefits of managing the capacity dynamically for systems with both invisible and visible queue sizes and provide guidance for computing relevant performance measures. Aksin and Harker (2000) consider methods for computing queuing performance measures such as the blocking probability of a call and the abandonment probability of customers that are not serviced immediately using multidimensional convolutions and Monte Carlo summation techniques. Duder and Rosenwein (2001) provide rules of thumb to quantify the cost of abandoning and compare it against those implied by standard queuing results. Furthermore, Koole and Mandelbaum (2002) and Aksin, Armony, and Mehrotra (2007) provide thorough surveys of recent call center research from queuing and operations points of view. The Bayesian point of view has received some recent attention in the call center literature with main focus on the modeling of

the arrival process and forecasting. Weinberg et al. (2007) introduce a non-homogeneous Poisson model (using the square-root transformation technique) with correlated intra-day and interday arrival rates. Soyer and Tarimcilar (2008) consider the analysis of the arrival process of a consumer electronics producer using a modulated Poisson process which is a function of both advertising specific characteristics and of time. Landon, Ruggeri, Tarimcilar, and Soyer (2010) discuss latent sources in call center arrival data via Bayesian hidden Markov models. Aktekin and Soyer (2011) introduce discrete time Bayesian state space models with Poisson measurements in forecasting inter-day and intra-day call center arrivals. Aktekin and Soyer (2012b) consider the Bayesian inference of the M=M=s þ M, the Erlang-A queuing model. Recently, Aktekin and Soyer (2012) model the call center abandonment using Bayesian duration models and discuss differences in customer behavior. The following is a quick summary of our paper. In Section 2, we introduce three mixture modeling strategies for the service process and consider their Bayesian inference. In Section 3, we use actual call center data to show how the proposed models can be implemented. In doing so, we discuss service differences between customer profiles and the implications of Bayesian inference. In addition, we show how the proposed models can be used in conjunction with the existing theory for staffing. In Section 4, we conclude with the summary of findings and directions for future research. 2. Mixture modeling of the service process and its bayesian inference The proposed models are motivated by the observation that different customer profiles might require different agent skill sets. There are three distinct service profiles in our dataset: NE, NW and PS. Fig. 1 shows the histograms of the three service profiles in the log-scale from which it can be observed that each profile exhibits distinct service behavior in addition to potential multimodality, skewness and excess kurtosis density behavior. To model such non-standard density behavior, a natural set of models to consider are the mixtures. In what follows, we introduce three classes of mixture models, parametric with known number of components, parametric with unknown number of components and semi-parametric and discuss their Bayesian analysis. 2.1. Parametric family of models: mixtures with known number of components The first class of mixture models we consider are referred to as finite mixture models with known number of components. Some of the attractive features of finite mixture models are their ability to capture the unobserved heterogeneity introduced by a latent categorical variable and their ability to estimate several density shapes. Diebolt and Robert (1994) points out that mixture models can be thought of as alternatives to non-parametric models and they are less restrictive as opposed to standard parametric distributions. If we let t represent the service completion time of a given service profile. A finite mixture model for t can be defined as

ðtjp1 ; . . . ; pk ; h1 ; . . . ; hk Þ 

k X

pj f ðtjhj Þ;

ð2:1Þ

j¼1

P where kj¼1 pj ¼ 1 (referred to as mixing weights), k > 0 is any integer, and f ðtjhj Þ for j ¼ 1; . . . ; k represent the densities of k different categories (also referred to as the components of the mixture). Since Brown et al. (2005) argues that service process in call centers roughly follows a lognormal distribution, we use normal mixture

711

T. Aktekin / European Journal of Operational Research 234 (2014) 709–719 NE Customers 200

600

Frequency

400

100

Frequency

40 30

−4

−2

0

2

4

6

8

0

0

0

10

200

50

20

Frequency

150

50

800

60

PS Customers

1000

NW Customers

−4

−2

log−service times

0

2

4

6

8

log−service times

−4

−2

0

2

4

6

8

log−service times

Fig. 1. Histograms for different types of service profiles.

components for the f ðtjhj Þs after transforming the data into the logscale and thus write the full mixture model as

ðtjp1 ; . . . ; pk ; l1 ; . . . ; lk ; r21 ; . . . ; r2k Þ 

k X pj f ðtjlj ; r2j Þ;

ð2:2Þ

j¼1

where f ðtjlj ; r2j Þ represents the normal density for the jth component with parameters lj and r2j . Depending on the context the mixture is being used for, the components of the mixture may or may not have any physical meaning. For the service process modeling, the latent categories can represent the inherent easiness of a job in a given profile or the competency level of a given call center agent. For instance, some issues might be easily handled in a short period of time and some will require additional time for a given service profile. Prior to handling the issue, it is not possible to categorize a given call as quick to handle vs. time consuming. We believe that such behavior is observed in Fig. 1 for all three profiles, which justifies the use of mixture models that can capture such latent heterogeneity in data in a straightforward manner. The maximum likelihood (ML) estimation of mixture models is not straightforward in most cases and it may not exist in others as discussed in Diebolt and Robert (1994). Recently, Fruhwirth-Schnatter (2006) points out that for certain cases the ML estimates may not exist, may be degenerate or may lie outside of the parameter space. In addition, Fruhwirth-Schnatter (2006) compare the performance of the ML estimates against Bayesian estimates with simulations and find that deviations from asymptotic normality are easily accounted for using the Bayesian approach, a feature ML estimates fails to account for. Thus, taking a Bayesian approach for model inference is the natural way to approach the problem. It is possible to develop a Markov chain Monte Carlo algorithm to estimate the  model parameters using a data augmentation step. Let H ¼ p1 ; . . . ; pk ; l1 ; . . . ; lk ; r21 ; . . . ; r2k . Furthermore assume that zij is a latent variable indicating which component the ith Pk observation t i belongs to, that is, zij 2 f0; 1g and j¼1 zij ¼ 1. z ¼ fz1 ; . . . ; zn g is referred to as the incomplete data. Given z and H, the data augmented density of the ith observation from (2.2) can be rewritten as

f ðt i jH; zÞ ¼

k Y

z

z

pj ij f ðti jlj ; r2j Þ ij ;

for i ¼ 1; . . . n;

ð2:3Þ

j¼1

where n represents the number of service times for a given profile. Thus, the likelihood term with the missing data structure can be written as

LðH; D; zÞ ¼

n Y k Y i¼1 j¼1

z

z

pj ij f ðti jlj ; r2j Þ ij :

ð2:4Þ

We assume that the mixing weights, ðp1 ; . . . ; pk Þ follow a Dirichlet prior with parameters ða1 ; . . . ; ak Þ. In addition, we assume the following flat but proper priors, lj  Nða; bÞ and r2  Gammaðc; dÞ j for all j. We discuss the values of the prior distribution parameters in our numerical example. The posterior samples for the model parameters can be obtained by implementing a Gibbs sampler with the data augmentation step which is standard in the literature. In order to preserve space, we will omit the details of its MCMC discussion and refer the reader to the algorithm discussed in Diebolt and Robert (1994). Additionally, the mixture models with known components can be extended in a straight forward manner to incorporate covariate information if it were to be available in service data. In our numerical example, we make use of such an extension by introducing the time of the day as a covariate and discuss its implications on the service process inference. 2.2. Parametric family of models: mixtures with unknown number of components One of the unattractive features of the mixture models discussed in the previous section is that they require that the number of latent categories are known with certainty. An alternative method discussed by Richardson and Green (1997) is the analysis of mixtures with unknown number of components. Such a method treats k, the number of components of the mixture, as an uncertain quantity and puts a prior and carries out its inference via its probability distribution. Since the estimation is not straightforward, from computational and practical perspectives, mixtures with known k can be preferred over mixtures with unknown k, especially if the goal is to update the service distribution in a quick and efficient manner for dynamic staffing. In our analysis of service times data, we consider both cases. Following Richardson and Green (1997), the proposed structure will be very similar to the one proposed in the previous section. The mixture model is given by (2.2) with the exception of k now being treated as a random variable. Thus, we write the full mixture model as ðtjk; p1 ; . . . ; pk ; l1 ; . . . ; lk ; r21 ; . . . ; r2k Þ. In addition, the following priors are assumed for the model parameters, starting with k as

k  Uð0; kmax Þ;

ð2:5Þ

where U represents the uniform distribution and kmax the maximum value k can take. In our numerical examples, we set kmax to 30. Once again the mixing weights, ðp1 ; . . . ; pk Þ follow a Dirichlet prior with parameters ða1 ; . . . ; ak Þ as before. The priors for the normal components will be given by

lj  Nða; bÞ and r2  Gammaða; bÞ; j

ð2:6Þ

712

T. Aktekin / European Journal of Operational Research 234 (2014) 709–719

where the values of a and b are set to values coming from the data and b is assumed to have a hyper-prior structure. a is set to the midpoint of the service times for a given profile as a ¼ 0:5fmaxðtÞþ minðtÞg. b is set to multiples of the squared service times length 2 as b ¼ C fmaxðtÞ  minðtÞg for C ¼ 2; 5; 10 in our numerical example. The hyper-prior for b is given by

b  Gammaðg; hÞ;

ð2:7Þ 2

where h is set to a small multiple of fmaxðtÞ  minðtÞg and g is selected such that a > 1 > g which reflects the similarity of the r2j s. We discuss the implications of these hyper-parameters in our numerical example. In estimating the joint posterior distribution of the model parameters, we implement the reversible jump Markov chain Monte Carlo (RJMCMC) algorithm as discussed in Richardson and Green (1997). The idea is based on a random sweep Metropolis– Hastings algorithm implemented on the state space of the model parameters, ðk; p1 ; . . . ; pk ; l1 ; . . . ; lk ; r21 ; . . . ; r2k Þ. The random sweeps move from a model with a given state to another model with another state for a given acceptance probability. The dimension of the state space is adjusted according to the proposed value of k. In summary, the RJMCMC algorithm updates the model parameters and the allocation of z first, then considers splitting one mixture component into two, or combining two into one (by changing k by 1) which results with the birth or death of an empty component. To implement the last few steps, one needs to generate a three variable random vector denoted by fu1 ; u2 ; u3 g. The RJMCMC uses the following Beta distributions

u1  Betað2; 2Þ; u2  Betað2; 2Þ and u3  Betað1; 1Þ:

ð2:8Þ

We adopt the same above structure in implementing the RJMCC in our analysis of different service profiles; see Richardson and Green (1997) for a detailed discussion of the algorithm. 2.3. Semi-parametric family of models: mixtures of Dirichlet processes The last class of models we consider are semi-parametric (or non-parametric) models using normal mixtures of Dirichlet priors (MDP). Two of the most commonly applied areas of the MDP models are density estimation and clustering. We make use of such models to estimate service densities of different service profiles and also to identify the latent classes. Mixture models with countably infinite number of components can be modeled using a Dirichlet process as the prior for the mixing proportions which is referred to as the mixtures of Dirichlet priors (MDP). Using the MDP setup avoids the requirement to determine the correct number of components in finite mixture models. However, the finite mixture models can still be preferred over the MDP models due to their simpler and less expensive computational requirements with a more parsimonious set of model parameters. We follow and summarize the MDP model as discussed by Escobar and West (1995). The model for the service times of a given profile, ti , can be summarized as follows

ðti jli ; Ri Þ  Nðli ; Ri Þ;

ð2:9Þ

for i ¼ 1; . . . ; n. In (2.9), each ti follows a normal distribution with parameters li and Ri . The model has as many parameters as the dimension of the data (n) as opposed to the mixture models presented by (2.1) where the number of parameters are way smaller. In addition, the parameter pair, ðli ; Ri Þ is assumed to follow a general distribution (G) as

ðli ; Ri Þ  G;

ð2:10Þ

Table 1 Summary of descriptive statistics for the service profiles (in seconds). Statistics

NE

NW

PS

Min 25th Median Mean 75th Max Size (n)

0.00 102.80 199.00 323.40 396.00 4171.00 1476

0.00 35.00 61.00 106.50 120.00 2498.00 4017

0.00 69.00 123.00 192.00 231.00 5300.00 20,374

where a is referred to as the concentration parameter with the following prior

ðaja; bÞ  Gammaða; bÞ;

ð2:12Þ

and the baseline distribution G0 is a normal-inverted Wishart as

G0  Nðljm1 ; ð1=k0 ÞRÞIWðm1 ; w1 Þ:

ð2:13Þ

Additional hyper-priors of the model can be summarized as

ðm1 jm2 ; s2 Þ  Nðm2 ; s2 Þ;

ð2:14Þ

ðk0 js1 ; s2 Þ  Gammaðs1 =2; s2 =2Þ;

ð2:15Þ

ðw1 jm2 ; w2 Þ  IWðm2 ; w2 Þ;

ð2:16Þ

where IW represents an inverted Wishart distribution (IW q ) with parameterization with mean 1=ðw2 Þðm2  q  1Þ. For the MCMC estimation of the model parameters, we make use of Algorithm 8 as discussed in Neal (2000) which uses Metropolis–Hastings updates with the conditional prior as the proposal distribution and a partial form of Gibbs sampling. Another strategy for inducing regularity and sparsity for the number of components is to put additional priors on the hyper parameters of the a parameter in (2.12). Since the purpose of this study is to shed light on service processes for call center practitioners and to provide practical guidelines for the inputs of call center queuing models, more complex extensions can be left for future studies. We comment on the hyper-prior value selection in our numerical example. 3. Analysis of different service profile times in call centers In order to assess the differences between different service profiles and illustrate the implementation of the proposed models, we use real call center data from an anonymous bank operation.1 The bank is operational 7 days of the week from 7AM in the morning to midnight and serves several different customer profiles. For illustrative purposes, we use service information for the month of January for three different profiles: NE, NW and PS. The histograms for the service times in the log-scale of these three profiles were given in Fig. 1. In addition, Table 1 shows a summary of descriptive statistics for the three service profiles (in seconds) which shows further evidence of distinct service profiles. 3.1. Prior selection and MCMC implementation implications In implementing the three proposed classes of models discussed in Section 2, we used the following priors. For the finite mixture models with known number of components, to avoid a possible label switching problem in the Gibbs sampler, we imposed constraints on the parameters and used slightly stronger priors. More specifically, for the mean parameters we introduced a hierarchy in the form of l1 < l2 < . . . < lk . To do so, we added positive normal priors on the difference of two consecutive means

with G following a Dirichlet prior as

ðGja; G0 Þ  DPðaG0 Þ;

ð2:11Þ

1 Available at Technion, Israel Institute of Technology, http://iew3.technion.ac.il/ serveng/callcenterdata/.

713

median 97.5%

median 97.5%

1.01

1.01 40000

20000

25000

1.4

40000

20000

25000

median 97.5%

shrink factor

1.3 1.2

shrink factor

8 6

30000

35000

40000

last iteration in chain

median 97.5%

1.000

median 97.5%

10

35000

1.0

0.992

2

0.994

1.1

4

shrink factor

30000

last iteration in chain

1.002

35000

0.998

30000

last iteration in chain

0.996

25000

1.00

1.00

1.000 20000

1.02

shrink factor

1.03

shrink factor

1.02

1.010 1.005

shrink factor

1.03

1.04

median 97.5%

1.04

1.015

T. Aktekin / European Journal of Operational Research 234 (2014) 709–719

20000

25000

30000

35000

40000

20000

25000

last iteration in chain

30000

35000

40000

20000

25000

last iteration in chain

Fig. 2. Brooks and Gelman plots for some of the pi s and

in ðl2  l1 Þ; ðl3  l2 Þ; . . . ðlk  lk1 Þ  N þ ð4; 10Þ with l1  Nð4; 10Þ. Having the constraint structure with slightly as

stronger priors helped us to avoid a possible label switching problem. We used flat priors for the precision parameters as in r2  Gammað0:001; 0:001Þ for all j. The Dirichlet priors for the j mixture weights were set at ðp1 ; . . . ; pk Þ  DRð1; . . . ; 1Þ which is common practice in Bayesian mixture models. For the finite mixture models with the unknown number of components, the prior for the number of components was set at k  Uð0; 30Þ, where the upper bound was set to 30 to be cautious. Furthermore, we followed the suggested practical values of Richardson and Green (1997) as discussed in Section 2.2. For the NE service profile, lj  Nð1:86; 167:4Þ; r2  Gammað2; bÞ; b  Gammað0:2; j 0:0059Þ, for the NW service profile, lj  Nð1:60; 154:4Þ; r2  Gamma ð2; bÞ; b  Gammað0:2; 0:0064Þ and for the PS service j profile, lj  Nð1:98; 173:7Þ; r2  Gammað2; bÞ; b  Gammað0:2; j 0:0057Þ were used. The mixture weights were assumed to be ðp1 ; . . . ; pk Þ  DRð1; . . . ; 1Þ for all three profiles. For the MDP model, a  Gð2; 1Þ; m1  Nð0; 1000Þ; k0  Gammað0:5; 50Þ; w1  IWð4; 0:5Þ for all three profiles. We note that the posterior results were not sensitive to the choice of priors. For assessing the sensitivity for the model with unknown number of components, we increased the precision of the priors for lj s and the shape parameter of the prior of b by a multiple of 2,5 and 10 and did not observe any major changes in the posterior inference. We will omit a detailed discussion about MCMC convergence for the models considered in the study to preserve space in the narrative. We only present the results for the finite mixture models of the NE type customers and note that similar results were obtained for all other parameters. In implementing the MCMC algorithms, three parallel chains with reasonably separated starting initials were used. The chains were ran for 20,000 observations for the finite mixture models with both known and unknown components

30000

35000

40000

last iteration in chain

li s.

and 30,000 for the MDP model as the burn-in period. After the burn-in period, 20,000 samples were collected for posterior inference with thinning intervals ranging from 5 to 10 units for all three models. A formal way of assessing convergence is due to the Brooks and Gelman plots and the shrink factor; see Brooks and Gelman (1998). If the shrink factor is around 1 then convergence is said to have been attained. The Brooks and Gelman plots are shown in Fig. 2 for the 20,000 collected observations after a burn-in period of 20,000 iterations where the shrink factor approaches 1 as the number of iterations increases for the finite mixture models of the NE type customers. The estimated shrink factors were between 1 and 1.03 for all parameters with upper confidence intervals ranging between 1.01 and 1.10. Thus, we can formally conclude that we did not encounter any convergence issues for any of the models considered in the sequel. 3.2. Posterior inference of service profiles First we discuss the inference for mixtures with known number of components. We considered mixture models with two to seven components for each service profile. Most higher order mixtures (more than five) did not provide reasonable fits to data. A practical approach for determining the number of components is to consider a measure of fit such as the deviance information criteria (DIC). Celeux, Forbes, Robert, and Titterington (2006) propose the use of a DIC measure for mixture models, referred to as DIC-II which makes use of the predictive density. Table 2 shows the estimates of DIC-II for various profiles for different values of the number of components, where a small value indicates a more adequate fit to data (bold values indicate the best fit models). Reasonable number of components were determined to be 4 for the NE service, 5 for the NW service (the DIC-II starts to drop after k ¼ 5) and 4 for the PS service. The posterior statistics for these

714

T. Aktekin / European Journal of Operational Research 234 (2014) 709–719

manner despite the use of practical measures such as the DIC or the marginal likelihood, none of which puts a probability distribution on k. Mixture models with unknown number of components naturally handle such a case but are computationally expensive. From a practical point of view, mixtures with fixed number of components can be preferred due to their computational simplicity to model the service times as part of a call center queuing system (especially for dynamic staffing purposes). Next, we discuss the posterior inference for the mixtures with unknown k. Table 6 and Fig. 3 show the posterior summary statistics and posterior distributions of the service times as a result of the RJMCMC algorithm. Fig. 4 shows the trace plots for the unknown ks. The results are close but not exactly the same with what was found with the previous mixture models in an informal manner. Even though the prior for k had a maximum of 30, none of the models assigned probabilities greater than zero to values of k more than six. The posterior results show strong support for k ¼ 3 for the NE customers with probability 0.77, k ¼ 4 for the NW customers with probability 0.50 and k = 3–4 for the PS customers with roughly equal probabilities 0.50 and 0.48 that are shown in bold in Table 6 (this can also be observed from the trace plot of the PS type service profile at the right panel of Fig. 4 where the chain alternates between the values of 3 and 4 with almost equal probabilities). These posterior probabilities allows us to quantify and hypothesize our uncertainty about the existence of latent categories for different service profiles in a natural and formal manner. Another important finding about the service behavior of different profiles is that there is strong evidence of mixture existence unlike the common belief that they are lognormally distributed for any of the service profiles. In fact, when k ¼ 1, the proposed mixture models are standard lognormal distributions and we can quantify the probability that they are lognormally distributed by computing the posterior probability that k ¼ 1. Table 6 shows that Pðk ¼ 1Þ was approximately equal to 0 for all three service profiles. Thus, we can formally argue that the common assumption about the lognormal service times in call center operations was violated for this specific dataset. It would be interesting to observe the behavior of service data from other call center operations to gain more general insights about the call center service times. Fig. 5 shows the posterior predictive fits of the mixture models on the

Table 2 DIC-II: mixture models with fixed k. k

2

3

4

5

NE NW PS

4709.5 11938.4 58230.6

4667.1 11872.3 57358.1

4620.8 11814.2 57027.2

4621.5 11796.7 57471.7

three profiles are given in Tables 3–5. The posterior statistics for all three suggests that there are at least four latent categories for each service profile. These can informally be thought as service jobs requiring small, medium, high and very high amount of times for each profile. The existence of such latent categories might be due to the level of expertise of the agents or the nature of the service request, none of which can be known prior to completing the service job. Such an observation further justifies the use of mixture models with latent categories. An observation that can be drawn for all service profiles is that there is at least one latent category (with a small probability) with a negative mean (in the log-scale an estimate very close to zero in service times). The existence of such a latent category can also be observed from the left tail of the histograms given in Fig. 1. They might represent jobs that were handled extremely quickly or jobs which resulted with a connection problem. More specifically, the mean service times (in seconds) for the latent categories identified for the NE service profile were computed to be 0.01, 2, 6, and 221, for the NW profile 0.01, 0.33, 44.7, 52.9 and 109.9 and for the PS profile were 0.01, 5.4, 79.04 and 181.27. In addition, one interesting observation regarding some of the precision parameters from Tables 3–5, is that some r2 j s have large posterior means. For instance, the posterior mean of the second precision for the NE profile is 780.80, implying a very small variation in its posterior estimate. This can be attributed to the small number of observations belonging to that specific category which could explain the reason for a small variation estimate, as also evidenced by the small posterior mean estimate of p2 (0.01). It is possible to make similar arguments for other service profiles with the large precisions. One of the main challenges of finite mixtures with known components is to determine the value of k in a more formal

Table 3 Posterior statistics for the NE service profile: mixture models with fixed k. k¼4

p1

p2

p3

p4

l1

l2

l3

l4

r2 1

r2 2

r2 3

r2 4

2.5% Mean 97.5%

0.00 0.00 0.00

0.00 0.01 0.01

0.03 0.05 0.07

0.91 0.93 0.94

5.20 4.60 4.20

0.67 0.69 0.71

1.50 1.80 2.30

5.34 5.40 5.42

2.99 49.80 253

311.70 780.80 1442

0.60 1.40 2.40

1.00 1.10 1.20

Table 4 Posterior statistics for the NW service profile: mixture models with fixed k. k¼5

p1

p2

p3

p4

p5

l1

l2

l3

l4

l5

r2 1

r2 1

r2 2

r2 3

r2 4

2.5% Mean 97.5%

0.00 0.00 0.00

0.00 0.00 0.00

0.16 0.22 0.27

0.46 0.54 0.64

0.13 0.23 0.31

4.64 4.60 4.57

1.40 1.10 8.40

3.80 3.87 3.91

3.87 3.97 4.06

4.50 4.70 4.90

117 410 888

0.00 1.09 4.01

4.80 6.60 9.00

0.50 0.60 0.70

1.60 2.20 3.10

Table 5 Posterior statistics for the PS service profile: mixture models with fixed k. k¼4

p1

p2

p3

p4

l1

l2

l3

l4

r2 1

r2 2

r2 3

r2 4

2.5% Mean 97.5%

0.00 0.00 0.00

0.04 0.04 0.04

0.24 0.33 0.40

0.55 0.61 0.70

4.61 4.60 4.59

1.62 1.70 1.78

4.31 4.37 4.41

5.10 5.20 5.30

108 183 275

1.20 1.50 1.70

3.10 3.50 4.30

1.40 1.50 1.60

715

T. Aktekin / European Journal of Operational Research 234 (2014) 709–719

with no concern in obtaining a parsimonious model. For instance, the model with the unknown number of components had produced roughly posterior means of 3 for the NE profile, 4 for the NW profile and 3–4 for the PS profile as shown in Fig. 3. Thus, trying to interpret the meanings for the mixture components will not be a useful exercise unlike the finite mixture models. Fig. 11 shows the posterior predictive density fits for the service profiles which once again shows support in favor of mixture behavior of the call center service times for different profiles. We also note here that, as an attractive feature of the Bayesian inference in MDP models, it would be straight forward to show the credibility intervals of the predictive density fits as was done for the all model parameters. To summarize, based on the inference obtained from all three models, we can fairly conclude that none of the three distinct service profiles in our call center data follow non-mixture behavior unlike what has been claimed in the call center queuing literature. All mixture models support the existence of several unknown components, namely the existence of heterogeneous service behavior of all three profiles. From an inference perspective, the posterior summaries from Tables 3–5 and 7 show distinct service behavior for the three service profiles under consideration, confirming our suspicions drawn from the initially presented summary statistics of Table 1. We also note that the predictive density fits from Figs. 5 and 11 for the models considered show similar behavior for the three service profiles with differences in the component specific posterior inferences.

Table 6 Posterior probabilities of k for different service profiles. Profile

Pðk ¼ 1Þ

Pðk ¼ 2Þ

Pðk ¼ 3Þ

Pðk ¼ 4Þ

Pðk ¼ 5Þ

Pðk ¼ 6Þ

NE NW PS

0.000 0.004 0.000

0.018 0.181 0.009

0.775 0.142 0.505

0.160 0.508 0.484

0.035 0.147 0.001

0.007 0.016 0.000

0.5

The main focus of the call center related research has been on staffing and queuing theoretic implications. Since our modeling of service times using mixture models are not Markovian, one

2

4

6

8

0.4 0.0

0.0

0.0

0.2

0.1

0.2

0.4

0.2

0.3

Density

0.8 0.6

Density

0.6 0.4

Density

3.3. Call center staffing implications

1.0

0.8

actual service times of the three service profiles, which is further support in favor of the mixture behavior of the call center service times. The last class of models considered are semi-parametric (or non-parametric) mixtures of Dirichlet processes (MDP) that are commonly used for density estimation. In estimating the MDP models, the DPpackage in R was used as discussed in Jara, Hanson, Quintana, Müller, and Rosner (2011). Unlike the finite mixture models with known components, they are computationally expensive since there are at least as many parameters as the dimension of the data. There were 1476 observations for the NE profile, 4017 for the NW and 20,374 for the PS. To ease the burden of computing, we worked with a random sample of 5000 observations for the PS data and used the same number of observations for NE and NW data. We used couple of different samples of 5000 for the PS data and the results were identical. A summary of the posterior means of the MDP model is shown in Fig. 7 and the posterior distributions of model parameters with credibility intervals are given by Figs. 7–10. The variable ‘‘No. cluster’’ is similar to k of the finite mixtures whose posterior densities are shown in Fig. 6. The estimated number of components were roughly 5 for the NE profile, 11 for the NW profile and 8 for the PS profile. It is not surprising that the semi-parametric model estimates more components than the parametric mixture models since the goal is density estimation

1

2

3

4

K

5

6

7

2

3

4

K

5

6

K

4

K

4

K 30000

35000

40000

45000

Iteration

50000

55000

2

1

2

2

3

3

4

K

6

5

5

6

8

7

6

Fig. 3. Posterior distributions for k for NE (left), NW (middle) and PS (right).

30000

35000

40000

45000

50000

55000

Iteration

Fig. 4. Trace plots for k for NE (left), NW (middle) and PS (right).

30000

35000

40000

45000

Iteration

50000

55000

T. Aktekin / European Journal of Operational Research 234 (2014) 709–719

0.4

716

0.4

0.4

Actual Data Predictive Fit

Actual Data Predictive Fit

−4

−2

0

2

4

6

0.3 0.2 0.0

0.0

0.0

0.1

0.1

0.1

0.2

0.2

0.3

0.3

Actual Data Predictive Fit

8

−4

−2

0

2

4

6

8

−4

−2

0

2

4

6

8

Fig. 5. Posterior predictive density fits of NE (left), NW (middle) and PS (right) for the mixture model with unknown k.

Density of ncluster

5

10

15

0.15 0.10 0.00

0.0

0.00

0.02

0.1

0.05

0.04

0.06

probability

probability

0.08

0.4 0.3 0.2

probability

Density of ncluster

0.10

Density of ncluster

20

10

20

values

30

40

50

5

10

15

values

20

25

30

35

values

Fig. 6. Posterior (k) cluster distributions of NE (left), NW (middle) and PS (right) for the MDP model.

Table 7 Posterior means for the MDP model of NE, NW and PS. Profile

m1

k0

w1

No. cluster

a

NE NW PS

1.202 2.136 2.926

0.024 0.028 0.032

0.594 1.316 0.696

5.330 11.200 8.420

0.802 1.532 1.152

research and requires the estimation of the service process in addition to the arrival process. According to Aldor-Noiman et al. (2009), the offered workload (R) can be estimated for a given time interval with Markovian arrival distribution and a general service distribution by

b ¼ ^k^t; R

k is the expected number of arrivals for a given time interval where ^ and ^t is the expected service time for a given time interval. The estimation of ^ k has been studied in detail in the call center literature. Our proposed models can assist in estimating ^t in (3.1) for different profiles which can in turn be used with the square-root staffing rule to determine optimal staffing via

ð3:2Þ

0.23

0.27 1.20

values

0.06

0.12

density

0.00

0.07 0.00 −6.64

Density of m1−PS_Service

0.17

0.21 0.14

density

0.07 0.04 0.00

density

qffiffiffiffi b; b þb R ^s ¼ R

Density of m1−NW_Service

0.11

0.15

needs a queuing model that allows general service times. One such queuing system is the M t =G=st þ G queue where the arrival times are modeled by a non-homogeneous Poisson process with the arrival rate function kt , the service times and the abandonment times are distributed according to a general density and there are st agents for each t time interval as discussed in Whitt (2007). For practical staffing estimates, a measure referred to as the offered workload of a queuing system is commonly used in call center

Density of m1−NE_Service

ð3:1Þ

−1.10

2.14

5.02

−1.59

values

Fig. 7. Posterior distributions of m1 ; NE (left), NW (middle) and PS (right) for the MDP model.

2.93

6.91

values

717

T. Aktekin / European Journal of Operational Research 234 (2014) 709–719

24.52 6.13

12.26

density

21.16 14.11

density 0.02

0.06

0.00

0.00

0.00 0.00

Density of k0

18.39

28.22

Density of k0

7.05

7.46

14.92

density

22.38

29.84

Density of k0

0.00

0.03

0.07

values

0.00

0.03

0.08

values

values

Fig. 8. Posterior distributions of k0 ; NE (left), NW (middle) and PS (right) for the MDP model.

1.79 0.90

density

0.62 0.42

density 0.60

1.16

0.45 0.00

0.00

0.50 0.00 0.17

0.36

values

Density of psi1−PS_Service

1.35

0.83

Density of psi1−NW_Service

0.21

0.99

density

1.49

1.99

Density of psi1−NE_Service

1.32

2.69

0.24

values

0.70

1.32

values

Fig. 9. Posterior distributions of w1 ; NE (left), NW (middle) and PS (right) for the MDP model.

0.80

1.71

0.62 0.41

density

0.21 0.00

0.16 0.00 0.11

Density of alpha

0.83

0.63 0.32

density

0.82 0.55 0.00

0.27

density

Density of alpha

0.48

1.09

Density of alpha

0.25

1.54

0.19

3.29

values

1.15

2.36

values

values

Fig. 10. Posterior distributions of a; NE (left), NW (middle) and PS (right) for the MDP model.

Density of NW_Service

Density of PS_Service

−2

0

2

values

4

6

8

density

0.2 0.0

0.0 −4

0.4

0.4

density

0.2

0.4 0.2 0.0

density

0.6

0.6

0.6

0.8

0.8

0.8

Density of NE_Service

−4

−2

0

2

4

6

8

−4

−2

0

values

Fig. 11. Posterior predictive density fits of NE (left), NW (middle) and PS (right) for the MDP model.

2

values

4

6

8

718

T. Aktekin / European Journal of Operational Research 234 (2014) 709–719

where b represents a measure which determines the operational regime of a call center operation and its estimation is discussed in detail in Aldor-Noiman et al. (2009). For the estimation of ^t, we can evaluate the predictive service time distribution for a given profile via

pðtjDÞ ¼

Z

pðtjHÞpðHjDÞdH;

also estimated a model by mixing the btime coefficient in the mean of the normal mixture and found that they were insignificant. The above structure can easily take into account other variables that can contribute to the predictive variability portion of the service process if they were available.

ð3:3Þ 4. Conclusion and directions for future work

where H is the set of model parameters for a given mixture model of a given service profile. (3.3) can be approximated using straight Monte Carlo as

pðtjDÞ 

M 1X pðtjHðjÞ Þ; M j¼1

ð3:4Þ

where HðjÞ represents the jth generated posterior sample set of model parameters. pðtjDÞ was used to obtain the predictive fits shown in Figs. 5 and 11. Thus, the expected service time for a given profile can be obtained by numerically integrating t in the region it was defined by using

EðtjDÞ ¼ ^t ¼

Z

pðtjDÞtdt:

ð3:5Þ

For our numerical example, using the results of the finite mixture with known k, the estimates of ^t for different customer profiles would be ^tNE ¼ 207:4; ^tNW ¼ 84:2; ^t PS ¼ 153:6. These can be coupled with the estimation of the arrival counts and be used to determine optimal staffing (^s) by plugging in estimates in (3.1) and (3.2). 3.4. Analysis of the effect of time on the service profile distributions The implementation of the square-root staffing rule (3.2) is generally implemented for discrete time intervals in call center applications. Therefore, it might be useful to be able to capture possible deterministic time varying effects (if any) in the service profile distributions. Whitt (2007) points out that fundamental call center distributions such as the service can have two kinds of variability, predictable and stochastic. By adding a term which takes into account the time of the day will be an example of predictable variability, whereas the mixture structure will be an example of stochastic variability. This can be achieved in several ways for finite normal mixtures. We discuss the case where we add a covariate component to the precision parameter of the normal densities in the mixture models. Since the estimation of the finite mixture model with known k is the most straightforward one, we estimated the model by splitting a given day into discrete intervals of length 30 minutes from 7AM until midnight and assumed that r2 ¼ exp fbtime ðTIMEÞg in (2.2) with the prior btime  Nð0; 1000Þ. The posterior results are shown in Table 8. The results suggest that for the NW service profile the time of the day is less significant compared to that of the NE customers, since the credibility interval includes the hypothesized value of btime ¼ 0. Such an outcome makes intuitive sense since NE customers are looking for stock related services and the time of the day might be more important for them than it is for potential customers that have the luxury of calling back at another time. We also note that having the covariate structure on the shape parameter did not affect the inference of the other mixture parameters. We Table 8 Posterior statistics for the time effects of NE, NW and PS. btime

Mean

SD

2.5%

97.5%

NE NW PS

0.0073 0.0018 0.0154

0.0022 0.0026 0.0013

0.0028 0.0024 0.0137

0.0117 0.0079 0.0187

In this paper, we focused on the modeling and the inference of the service times in call centers which was not previously considered in detail in the literature. In general, the service times are assumed to be exponentially distributed in conventional call center queuing models. Recently, Brown et al. (2005) argued that the call center service times roughly follow a lognormal distribution. Our aim was to investigate the service times behavior of different service profiles using models that were capable of capturing nonstandard behavior such as multi-modality, skewness and excess kurtosis. In doing so, we used real call center service data and considered three distinct customer profiles; stock exchange, new and regular. Furthermore, we proposed and made use of three classes of Bayesian parametric and semi-parametric mixture models. The first two are parametric approaches using mixtures with known and unknown number of components. The third one is a semiparametric approach for mixtures of Dirichlet processes. We considered their posterior inference using Markov chain Monte Carlo (MCMC) methods such as the Gibbs sampler and the reversible jump algorithms. To the best of our knowledge, such Bayesian mixture models have not been considered in the analysis of call center service processes and represent the novelty of our approach in the call center literature. Our findings showed strong support in favor of distinct mixture behavior for all three profiles with at least four latent groups for each. Using the inference obtained for the mixture models with unknown number of components, we were able to show that the posterior probability that there was only one component (i.e. a standard lognormal distribution) was equal to zero for all three profiles. In addition, we showed how covariate effects such as the time of the day can be taken into account in estimating the call center service times and found that it was more important for stock exchange service profile than it was for the potential customers. There are many advantages of the proposed Bayesian mixture models. For instance, being able to quickly update service uncertainty as new service information becomes available to implement staffing policies is one area where the Bayesian approach would be the natural choice. In addition, modeling the inherent uncertainty of the system processes such as the service input is of concern to call center practitioners with implications in call center operating characteristic and service level estimation as emphasized by Mandelbaum and Zeltyn (2007, chap.). We also discussed how the proposed models can be used in conjunction with the existing theory in call center research to determine optimal staffing in queuing models with general service times such as the M t =G=st þ G queue. In summary, if the goal is pure density estimation then semi-parametric models can be preferred. If the goal is to efficiently and quickly update uncertainties as new service data becomes available then mixtures with known number of components can be used. For general staffing, mixtures with unknown k might be preferred, for at least determining the value of k supported most by data which can be used as a fixed value in standard mixtures as part of a dynamic staffing scheme. It is possible to further extend the current work in call center research from the Bayesian inferential point of view. One such area would be to investigate the existence of mixture behavior in other call center operations for staffing purposes and the effects of other covariates such as the expected (or announced) waiting time in the

T. Aktekin / European Journal of Operational Research 234 (2014) 709–719

queue. Also, from a queuing perspective, work using Bayesian inference is scarce with the exception of a few. It is possible to fully analyze the M t =G=st þ G queue from a Bayesian perspective that can provide models capable of exhibiting general service behavior. Additionally, estimating a time dependent and uncertain offered load from a Bayesian perspective and studying its implications in dynamic staffing would be of interest to call center practitioners. As pointed out by one of the referees, another possibility is to investigate the use of multivariate mixtures using all three profiles jointly. For the sake of offering practical guidelines for call center managers as part of the service input of a queuing model, such an extension can be left for a future study since their estimation is computationally more expensive and complex compared to models considered in this study. Acknowledgment The author would like to thank the three anonymous referees whose suggestions significantly improved the study. References Aguir, M. S., Aksin, O. Z., Karaesmen, F., & Dallery, Y. (2008). On the interaction between retrials and sizing of call centers. European Journal of Operational Research, 191(2), 398–408. Aksin, O. Z., Armony, M., & Mehrotra, V. (2007). The modern call center: A multidisciplinary perspective on operations management research. Production and Operations Management, 16(6), 665–688. Aksin, O. Z., & Harker, P. T. (2000). Computing performance measures in a multiclass multi-resource processor-shared loss system. European Journal of Operational Research, 123(1), 464–483. Aksin, O. Z., & Harker, P. T. (2003). Capacity sizing in the presence of a common shared resource: Dimensioning an inbound call center. European Journal of Operational Research, 147(3), 464–483. Aktekin, T., & Soyer, R. (2011). Call center arrival modeling: A Bayesian state-space approach. Naval Research Logistics, 58(1), 28–42. Aktekin, T., & Soyer, R. (2012a). Bayesian analysis of abandonment in call center operations. Applied Stochastic Models in Business and Industry. Advance online publication. doi: 10.1002/asmb.1949. Aktekin, T., & Soyer, R. (2012b). Bayesian analysis of queues with impatient customers: Applications to call centers. Naval Research Logistics, 59(6), 441–456. Aldor-Noiman, S., Feigin, P. D., & Mandelbaum, A. (2009). Workload forecasting for a call center: Methodology and a case study. Annals of Applied Statistics, 3(4), 1403–1447. Avramidis, A. N., Chan, W., Gendreau, M., L’Ecuyer, P., & Pisacane, O. (2010). Optimizing daily agent scheduling in a multiskill call center. European Journal of Operational Research, 200(3), 822–832.

719

Batta, R., Berman, O., & Wang, Q. (2007). Balancing staffing and switching costs in a service center with flexible servers. European Journal of Operational Research, 177(2), 924–938. Brooks, S. P., & Gelman, A. (1998). General methods for monitoring convergence of iterative simulations. Journal of Computational and Graphical Statistics, 7(4), 434–455. Brown, L., Gans, N., Mandelbaum, A., Sakov, A., Shen, H., Zeltyn, S., et al. (2005). Statistical analysis of a telephone call center: A queuing-science perspective. Journal of the American Statistical Association, 100(469), 36–50. Celeux, G., Forbes, F., Robert, C. P., & Titterington, D. M. (2006). Deviance information criteria for missing data models. Bayesian Analysis, 1(4), 651–673. Diebolt, J., & Robert, C. P. (1994). Estimation of finite mixture distributions through Bayesian sampling. Journal of the Royal Statistical Society. Series B, 56(2), 363–375. Duder, J. C., & Rosenwein, M. B. (2001). Towards zero ‘‘abandonments’’ in call center performance. European Journal of Operational Research, 135(1), 50–56. Escobar, M. D., & West, M. (1995). Bayesian density estimation and inference using mixtures. Journal of the American Statistical Association, 90(430), 577–588. Fruhwirth-Schnatter, S. (2006). Finite mixture and markov switching models. New York, NY: Springer. Garnett, O., Mandelbaum, A., & Reiman, M. (2002). Designing a call center with impatient customers. Manufacturing and Service Operations Management, 4(3), 208–227. Hoffman, K. L., & Harris, C. M. (1986). Estimation of a caller retrial rate for a telephone information system. European Journal of Operational Research, 27(2), 207–214. Jara, A., Hanson, T., Quintana, F., Müller, P., & Rosner, G. (2011). DPpackage: Bayesian semi- and nonparametric modeling in R. Journal of Statistical Software, 40(5), 1–30. Koole, G., & Mandelbaum, A. (2002). Queuing models of call centers: An introduction. Annals of Operations Research, 113(1-4), 41–59. Landon, J., Ruggeri, F., Tarimcilar, M., & Soyer, R. (2010). Modeling latent sources in call center arrival data. European Journal of Operational Research, 204(3), 597–603. Mandelbaum, A., & Zeltyn, S. (2007). Advances in service innovations. In Service engineering in action: The Palm/Erlang–A queue, with applications to call centers (pp. 17–48). New York, NY: Springer. Neal, R. M. (2000). Markov chain sampling methods for Dirichlet process mixture models. Journal of Computational and Graphical Statistics, 9(2), 249–265. Richardson, S., & Green, P. J. (1997). On bayesian analysis of mixtures with an unknown number of components. Journal of the Royal Statistical Society, Series B, 59(4), 731–792. Shen, H., & Huang, J. Z. (2008). Interday forecasting and intraday updating of call center arrivals. Manufacturing and Service Operations Management, 10(3), 391–410. So, K. C., & Tang, C. S. (1996). On managing operating capacity to reduce congestion in service systems. European Journal of Operational Research, 92(1), 83–98. Soyer, R., & Tarimcilar, M. (2008). Modeling and analysis of call center arrival data: A Bayesian approach. Management Science, 54(2), 226–278. Weinberg, J., D.Brown, L., & Stroud, J. R. (2007). Bayesian forecasting of an inhomogeneous Poisson process with applications to call center data. Journal of the American Statistical Association, 102(480), 1185–1198. Whitt, W. (2007). What you should know about queueing models to set staffing requirements in service systems. Naval Research Logistics, 54(5), 476–484.