Forecasting emergency medical service call arrival rates\thanksrefTitl
Abstract
We introduce a new method for forecasting emergency call arrival rates that combines integervalued time series models with a dynamic latent factor structure. Covariate information is captured via simple constraints on the factor loadings. We directly model the countvalued arrivals per hour, rather than using an artificial assumption of normality. This is crucial for the emergency medical service context, in which the volume of calls may be very low. Smoothing splines are used in estimating the factor levels and loadings to improve longterm forecasts. We impose time series structure at the hourly level, rather than at the daily level, capturing the finescale dependence in addition to the longterm structure.
Our analysis considers all emergency priority calls received by Toronto EMS between January 2007 and December 2008 for which an ambulance was dispatched. Empirical results demonstrate significantly reduced error in forecasting call arrival volume. To quantify the impact of reduced forecast errors, we design a queueing model simulation that approximates the dynamics of an ambulance system. The results show better performance as the forecasting method improves. This notion of quantifying the operational impact of improved statistical procedures may be of independent interest.
10.1214/10AOAS442 \volume5 \issue2B 2011 \firstpage1379 \lastpage1406
Forecasting emergency medical service call arrival rates \runtitleForecasting EMS call arrival rates \thankstextTITLSupported in part by NSF Grant CMMI0926814.
A]\fnmsDavid S. \snmMatteson\correflabel=e1]dm484@cornell.edu, A]\fnmsMathew W. \snmMcLean, A]\fnmsDawn B. \snmWoodard and B]\fnmsShane G. \snmHenderson
Ambulance planning \kwddynamic factor model \kwdnonhomogeneous Poisson process \kwdintegervalued time series \kwdsmoothing splines.
1 Introduction
Considerable attention has been paid to the problem of how to best deploy ambulances within a municipality to minimize their response times to emergency calls while keeping costs low. Sophisticated operations research models have been developed to address issues such as the optimal number of ambulances, where to place bases, and how to move ambulances in real time via systemstatus management [Swersey (1994); Goldberg (2004); Henderson (2009)]. However, methods for estimating the inputs to these models, such as travel times on road networks and call arrival rates, are ad hoc. Use of inaccurate parameter estimates in these models can result in poor deployment decisions, leading to low performance and diminished user confidence in the software. We introduce methods for estimating the demand for ambulances, that is, the total number of emergency calls per period, that are highly accurate, straightforward to implement, and have the potential to simultaneously lower operating costs while improving response times.
Current practice for forecasting call arrivals is often rudimentary. For instance, to estimate the call arrival rate in a small region over a specific time period, for example, next Monday from 8–9 a.m., simple estimators have been constructed by averaging the number of calls received in the corresponding period in four previous weeks: the immediately previous two weeks and the current and previous weeks of the previous year. Averages of so few data points can produce highly noisy estimates, with resultant cost and efficiency implications. Excessively large estimates lead to overstaffing and unnecessarily high costs, while low estimates lead to understaffing and slow response times. Setzler, Saydam and Park (2009) document an emergency medical service (EMS) agency which extends this simple moving average to twenty previous observations: the previous four weeks from the previous five years. A more formal time series approach is able to account for possible differences from week to week and allows inclusion of neighboring hours in the estimate.
We generate improved forecasts of the callarrival volume by introducing an integervalued time series model with a dynamic latent factor structure for the hourly call arrival rate. Dayofweek and weekofyear effects are included via simple constraints on the factor loadings. The factor structure allows for a significant reduction in the number of model parameters. Further, it provides a systematic approach to modeling the diurnal pattern observed in intraday counts. Smoothing splines are used in estimating the factor levels and loadings. This may introduce a small bias in some periods, but it offers a significant reduction in longhorizon outofsample forecasterror variance. This is combined with integervalued time series models to capture residual dependence and to provide adaptive shortterm forecasts. Our empirical results demonstrate significantly reduced error in forecasting hourly callarrival volume.
Few studies have focused specifically on EMS call arrival rates, and of those that have proposed methods for time series modeling, most have been based on Gaussian linear models. Even with a continuity correction, this method is highly inaccurate when the call arrival rate is low, which is typical of EMS calls at the hourly level. Further, it conflicts with the Poisson distribution assumption used in operations research methods for optimizing staffing levels. For example, Channouf et al. (2007) forecast EMS demand by modeling the daily call arrival rate as Gaussian, with fixed dayofweek, monthofyear, special day effects and fixed daymonth interactions. They also consider a Gaussian autoregressive movingaverage (ARMA) model with seasonality and special day effects. Hourly rates are later estimated either by adding hourofday effects or assigning a multinomial distribution to the hourly volumes, conditional on the daily call volume estimates.
Setzler, Saydam and Park (2009) provide a comparative study of EMS call volume predictions using an artificial neural network (ANN). They forecast at various temporal and spatial granularities with mixed results. Their approach offered a significant improvement at low spatial granularity, even at the hourly level. At a high spatial granularity, the mean square forecast error (MSFE) of their approach did not improve over simple binning methods at a temporal granularity of three hours or less.
Methods for the closely related problem of forecasting call center demand have received much more study. Bianchi, Jarrett and Choudary Hanumara (1998) and Andrews and Cunningham (1995) use ARMA models to improve forecasts for daily call volumes in a retail company call center and a telemarketing center, respectively. A dynamic harmonic regression model for hourly call center demand is shown in Tych et al. (2002) to outperform seasonal ARMA models. Their approach accounts for possible nonstationary periodicity in a time series. The major drawback common to these studies is that the integervalued observations are assumed to have a continuous distribution, which is problematic during periods with low arrival rates.
The standard industry assumption is that hourly callarrival volume has a Poisson distribution. The Palm–Khintchine theorem—stating that the superposition of a number of independent point processes is approximately Poisson—provides a theoretical basis for this assumption [see, e.g., Whitt (2002)]. Brown et al. (2005) provide a comprehensive analysis of operational data from a bank call center and thoroughly discuss validating classical queueing theory, including this theorem. Henderson (2005) states that we can expect the theorem to hold for typical EMS data because there are a large number of callers who can call at any time and each has a very low probability of doing so.
Weinberg, Brown and Stroud (2007) use Bayesian techniques to fit a nonhomogeneous Poisson process model for call arrivals to a commercial bank’s call center. This approach has the advantage that forecast distributions for the rates and counts may be easily obtained. They incorporate smoothness in the withinday pattern. They implement a variance stabilizing transformation to obtain approximate normality. This approximation is most appropriate for a Poisson process with high arrival rates, and would not be appropriate for our application in which very low counts are observed in many time periods.
Shen and Huang (2008b) apply the same variance stabilizing transformation and achieve better performance than Weinberg, Brown and Stroud (2007). They use a singular value decomposition (SVD) to reduce the number of parameters in modeling arrival rates. Their approach is used for intraday updating and forecasts up to one day ahead.
Shen and Huang (2008a) propose a dynamic factor model for 15minute call arrivals to a bank call center. They assume that call arrivals are a Cox process. A Cox process [cf. Cox and Isham (1980)] is a Poisson process with a stochastic intensity, that is, a doubly stochastic Poisson process. The factor structure reduces the number of parameters by explaining the variability in the call arrival rate with a small number of unobserved variables. Estimation proceeds by iterating between an SVD and fitting Poisson generalized linear models to successively estimate the factors and their respective loadings. The intensity functions are assumed to be serially dependent. Forecasts are made by fitting a simple autoregressive time series model to the factor score series.
We assume that the hourly EMS callarrival volume has a Poisson distribution. This allows parsimonious modeling of periods with small counts, conforms with the standard industry assumption, and avoids use of variance stabilizing transformations. We assume the intensity function is a random process and that it can be forecast using previous observations. This has an interpretation very similar to a Cox process, but is not equivalent since the random intensity is allowed to depend on not only its own history, but also on previous observations. We partition the random intensity function into stationary and nonstationary components.
Section 2 describes the general problem and our data set. Section 3 presents the proposed methodology. We consider a dynamic latent factor structure to model the nonstationary pattern in intraday call arrivals and greatly reduce the number of parameters. We include dayofweek and weekofyear covariates via simple constraints on the factor loadings of the nonstationary pattern. Smoothing splines are easily incorporated into estimation of the proposed model to impose a smooth evolution in the factor levels and loadings, leading to improved longhorizon forecast performance. We combine the factor model with stationary integervalued time series models to capture the remaining serial dependence in the intensity process. This is shown to further improve shortterm forecast performance of our approach. A simple iterative algorithm for estimating the proposed model is presented. It can be implemented largely through existing software. Section 4 assesses the performance of our approach using statistical metrics and a queueing model simulation. Section 5 gives our concluding remarks.
2 Notation and data description
We assume that over short, equallength time intervals, for example, one hour periods, the latent call arrival intensity function can be well approximated as being constant, and that all data have been aggregated in time accordingly. We suppose aggregated call arrivals follow a nonhomogeneous counting process , with discrete time index . Underlying this is a latent, realvalued, nonnegative intensity process . We further assume that conditional on , has a Poisson distribution with mean .
As shown in Figure 1, the pattern of call arrivals over the course of a typical day has a distinct shape. After quickly increasing in the late morning, it peaks in the early afternoon, then slowly falls until it troughs between 5 and 6 a.m. See Section 4 for more discussion. In our analysis, we consider an arrival process that has been repeatedly observed over a particular time span, specifically, a 24 hour day. Let
denote the sequence of call arrival counts, observed over time period , which corresponds onetoone with the th subperiod of the th day, so that . Our baseline approach is to model the arrival intensity for the distinct shape of intraday call arrivals using a small number of smooth curves.
We consider two disjoint information sets for predictive conditioning. Let denote the field generated by , and let denote any available deterministic covariate information about each observation. We incorporate calendar information such as dayofweek and weekofyear in our analysis. We define as the conditional expectation of given and . We defined this above as the mean of . In our model these coincide; however, this mean may not be the same as the conditional expectation since may depend on other unobserved random variables. Let denote the conditional mean of given only the covariates . Let
(1) 
in which is referred to as the conditional intensity inflation rate (CIIR). By construction,
The CIIR process is intended to model any remaining serial dependence in call arrival counts after accounting for available covariates. In the EMS context we hypothesize that this dependence is due to sporadic events such as inclement weather or unusual traffic patterns. Since information regarding these events may not be available or predictable in general, we argue that an approach such as ours which explicitly models the remaining serial dependence will lead to improved shortterm forecast accuracy. In Section 3 we consider a dynamic latent factor model estimated with smoothing splines for modeling , various time series models for modeling , and finally a conditional likelihood algorithm for estimating the latent intensity process via estimation of given .
The call arrival data used consists of all emergency priority calls received by Toronto EMS between January 1, 2007 and December 31, 2008 for which an ambulance was dispatched. This includes some calls not requiring lightsandsirens response, but does not include scheduled patient transfers. We include only the first call arrival time in our analysis when multiple calls are received for the same event. The data were processed to exclude calls with no reported location. These removals totaled less than 1% of the data.
Many calls resulted in multiple ambulances being dispatched. Exploratory analysis revealed that the number of ambulances deployed for a single emergency did not depend on the day of the week, the week of the year, or exhibit any serial dependence. However, such instances were slightly more prevalent in the morning hours. Our analysis of hourly ambulance demand defines an event as a call arrival if one or more ambulances are deployed.
We removed seven days from the analysis because there were large gaps, over at least two consecutive hours, in which no emergency calls were received. These days most likely resulted from malfunctions in the computeraided dispatch system which led to failures in recording calls for extended periods. Strictly speaking, it is not necessary to remove the entire days; however, we did so since it had a negligible impact on our results and it greatly simplified outofsample forecast comparisons and implementation of the simulation studies in Section 4.
Finally, we gave special consideration to holidays. We found that the intraday pattern on New Year’s Eve and Day was fundamentally different from the rest of the year and removed these days from our analysis. This finding is similar to the conclusions of Channouf et al. (2007) who found that New Year’s Day and the dates of the Calgary Stampede were the only days requiring special consideration in their methodology when applied to the city of Calgary. In practice, staffing decisions for holidays require special planning and consideration of many additional variables.
3 Modeling
Factor models provide a parsimonious representation of high dimensional data in many applied sciences, for example, econometrics [cf. Stock and Watson (2002)]. We combine a dynamic latent factor model with integervalued time series models. We include covariates via simple constraints on the factor loadings. We estimate the model using smoothing splines to impose smooth evolution in the factor levels and loadings. The factor model provides a parsimonious representation of the nonstationary pattern in intraday call arrivals, while the time series models capture the remaining serial dependence in the arrival rate process.
3.1 Dynamic latent factor model
For notational simplicity, assume consecutive observations per day are available for consecutive days with no omissions in the record. Let denote the matrix of observed counts for each day over each subperiod . Let , and let denote the corresponding latent nonstationary intensity matrix. To reduce the dimension of the intensity matrix , we introduce a factor model.
We assume that the intraday pattern of expected hourly call arrivals on the log scale can be well approximated by a linear combination of (a small number) factors or functions, denoted by for . The factors are orthogonal length vectors. The intraday arrival rate model over a particular day is given by
(2) 
Each of the factors varies as a function over the periods within a day, but they are constant from one day to the next. Daytoday changes are modeled by allowing the various factor loadings to vary across days. When is much smaller than either or , the dimensionality of the general problem is greatly reduced. In practice, must be chosen by the practitioner; we provide some discussion on choosing in Section 4.
In matrix form we have
(3) 
in which denotes the matrix of underlying factors and denotes the corresponding matrix of factor loadings, both of which are assumed to have full column rank. Although other link functions are available, the componentwise log transformation implies a multiplicative structure among the common factors and ensures a positive estimate of each hourly intensity . Since neither nor are observable, the expression (3) is not identifiable. We further require to alleviate this ambiguity and we iteratively estimate and .
3.2 Factor modeling with covariates via constraints
To further reduce the dimensionality, we impose a set of constraints on the factor loading matrix . Let denote a full rank matrix () of given constraints (we discuss later what these should be for EMS). Let denote an matrix of unconstrained factor loadings. These unconstrained loadings linearly combine to constitute the constrained factor loadings , such that . Our factor model may now be written as
A considerable reduction in dimensionality occurs when is much smaller than .
Constraints to assure identifiability are standard in factor analysis. The constraints we now consider incorporate auxiliary information about the rows and columns of the observation matrix to simplify estimation and to improve outofsample predictions. Similar constraints have been used in Takane and Hunter (2001), Tsai and Tsay (2010) and Matteson and Tsay (2011).
For example, the rows of might consist of incidence vectors for particular days of the week, or special days which might require unique loadings on the common factors. We may choose to constrain all weekdays to have identical factor loadings and similarly constrain weekend days. However, this approach is much more general than simple equality constraints, as demonstrated below.
The intraday pattern of hourly call arrivals varies from one day to the next, although the same general shape is maintained. As seen in Figure 1, different days of the week exhibit distinct patterns. We do not observe large changes from one week to the next, but there are significant changes over the course of the year. We allow loadings to slowly vary from week to week. Both of these features are incorporated into the factor loadings by specifying appropriate constraints . Let
(4) 
in which the first term corresponds to dayofweek effects and the second to smoothly varying weekofyear effects. is a matrix in which each row is an incidence vector for the dayofweek. Similarly, is a matrix in which each row is an incidence vector for the weekofyear. (We use a 53 week year since the first and last weeks may have fewer than 7 days.) The matrix contains unconstrained factor loadings for the dayofweek and is a matrix of factor loadings for the weekofyear.
3.3 Factor model estimation via smoothing splines
We assume that as the nonstationary intensity process varies over the hours of each day , it does so smoothly. If each of the common factors varies smoothly over subperiods , then the smoothness of is guaranteed for each day. Increasing the number of factors reduces possible discontinuities between the end of one day and the beginning of the next. To incorporate smoothness into the model (2), we use Generalized Additive Models (GAMs) in the estimation of the common factors . GAMs extend generalized linear models, allowing for more complicated relationships between the response and predictors by modeling some predictors nonparametrically [see, e.g., Hastie and Tibshirani (1990); Wood (2006)]. GAMs have been successfully used for countvalued data in the study of fish populations [cf. Borchers et al. (1997); Daskalov (1999)]. The factors are a smooth function of the intraday time index covariate . The loadings are defined as before. If the loadings were known covariates, equation (2) would be a varying coefficient model [cf. Hastie and Tibshirani (1993)].
There are several excellent libraries available in the statistical package R [R Development Core Team (2009)] for fitting GAMs, thus making them quite easy to implement. We used the gam function from the mgcv library [Wood (2008)] extensively. Other popular libraries include the gam package [Hastie (2009)] and the gss package [Gu (2010)]. See Wood [(2006), Section 5.6] for an introduction to GAM estimation using R.
In estimation of the model (2) via the gam function, we have used thin plate regression splines with a tendimensional basis, the Poisson family, and the loglink function. Thin plate regression splines are a low rank, isotropic smoother with many desirable properties. For example, no decisions on the placement of knots is needed. They are an optimal approximation to thin plate splines and, with the use of Lanczos iteration, they can be fit quickly even for large data sets [cf. Wood (2003)].
When the factors are treated as a fixed covariate, the factor model can again be interpreted as a varying coefficient model. Given the calendar covariates , let
in which is a piecewise constant function of the dayofweek, and is a smoothly varying function over the weekofyear. We may again proceed with estimation via the gam function in R. Dayofweek covariates are simply added to the linear predictor as indicator variables. These represent a level shift in the daily loadings on each of the factors . In our application it is appropriate to assume a smooth transition between the last week of one year and the first week of the next. To ensure this in estimation of , we use a cyclic cubic regression spline for the basis [cf. Wood (2006), Section 5.1]. Iterative estimation of , and via , for a given number of factors is discussed in Section 3.5.
We allow the degree of smoothness for the factors and the loadings function to be automatically estimated by generalized cross validation (GVC). We expect short term serial dependence in the residuals for our application. For smoothing methods in general, if autocorrelation between the residuals is ignored, automatic smoothing parameter selection may break down [see, e.g., Opsomer, Wang and Yang (2001)]. The proposed factor model may be susceptible to this if the number of days included is not sufficiently large compared to the number of smooth factors and loadings, or if the residuals are longrange dependent. We use what is referred to as a performance iteration [cf. Gu (1992)] versus an outer iteration strategy which requires repeated estimation for many trial sets of the smoothing parameters. The performance iteration strategy is much more computationally efficient for use in the proposed algorithm, but convergence is not guaranteed, in general. In particular, cycling between pairs of smoothing parameters and coefficient estimates may occur [cf. Wood (2006), Section 4.5], especially when the number of factors is large.
3.4 Adaptive forecasting with time series models
Let denote the multiplicative residual in period implied by the fitted values from a factor model estimated as described in the previous sections. Time series plots of this residual process appear stationary, but exhibit some serial dependence. In this section we consider time series models for the latent CIIR process to account for this dependence.
To investigate the nature of the serial dependence, we study the bivariate relationship between the process versus several lagged values of the process . Scatterplots reveal a roughly linear relationship. Residual autocorrelation and partial autocorrelation plots for one of the factor models fit in Section 4 are given in Figure 5(b) and (c). These quantify the strength of the linear relationship as the lag increases. It appears to persist for many periods, with an approximately geometric rate of decay as the lag increases.
To explain this serial dependence, we first consider a generalized autoregressive linear model, defined by the recursion
(6) 
To ensure positivity, we restrict and . When is constant, the resulting model for is an (IntGARCH) model [e.g., Ferland, Latour and Oraichi (2006)]. It is worth noting some properties of this model for the constant case. To ensure the stationarity of , we further require that . This sum determines the persistence of the process, with larger values of leading to more adaptability. When this stationarity condition is satisfied, and has reached its stationary distribution, the expectation of given is
To ensure for the fitted model, we may parameterize . This constraint is simple enough to enforce for the model (6) and we do so. However, additional parameter constraints such as this may make numerical estimation intractable in more complicated models and they are not enforced by us in the models outlined below.
When is a nonstationary process, the conditional intensity
is also nonstationary. Since , we interpret as the stationary multiplicative deviation, or inflation rate, between and . The process is mean reverting to the process. Let
denote the multiplicative standardized residual process given an estimated CIIR process . If a fitted model defined by (6) sufficiently explains the observed linear dependence in , then an autocorrelation plot of should be statistically insignificant for all lags . As a preview, the standardized residual autocorrelation plot for one such model fit in Section 4 is given in Figure 5(d). The serial correlation appears to have been adequately removed.
Next, we formulate three different nonlinear generalizations of (6) that may better characterize the serial dependence, and possibly lead to improved forecasts. The first is an exponential autoregressive model defined as
(7) 
in which . Exponential autoregressive models are attractive in application because of their thresholdlike behavior. For large , the functional coefficient for is approximately , and for small it is approximately . Additionally, the transition between these regimes remains smooth. As in Fokianos, Rahbek and Tjøstheim (2009), for one can verify the process has a stationarity version when is constant.
We also consider a piecewise linear threshold model
(8) 
in which is an indicator variable and the threshold boundaries satisfy . To ensure positivity of , we assume , , and . Additionally, we take and , such that when is outside the range the CIIR process is more adaptive, that is, puts more weight on and less on . When is constant, the process has a stationary version under the restriction ; see Woodard, Matteson and Henderson (2010). In practice, the threshold boundaries and are fixed during estimation, and may be adjusted as necessary after further exploratory analysis. We chose and , that is, thresholds at 15% above and below 1.
Finally, we consider a model with regime switching at deterministic times, letting
(9) 
This model is appropriate assuming the residual process has two distinct regimes for different periods of the day. For example, one regime could be for normal workday hours with the other regime being for the evening and early morning hours. No stationarity is possible for this model. A drawback of this model is that the process has jumps at and . As was the case for and in (8), and are fixed during estimation. After exploratory analysis, we chose a.m. and p.m.
3.5 Estimation algorithm
The estimation procedure below begins with an iterative algorithm for estimating the factor model from Sections 3.1–3.3 through repeated use of the gam function from the mgcv library in R. Any serial dependence is ignored during estimation of for simplicity. Given estimates for the factor model , conditional maximum likelihood is used to estimate the conditional intensity via one of the time series models given in (6)–(9) for the CIIR process .

Initialization:

[(a)]

Fix and .

Choose some and define .

Apply a singular value decomposition (SVD) to find .

[(iii)]

Let denote the first columns of the left singular matrix .

Let denote the first columns of the right singular matrix .

Let denote the upperleft submatrix of , the diagonal matrix of singular values.


Assign and .

No smoothing is performed and the constraints are omitted in initialization.


Update:

[(a)]

Fit the Poisson GAM model described in Section 3.3 with and as fixed covariates.

Assign as the estimated parameter values from this fit and let .


Fit the Poisson GAM model described in Section 3.3 with as a fixed covariate.

Assign as the estimated parameter values from this fit.


Apply an SVD to find .

[(iii)]

Assign .

Assign .

Assign .


Let .


Repeat the update steps recursively until convergence.
Convergence is reached when the relative change in is sufficiently small. After convergence we can recover from the rows of the final estimate of . These values are then treated as fixed constants during estimation of . We use conditional maximum likelihood to estimate the parameters associated with a time series model for . The recursion defined by (6)–(9) requires initialization by choosing a value for ; the estimates are conditional on the chosen initialization.
We may always specify the joint distribution of the observations as an iterated product of successive conditional distributions for given as
We follow the standard convention of fixing in estimation. For large sample sizes the practical impact of this decision is negligible. We may therefore write the log likelihood function as the sum of iterated conditional log likelihood functions. The conditional distribution for the observations is assumed to be Poisson with mean
For uninterrupted observations over periods we define the log likelihood function as
This recursion requires an initial value for . For simplicity, we use its expected value, . When there are gaps in the observation record, equation (3.5) is calculated over every contiguous block of observations. This requires reinitialization of at the beginning of each block. The log likelihood for the blocks are then added together to form the entire log likelihood. The maximum likelihood estimate is the of this quantity, subject to the constraints given in Section 3.4. Finally, is estimated by the respective recursion given by equations (6)–(9) with parameters replaced by their estimates, again with reinitialization of at the beginning of each contiguous block of observations. Blocks were large enough in our application that the effect of reinitialization was negligible.
4 Empirical analysis
Using the data described in Section 2, we perform the following analysis: (a) we define various statistical goodnessoffit metrics suitable for the proposed models; based on insample performance, these metrics are used to determine the number of factors for use in the dynamic factor models. (b) We compare the outofsample forecast performance for the factor model in (3), the factor model with constraints in (4), and the factor model with constraints and smoothing splines in (3.3). These comparisons help ascertain the improvement from each refinement and validate the proposed selection methods for . (c) For the latter factor model, we compare the outofsample forecast performance with the addition of the CIIR process, via use of the various time series models defined in Section 3.4. (d) We quantify the practical impact of these successive statistical improvements with a queueing application constructed to approximate ambulance operations.
4.1 Interpreting the fitted model
The mean number of calls was approximately 24 per hour for 2007 and 2008, and no increasing or decreasing linear trend in time was detected during this period. We partition the observations by year into two data sets referred to as 2007 and 2008, respectively. Each year is first regarded as a training set, and each model is fit individually to each year. The opposite year is subsequently used as a test set to evaluate the outofsample performance of each fitted model. To account for missing days, we reinitialize the CIIR process in the first period following each day of missing data. This was necessary at most five times per year including the first day of the year.
We found the factor model fit with constraints, smoothing splines, and factors to be the most appropriate of the factor models considered. The estimated factors for 2008 are shown in Figure 2(a). Each of the four factors varies smoothly over the hours of the day via use of smoothing splines. The first factor is strictly positive and the least variable. It appears to capture the mean diurnal pattern. The factor appears to isolate the dominant relative differences between weekdays and weekend days. The defining feature of and is the large increase late in the day, corresponding closely to the relative increase observed on Friday and Saturday evenings. However, decreases in the morning, while increases in the morning and decreases in the late afternoon. As increases, additional factors become increasingly more variable over the hours of the day. Too many factors result in overfitting the model, as the extra factors capture noise.
The corresponding daily factor loadings for the first four weeks of 2008 are shown in Figure 2(b). The loadings are shown to simplify comparisons. The much higher loadings on confirm its interpretation as capturing the mean. The peaks on Fridays coincide with Friday having the highest average number of calls, as seen in Figure 1. Weekdays get a positive loading on , while weekend days get negative loading. Loadings on are lowest on Sundays and Mondays and loadings on are largest on Fridays and Saturdays. As increases, the loadings on additional factors become increasingly close to zero. This partially mitigates the overfitting described above. Factors with loadings close to zero have less impact on the fitted values . Nevertheless, they can still reduce outofsample forecast performance.
The daily factor loadings for all of 2008 are shown in Figure 2(c). The relative magnitude of each loading vector with respect to dayofweek is constant. This results from use of the constraint matrix in (4). As the loadings vary over the days of the week, they also vary smoothly over the course of the year, via use of the constraint matrix and the use of cyclic smoothing splines in estimation of in (4). The loadings on show how the expected number of calls per day varies over the year. The week to week variability in the other loadings influences how the days of the week change relative to each other over the year. Figure 3 shows the estimated intensity process for every day in 2008, shaded by dayofweek. The curves vary smoothly over the hours of the day. The fit for each day of the week keeps the same relative shape, but it varies smoothly over the weeks of the year.
Section 3.4 described incorporating time series models to improve the shortterm forecasts of a factor model. The models capture the observed serial dependence in the multiplicative residuals from a fitted factor model; see Figure 5. Parameter estimates and approximate standard errors for the IntGARCH model are given in Supplemental material (Table 1). A fitted factor model using constraints, smoothing splines and , as well as the factor model including a fitted model , are also shown in Figure 6(a), with the observed call arrivals per hour for Weeks 8 and 9 of 2007. The process is mean reverting about the process. They are typically close to each other, but when they differ by a larger amount, they tend to differ for several hours at a time. The corresponding fitted CIIR process is shown in Figure 6(b). This clearly illustrates the dependence and persistence exhibited in Figure 6(a). The CIIR process ranges between 6% during this period. With a mean of 24 calls per hour, this range corresponds to varying about by about 1.5 expected calls per hour.
4.2 Goodness of fit and model selection
To evaluate the fitted values and forecasts of the proposed models, three types of residuals are computed: multiplicative, Pearson and Anscombe. Their respective formulas for the Poisson distribution are given by
We refer to the root mean square error (RMSE) of each metric as RMSME, RMSPE and RMSAE, respectively. The multiplicative residual is defined as before and is a natural choice given the definition for the CIIR. Since the variance of a Poisson random variable is equal to its mean, the Pearson residual is quite standard. However, the Pearson residual can be quite skewed for the Poisson distribution [cf. McCullagh and Nelder (1989), Section 2.4]. The Anscombe residual is derived as a transformation that makes the distribution of the residuals as close to Gaussian as possible while suitably scaling to stabilize the variance. See Pierce and Schafer (1986) for further discussion of residuals for generalized linear models. While the three methods always yielded the same conclusion, we found use of the Anscombe residuals gave a more robust assessment of model accuracy and simplified paired comparisons between the residuals of competing models.
The three RMSE metrics were used for both in and outofsample model comparisons. For insample comparisons of the factor models, we also computed the deviance of each fitted model . As a goodnessoffit metric, deviance is derived from the logarithm of a ratio of likelihoods. For a log likelihood function , it is defined as
in general. For a fitted factor model, ignoring serial dependence, the deviance corresponding to a Poisson distribution is
in which the first term is zero if .
We compare the fitted models’ relative reduction in deviance and RMSE as we increase the number of factors . Figure 4 shows these results for factor models fit to 2007 data with constraints and smoothing splines. The results for other models and for 2008 were very similar. This plot may be interpreted similarly to a scree plot in PCA by identifying the point at which performance tapers off and the marginal improvement from additional factors is negligible. Under each scenario we consistently selected factors through this graphical criterion. To further justify this as a factor selection strategy, we also consider the impact the number of factors has on outofsample performance for each of the proposed models below. This approach is straightforward, but it does not fully account for the uncertainty on the number of factors. Bayesian estimation would require specialized computation, but it may improve model assessment [see, e.g., Lopes and West (2004)].
4.3 Outofsample forecast performance
Outofsample comparisons were made by fitting models to the 2007 training set and forecasting on the 2008 test set, and vice versa. To make predictions comparable from one year to the next, we align corresponding calendar weeks of the year, not days of the year. This ensures that estimates for Sundays are appropriately compared to Sundays, etc.
The first model considered was the simple prediction (SP) method. This simple moving average involving four observations was defined in theIntroduction. Next, the forecasts of various factor models (FM) were considered. For we evaluated the forecasts from the FM in (3), the FM with constraints in (4), and the FM with constraints and smoothing splines in (3.3). Finally, for the latter FM, with , we calculate the implied fit from the training set with the inclusion of the CIIR process via the various time series models defined in Section 3.4. We compute the forecast RMSE of each model for the three residual types, for both years.
The forecast results are shown in Table 1. The basic FMs did slightly worse than the SP both years. With only one year of observations, these FMs tend to overfit the training set data, even with a small number of factors. The FMs with constraints give a very significant improvement over the previous models. The forecast RMSE is lowest at for the 2007 test set, and at for the 2008 test set. There was also a very large decrease between and . The FMs with constraints and smoothing splines offered an additional improvement. The forecast RMSE is lowest at for both test sets. With the addition of the IntGARCH model for the CIIR process to this model, the RMSE improved again. Application of the nonlinear time series models instead offered only a slight improvement over the IntGARCH model.
With only one year of training data, each FM begins to overfit with factors. Results were largely consistent regardless of the residual used, but the Anscombe residuals were the least skewed and allowed the simplest pairwise comparisons. Although the FMs with constraints had superior insample performance, the use of smoothing splines reduced the tendency to overfit and resulted in improved forecast performance. The CIIR process offered improvements in fit over FMs alone.
We also fit each of the nonlinear time series models discussed in Section 3.4 using a FM with . The regime switching model had the best performance. It had the lowest RMSE for both test sets. The exponential autoregressive and the piecewise linear threshold models performed similarly to the IntGARCH model for both test sets. Although the nonlinear models consistently performed better insample, their outofsample performance was similar to the IntGARCH model.
4.4 Queueing model simulation to approximate ambulance operations
To comprehensively improve ambulance operations, it would be advantageous to simultaneously model the service duration of dispatched ambulances in addition to the demand for ambulance service. Unfortunately, such information was not available. We are currently working with Toronto EMS to use our improved estimates of call arrival rates to improve staffing in their dispatch call center. Extending our approach to a spatialtemporal forecasting model will likely be used to help determine when and where to deploy ambulances.
We present a simulation study that uses a simple queueing system to quantify the impact that improved forecasts have on staffing decisions and relative operating costs, for the Toronto data. The queueing model is a simplification of ambulance operations that ignores the spatial component. Similar queueing models have been used frequently in EMS modeling [see Swersey (1994), page 173]. This goodnessoffit measure facilitates model comparisons and a similar approach may be useful in other contexts.
We use the terminology employed in the call center and queueing theory literature throughout the section; for our application, servers are a proxy for ambulances, callers or customers are those requiring EMS, and a server completing service is equated to an ambulance completing transport of a person to a hospital, etc. As before, let denote the observed number of call arrivals during hour . Our experiment examines the behavior of a simple queueing system. The arrival rate in time period is . During this period, let denote the number of servers at hand. For simplicity, we assume that the service rate for each server is the same, and constant over time. Furthermore, intrahour arrivals occur according to a Poisson process with rate , and service times of callers are independent and exponentially distributed with rate .
As in Section 4.3, models are calibrated on one year of observations and forecasts for are made for the other year. Each model’s forecasts are then used to determine corresponding staffing levels for the system.
To facilitate comparisons of shortterm forecasts, we assume that the number of servers can be changed instantaneously at the beginning of each period. In practice, it is possible to adjust the number of ambulances in real time, but not to the degree that we assume here.
Each call has an associated arrival time and service time. When a call arrives, the caller goes immediately into service if a server is available, otherwise it is added to the end of the queue. A common goal in EMS is to ensure that a certain proportion of calls are reached by an ambulance within a prespecified amount of time. We approximate this goal by instead aiming to answer a proportion, , of calls immediately; this is a standard approximation in queueing applications in many areas including EMS [Kolesar and Green (1998)]. For each call arrival, we note whether or not the caller was served immediately. As servers complete service, they immediately begin serving the first caller waiting in the queue, otherwise they await new arrivals if the queue is currently empty. One simulation replication of the queueing system simulates all calls in the test year.
To implement the queueing system simulation, it is first necessary to simulate arrival and service times for each caller in the forecast period. We use the observed number of calls for each hour as the number of arrivals to the system in period . Since arrivals to the system are assumed to follow a Poisson process, we determine the call arrival times using the wellknown result that, conditional on the total number of arrivals in the period , the arrival times have the same distribution as the order statistics of independent Uniform() random variables. We exploit this relationship to generate the intrahour arrival times given the observed arrival volume . The service times for each call are generated independently with an Exponential() distribution.
The final input is the initial state of the queue within the system. We generate an initial number of callers in the queue as Poisson(), then independently generate corresponding Exponential() residual service times for each of these callers. This initialization is motivated through an infiniteserver model; see, for example, Kolesar and Green (1998). Whenever there is a missing day, in either the test set or corresponding training set period, we similarly reinitialize the state of the queue but with replaced by the number of calls observed in the first period following the missing period. These initializations are common across the different forecasting methods to allow direct comparisons.
To evaluate forecast performance, we define a cost function and an appropriate method for determining server levels from arrival rate estimates. Let denote the number of callers served immediately in period . The hourly cost function is given by
in which
is the targeted proportion of calls served immediately, and is the cost of not immediately serving a customer, relative to the cost of staffing one server for one hour. The total cost, with respect to the hourly server cost, for the entire forecast period is
This approach, where penalties for poor service are balanced against staffing costs, is frequently used; see, e.g., Andrews and Parsons (1993), Harrison, Zeevi and Shum (2005).
At time , the number of call arrivals and the number served immediately are random variables, denoted as and , respectively. A natural objective is to choose staffing levels that minimize the hourly expected cost as
(11) 
in which is assumed to have a Poisson distribution with mean equal to the arrival rate forecast . The staffing levels are then a function of arrival rate forecasts, . We approximate this expectation numerically by randomly generating independent realizations as . Then, for each we simulate one independent realization of . For a fixed value of the expectation is approximated by . We found that provided adequate accuracy.
Independent realizations of require running the queueing system forward one hour, but this is very computationally intensive. To approximate , we use a Binomial distribution. Let . The function gives the steady state probability that a customer is served immediately for a queueing system with a constant arrival rate, server level and service rate, and , respectively. Derivation of this function is available in any standard text on queueing theory [e.g., Gross and Harris (1998), Chapter 2].
Let denote the long run proportion of time such a system contains customers and let . Then
(12) 
When , the arrival rate is faster than the net service rate, and the system is unstable; the long run probability that a customer is served immediately is zero. The binomial approximation greatly reduces the computational costs and provides reasonable results, though it tends to underestimate the true variability of due to the positive correlation in successive caller delays.
A final deliberation is needed on the removal of servers when decreases. In our implementation, idle servers were removed first, and, if necessary, busy servers were dropped in ascending order with respect to remaining service time. We also considered random selection of servers to be dropped. Doing so produced highly variable results, and is under further study. To further simplify the implementation, if it was necessary to drop a busy server, it was simply discarded, along with any remaining service time for that caller. The effect of this simplification depends on the service rate ; our results did not appear to be sensitive to this simplification.
Simulation of the queueing system is now rather straightforward. On each iteration , we note whether each caller was served immediately or not. Forecast performance is assessed by examining the total cost over the test period. For both years, we performed 100 simulations over the test year for each forecast method. To demonstrate the robustness of this methodology, we performed the experiment for several different values of the queuing system’s parameters. Specifically, all combinations of and were considered, after consultation with EMS experts.
Results for the mean hourly cost over the 100 simulations for each forecasting method, for each test year, are summarized in Figure 7. We see that the mean hourly cost is lowest for the FM w/ IntGARCH, followed by the FM only, and finally by SP. All pairwise differences in mean were highly significant; the smallest ratio was 80. In fact, this ordering in performance held for almost every iteration of the queueing system, not just on average.
The mean percentage of callers served immediately can be found in Figure 8. The total number of server hours used was also recorded for each model for each set of parameter values. A table containing the values of all these quantities can be found in the online supplemental material. Both mean percentage served immediately and mean hourly cost increase with . For each test year, for each level of , differed by between one and three thousand serverhours, for the different models.
5 Conclusions
Our analysis was motivated by a data set provided by Toronto EMS. The proposed forecasting method allows parsimonious modeling of the dependent and nonstationary countvalued EMS call arrival process. Our method is straightforward to implement and demonstrates substantial improvements in forecast performance relative to simpler forecasting methods. We measured the impact of our successive refinements to the model, showing the merit of factor model estimation with covariates and smoothing splines. The factor model was able to capture the nonstationary behavior exhibited in call arrivals. Introduction of the CIIR process allowed adaptive forecasts of deviations from this diurnal pattern.
Assessing the impact that different arrival rate forecasts can have on call centers and related applications has received very little attention in the literature. Our databased simulation approach is straightforward to implement, and was able to clearly distinguish the effectiveness of each forecasting method. The simulation results coincide with the outofsample RMSE analysis in Section 4.3 and provide a practical measure of forecast performance. Relative operating cost is a natural metric for measuring call arrival rate forecasts, and our implementation may easily be extended to many customized cost functions and a wide variety of applications.
Ultimately, we seek to strengthen emergency medical service by improving upon relevant statistical methodology. Future work will consider inclusion of additional covariates and study of other nonlinear time series models. Bayesian methods which directly model countvalued observations have desirable properties for inference and many applications, and are under study. Spatial and spatial–temporal analysis of call arrivals will also offer new benefits to EMS.
Acknowledgments
The authors sincerely thank Toronto EMS for sharing their data, in particular, Mr. Dave Lyons for his comments and support.
[id=suppA] \snameSupplement A \stitleAdditional tables \slink[doi]10.1214/10AOAS442SUPPA \slink[url]http://lib.stat.cmu.edu/aoas/442/supplementA.pdf \sdatatype.pdf \sdescriptionTables 1 and 2.
[id=suppB] \snameSupplement B \stitleEstimation algorithms \slink[doi]10.1214/10AOAS442SUPPB \slink[url]http://lib.stat.cmu.edu/aoas/442/supplementB.R \sdatatype.R \sdescriptionR code for estimating the models in Section 3 and for calculating the RMSE metrics in Section 4.
[id=suppC] \snameSupplement C \stitleSimulation algorithms \slink[doi]10.1214/10AOAS442SUPPC \slink[url]http://lib.stat.cmu.edu/aoas/442/supplementC.R \sdatatype.R \sdescriptionR code for implementing the queueing model simulation in Section 4.4.
References
 Andrews and Cunningham (1995) {barticle}[auto:STB—20101118—09:18:59] \bauthor\bsnmAndrews, \bfnmB.\binitsB. \AND\bauthor\bsnmCunningham, \bfnmS.\binitsS. (\byear1995). \btitleLL Bean improves callcenter forecasting. \bjournalInterfaces \bvolume25 \bpages1–13. \endbibitem
 Andrews and Parsons (1993) {barticle}[auto:STB—20101118—09:18:59] \bauthor\bsnmAndrews, \bfnmB.\binitsB. \AND\bauthor\bsnmParsons, \bfnmH.\binitsH. (\byear1993). \btitleEstablishing telephoneagent staffing levels through economic optimization. \bjournalInterfaces \bvolume23 \bpages14–20. \endbibitem
 Bianchi, Jarrett and Choudary Hanumara (1998) {barticle}[auto:STB—20101118—09:18:59] \bauthor\bsnmBianchi, \bfnmL.\binitsL., \bauthor\bsnmJarrett, \bfnmJ.\binitsJ. \AND\bauthor\bsnmChoudary Hanumara, \bfnmR.\binitsR. (\byear1998). \btitleImproving forecasting for telemarketing centers by ARIMA modeling with intervention. \bjournalInternat. J. Forecasting \bvolume14 \bpages497–504. \endbibitem
 Borchers et al. (1997) {barticle}[auto:STB—20101118—09:18:59] \bauthor\bsnmBorchers, \bfnmD.\binitsD., \bauthor\bsnmBuckland, \bfnmS.\binitsS., \bauthor\bsnmPriede, \bfnmI.\binitsI. \AND\bauthor\bsnmAhmadi, \bfnmS.\binitsS. (\byear1997). \btitleImproving the precision of the daily egg production method using generalized additive models. \bjournalCanad. J. Fisheries and Aquatic Sciences \bvolume54 \bpages2727–2742. \endbibitem
 Brown et al. (2005) {barticle}[mr] \bauthor\bsnmBrown, \bfnmLawrence\binitsL., \bauthor\bsnmGans, \bfnmNoah\binitsN., \bauthor\bsnmMandelbaum, \bfnmAvishai\binitsA., \bauthor\bsnmSakov, \bfnmAnat\binitsA., \bauthor\bsnmShen, \bfnmHaipeng\binitsH., \bauthor\bsnmZeltyn, \bfnmSergey\binitsS. \AND\bauthor\bsnmZhao, \bfnmLinda\binitsL. (\byear2005). \btitleStatistical analysis of a telephone call center: A queueingscience perspective. \bjournalJ. Amer. Statist. Assoc. \bvolume100 \bpages36–50. \biddoi=10.1198/016214504000001808, mr=2166068 \endbibitem
 Channouf et al. (2007) {barticle}[auto:STB—20101118—09:18:59] \bauthor\bsnmChannouf, \bfnmN.\binitsN., \bauthor\bsnmL’Ecuyer, \bfnmP.\binitsP., \bauthor\bsnmIngolfsson, \bfnmA.\binitsA. \AND\bauthor\bsnmAvramidis, \bfnmA.\binitsA. (\byear2007). \btitleThe application of forecasting techniques to modeling emergency medical system calls in Calgary, Alberta. \bjournalHealth Care Management Science \bvolume10 \bpages25–45. \endbibitem
 Cox and Isham (1980) {bbook}[mr] \bauthor\bsnmCox, \bfnmDavid Roxbee\binitsD. R. \AND\bauthor\bsnmIsham, \bfnmValerie\binitsV. (\byear1980). \btitlePoint Processes. \bpublisherChapman & Hall, \baddressLondon. \bidmr=0598033 \endbibitem
 Daskalov (1999) {barticle}[auto:STB—20101118—09:18:59] \bauthor\bsnmDaskalov, \bfnmG.\binitsG. (\byear1999). \btitleRelating fish recruitment to stock biomass and physical environment in the Black Sea using generalized additive models. \bjournalFisheries Research \bvolume41 \bpages1–23. \endbibitem
 Ferland, Latour and Oraichi (2006) {barticle}[mr] \bauthor\bsnmFerland, \bfnmRené\binitsR., \bauthor\bsnmLatour, \bfnmAlain\binitsA. \AND\bauthor\bsnmOraichi, \bfnmDriss\binitsD. (\byear2006). \btitleIntegervalued GARCH process. \bjournalJ. Time Ser. Anal. \bvolume27 \bpages923–942. \biddoi=10.1111/j.14679892.2006.00496.x, mr=2328548 \endbibitem
 Fokianos, Rahbek and Tjøstheim (2009) {barticle}[mr] \bauthor\bsnmFokianos, \bfnmKonstantinos\binitsK., \bauthor\bsnmRahbek, \bfnmAnders\binitsA. \AND\bauthor\bsnmTjøstheim, \bfnmDag\binitsD. (\byear2009). \btitlePoisson autoregression. \bjournalJ. Amer. Statist. Assoc. \bvolume104 \bpages1430–1439. \biddoi=10.1198/jasa.2009.tm08270, mr=2596998 \endbibitem
 Goldberg (2004) {barticle}[auto:STB—20101118—09:18:59] \bauthor\bsnmGoldberg, \bfnmJ. B.\binitsJ. B. (\byear2004). \btitleOperations research models for the deployment of emergency services vehicles. \bjournalEMS Management J. \bvolume1 \bpages20–39. \endbibitem
 Gross and Harris (1998) {bbook}[mr] \bauthor\bsnmGross, \bfnmDonald\binitsD. \AND\bauthor\bsnmHarris, \bfnmCarl M.\binitsC. M. (\byear1998). \btitleFundamentals of Queueing Theory, \bedition3rd ed. \bpublisherWiley, \baddressNew York. \bidmr=1600527 \endbibitem
 Gu (1992) {barticle}[auto:STB—20101118—09:18:59] \bauthor\bsnmGu, \bfnmC.\binitsC. (\byear1992). \btitleCrossvalidating nonGaussian data. \bjournalJ. Comput. Graph. Statist. \bvolume1 \bpages169–179. \endbibitem
 Gu (2010) {bmisc}[auto:STB—20101118—09:18:59] \bauthor\bsnmGu, \bfnmC.\binitsC. (\byear2010). \bhowpublishedgss: General smoothing splines. R Package Version 1.13. \endbibitem
 Harrison, Zeevi and Shum (2005) {barticle}[auto:STB—20101118—09:18:59] \bauthor\bsnmHarrison, \bfnmJ.\binitsJ., \bauthor\bsnmZeevi, \bfnmA.\binitsA. \AND\bauthor\bsnmShum, \bfnmS.\binitsS. (\byear2005). \btitleA method for staffing large call centers based on stochastic fluid models. \bjournalManufacturing & Service Operations Management \bvolume7 \bpages20–36. \endbibitem
 Hastie (2009) {bmisc}[auto:STB—20101118—09:18:59] \bauthor\bsnmHastie, \bfnmT.\binitsT. (\byear2009). \bhowpublishedgam: Generalized additive models. R Package Version 1.01. \endbibitem
 Hastie and Tibshirani (1990) {bbook}[mr] \bauthor\bsnmHastie, \bfnmT. J.\binitsT. J. \AND\bauthor\bsnmTibshirani, \bfnmR. J.\binitsR. J. (\byear1990). \btitleGeneralized Additive Models. \bseriesMonographs on Statistics and Applied Probability \bvolume43. \bpublisherChapman & Hall, \baddressLondon. \bidmr=1082147 \endbibitem
 Hastie and Tibshirani (1993) {barticle}[mr] \bauthor\bsnmHastie, \bfnmTrevor\binitsT. \AND\bauthor\bsnmTibshirani, \bfnmRobert\binitsR. (\byear1993). \btitleVaryingcoefficient models. \bjournalJ. Roy. Statist. Soc. Ser. B \bvolume55 \bpages757–796. \bidmr=1229881 \endbibitem
 Henderson (2005) {binproceedings}[auto:STB—20101118—09:18:59] \bauthor\bsnmHenderson, \bfnmS. G.\binitsS. G. (\byear2005). \btitleShould we model dependence and nonstationarity, and if so how? In \bbooktitleProceedings of the 37th Conference on Winter Simulation (\beditorM. Kuhl, N. Steiger, F. Armstrong and J. Joines, eds.) \bpages120–129. \bpublisherIEEE, \baddressPiscataway, NJ. \endbibitem
 Henderson (2009) {bincollection}[auto:STB—20101118—09:18:59] \bauthor\bsnmHenderson, \bfnmS. G.\binitsS. G. (\byear2009). \btitleOperations research tools for addressing current challenges in emergency medical services. In \bbooktitleWiley Encyclopedia of Operations Research and Management Science (\beditorJ. J. Cochran, ed.). \bpublisherWiley, \baddressNew York. \endbibitem
 Kolesar and Green (1998) {barticle}[auto:STB—20101118—09:18:59] \bauthor\bsnmKolesar, \bfnmP.\binitsP. \AND\bauthor\bsnmGreen, \bfnmL.\binitsL. (\byear1998). \btitleInsights on service system design from a normal approximation to Erlang’s delay formula. \bjournalProduction and Operations Management \bvolume7 \bpages282–293. \endbibitem
 Lopes and West (2004) {barticle}[mr] \bauthor\bsnmLopes, \bfnmHedibert Freitas\binitsH. F. \AND\bauthor\bsnmWest, \bfnmMike\binitsM. (\byear2004). \btitleBayesian model assessment in factor analysis. \bjournalStatist. Sinica \bvolume14 \bpages41–67. \bidmr=2036762 \endbibitem
 Matteson and Tsay (2011) {bmisc}[auto:STB—20101118—09:18:59] \bauthor\bsnmMatteson, \bfnmD. S.\binitsD. S. \AND\bauthor\bsnmTsay, \bfnmR. S.\binitsR. S. (\byear2011). \bhowpublishedConstrained independent component analysis. Unpublished manuscript. \endbibitem
 McCullagh and Nelder (1989) {bbook}[mr] \bauthor\bsnmMcCullagh, \bfnmP.\binitsP. \AND\bauthor\bsnmNelder, \bfnmJ. A.\binitsJ. A. (\byear1989). \btitleGeneralized Linear Models, \bedition2nd ed. \bpublisherChapman & Hall, \baddressLondon. \endbibitem
 Opsomer, Wang and Yang (2001) {barticle}[mr] \bauthor\bsnmOpsomer, \bfnmJean\binitsJ., \bauthor\bsnmWang, \bfnmYuedong\binitsY. \AND\bauthor\bsnmYang, \bfnmYuhong\binitsY. (\byear2001). \btitleNonparametric regression with correlated errors. \bjournalStatist. Sci. \bvolume16 \bpages134–153. \biddoi=10.1214/ss/1009213287, mr=1861070 \endbibitem
 Pierce and Schafer (1986) {barticle}[mr] \bauthor\bsnmPierce, \bfnmDonald A.\binitsD. A. \AND\bauthor\bsnmSchafer, \bfnmDaniel W.\binitsD. W. (\byear1986). \btitleResiduals in generalized linear models. \bjournalJ. Amer. Statist. Assoc. \bvolume81 \bpages977–986. \bidmr=0867620 \endbibitem
 R Development Core Team (2009) {bmisc}[auto:STB—20101118—09:18:59] \borganizationR Development Core Team (\byear2009). \bhowpublishedR: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. \endbibitem
 Setzler, Saydam and Park (2009) {barticle}[auto:STB—20101118—09:18:59] \bauthor\bsnmSetzler, \bfnmH.\binitsH., \bauthor\bsnmSaydam, \bfnmC.\binitsC. \AND\bauthor\bsnmPark, \bfnmS.\binitsS. (\byear2009). \btitleEMS call volume predictions: A comparative study. \bjournalComput. Oper. Res. \bvolume36 \bpages1843–1851. \endbibitem
 Shen and Huang (2008a) {barticle}[mr] \bauthor\bsnmShen, \bfnmHaipeng\binitsH. \AND\bauthor\bsnmHuang, \bfnmJianhua Z.\binitsJ. Z. (\byear2008a). \btitleForecasting time series of inhomogeneous Poisson processes with application to call center workforce management. \bjournalAnn. Appl. Statist. \bvolume2 \bpages601–623. \biddoi=10.1214/08AOAS164, mr=2524348 \endbibitem
 Shen and Huang (2008b) {barticle}[auto:STB—20101118—09:18:59] \bauthor\bsnmShen, \bfnmH.\binitsH. \AND\bauthor\bsnmHuang, \bfnmJ.\binitsJ. (\byear2008b). \btitleInterday forecasting and intraday updating of call center arrivals. \bjournalManufacturing and Service Operations Management \bvolume10 \bpages391–410. \endbibitem
 Stock and Watson (2002) {barticle}[mr] \bauthor\bsnmStock, \bfnmJames H.\binitsJ. H. \AND\bauthor\bsnmWatson, \bfnmMark W.\binitsM. W. (\byear2002). \btitleMacroeconomic forecasting using diffusion indexes. \bjournalJ. Bus. Econom. Statist. \bvolume20 \bpages147–162. \biddoi=10.1198/073500102317351921, mr=1963257 \endbibitem
 Swersey (1994) {bincollection}[auto:STB—20101118—09:18:59] \bauthor\bsnmSwersey, \bfnmA.\binitsA. (\byear1994). \btitleThe deployment of police, fire, and emergency medical units. In \bbooktitleHandbooks in Operations Research and Management Science \bvolume6 \bnote151–200. \baddressNorthHolland, Amsterdam. \endbibitem
 Takane and Hunter (2001) {barticle}[mr] \bauthor\bsnmTakane, \bfnmYoshio\binitsY. \AND\bauthor\bsnmHunter, \bfnmMichael A.\binitsM. A. (\byear2001). \btitleConstrained principal component analysis: A comprehensive theory. \bjournalAppl. Algebra Engrg. Comm. Comput. \bvolume12 \bpages391–419. \biddoi=10.1007/s002000100081, mr=1864610 \endbibitem
 Tsai and Tsay (2010) {barticle}[auto:STB—20101118—09:18:59] \bauthor\bsnmTsai, \bfnmH.\binitsH. \AND\bauthor\bsnmTsay, \bfnmR. S.\binitsR. S. (\byear2010). \btitleConstrained factor models. \bjournalJ. Amer. Statist. Assoc. \bvolume105 \bpages1593–1605. \endbibitem
 Tych et al. (2002) {barticle}[auto:STB—20101118—09:18:59] \bauthor\bsnmTych, \bfnmW.\binitsW., \bauthor\bsnmPedregal, \bfnmD.\binitsD., \bauthor\bsnmYoung, \bfnmP.\binitsP. \AND\bauthor\bsnmDavies, \bfnmJ.\binitsJ. (\byear2002). \btitleAn unobserved component model for multirate forecasting of telephone call demand: The design of a forecasting support system. \bjournalInternat. J. Forecasting \bvolume18 \bpages673–695. \endbibitem
 Weinberg, Brown and Stroud (2007) {barticle}[mr] \bauthor\bsnmWeinberg, \bfnmJonathan\binitsJ., \bauthor\bsnmBrown, \bfnmLawrence D.\binitsL. D. \AND\bauthor\bsnmStroud, \bfnmJonathan R.\binitsJ. R. (\byear2007). \btitleBayesian forecasting of an inhomogeneous Poisson process with applications to call center data. \bjournalJ. Amer. Statist. Assoc. \bvolume102 \bpages1185–1198. \biddoi=10.1198/016214506000001455, mr=2412542 \endbibitem
 Whitt (2002) {bbook}[mr] \bauthor\bsnmWhitt, \bfnmWard\binitsW. (\byear2002). \btitleStochasticProcess Limits. \bpublisherSpringer, \baddressNew York. \bidmr=1876437 \endbibitem
 Wood (2003) {barticle}[mr] \bauthor\bsnmWood, \bfnmSimon N.\binitsS. N. (\byear2003). \btitleThin plate regression splines. \bjournalJ. R. Stat. Soc. Ser. B Stat. Methodol. \bvolume65 \bpages95–114. \biddoi=10.1111/14679868.00374, mr=1959095 \endbibitem
 Wood (2006) {bbook}[mr] \bauthor\bsnmWood, \bfnmSimon N.\binitsS. N. (\byear2006). \btitleGeneralized Additive Models: An Introduction with . \bpublisherChapman & Hall, \baddressBoca Raton, FL. \bidmr=2206355 \endbibitem
 Wood (2008) {bmisc}[auto:STB—20101118—09:18:59] \bauthor\bsnmWood, \bfnmS.\binitsS. (\byear2008). \bhowpublishedmgcv: GAMs with GCV/AIC/REML smoothness estimation and GAMMs by PQL. R Package Version 1.61. \endbibitem
 Woodard, Matteson and Henderson (2010) {bmisc}[auto:STB—20101118—09:18:59] \bauthor\bsnmWoodard, \bfnmD. B.\binitsD. B., \bauthor\bsnmMatteson, \bfnmD. S.\binitsD. S. \AND\bauthor\bsnmHenderson, \bfnmS. G.\binitsS. G. (\byear2010). \bhowpublishedStationarity of countvalued and nonlinear time series models. Unpublished manuscript. \endbibitem