Proceedings of the 2021 Winter Simulation Conference
S. Kim, B. Feng, K. Smith, S. Masoud, Z. Zheng, C. Szabo, and M. Loper, eds.

SHORT-TERM ADAPTIVE EMERGENCY CALL VOLUME PREDICTION
Elioth Sanabria
Henry Lam
Enrique Lelo de Larrea
Jay Sethuraman

Sevin Mohammadi
Audrey Olivier
Andrew W. Smyth

Department of IE&OR
Columbia University
500 West 120th Street
New York, NY 10027, USA

Department of Civil Eng. & Eng. Mechanics
Columbia University
500 West 120th Street
New York, NY 10027, USA

Edward M. Dolan
Nicholas E. Johnson
Timothy R. Kepler
Afsan Quayyum
Kathleen S. Thomson
Bureau of Management Analysis and PlanningFire Department, City of New York
9 MetroTech Center
Brooklyn, NY 11201, USA

ABSTRACT
Sudden periods of extreme and persistent changes in the distribution of medical emergencies can trigger
resource planning inefficiencies for Emergency Medical Services, causing delayed responses and increased
waiting times. Predicting such changes and reacting adaptively can alleviate these adversarial impacts. In
this paper, we propose a simple framework to enhance historically calibrated call volume models, the lattera focus of study in the arrival estimation literature, to give more accurate short-term prediction by refitting
their residuals into time series. We discuss some justification of our framework from the perspective of
doubly stochastic Poisson processes. We illustrate our methodology in predicting the hourly call volume to
the 911 call center during the Covid-19 pandemic in NYC, showing how it could improve the performanceof baseline historical estimators by close to 50% measured by the out-of-sample prediction error for the
next hour.
1

INTRODUCTION

During the Covid-19 pandemic, Emergency Medical Services (EMS) suffered episodes of unusual increased
demand that were severe both in magnitude and duration. EMS worldwide were overburdened to levels
beyond maximum capacity during the peak periods of the Covid-19 waves. For example, during the worst
months of the early Covid-19 pandemic in NYC, the call volume to the 911 ambulance call center was
significantly above the historical average (see Figure 1).

Sanabria, Dolan, Johnson, Kepler, Lam, Lelo de Larrea, Mohammadi, Olivier, Quayyum, Sethuraman,
Smyth, and Thomson

450

Historical Average λ(t)400Realized Call Volume

Call Volume

350
300
250
200
150
100
2020-03-21

2020-03-25

2020-03-292020-04-01

2020-04-05

2020-04-09

2020-04-13

2020-04-17

2020-04-21

Date

Figure 1: Realized hourly call volume (red line) vs baseline historical average λ (t) for the month withhighest call volume during 2020.
During unusual periods such as those encountered in Figure 1, using solely a historical estimationof the demand could suffer poor performances. On the other hand, much of the call volume or arrival
process literature, at least within the operations research community, has focused on historical models.
Typically, an arrival process is modeled by means of specifying an intensity function (either parametricor non-parametric) of a Non-homogeneous Poisson process, and fitted to historical data using Maximum
Likelihood (ML). These include, for instance, piece-wise constant intensities in Morgan et al. (2016) andlinear in Glynn and Zheng (2019), and splines in Morgan et al. (2019). Taylor and Letham (2018) study
a “prophet” forecasting paradigm to estimate intensity functions by combining many popular models usedwidely by practitioners. Also related is Kim and Whitt (2014) that test the Poisson process assumption ofcall arrival processes with piece-wise constant or linear intensities, by dividing time into subintervals andtesting the conditional order statistic distribution using the Kolmogorov-Smirnov test. Moreover, in thespirit of time series modeling, Ibrahim and L’Ecuyer (2013) study linear models for call center volume
prediction, and Matteson et al. (2011) use a latent variable time series model to forecast based on seasonalcovariates of medical emergency arrivals. Both of these works are similar to our approach, although we
focus on studying the effect of the estimation window length for capturing the changes of distribution, as
well as capturing the conditional changes of the arrival counts distribution. Finally, for a survey on callcenter arrival modeling, see Ibrahim et al. (2016).
To handle and react to unusual patterns that deviate from the historical, in this paper we focus onshort-term prediction that corrects for the latter. In particular, we consider a simple framework to correct agiven historical estimator (which could be built from any approaches described above) in order to adapt torecent changes in the call arrival pattern. More specifically, our framework models the residuals between
the recent observed volume and the historical model prediction via an additional layer of time series, which
gives a short-term adjustment on the historical model based on the most recent data. This allows us to makeshort-term forecast that is better than using a standalone historical model or a pure autoregression-type timeseries. Our framework resembles time series with seasonality (e.g., Brockwell and Davis (2016) Section
1.4), but the seasonality component can be flexibly modeled by any established historical call volume
models in the arrival estimation literature. We discuss how our framework connects to the use of doubly
stochastic Poisson processes.

Sanabria, Dolan, Johnson, Kepler, Lam, Lelo de Larrea, Mohammadi, Olivier, Quayyum, Sethuraman,
Smyth, and Thomson
To give some motivational illustration, suppose we simply use the historical average of the same time of
the year as the historical model. Figure 2 shows that the higher call volume days of the Covid-19 pandemic
in NYC bore a sizable residual (i.e., gap) between the observed volume and the historical average very
early in the morning (the figure shows that after 3 a.m. the call volume was already outside the historicalconfidence interval). If we capture the trend of the residual as a time series as in our proposed framework
early in the morning (for example at 6a.m. in Figure 2, when the residual gap was already sizable), thehours of higher volume later in the day (after 10 a.m.) could have been forecasted, and action could havebeen preemptively taken to deal with the increased demand. Indeed, the primary motivation for studying
our approach is to enhance downstream decision-making tasks including adaptive staffing and call shedding
to more accurately react to the real-time call volume trend.

350

Call Volume

300

250

200

150

Realized Call Volume Nt
Historical Average Volume lt

100
03-29 00

03-29 03

03-29 06

03-29 09

03-29 1203-29 15

03-29 18

03-29 2103-30 00

Date

Figure 2: Realized hourly call volume vs historical average estimated using hourly arrivals (by hour of theweek with data from 2018-2019), for a day with previously unseen high call volume during the Covid-19
crisis.
We close this introduction by discussing a small number of works on short-term prediction for callvolumes to our knowledge. L’Ecuyer et al. (2018) specify a doubly stochastic intensity function tocharacterize short lived call bursts in a call center (depending on 3 random parameters, denoting theamplitude, decay and duration of the burst). Along the same vein, Oreshkin et al. (2016) study a doubly
stochastic intensity, specified as a deterministic intensity multiplied by a time-varying stochastic factor.
These approaches have in common the need of specifying an underlying doubly stochastic intensity tocapture the periods when the call distribution changes. Our fitting approach can also be viewed as modelinga doubly stochastic process, but instead of fitting an explicit model for the intensity, we formulate adata-driven time series estimation of the prediction error of a given historical prediction model which isthen related to an implicit underlying doubly stochastic intensity.
The rest of this paper is as follows. Section 2 presents our algorithm and analysis. Section 3 presents anumerical example of our methodology to predict 911 call volume during the Covid-19 pandemic in NYC.
Section 4 gives some concluding remarks.
2

METHODOLOGY

We consider discrete time intervals of size ∆, say 1 hour. Time is indexed as ∆t = t for t ∈ Z. We set thefollowing notation:

Sanabria, Dolan, Johnson, Kepler, Lam, Lelo de Larrea, Mohammadi, Olivier, Quayyum, Sethuraman,
Smyth, and Thomson••••

Nt := Random arrival count in the time interval [t − 1,t). Think of this as the arrival counts for acall center.
lt := Historical estimator for ENt . This estimator can be a seasonal historical
average of the arrivals
Rtor be associated to a calibrated intensity function λ (u) such that lt = t−1
λ (u)du.
εt = Nt − lt := Residual (error) between the random unobserved volume Nt and the historicallycalibrated expected volume lt .
F(t−k):t := Information from time t − k to t, namely {N j }tj=t−k . This is used to represent the recently
observed arrival counts.

We assume the historical model or estimator {lt } is given and is independent of the most recent
observations. We aim to model the series {εt } which is the residual between recent volume Nt andthe historical model lt . Then, suppose we want to predict Nt+1 for the next hour t + 1, we would firstcompute E(ε̂t+1 |F(t−k):t ), i.e., the expected residual conditional on the most recent information, and then the
short-term estimator for Nt+1 is ŝt+1 := lt+1 + E(ε̂t+1 |F(t−k):t ). In other words, ŝt+1 corrects the historicalestimator lt+1 shifting it by the amount implied by the expected residual E(ε̂t+1 |F(t−k):t ). We summarize
our proposed procedure in Algorithm 1, where we use an AR-GARCH (Brockwell and Davis (2016)) tofit the residuals. The AR component is used to model a time-varying mean of the residual εt , while the
GARCH is used to model the variance. We could choose other time series models, but the one in Algorithm
1 appears to work reasonably in our experimental tests, and we can also make some intuitive (though not
rigorous) argument to support its use, which we discuss next.
Algorithm 1: Short-term estimator with estimation window of size k at time t.
Inputs: Estimation window k, historical arrival counts {N j }tj=t−k , given historical model {l j }tj=t−k .
1. Compute residual series {ε j := N j − l j }tj=t−k .
p
2 + r β σ2 , η =ϕi εt−i + Zt σt , σt2 = α + ∑qj=1 α j ηt−
2. Fit an AR-GARCH model εt = ϕ + ∑i=1
∑s=1 s t−s tjZt σt , by Maximum Likelihood.
p
3. Output the short-term predictor ŝt+1 = lt+1 + E(ε̂t+1 |F(t−k):t ) = lt+1 + ϕ̂ + ∑i=1
ϕ̂i εt+1−i .
2.1 Intuitive Justification of AR-GARCH via Doubly Stochastic Poisson Processes
We provide some intuitive justification to using Algorithm 1, based on the perspective that the arrivals
naturally arise from a doubly stochastic Poisson process (DSPP), a model that is widely used in the arrivalmodeling literature. First, supposing the arrivals follow a non-homogeneous Poisson Process (NHPP)
with deterministic
intensity λ (u), then the arrival count in time interval [t − 1,t] is Poisson with rateRt
lt = ENt = t−1
λ (u)du. In this setting, the residual series εt = Nt − lt is a sequence of independent meanzero random variables (with variance lt ). If we use Algorithm 1, then when k increases, the AR-GARCH
model would simply converge to a constant 0 mean for the residual and time dependent variance lt . In
this case, our approach does not do any enhancement to the historical model, but we also don’t harm itmaterially either.
For a general DSPP, the intensity function λ is itself random, and given the intensity realization
the arrivals are Poisson like in NHPP. In this case, the ultimate arrival counts possess both the intensity
randomness and the arrival randomness. Given a realization of the intensity, suppose it remains fixed forsome amount of time (reducing the intensity randomness in the short-term) and that the realized intensity
can be approximated by an autoregressive recursion (more on this later), then, the residual is composed of adeterministic autoregressive recursion (from the intensity) plus the randomness of the arrivals, leading to an
autoregressive process (AR) for the residual. Likewise, as the variance of the Poisson is equal to its mean,
the variance also has an autoregresive structure in the short-term, captured by the GARCH specification in
our algorithm. Thus, both the mean and variance of the residual depend in some sense on the observationsin the recent past periods. Intuitively, estimating the mean of the residual (AR part), requires accounting(controlling) for the time-varying variance of the arrivals (GARCH) whenever present.

Sanabria, Dolan, Johnson, Kepler, Lam, Lelo de Larrea, Mohammadi, Olivier, Quayyum, Sethuraman,
Smyth, and Thomson
2.2 Regime-Switching Models: A Simple Example
Consider the following regime-switching model with two regimes: Arrivals are generated by a doubly
stochastic Poisson intensity λ (u) = `t λ1 (u) + (1 − `t )λ2 (u) where `t ∈ {0, 1} is an indicator fixed for adeterministic amount of time T , and after that, it is sampled as a Bernoulli random variable with success rate p.
In this example, `t represents an unobservable random mechanism that selects one of the two intensities λ1 , λ2
and fixes either for a deterministic amount of time T . Moreover, p % 1, meaning that with high probability
only the intensity λ1 (u) is observed most of the time. Calibrating a (discrete) non-homogeneous process onRR
this doubly stochastic process gives lt+1 = ENt+1 = tt+1 [pλ1 (u) + (1 − p)λ2 (u)]du ≈ tt+1 λ1 (u)du. In thiscase, the distribution of the residual εt+1 conditional on `t+1 can be written as mixture of two independent
normals as:
h
p i


d
0 √
st+1 ,(1)
εt+1 |`t+1 = `t+1 (lt+1 − lt+1 ) + Zt+1 lt+1 + (1 − `t+1 ) (st+1 − lt+1 ) + Zt+10
where Zt+1 , Zt+1
are independent standard normals. st+1 := tt+1 λ2 (u)du is the Poisson mean when `t+1 = 0(in this case the distribution is centered at st+1 − lt+1 as opposed to 0 when `t+1 = 1). This expression
comes from a CLT approximation of the Poisson distribution to the normal assuming the rates lt+1 , st+1
are large enough.
Without loss of generality the residual can be written as a normal with random mean µt+1 and random
2 . Where the randomness is conditional on `variance σt+1
t+1 :

R

d

εt+1 = µt+1 (`t+1 ) + Zt+1 σt+1 (`t+1 ).

(2)

While µt+1 , σt+1 are random (depending on `t+1 ), assume It+1 has been fixed at least for k amount of time(with k << T ), the path defined by µt−k , σt−k · · · , µt−1 , σt−1 , µt , σt has the same intensity selected, that is,`t−k = · · · = `t−1 = `t = i. Thus, if µt+1 , σt+1 can be written as functions f , g of the previous k observationsµt−k , σt−k · · · , µt , σt , then, the distribution of εt+1 can be estimated at time t and be used to correct the biasof the deterministic predictor lt+1 . As:dεt+1 = µt+1 (i) + Zt+1 σt+1 (i) = f (µt−k , σt−k · · · , µt , σt ) + Zt+1 g(µt−k , σt−k · · · , µt , σt ),(3)which in turn can be used to correct the prediction lt+1 as E(εt+1 |F(t−k):t ) is equal to E(µt+1 |F(t−k):t ) =f (µt−k , σt−k · · · , µt , σt ), that is, the expected bias (residual) of the estimator lt+1 conditional on F(t−k):t .
This follows from:lt+1 + E(µt+1 |F(t−k):t ) = lt+1 + E(εt+1 |F(t−k):t ) = lt+1 + E(Nt+1 − lt+1 |F(t−k):t ) = E(Nt+1 |F(t−k):t ).

(4)

The first equality follows from the above discussion, the second one follows by definition, while for the last
one note that lt+1 is constant and Ft measurable (known). Then, denote our corrected short-term estimatorby st+1 = lt+1 + E(µt+1 |F(t−k):t ).
We summarize the intuition of this section in the following points:•••

The residual εt+1 = Nt+1 − lt+1 can be written as a normal random variable with random (time2 dependent on a unobservable random variable `dependent) mean and variance µt+1 , σt+1
t+1 thatis persistent over time (meaning that the unobservable variable remains constant for some periodof time).
If µt+1 , σt+1 can be written as functions f , g of k previous observations µt−k , σt−k · · · , µt , σt . Then,the distribution of the residual εt+1 can be estimated at time t as µt−k , σt−k · · · , µt , σt and µt+1 , σt+1
share the same unobservable driving random variable.
With an estimator of µt+1 at time t we get a conditional prediction of the unobserved volume
E(Nt+1 |F(t−k):t ) given by lt+1 + E(µt+1 |F(t−k):t ).

Sanabria, Dolan, Johnson, Kepler, Lam, Lelo de Larrea, Mohammadi, Olivier, Quayyum, Sethuraman,
Smyth, and Thomson
In this paper we take the functional form of f , g as a time series model with time-varying mean/variancep
ϕi εt−i + Zt σt(modeling the variance as a GARCH model, see Bollerslev (1986)) given by εt = ϕ + ∑i=1
p
in Algorithm 1. Letting µt := ϕ + ∑i=1 ϕi εt−i leads to an equivalent expression to Equation (2) with
εt = µt + Zt σt . In this case, the predicted correction is given by E(εt+1 |F(t−k):t ) = E(µt+1 |F(t−k):t ) =p
ϕ + ∑i=1ϕi εt+1−i .
In the next subsection we generalize the unobservable mechanism to change distributions over time
beyond two distributions and a fixed amount of time to a case with countably many distributions andrandom time for switching them. We show that the time series model described in Algorithm 1, while
simple, is flexible enough to accommodate cases of interest for changes of regime.
2.3 General Regime-Switching Models
For a generalization of the previous case, considerR a countable collection of γ-periodic, deterministictintensities {λi (u)}∞i=1 . As time is discrete, let st (i) := t−1 λi (u)du, the periodicity of λi implies st (i) = st+γ (i)for all t. Let I be an independent discrete random variable with support {1, . . . , ∞} sampled with probabilities{pi }∞
i=1 denoting which regime is sampled when there is a regime change. Denote τ as a geometric randomvariable with mean Eτ. Let `t be a stochastic process written as:`t = `t−1 (1 − B) + IB,(5)

where B is a Bernoulli random variable with success probability 1/Eτ. `t can be interpreted as the indicator
to select one of the intensities {λi (u)}∞i=1 that will be set for a geometric time until `t changes (switches)to another intensity (once B is 1 and I is sampled). That is, the arrivals are distributed Poisson conditionalon the intensity selected by `t (for a geometric amount of time) as:
Nt ∼ Poisson(st (`t )|`t ).

(6)

Likewise, denote lt := ENt = E[st (`t )] = ∑∞i=1 pi st (i) and εt+1 = Nt+1 − lt+1 . As in the previous section thepdresidual can be written as εt+1 = µt+1 + Zt+1 σt+1 with µt+1 = st+1 (`t+1 ) − lt+1 and σt+1 = st+1 (`t+1 ).
In the previous section and in Algorithm 1 we propose an AR-GARCH F(t−k):t measurable approximation of εt+1 . In general, for such approximations to be correct further assumptions are needed on themean of the geometric time Eτ, the functional form of the intensities {λi (u)}∞i=1 and how concentrated
are the sampling probabilities {pi }∞
.
Next,
we
discuss
some
cases
of
interest
when the mean µt+1 andi=1variance σt+1 can be written as time series model and estimated using our procedure in Algorithm 1. The
first, is the most common way of thinking of shifts in the trend of the residual. The second shows that
continuous and periodic intensity rates have an autoregressive structure representation that leads to a timeseries model for the residual. And in the last example a totally general change in the distribution of thearrival counts (intensity) is considered where the estimation price is that the regime changes should be
longer.
Shifted Intensities: For example, the time series model in Algorithm 1 corrects the bias of the historicalestimator lt if the intensities are just “shifted” from some baseline
intensity, that is λi (u) = ci + λ (u) whereRt
ci is the shift from the baseline intensity. Here, lt = Eci + t−1 λ (u)du. In this case, µt+1 is constant andequal to (c`t − Eci )∆ while `t remains fixed (that is, for the geometric amount of time τ). Conditional thatthe intensity selected remains the same, that is, `t+1 = `t = · · · . We have µt+1 could be written as a AR(1)ϕprocess with ϕ + ϕ1 εt such that µt = 1−ϕ= (c`t − Eci )∆ = µt+1 (recalling the mean of an AR process is1
constant, then taking expected values on both sides of εt+1 = ϕ + ϕ1 εt yields µt+1 = ϕ + ϕ1 µt , as µt+1 = µtϕimplies µt = µt+1 = 1−ϕ).
1
To see why (c`t − Eci )∆ = µt+1 , conditional `t+1 note that E(Nt+1 |`t+1 ) = ∆c`t+1 + tt+1 λ (u)du, addiR
tionally, E(εt+1 |`t+1 ) = µt+1 = E(Nt+1 − lt+1 |`t+1 ) = E(Nt+1 |`t+1 ) − [∆Eci + tt+1 λ (u)du]. Putting the two
expressions together yields µt+1 = (c`t+1 − Eci )∆, which as long as the intensity does not change, that is,`t+1 = `t = · · · , the means µt+1 = µt will remain constant (for a geometric amount of time).
RSanabria, Dolan, Johnson, Kepler, Lam, Lelo de Larrea, Mohammadi, Olivier, Quayyum, Sethuraman,
Smyth, and ThomsonqR
On the other hand, the variance is time dependent. That is σt+1 = c`t+1 + tt+1 λ (u)du, thus the needto specify the GARCH component of the autoregressive process (with a constant intercept α = c`t+1 and aRseasonal time-varying term tt+1 λ (u)du). Estimating the model with a rolling window of k observations(as seen later on a condition for k << Eτ, such that the intensity indicator `t remains constant with highprobability for the estimation window) allows to estimate the mean and variance of the residual using
Algorithm 1 for large enough k.
Continuous intensities/rates: As another key use case of the generality of our approach, here weshow the intuition that continuous intensities/rates (in a closed interval) can be written as autorregresive
processes specified by our Algorithm 1: By Fourier analysis it is a well known fact that any continuousfunction in a compact interval can be approximated by a linear combination of sines and cosines. Then,as the rates considered are periodic, they can be approximated as a Fourier series. Next, it can be shownthat a sine function can be written recursively as sin(t + 1) = 2 cos(∆) sin(t) − sin(t − 1) (recalling timeis indexed t = ∆t. The identity follows from elementary trigonometric identities, and the same can bedone for cosines). For example, letting the rates be st (i) = sin(t) + ci , we have that they follow an AR(2)recursion as st+1 (i) = φ + φ1 st (i) − st−1 (i) with φ = 2ci (1 − cos ∆) and φ1 = 2 cos ∆. As this procedure canbe done for sines and cosines with arbitrary amplitudes, by Fourier analysis this can be generalized toapproximate any continuous and periodic rate function (including lt ), that is, as every sine and cosine hasan AR representation, a linear combination of them (the Fourier series) is also an AR model. In summary,continuous and periodic rate functions have an associated AR representation.
pdRecall the residual is equal to εt+1 = Nt+1 − lt+1 = (st+1 (`t+1 ) − lt+1 ) + Zt+1 st+1 (`t+1 ). As discussed
in the previous paragraph st+1 (`t+1 ) has an autoregressive recursion (and so does st+1 (`t+1 ) − lt+1 , as adifference of periodic functions is also periodic), the residual follows an AR-GARCH functional form as
in Algorithm 1 while `t+1 = `t−k , i.e. the regime has not switched. Intuitively, the mean of the residual
st+1 (`t+1 ) − lt+1 has a different autoregressive representation than its variance st+1 (`t+1 ), which reflects theneed to specify the GARCH component when fitting the time series model in addition to the AR component
for the mean.
General Intensities and Seasonality: Exploiting the γ-periodicity of the intensities leads to anotherinteresting case: Recall that the point of doing prediction in this setting amounts to find F(t−k):t measurable
approximation/estimation of µt+1 and σt+1 while the regime (intensity) does not change. As the intensities
are γ periodic, µt+1 and σt+1 could be estimated by looking at the observations {εt+1−γk }Kk=1 , that, as longas the indicator for intensities does not change, these are i.i.d. observations with the same distribution as
εt+1 which is the same as the NHPP case except in the sense that the estimation is done in a window ofsize K instead of taking all history. Intuitively, this could be stated as: if the periodicity of the intensity
is daily and there has not been a change in the distribution of the arrivals, then, a good estimator of theprediction error (residual) for the next hour is the average error of the same hour yesterday and the days
before (while the intensity has not changed). This can be stated as a seasonal time series model by:
εt+1 = α + β εt+1−γk + Zσt+1 .

(7)

As long as the regime has not changed interval considered (of size γK), the observations are i.i.d and themean and variance of the error are constant (modulo γ). In practice, the price to pay for the rolling window
estimation is that the time τ for switching the intensity distribution should be in average longer, that is,γK << Eτ (more on this later). As long as that is the case, estimating in a rolling window of size Kγ yields
consistent estimates for the residual εt+1 . Note that in this case there are no restrictions on the functional
form of the intensities nor their sampling probabilities, moreover, this regression can be estimated using
our procedure in Algorithm 1 with only an AR specification of the model (modulo γ).
Limitations: The success on the convergence of the parameters of the window estimation in Algorithm
1 (of size γK or k) hinges on the assumption that the changes of regime are not very frequent, that is,γK << Eτ or k << Eτ, such that the distribution of the arrivals remains the same in the estimation andforecasting window, and the mean and variance of the residual (modeled as time series) corresponds to a

Sanabria, Dolan, Johnson, Kepler, Lam, Lelo de Larrea, Mohammadi, Olivier, Quayyum, Sethuraman,
Smyth, and Thomson
single regime. If Eτ is small, the changes of arrival distribution are frequent and the window estimation
will intuitively fail as the likelihood of having a regime change in the estimation window is equal to1 − (1 − Eτ −1 )k (one minus the probability there are no regime changes in k hours). For a window of size
k to have only one regime (with probability a) it should be true that Eτ is at least 1−a1 1/k (this comes fromsolving (1 − Eτ −1 )k ≥ a). For example, for a window of size k = 100 hours for the estimation to have only
one regime (w.p. a = 95%), Eτ should be at least 2000 hours in average.
Nonetheless, the case when the regimes change too often highlights a key limitation of the window
estimation method and the need to use different methods in this case (or restricting the assumptions on thechanges of the intensities when the regime changes). In general, having larger windows is better for theasymptotic behavior of the parameters of the time series model with the caveat the estimation is correct(only has one regime) with probability (1 − Eτ −1 )k , which is a decreasing function on k and dependent on
Eτ.
3

NUMERICAL EXAMPLESIn this section we start by calibrating a piecewise constant intensity λ (t) as described in Section 2 for thecall arrivals to the FDNY’s 911 EMS call center with historical data from January 2018 to December 2019.
The calibrated intensity is plotted in Figure 3. We use this calibrated intensity as the baseline predictionR
lt+1 = tt+1 λ (u)du = ENt+1 . This calibrated intensity can be simply seen as the historical average weekly
incidents for every hour of the week.

300

λ(t)

250
200
150
100

0

25

50

75

100

125

150

175

t
Figure 3: Example of a calibrated weekly arrival intensity λ (t) for the 168 hours of the week with confidence
intervals.
For testing our model vs the baseline prediction we use the hourly call arrivals from January to December
2020. In Figure 1 we plot a period where the call volume was higher than previously seen. In particular,
the dates around April 1st showed historical highs in total call volumes at the FDNY 911 call center.
We predict the hourly call volume using Algorithm 1 using an estimation window of k = 150 hours.
That is, to predict the call volume next j-th hour t + j at time t, the previous k = 150, 200, 250 hours
are used to estimate the parameters of the AR-GARCH model, which are used to compute the short-term
prediction ŝt+ j = lt+ j + E(ε̂t+ j |F(t−k):t ). We measure the performance by estimating the average absoluteSanabria, Dolan, Johnson, Kepler, Lam, Lelo de Larrea, Mohammadi, Olivier, Quayyum, Sethuraman,
Smyth, and Thomson
error, that is, E(|Nt+ j − ŝt+ j |) out of sample. As a reference benchmark we compare with the mean absolute
error of the historical prediction, that is, E(|Nt − lt |) for the same testing period. In Table (1) we report
the performance of the model. For the AR-GARCH model we use order (p = 2, q = 1, r = 1) in reference
to Algorithm 1 (this order of the model was the most numerically stable in our data. Nonetheless, testing
different order combinations is encouraged on a case by case basis).
Our short-term predictor ŝt+ j improves the performance of the long-term estimation lt+ j close to 50%
every hour when predicting the next hour ( j = 1). As the forecasting horizon j increases, the performance
deteriorates but plateaus around a mean absolute error of 20 calls per hour, which is still a 25% improvement
over the historical estimator lt . The plateau phenomenon is expected due to the mean reverting nature of
ARIMA models as time goes to infinity, making the prediction of the residual time series constant forlong prediction horizons, thus, the emphasis on the short-term prediction use of the model and the need to
re-calibrate the parameters on a rolling window. It can also be seen that changing the prediction window
k does not give a significant performance improvement (lower values of k created numerical instabilities
in the maximum likelihood estimation of the GARCH component and are omitted).
Table 1: Out-of-sample mean absolute error E(|Nt+ j − ŝt+ j |) of our short-term estimator ŝt+ j in Algorithm
1 for different estimation windows k and prediction horizons j. As a benchmark, the mean absolute errorof the historical predictor E(|Nt − lt |) is equal to 27.48 calls per hour.

k = 150k = 200
k = 250

j=1
14.63
14.64
14.64

Hours predicted into the future t + j
j = 2 j = 5 j = 10 j = 24 j = 36
16.04 18.56 19.52
18.79
19.96
16.04 18.59 19.51
19.47
20.08
16.04 18.58 19.56
20.01
20.46

In Figure 4 we plot our short-term predictor ŝt+1 against the long-term estimator lt+1 for the week withhighest observed volume in 2020. In Figure 5 the distribution of the error of the long-term estimator lt+1
vs our predictor ŝt+1 is plotted. Note that with our method the error is centered at 0 with less variance andbias w.r.t. the long-term estimator, which is a desirable property.

450
400

Call Volume

350

Realized
ŝt+1
lt+1

300
250
200
150
100
2020-03-29 2020-03-30 2020-03-31 2020-04-01 2020-04-02 2020-04-03 2020-04-04 2020-04-05 2020-04-06 2020-04-07

Figure 4: Comparison between our proposed short-term prediction ŝt+1 vs long-term prediction lt+1 , actual
call arrival is the red line.

Sanabria, Dolan, Johnson, Kepler, Lam, Lelo de Larrea, Mohammadi, Olivier, Quayyum, Sethuraman,
Smyth, and Thomson

700

Nt+1 − lt+1
Nt+1 − ŝt+1

600
500
400
300
200
100
0

−100

−50

0

50

100

150

200

250

Figure 5: Histogram of the long-term prediction error (Nt+1 − lt+1 with mean -12.40 and standard deviation
34.23) compared with the error of our method (Nt+1 − ŝt+1 with mean -0.04 and standard deviation 18.84).
4

CONCLUDING REMARKSIn this paper we characterized a model and an algorithm to improve a long-term historical prediction ofarrivals, that is both tractable and easy to implement. This model can be used in multiple settings to forecast
demand near real-time. Moreover, the model provides a complete distribution of the prediction error thatcan be used to build confidence intervals or do robust allocation of resources based on the distribution of
the prediction error Nt+1 − lt+1 = εt+1 = µt+1 + σt+1 Zt+1 .
An important use case of the model is the prediction of call arrivals to an EMS call center, wheresuch model could direct efforts in mitigating the effects of increases in the call volume and ensuring thereliability of operations during such events (by either increasing staffing or changing the processing of calls
adaptively relative to the demand).
Further research effort is needed in finding F(t−k):t measurable approximations of µt+1 , σt+1 andcharacterizing the convergence of data-driven procedures in their calibration relative to the expected regimechange time Eτ, the assortment of intensities {λi } and their sampling probabilities beyond the time seriesmodel presented.
ACKNOWLEDGMENTS
We are grateful to the members of the Bureau of Management Analysis and Planning at the FDNY for the
many enlightening conversations and for providing their expert knowledge and data on the EMS system.
We appreciate the guidance of Chief Denise Werner and Chief Ian Swords. We gratefully acknowledge
support from Google and the Tides Foundation under the grant “EMS Resource Deployment Modeling”
and the Columbia University Urban Technology Pilot Award.
REFERENCES
Bollerslev, T. 1986. “Generalized autoregressive conditional heteroskedasticity”. Journal of econometrics 31(3):307–327.
Brockwell, P. J., and R. A. Davis. 2016. Introduction to time series and forecasting. Springer.
Glynn, P. W., and Z. Zheng. 2019. “Estimation and inference for non-stationary arrival models with a linear trend”. In 2019
Winter Simulation Conference (WSC), 3764–3773. IEEE.

Sanabria, Dolan, Johnson, Kepler, Lam, Lelo de Larrea, Mohammadi, Olivier, Quayyum, Sethuraman,
Smyth, and Thomson
Ibrahim, R., and P. L’Ecuyer. 2013. “Forecasting call center arrivals: Fixed-effects, mixed-effects, and bivariate models”.
Manufacturing & Service Operations Management 15(1):72–85.
Ibrahim, R., H. Ye, P. L’Ecuyer, and H. Shen. 2016. “Modeling and forecasting call center arrivals: A literature survey and acase study”. International Journal of Forecasting 32(3):865–874.
Kim, S.-H., and W. Whitt. 2014. “Are call center and hospital arrivals well modeled by nonhomogeneous Poisson processes?”.
Manufacturing & Service Operations Management 16(3):464–480.
L’Ecuyer, P., K. Gustavsson, and L. Olsson. 2018. “Modeling bursts in the arrival process to an emergency call center”. In
2018 Winter Simulation Conference (WSC), 525–536. IEEE.
Matteson, D. S., M. W. McLean, D. B. Woodard, S. G. Henderson et al. 2011. “Forecasting emergency medical service callarrival rates”. Annals of Applied Statistics 5(2B):1379–1406.
Morgan, L. E., B. L. Nelson, A. C. Titman, and D. J. Worthington. 2019. “A spline-based method for modelling and generatinga nonhomogeneous Poisson process”. In 2019 Winter Simulation Conference (WSC), 356–367. IEEE.
Morgan, L. E., A. C. Titman, D. J. Worthington, and B. L. Nelson. 2016. “Input uncertainty quantification for simulation
models with piecewise-constant non-stationary Poisson arrival processes”. In 2016 Winter Simulation Conference (WSC),
370–381. IEEE.
Oreshkin, B. N., N. Réegnard, and P. L’Ecuyer. 2016. “Rate-based daily arrival process models with application to call centers”.
Operations Research 64(2):510–527.
Taylor, S. J., and B. Letham. 2018. “Forecasting at scale”. The American Statistician 72(1):37–45.

AUTHOR BIOGRAPHIESELIOTH SANABRIA is a Ph.D. student in Operations Research at Columbia University. He is primarily interested in theinterplay between simulation, machine learning and optimization from a probabilistic point of view, as well as their broad
applications in healthcare to improve patients outcomes. His email address is m.elioth@columbia.edu.
EDWARD DOLAN is the Deputy Commissioner for Strategic Initiatives and Policy at the Fire Department of the City of
New York (FDNY). His email address is edward.dolan@fdny.nyc.gov.
NICHOLAS JOHNSON is the Director of Operations Research for the FDNY. Prior to joining the FDNY, he was a Postdoctoral
Associate at NYU’s Marron Institute of Urban Management where his research focused on modeling real-time urban populations
using mobility data. He earned a Ph.D. in Urban Science from the University of Warwick in the United Kingdom and holds a
Master’s degree from NYU’s Interactive Telecommunications Program where he used physical computing and interaction design to
explore the impact and pervasiveness of waste streams in urban environments. His email address is nicholas.johnson@fdny.nyc.gov.
TIMOTHY KEPLER is the Director of Data Quality for the FDNY. He has a Master’s degree in Public Administration fromthe City University of New York at Baruch College. For the past five years he has worked as an analyst for the FDNY, working
intensively with data from computer aided dispatch systems. His email address is timothy.kepler@fdny.nyc.gov.
HENRY LAM is an Associate Professor in the Department of Industrial Engineering and Operations Research at Columbia
University. His research focuses on simulation and optimization under uncertainty. His email address is khl2114@columbia.edu.
ENRIQUE LELO DE LARREA is a Ph.D. student in Operations Research at Columbia University. Before joining Columbia,
he worked as a credit risk analyst at BBVA Mexico. He holds a double Bachelor’s degree in Applied Mathematics and Actuarial
Science from ITAM and a Master’s degree in Operations Research from Columbia. His research interests include stochastic
simulation, applied probability, and financial engineering. His email address is enrique.lelodelarrea@columbia.edu.
SEVIN MOHAMMADI is currently pursuing a Ph.D. degree in Civil Engineering and Engineering Mechanics at Columbia
University. Prior to her current position, she obtained a Master’s degree in Civil Engineering, majoring in Transportation,from the University of Tennessee Knoxville. Her research interests include data-driven and machine learning approaches inTransportation Engineering. Her email address is sm4894@columbia.edu.
AUDREY OLIVIER is currently an Associate Research Scientist in Civil Engineering at Columbia University, and will bejoining the University of Southern California as an Assistant Professor in fall 2021. She holds a Ph.D. in Civil Engineeringand Engineering Mechanics from Columbia University and a Diplôme d’Ingénieur from Ecole Centrale de Nantes, France. Her
research interests revolve around probabilistic data analytics and physics-based modeling for civil engineering applications. Her
email address is audreyol@usc.edu.

Sanabria, Dolan, Johnson, Kepler, Lam, Lelo de Larrea, Mohammadi, Olivier, Quayyum, Sethuraman,
Smyth, and Thomson
AFSAN QUAYYUM is a Data Scientist at the Bureau of Management Analysis and Planning at the FDNY. He has an M.S.
in Mathematics and a B.S. in Mathematics with a minor concentration in Applied Physics from New York University. His
work focuses on implementing statistical learning techniques for inference and prediction to help EMS and Fire Operations.
His email address is afsan.quayyum@fdny.nyc.gov.
JAY SETHURAMAN is a Professor of Industrial Engineering and Operations Research at Columbia University. Currently, heserves as the chair of the IEOR department at Columbia. His research interests are in discrete optimization and applications,game theory, mechanism design, and applied probability. His email address is jay@ieor.columbia.edu.
ANDREW SMYTH is the Carleton Professor of Civil Engineering & Engineering Mechanics and also serves as the Co-Chair of
the Smart Cities Center of the Data Science Institute at Columbia University. His research focuses on infrastructure monitoring,
dynamic system identification and modeling. He received his Ph.D. in Civil Engineering from the University of Southern
California as well as an M.S. in Electrical Engineering, an M.S. from Rice University, and a Sc.B. and A.B. from Brown
University. His email address is smyth@civil.columbia.edu.
KATHLEEN THOMSON is the Assistant Commissioner for the Bureau of Management Analysis and Planning at the FDNY.
Her email address is kat.thomson@fdny.nyc.gov.