16 Call Center Model

1 Introduction

1 Outline

Call centers are important and rapidly developing commercial activities. Call centers serve customers by phone, by fax and by Internet. There are different call centers depending on their objectives and environments. However, investigating the call centers one encounters three problems:

• modeling;
  queuing systems are common tools while modeling call centers,
• optimization;
  one needs methods of stochastic and global optimization are to optimize these systems,
• prediction; parameters, which are needed for modeling and optimization, depends on predictions, predictions involve both the observed data and the expert knowledge.

2 Assumptions, Notations, and Objectives

One models call centers as queuing systems. The modeling system includes features for call rate predictions and service optimization. One assumes that

• incoming calls are united into one stream with the call rate $\lambda$, a call rate is the average number of calls in a time unit,
• there are $m$ servers with the same service rate $\mu$, the service rate is the average number of calls that can be served in a time unit,
• each server can serve any call, one at a time,
• calls are Poisson with the rate $\lambda$,
• service times are exponential with the rate $\mu > 1/m\ \lambda$, in such a case
  $$\begin{eqnarray} \lambda=1/T_{\lambda}, \\ \mu=1/T_{\mu}, \end{eqnarray}$$
  
  
  
  
  where $T_{\lambda}$ is an average time between the calls and $T_{\mu}$ is an average service time,
• arriving calls enter the first available server, if all the servers are busy then the call takes a free waiting place,
• there are $r$ waiting places, that means that a call waits, if there are not more than $r$ other calls waiting, otherwise, the call disappears,
• waiting places are common for all the calls
• the system is stationary, meaning that the parameters $\mu$ and $\lambda$ are constant,
• $c_m$ is the server running cost ( $ per time unit),
• $c_t$ is the customer time cost ( $ per time unit),
• $c_p$ is the lost customer cost ( $ per lost call),
• the optimal number of servers $m=m(c,r)$ minimizes the total cost $C(c,m)$ per time unit, including the server running cost, the customer time cost, and the lost customer cost at fixed parameters $r$ and $c=(c_s,c_t,c_p)$,
• results are presented as the waiting time distribution functions
  $F_{\tau}^{m,r}(t)= P^{m,r}\{\tau <t\}$, where $P^{m,r}\{\tau <t\}$ is the probability that the waiting time $\tau$ will be less than $t$ at fixed number of servers $m$ and waiting places $r$,
• a family of functions $F_{\tau}^{c,r}(t)$ is defined for different parameters $c,r$, these functions are presented in a simple format, assuming that the number of servers $m=m(c,r)$ is obtained by minimizing the total cost $C(c,m,r)$.

2 Calculation of Stationary Probabilities

Under these assumptions [Tijms, 1994], the probability of $k$ calls in the system, including both, waiting calls and calls in services

$$\begin{eqnarray} P_{k}(\mu,\lambda)&=&\cases{\frac{\rho^k}{k!} P_0(\mu,\lambda),... ...rac{\rho^k}{m! m^{k-m}} P_0(\mu,\lambda),&if $m < k \le m+r$,\cr} \end{eqnarray}$$

where $\rho=\lambda/\mu$. The no-calls probability

$$\begin{eqnarray} P_0(\mu,\lambda)=(\sum_{k=0}^{m-1} \frac {\rho^k}{ k!} +\frac {m^m}{m!}\sum_{k=m}^{m+r} (\frac{\rho}{m})^k)^{-1} \end{eqnarray}$$

follows from the condition

$$\begin{eqnarray} \sum_{k=0}^{m+r} P_k(\mu,\lambda)=1. \end{eqnarray}$$

The probability of losing a call at given numbers $m,r,\mu,\lambda$ is

$$\begin{eqnarray} P_{m+r}(\mu,\lambda)= \frac {\rho^{m+r}}{ m! m^r} P_0(\mu,\lambda). \end{eqnarray}$$

The average waiting time

$$\begin{eqnarray} T(m,r,\mu,\lambda)= \sum_{l=0}^{r-1} P_{m+l}(\mu,\lambda)\ t_l,\ r >0, \end{eqnarray}$$

where $P_k(\mu,\lambda)$ is defined by expression (16.3) and $t_l=\frac {l+1}{m \mu}$. The waiting time distribution functions

$$\begin{eqnarray} F_{\tau}^{m,r}(t)= P^{m,r}(\mu,\lambda)\{\tau<t\}=1-\Psi_{\tau}... ...^{m,r}(t)= \sum_{k=m}^{m+r-1} P_{k}(\mu,\lambda)\ P_k\{\tau >t\}, \end{eqnarray}$$

Here $P_k\{\tau >t\}$ denotes the probability that waiting time $\tau$ will exceed some fixed $t$ under condition that there are $k$ calls

$$\begin{eqnarray} P_k\{\tau >t\}= \sum_{s=0}^{k-m} q_s(t),\ k \ge m. \end{eqnarray}$$

Here $q_s(t)$ is the probability that $s$ calls will be served during the time $t$ and

$$\begin{eqnarray} q_s(t)=e^{-m \mu t} \frac {(m \mu t)^s}{s!},\ s=0,1,...,r. \end{eqnarray}$$

3 Asymptotic Expressions

If the number $r$ of waiting places is very large, one considers asymptotic expressions ( $r \rightarrow \infty,\ \rho < m$) as a reasonable approximation. The asymptotic values are denoted by the index $\infty$. They are used to test the software designed for finite $r$. The no-calls probability

$$\begin{eqnarray} P_0^{\infty}(\mu,\lambda )=(\sum_{k=0}^{m-1} \frac {\rho^k}{ k!} +\frac {\rho^m}{(m-1)! (m-\rho)})^{-1}. \end{eqnarray}$$

The average waiting time

$$\begin{eqnarray} T(m,\infty,\mu,\lambda)=\frac {\rho^m}{m!m\mu (1-\rho / m)^2} P_0^{\infty}(\mu,\lambda) . \end{eqnarray}$$

The waiting time distribution functions

$$\begin{eqnarray} F_{m,\infty}(t)=P_{m,\infty}(\mu,\lambda)\{\tau<t\} \nonumber\\... ...{1-\rho /m} e^{-(m \mu - \lambda) t} P_{m}^{\infty}(\mu,\lambda). \end{eqnarray}$$

Here

$$\begin{eqnarray} P_m^{\infty}(\mu,\lambda)=\frac { \rho^m}{m!} P_0^{\infty}(\mu,\lambda) . \end{eqnarray}$$

The asymptotic probability of losing a call is zero.

4 "Surrogate" Services

An incoming call gets a regular service by an agent, if the actual or estimated waiting time is less then $\tau_w$ sec. Otherwise, a call gets a "surrogate" service by voice-mail. For simplicity, one replaces this service system by a system with $r=\tau_w / \tau_s$, where the expected service time is $\tau_s=1 / m \mu$. Then one can use all the expressions for the queuing system with $r$ waiting places. Here, the call served by voice-mail is considered as the "lost" call. One considers calls of different type separately, because they are served by different agents.

5 Call Rate Estimate

The lost calls are not registered, as usual. Therefore, it is difficult to estimate the call rate $\lambda$ directly, if $r$ is limited. Then one uses least square estimates. The square deviations $\Delta (\mu,\lambda)$ between stationary probabilities $P_{k}(\mu,\lambda)$ and their estimates $P_{k}^0$ are minimized.

$$\begin{eqnarray} \Delta (\mu,\lambda)=\sum_{k=0}^{m+r} (P_{k}(\mu,\lambda)-P_{k}^0)^2 . \end{eqnarray}$$

The stationary probabilities that there are $k$ calls in the system are defined by expression

$$\begin{eqnarray} P_{k}(\mu,\lambda)&=&\cases{\frac {\rho^k}{k!} P_0(\mu,\lambda)... ...ho^k}{m! m^{k-m}} P_0(\mu,\lambda) , & if $m \le k \le m+r$.\cr}. \end{eqnarray}$$

The least squares estimate of the call rate

$$\begin{eqnarray} \lambda^o= \arg \min_{\lambda} \Delta (\mu,\lambda). \end{eqnarray}$$

The estimates $P_{k}^0$ are obtained by counting the numbers of waiting calls at different time moments. Additional errors are expected, if the average number of waiting calls in a time unit is counted, instead of moment numbers.

6 Optimization of Number of Servers

The total cost of a service system

$$\begin{eqnarray} C(c,m)=c_m m + c_t T(m,r,\mu,\lambda) + c_p P_{m+r}(\mu,\lambda). \end{eqnarray}$$

One optimizes the server number $m$ by simple comparison of different $m$ within reasonable bounds $m_{min} \le m \le m_{max},\ m_{min} > \lambda / \mu$

$$\begin{eqnarray} m(c,r)=arg \min_{m_{min} \le m \le m_{max}} C(c,m,r). \end{eqnarray}$$

7 Monte Carlo Simulation (MCS)

The explicit steady-state solutions (16.6)-(16.13) are simple and exact, under the assumptions. However, there is no general explicit solution, if the call rate $\lambda$, and/or the service rate $\mu$ depends on time. That often happens in real service systems. The configuration and operating rules of real service systems are too complicated for the exact solution, as usual.

Then the Monte Carlo Simulation (MCS) is needed. We start the discussion of MCS from the steady-state systems keeping all the assumptions described in Section 1.2. That helps to test both the analytical and the Monte Carlo solutions by comparing the results. Later, MCS is extended to more complicated configurations and to time-dependent cases.

1 Event Generation

The events in a queuing system are the moments when a call arrives, when a call enters a server, and when a call leaves the system. To generate events, one needs two types of random number generators. The first type generates times until the next call. The second type generates service times.

Denote by $F_a(t)=P_a\{\tau < t\}$ the distribution function, where $P_a\{\tau < t\}$ is the probability that a random time $\tau$ will be less then $t$. Here $a$ is the expected value of $\tau$. Denote by $\xi \in [0,1]$ the random variable uniformly distributed between zero and one. Then

$$\begin{eqnarray} \tau=F^{-1}(\xi). \end{eqnarray}$$

In exponential cases

$$\begin{eqnarray} F_a(t)=1-e^{-1/a t},\\ \tau=-a ln (1-\xi). \end{eqnarray}$$

Here $a=1/\lambda_s$, if times until the next call are generated.
$a=1/(m_s\mu_s)$, if service times are generated.
The model defines moments when a call enters a server and leaves the system. These moments depend on the specific structure of the system. The model is designed trying to represent the actual operations as realistically as possible. However, the are limits that depend on the available computing power.

Here is a description of an algorithm of modeling and optimization of call centers:

• denote by $t_n$ the "arriving" time of the $n$th call, then $t_{n+1}=t_n+\xi_n(\lambda),\ \xi(\lambda)= - 1/\lambda\ log(1-\eta)$, where $\eta$ is a random number, uniform in $[0,1]$, and $\lambda$ is the average number of calls in a time unit, the moment of the first call moment $t_1=0$ is the start of MCS,
• denote by $\tau_n$ the "departing" time of $n$th call, then $\tau_{n+1}=\tau_n+\xi_n(\mu_n),\ \xi(\mu_n) = - 1/\mu_n\ log(1-\eta)$, where $\mu_n$ is the average number of calls which can be served in a time unit by all $1 \le m_n \le m$ servers operating at the moment $\tau_n$, thus $\mu_n= m_n \ \mu$, and $\mu$ is the average number of calls, served by single server in a unit time,
• count all calls in the system at the moment $t_{n+1}$ including those in $m$ servers and $r$ waiting places by this expression
  $$\begin{eqnarray} I_{n+1}&=&\cases {I_n-1, &if $t_{n+1} > \tau_n$\ \cr \min(I_n+1, m+r) , & if $t_{n+1} < \tau_n$,\cr} \end{eqnarray}$$
  
  
  
  
  where $I_n$ is the number of calls in the system at the moment $t_n$,
• define the number of servers $m_n$ operating at the moment $\tau_n$ by this condition
  $$\begin{eqnarray} m_n=\min (I_l,m),\ t_l \le \tau_n \le t_{l+1}, \end{eqnarray}$$
  
  
  
  

• count calls lost up to the moment $t_{n+1}$ by the expression
  $$\begin{eqnarray} L_{n+1}&=&\cases {L_n+1, &if $I_n+1 > m+r, t_{n+1}<\tau_n$\ \cr L_n , & otherwise,\cr} \end{eqnarray}$$
  
  
  
  
  where $L_n$ is the number of calls lost up to the the moment $t_n$,
• count the waiting time by this counter
  $$\begin{eqnarray} T_{n+1} =T_n+\delta_{n},\\ \delta_{n}=\tau_{n}-t_{n}, \end{eqnarray}$$
  
  
  
  

• count the number of calls $q^k$ that waited less than $t^k,\ t^{k+1} > t^k,\ k=1,...,K$
  $$\begin{eqnarray} q_{n+1}^k&=&\cases {\min(q_n^k+1,m+r), &if $\delta_n < t^k $\ \cr q_{n}^k , & if $\delta_n \ge t^k$,\cr} \end{eqnarray}$$
  
  
  
  

• define the distribution of the call numbers $q^k$ that waited less than $t^k,\ t^{k+1} > t^k,\ k=1,...,K$
  $$\begin{eqnarray} Q(t^k)=1/N\ q_{N}^k, \end{eqnarray}$$
  
  
  
  

• define the optimal number of servers $m$ by simple comparison of different $m$ within some reasonable bounds $m_{min} \le m \le m_{max}$
  $$\begin{eqnarray} m(c,r)=arg \min_{m_{min} \le m \le m_{max}} C(c,m,r), \end{eqnarray}$$
  
  
  
  

• compare the results with the corresponding steady-state explicit solutions, including the average waiting time, the average of lost calls and the waiting time distribution functions.

2 Monte Carlo Errors

Under independence conditions, the standard deviation $\sigma(K)$ of the results $\theta(K)$ of a statistical model depends on the number of observations $K$ as

$$\begin{eqnarray} \sigma(K)=\sigma(1) / \sqrt K. \end{eqnarray}$$

Here $\sigma(1)$ is the "initial" standard deviation

$$\begin{eqnarray} \sigma(1)=\sqrt {E (\theta(1)-\Theta(1))^2}, \end{eqnarray}$$

where $\Theta(1)=E \theta(1)$ and $E$ is the expectation symbol.

In Monte Carlo simulation, $\Theta(1)$ cannot be estimated directly using the results $\theta(1)$ obtained by the 1-th observation. Many observations are needed, to reach a stationary state (one observation means one call). Therefore, we repeat the Monte Carlo simulation $L$ times doing $K$ observations each time. Then the error variance can be estimated as

$$\begin{eqnarray} \sigma_L^2(K)=1/(L-1) \sum_{l=1}^L (\theta_l(K)-1/L\sum_{j=1}^L\theta_j(K))^2. \end{eqnarray}$$

Expression (16.33) is simpler, if the exact solution $\Theta=\lim_{K \rightarrow \infty} \theta(K)$ is known . Then

$$\begin{eqnarray} \sigma_L^2(K)=1/(L-1) \sum_{l=1}^L (\theta_l(K)-\Theta)^2 \end{eqnarray}$$

Therefore, to test a Monte Carlo procedure, one considers a model with exact solution, first. Then one applies MCS to more complicated models.

3 Stopping Monte Carlo

The Monte Carlo errors are important to define stopping rules. The traditional idea is that the computing errors should not be greater then the data gathering errors. In call centers, data errors depend on the errors of call rate predictions, as usual. Thus, the stopping rule is:
- stop at the first observation $k \ge K$ satisfying this inequality:

$$\begin{eqnarray} \sigma(K) < \sigma. \end{eqnarray}$$

Here $\sigma$ is a standard deviation of call rate predictions.

8 Common Waiting

The common waiting system is an example were one applies approximate models, including the Monte Carlo. By common waiting we define a system that serves calls of different type. One represents them as a call-vector $\lambda=(\lambda_1,...,\lambda_l)$. Different calls are served by different servers. However, there are $r(c)$ of common waiting places. There are no simple formulas, such as (16.16) and (16.15). Therefore, we consider two approximate solutions: the analytical "reservation" model, and the Monte Carlo one.

1 Analytical Approximation: Reservation Model

One estimates $\lambda=(\lambda_1,...,\lambda_l)$ by minimizing the square deviations

$$\begin{eqnarray} \Delta_s (\mu,\lambda, r_s)=\sum_{k_s=0}^{m_s+r_s} (P_{k_s}(\mu_s,\lambda_s, r_s)-P_{k_s}^0)^2. \end{eqnarray}$$

Here $P_{k_s}(\mu_s,\lambda_s, r_s)$ is the stationary probability that there are $k_s$ calls in the $s$-th server. This probability is defined by expressions similar to (16.3), assuming that there are $r_s$ waiting places reserved for the calls $\lambda_s$. $r=(r_1,...,r_l)$ is the reservation vector and $\sum_s r_s=r(c)$.

$$\begin{eqnarray} P_{k_s}(\mu_s,\lambda_s, r_s)&=&\cases{\frac {\rho^k_s}{k_s!} P... ...,\lambda_s, r_s) , & if $m_s \le k_s \le m_s+r_s$.\cr}. \nonumber \end{eqnarray}$$

The estimates $P_{k_s}^0$ are obtained by counting the numbers of waiting $s$-calls at different time moments. The least square estimate of the call-vector $\lambda$ is as follows

$$\begin{eqnarray} \lambda^o= \arg \min_{\lambda} \min_r (\sum_s \Delta (\mu_s,\lambda_s, r_s)). \end{eqnarray}$$

The reservation model is simple and clear. However, one must test the reservation assumption by the Monte Carlo model.

2 Statistical Approximation: Monte Carlo Model

3 Call Rate Estimate

The statistical model can be used to estimate the call rates $\lambda$ in the same way as the analytical one. One estimates $\lambda=(\lambda_1,...,\lambda_l)$ by minimizing the square deviations $\Delta (\mu,\lambda,r)$ between the probabilities and their estimates

$$\begin{eqnarray} \Delta (\mu,\lambda,r)=\sum_s \sum_{k_s=0}^{m_s+r} (P_{k_s}(\mu,\lambda, r)-P_{k_s}^0)^2 . \end{eqnarray}$$

Here $P_{k_s}(\mu,\lambda, r)$ is the probabilitythat there are $k_s$ calls in the $s$-th server. The estimates $P_{k_s}^0$ are obtained by counting the numbers of waiting $s$-calls at different time moments.

The least squares estimate of the call-vector $\lambda$

$$\begin{eqnarray} \lambda^o= \arg \min_{\lambda} \Delta (\mu,\lambda, r)). \end{eqnarray}$$

The statistical model needs considerable computing power.

4 Testing Analytical Approximation

The analytical approximation (see Chapter 8.1) is simpler. However, it is based on the reservation assumption. The statistical model may represent most of the important factors. However, it takes a long time to "filter out" the random deviations. In the on-line operations, this time is too long, as usual.

To test analytical approximation, one compares estimates of call rates $\lambda=(\lambda_1,...,\lambda_l)$ obtained by both the statistical and the reservation models. A small deviation means that the reservation model is acceptable. The large deviation means that the analytical approximation should be improved.

9 Time-Dependant Cases

Consider the MCS when the call rates $\lambda=\lambda(t)$, the number of servers $m=m(t)$, and the number of waiting places $r=r(t)$ depend on time. Therefore, all the results should be represented as functions of time. One modifies algorithms, similar to those in the previous chapter, by including the time factor:

• denote by $t_n$ the "arriving" time of the $n$th call, then $t_{n+1}=t_n+\xi_n(\lambda(t)),\ \xi(\lambda(t))= - 1/\lambda(t)\ log(1-\eta)$, where $\eta$ is a random number uniform in $[0,1]$ and $\lambda$ is the average number of calls arriving in a time unit that includes the moment $t$,
• denote by $\tau_n$ the "departing" time of $n$th call, then $\tau_{n+1}=\tau_n+\xi_n(\mu_n),\ \xi(\mu_n) = - 1/\mu_n\ log(1-\eta)$ where $\mu_n$ is the average number of calls which can be served in a time unit by all the $1 \le m_n \le m(\tau_n)$ servers operating at the moment $\tau_n$, thus $\mu_n= m_n \ \mu$ and $\mu$ is the average number of calls served by single server in unit time,
• count the number of all the calls in the system at the moment $t_{n+1}$, including those in $m(t_{n+1})$ servers and $r(t_{n+1})$ waiting places
  $$\begin{eqnarray} I_{n+1}&=&\cases {I_n-1, &if $t_{n+1} > \tau_n$\ \cr \min(I_n+1, m(t_{n+1})+r(t_{n+1})) , & if $t_{n+1} < \tau_n$,\cr} \nonumber \end{eqnarray}$$
  
  
  
  
  where $I_n$ is the number of calls in the system at the moment $t_n$,
• count the number of servers $m_n$ operating at the moment $\tau_n$
  $$\begin{eqnarray} m_n=\min (I_l,m(\tau_n)),\ t_l \le \tau_n \le t_{l+1}, \end{eqnarray}$$
  
  
  
  

• count the number of calls lost up to the moment $t_{n+1}$
  $$\begin{eqnarray} L_{n+1}&=&\cases {L_n+1, &if $I_n+1 > m(t_{n+1})+r(t_{n+1}), t_{n+1}<\tau_n$\ \cr L_n , & otherwise,\cr} \nonumber \end{eqnarray}$$
  
  
  
  
  where $L_n$ is the number of calls lost up to the the moment $t_n$,
• define the waiting time
  $$\begin{eqnarray} T_{n+1} =T_n+\delta_{n},\\ \delta_{n}=\tau_{n}-t_{n}, \end{eqnarray}$$
  
  
  
  

• count the number of calls s $q^k$ that waited less than $t^k,\ t^{k+1} > t^k,\ k=1,...,K$
  $$\begin{eqnarray} q_{n+1}^k&=&\cases {\min(q_n^k+1,m(t_{n+1})+r(t_{n+1})), &if $\... ...a_n < t^k $\ \cr q_{n}^k , & if $\delta_n \ge t^k$,\cr} \nonumber \end{eqnarray}$$
  
  
  
  

• define the distribution of the call numbers $q^k$ that waited less than $t^k,\ t^{k+1} > t^k,\ k=1,...,K$
  $$\begin{eqnarray} Q(t^k)=1/N\ q_{N}^k, \end{eqnarray}$$
  
  
  
  

• optimize the number of servers $m(t)$ at different times $t$ by solving the corresponding stochastic scheduling problem
• define the difference between the steady-state and dynamic solutions by comparing the average results with the corresponding steady-state explicit solutions, including the waiting time distribution functions $F_{\tau}^{m,r}(t)= P^{m,r}\{\tau <t\}$, where $P^{m,r}\{\tau <t\}$ is the probability that the waiting time $\tau$ will be less than $t$,
• define the iteration number $N$ by Monte Carlo experimentation.

1 Simple Example

Consider this simple example

$$\begin{eqnarray} \lambda(t)=a+b\ sin(\phi+ 2 \pi \omega t). \end{eqnarray}$$

Here $\omega=1/24$ and $\phi = 8/24$. $t$ is caunted in hours. $a$ is average number of calls in 24 hours. $b$ is amplitude of daily variation of calls, $r(t)=r$.

$$\begin{eqnarray} m(t)&=&\cases {m^0, &if $sin(\phi+ 2 \pi \omega t ) < 0$, \cr m^1 , & if $sin(\phi+ 2 \pi \omega t) > 0$. \cr} \end{eqnarray}$$

One optimizes the numbers $m^0$ and $m^1$ of two different servers considering two periods. This is a short description of the optimization problem:

• there are two optimization parameters $m^0$ and $m^1$ restricted by bounds $m_{min} \le m^0 \le m_{max}^0$ and $m_{min}^1 \le m^1 \le m_{max}^1$,
• one may obtain the optimal pair by direct comparison of all the feasible pairs
  $$\begin{eqnarray} m*=(m^{*0},m^{*1})=arg \min_{m^0,m^1} C(c^0,c^1,m^0,m^1,r), \end{eqnarray}$$
  
  
  
  

• the total system cost is defined as
  $$\begin{eqnarray} C(c^0,c^1,m^0,m^1,r)=\nonumber \\ c^0 m^0 + c^1 m^1 + c_t T(m^... ...\omega) +\nonumber \\ c_p P_{m^0,m^1,r}(\mu, a, b,\phi, \omega). \end{eqnarray}$$
  
  
  
  

10 Call Rate Predictions

1 Introduction

Consider ideas and models that supplement the ones described in chapter 15. The call rate depends on many factors (see Figures 16.1 and 16.2). Therefore, one applies the extended version of the Auto-Regression-Moving-Average (ARMA) model and software (see Chapter 15).

Figure 16.1: A fragment of the the DATAFILE 'call.data' including the number, the date, the call rate (repeated twice) , the real-valued external factor (repeated twice) and the indicators of eight Boolean external factors.

Figure 16.2: Eight Boolean external factors related to call rate.

The results show that one does not improve predictions by including the external factors directly into the ARMA model. A reason is that in the ARMA framework it is difficult to estimate the delay time (see expression (15.22)) and the duration of impacts Special Events (SE). The additional difficulty is that most of SE are rare (see Figure 16.1). Therefore, columns of date file are filled mostly by zeros. Most of the SE indicated in Figure 16.2 are predictable. For example, one can predict factor 6 ( public holidays, traditional celebrations) and factor 7 (last day for ordering) exactly. Other SE can be predicted approximately. The knowledge of future values of external factors helps to predict the main one. However, the present software version cannot use this possibility yet.

Here, we consider call rates $\lambda(t)$ as a sum of two stochastic functions

$$\begin{eqnarray} \lambda(t)=z(t)+\nu(t). \end{eqnarray}$$

Often the rate $\lambda(t)$ has several components $\lambda(t)=(\lambda_1(t),...,\lambda_l(t))$ corresponding to different types of calls. In expression (16.48) $\nu(t)$ denotes a "stationary" component which is described by the ARMA model. A "non-stationary" component is denoted by $z(t)$. This way we separate the stationary part from the non-stationary one.

The theoretical analysis of non-stationary stochastic functions is difficult. Therefore, the separation of non-stationary component is important. The non-stationary component is defined by a local expert, as usual. The expert applies his knowledge while using the previous data and making the future predictions. Thus, we call $z(t)$ as the "expert" component and $\nu(t)$ as the "statistical" one. The estimate of the statistical component $\nu(t)$ is investigated in Chapter 15, considering the ARMA model. An alternative is "Scale" models. They predict the expert component $z(t)$ by estimated scales. The scales express the differences between different SE and different times.

2 Call Rate Prediction by Scale Models

The ARMA models considered in Chapter 11 reflect non-stationarity by eliminating unstable parameters applying the structural stabilization techniques. This way some data and some model parameters $a_i,\ b_j$ are eliminated. The estimates of the remaining ones are obtained. Therefore, one may regard ARMA parameters $a_i,\ b_j$ as some scales, too. These scales reflect the influence of the corresponding data sets by minimization of prediction errors. No expert knowledge is involved.

Now we shall consider the scale models, where the expert opinion is involved. This is done by choosing data sets and scales reflecting the expert opinion. Two versions of scales models are considered.

In the first model it is supposed that the prediction $z_i$ is a product of the present call rate $z_{i-1}$ and some scale $s_i$

$$\begin{eqnarray} z_i=z_{i-1} \ s_i. \end{eqnarray}$$

Here $z_i$ is a predicted call rate, called a "prediction." $z_{i-1}$ is an observed call rate, called a "data set." The index $i-1$ defines the data set used to predict $z_i$. The parameter $s_i$ is the "scale" to predict $z_i$ by the data set $z_{i-1}$. Expression (16.49) is a time scale model. In this model the scale is estimated as

$$\begin{eqnarray} s_i=\frac {z_{p(i)}}{z_{p(i-1)}}. \end{eqnarray}$$

In (16.50) the subscript $p(i)$ defines a period which precedes $i$. The subscript $p(i-1)$ denotes a period preceding $i-1$. The term "preceding" means the nearest previous period of the same type. For example, Saturdays, Sundays, holidays, Christmas weeks, sporting events, days of marketing messages, e.t.c., all are different types.

In the second model

$$\begin{eqnarray} z_i=z_{p(i)} \ s_{p_i}. \end{eqnarray}$$

Here the scale is estimated as

$$\begin{eqnarray} s_{p(i)}=\frac {z_{i-1}}{z_{p(i-1)}}. \end{eqnarray}$$

Both expressions(16.49) and (16.51) predicts the same call rates. That means that one obtains the same results by the event scales (16.50) and the time scales (16.52). Only the interpretation differs.

In the vector case, expressions (16.49) and (16.51) predict different call rate graphs. One represents these graphs as vectors $z_i=(z_{ij}, \ j=1,...,J)$. The components $z_{ij}$ denote call rates of different parts $j$ of a period $i$ For example, hours if we predict next day call rates.

The assumption of the first model (16.49) is that the next day graph $z_i=(z_{ij}, \ j=1,...,J)$ is equal to the present one $z_{i-1}=(z_{ij},\ j=1,...,J)$ multiplied by scales $s_i$

$$\begin{eqnarray} s_i=\frac {Z_{p(i)}}{Z_{p(i-1)}}. \end{eqnarray}$$

Here $Z_i=1/J \sum_{j=1}^J z_{ij}$ denotes the average call rates of the period $i$. Expression (16.53) means that the shapes of the next and the present graphs remains the same. Only the scales differ.

The second model (16.51) assumes that the next day graph is equal to the graph of the preceding period $z_{p(i)}=(z_{p(i),j},\ j=1,...,J)$ multiplied by scales

$$\begin{eqnarray} s_{p(i)}=\frac{Z_{i-1}}{Z_{p(i-1})}. \end{eqnarray}$$

This means that the shapes of the next graph and the preceding one are the same. Only the scales differ. The numerator $Z_{p(i)}$ of scales $s_i$ of the first model depends on the Special Events (SE) of the next period $i$. We call the first model (16.49) as the event scale model.

The numerator $Z_{p(i)}$ of scales $s_i$ of the first model depends on the Special Events (SE) of the next period $i$. By definition, the preceding period $p(i)$ depends on the type of next one. The period type is defined by the Special Event (SE) which is active in the period. Examples of SE are Saturdays, Sundays, holidays, Christmas weeks, sporting events, days of marketing messages e.t.c..

The numerator $Z_{i-1}$ of scales $s_{p(i)}$ of the second model depends on the average call rate now $i-1$. Therefore, we call the second model (16.51) as the time scale model.

Using the event scale model (16.49), the next day event SE defines the next day scale. The shape of the next day graph remains the same as today. Therefore, this model reflects the changes of the graph shapes without delay.

Using the time scale model (16.51), the shape of the next day graph is the same as the shape of the preceding one. Here the graph changes, but with some delay. The interval between the next period and the previous period of the same type defines the delay time.

3 Expert Model, Event Scale Version

1 Definition of Special Events

Here the call rate in the next period (16.49) is the same as in the present one multiplied by the scale $s_i$. The scale depends on the Special Event (SE) which is expected to be active in the next period. The periods where no special events are active, we consider as periods with "neutral" special events. Scales $s_i$ determine the impact of SE to the call rate $z_i$ of next period.

The scale $s_i$ reflecting the impact of a special event is estimated using empirical data by expressions (16.50), (16.53) and/or directly by expert opinions. The impact of some special events, such as Saturdays, Sundays, holidays, Christmas week e.t.c.. is well defined and instant. One calls them Instant Special Events (ISE).

The impact of some others, such as marketing messages, is delayed. Call rates react to these events with some delay. Therefore, one calls them Delay Special Events (DSE). The starts and the durations of DSE are fixed. The delay times $d$ and the durations $\tau$ of DSE impacts are estimated using the observed data and expert knowledge. In some periods several DSE may be active simultaneously. One calls these Multiple Special Event (MSE).

The starts and the durations of ISE and their impacts are equal to starts and durations of the corresponding periods. Beginnings and durations of DSE impacts are not necessarily equal to the beginnings and durations of the periods. Therefore, a special sorts of DSE appear, called Partial Special Events (PSE). In PSE the impact of a special event is active only during some part of the period. The delay times $d$, the durations $\tau$ of the impacts of DSE, and the scales $s$ defining the impact of all the special events, are estimated using the available data and the expert knowledge.

2 Procedures

We consider three expert procedures and the method of least squares, while defining estimates of scales $s_i$. The first procedure is the empirical one. The scales $s$ describing the impact of all sorts of events, including the partial and the multiple ones, are estimated by expressions (16.50) or (16.53). These expressions define the relation of the call rates in the pair of periods $p(i-1)$ and $p(i)$. A period $p(i-1)$ precedes the present period $i-1$. Denote by $p(i)$ a period that precedes the next one, denoted $i$. The numerator $z_{p(i)}$ is the call rate of the period preceding the next one. The denominator $z_{p(i-1)}$ is the call rate of the period preceding the present one.

Using the empirical procedure, partial and multiple special events are considered as different SE. Here scales $s$ are defined separately for each of them. Therefore, the empirical procedure is convenient, if the data is available for all types of periods, including the partial and multiple ones. In such a case, estimate of the delays $d$ and durations $\tau$ of DSE is not needed. It is supposed that an identical situationis in the available data. One needs large data sets for that.

The second direct procedure is the subjective one. All unknown parameters of the expert model, such as the scales $s$, the delay times $d$ and the duration times $\tau$, are defined by the expert opinion. This is a reasonable way to start a new system, when no data is available.

The third direct procedure is a mixture of the empirical and the subjective ones. We use empirical estimates, if the corresponding data is available. Otherwise, one uses subjective estimates. Using the method of least squares, we search for such scales $s$, such delay and duration times $d,t$ that minimize the sum of squared deviations. The deviations are differences between the call rates predicted by the expert model and the observed call rates $z_{ij}$. All four procedures are applied to different times: hours, days, weeks and seasons. First, we consider the scalar case.

3 Scalar Prediction

In the scalar case, one predicts a single number, the average call rate of the next period. In the formal terms, one considers the call rate function $z_i$ as a sample of some non-stationary stochastic process. Four different periods; hours, days, weeks, and seasons, are considered independently. Let us start by describing hours.

1 Hourly Prediction

Using the event scale model, a call rate $z_i$ of the next hour $i$ is expressed this way

$$\begin{eqnarray} z_i=z_{i-1}\ s_i. \end{eqnarray}$$

Here $z_{i-1}$ is the present call rate. $s_i$ is a scale of the next hour $i$. The hourly scales $s_i$ differ depending on the impact of SE felt at the day $i$. For example, working days, weekends, Christmas days, marketing messages days, e.t.c.

$$\begin{eqnarray} s_i=\frac {z_{p(i)} }{ z_{p(i-1)}} . \end{eqnarray}$$

Here $p(i)$ is a subscript of the hour preceding the next one. $p(i-1)$ is a subscript of the hour preceding the present one. Following examples illustrate what one means by the term "preceding."

2 Examples of Multiple Special Events (MSE)

Suppose that the impact of the special event "marketing messages send" starts at the next hour $i$, that is the second one of a working day. Assume that the present hour $i-1$ is the first one of a working day. Then

$$\begin{eqnarray} s_i=\frac {z_{p(i)} }{ z_{p(i-1)}}. \end{eqnarray}$$

Here $p(i)$ is the subscript of the most recent previous second hour of a working day with the same impact SE. $p(i-1)$ is the subscript of the most recent previous first hour of a working day with no SE. In this example, we consider a multiple SE as two single SE. The first single SE means the first working hour without the impact of marketing messages. The second single SE defines the second working hour with the impact of marketing messages.

Another way is to consider the impacts of those two SE separately. Then one defines the scale $s_i$ as a product of two scales

$$\begin{eqnarray} s_i=s^1_i \ s^2_i . \end{eqnarray}$$

Here

$$\begin{eqnarray} s^1_i=\frac {z_{p^1(i)} }{ z_{p^1(i-1)}}, \end{eqnarray}$$

where $p^1(i)$ is the index of the first hour and $p^1(i-1)$ is the index of the second hour. Both indices denote the previous working day.

$$\begin{eqnarray} s^2_i=\frac {z_{p^2(i)} }{ z_{p^2(i-1)}}, \end{eqnarray}$$

where $p^2(i)$ is the index of the most recent hour with "marketing" and $p^1(i-1)$ is the index the most recent hour with no "marketing." Expression (16.58) needs less data. However, it is based on the assumption that scales are multiplicative. Expression (16.58) involves expert opinion by assuming multiplicativity of scales. Expression (16.57) is an example of empirical approach; multiple SE are considered as different special events.

3 Examples of Partial Special Events (PSE)

The only difference from the previous example is that the impact of "marketing messages" starts at the 20-th minute of the next hour $i$. Supose that PSE is a special sort of SE with the scale

$$\begin{eqnarray} s_i=\frac {z_{p(i)} }{ z_{p(i-1)}} . \end{eqnarray}$$

Here $p(i)$ is the subscript of the most recent previous day when the "marketing" starts at the 20-th minute of second working hour. $p(i-1)$ is the subscript of the most recent previous first hour of a working day with no marketing. Another way is to consider the impact of PSE by setting the scale $s_i$ of marketing messages, in proportion to the duration of their impact

$$\begin{eqnarray} s^{20}_i=2/3\ s_i , \end{eqnarray}$$

where

$$\begin{eqnarray} s_i=\frac {z_{p(i)} }{ z_{p(i-1)}}. \end{eqnarray}$$

Here $p(i)$ is the subscript of the most recent previous day when the "marketing" starts at the second working hour. $p(i-1)$ is the subscript of the most recent previous first hour of a working day with no marketing.

Expression (16.62) needs less data. However, it is based on the assumption that scales are proportional to the impact duration. Expression (16.62) involves expert opinion by assuming the proportionality of scales. Expression (16.61) is an example of the empirical approach where a PSE is considered as a different special event.

4 Method of Least Squares

Denote by $\psi_i(s,d,\tau)$ the call rate prediction by an expert model with fixed parameters $s,d,\tau$. One predicts call rates of the next period $i$, using data up to the period $i-1$. The parameters of the model include scales $s$ and delay and duration times $d,\tau$. We minimize the sum

$$\begin{eqnarray} \min_{s,d,\tau} \sum_{i=T_0}^T(z_i- \psi_i(s,d,\tau))^2. \end{eqnarray}$$

Here $T$ is the end and $T_0$ is the beginning of the "learning" set.

One minimizes (16.64) by various global and local optimization methods. That means that estimates of parameters of expert scale models are obtained in the same way as estimates of ARMA parameters (see Section 15). The difference of (16.64), is that parameters $s,d,\tau$ of the expert model EM are optimized, instead of ARMA parameters $b$. The optimization problem (16.64) is very difficult. Here, the number of variables is great and the objective function is multi-modal.

5 Vector Prediction, Event Scale Version

One often predicts the graph of call rates of the next perid. The graph is represented as call rates of different parts of the next period. These call rates are components of some vector. That justifies the term "vector prediction".

The only difference from the scalar prediction is that we predict not a single call rate $z_i$ of the next period $i$, but some vector $z_i=(z_{ij}, \ j=1,...,J)$. Here the component $z_{ij}$ definines the call rate of the part $j$ of the period $i$. First, we consider daily predictions, where index $i$ denotes days and the index $j$ denotes hours.

1 Daily Vector Prediction

We predict the hourly call rates for the next day. The predicted call rate $z_{ij}$ of the hour $j$ of the day $i$ is

$$\begin{eqnarray} z_{ij}=z_{i-1,j}\ s_{i}. \end{eqnarray}$$

Here $z_{i-1,j}$ is the call rate of the hour $j$ of the present day $i-1$. $s_{i}$ is a scale of the next day $i$. The scales $s_{i}$ depend on the impact of SE felt at the day $i$. For example, working day, Saturday, Sunday, Christmas day, marketing day, e.t.c.

$$\begin{eqnarray} s_{i}=\frac {Z_{p(i)} }{ Z_{p(i-1)}} . \end{eqnarray}$$

Denote by $p(i)$ the subscript of the nearest previous day, similar to the next one. Denote by $p(i-1)$ the subscript of the nearest previous day, similar to the present one. The meaning of "nearest previous day similar to" is illustrated by following examples.

2 Examples of MSE

Suppose that the impact of the SE "marketing messages send" starts at the next day $i$. That is Tuesday. The present day $i-1$ is Monday with no additional special events. Then the scale

$$\begin{eqnarray} s_i=\frac {Z_{p(i)} }{ Z_{p(i-1)}}. \end{eqnarray}$$

Here $p(i)$ is the subscript of the most recent previous Tuesday with the same SE. $p(i-1)$ is the subscript of the most recent previous Monday with no SE. In this example we considered a multiple SE as two single SE: Monday without marketing, and Tuesday with marketing.

Another way is to consider the impacts of those two SE by defining the scale $s_i$ as a product of two scales.

$$\begin{eqnarray} s_i=s^1_i \ s^2_i \end{eqnarray}$$

Here

$$\begin{eqnarray} s^1_i=\frac {Z_{p^1(i)} }{ Z_{p^1(i-1)}}, \end{eqnarray}$$

where $p^1(i)$ is the index of Monday, $p^1(i-1)$ is the index of Tuesday of the previous week, and

$$\begin{eqnarray} s^2_i=\frac {Z_{p^2(i)} }{Zz_{p^2(i-1)}}. \end{eqnarray}$$

$p^2(i)$ is the index of the most recent day with "marketing." Index $p^1(i-1)$ denotes the most recent day with no "marketing." Expression (16.67) is an example of empirical approach. Here multiple SE are regarded as different special events. Expression (16.68) needs less data but is based on an expert assumption that scales are multiplicative.

3 Examples of PSE

One feels the impact of the special event, called "marketing," next day $i$ on the sixth working hour. That is the only difference from the previous example. Let us to consider PSE as another SE. Then

$$\begin{eqnarray} s_i=\frac {Z_{p(i)} }{ Z_{p(i-1)}}. \end{eqnarray}$$

Here $p(i)$ is the subscript of the most recent previous week when the "marketing" starts Tuesday at the sixth working hour. $p(i-1)$ is the subscript of the most recent previous Monday with no marketing.

Another way is to regard PSE as a reduced SE. Then, the scale $s_i$ is proportional to the impact duration

$$\begin{eqnarray} s^{6}_i=1/4 \ s_i . \end{eqnarray}$$

Here

$$\begin{eqnarray} s_i=\frac {Z_{p(i)} }{Zz_{p(i-1)}}, \end{eqnarray}$$

$p(i)$ is the subscript of the most recent previous week, such that "marketing" starts Tuesday on the sixth hour, and $p(i-1)$ is the subscript of the most recent previous Monday with no marketing.

Expression (16.71) is an example of empirical approach when a PSE is regarded as a different special event. Expression (16.72) needs less data. However, it is based on an expert assumption that scales are proportional to the impact duration.

6 Method of Least Squares

Denote by $\psi_{ij}(s,d,\tau)$ the call rate prediction by an expert model with fixed parameters $s,d,\tau$. One predicts call rates of the part $j$ of the next period $i$ using data up to the period $i-1$. The parameters of the expert model include scales $s$ and delay and duration times $d,\tau$. All special events, including DSE and PSE, are regarded. We minimize the sum

$$\begin{eqnarray} \min_{s,d,\tau} \sum_{i=T_0}^T \sum_{j=1}^J (z_{ij} - \psi_{ij}(s,d,\tau))^2. \end{eqnarray}$$

Here $T$ is the end and $T_0$ is the beginning of the "learning" set.

Define the expression (16.74) as a function 'fi' in the file 'fi.C'. Then one minimizes (16.74) by various global and local optimization methods. That means that estimates of parameters of expert scale models are obtained in the same way as estimates of ARMA parameters (see Section 15).

4 Time Scale Version, Vector Prediction

Scalar predictions by time scales are identical to those using event scales. Vector predictions by time scale models (16.51) are different from those using event scale models (16.49). Here the shape $z_I=(z_{ij},\ j=1,...,J)$ of the next day graph remains the same as that on the preceding day $p(i)$. Only scales $s_{i-1}$ change. This follows from (16.52). Scales $s_{i-1}$ of the present period $i-1$ reflect changes of average call rates between the present period and the preceding one. The scales are estimated using empirical data and/or expert opinion. The definitions and examples of all the special events are the same as in the Event Scale case (see Section 10.3).

1 Procedures

We shall consider four procedures to obtain estimates: empirical, subjective, mixed, and least squares. The estimates are similar to those in the Event Scale Model (see Section 10.3). Only scales and preceding periods are different. One predicts a vector $z_i=(z_{ij}, \ j=1,...,J)$. Here $z_{ij}$ means the call rate of the part $j$ of the period $i$. First, we consider daily predictions. Then the index $i$ denotes days, and the index $j$ denotes hours.

2 Daily Vector Prediction

We predict the hourly call rates of the next day. The predicted call rate $z_{ij}$ of the hour $j$ of the next day $i$ is

$$\begin{eqnarray} z_{ij}=z_{p(i)j}\ s_{i-1}. \end{eqnarray}$$

Here $z_{p(i)j}$ is the call rate of the hour $j$ of the day preceding the next one. $s_{i-1}$ is a scale of today $i-1$. Scales $s_{i-1}$ depend on the impact of SE felt at the day $i$. For example, working day, Saturday, Sunday, Christmas day, marketing day, e.t.c.

$$\begin{eqnarray} s_{i}=\frac {Z_{i-1} }{ Z_{p(i-1)}}. \end{eqnarray}$$

Here $p(i-1)$ is a subscript of the day preceding the present one. The term "preceding" is illustrated by the following examples.

1 Examples of MSE

Suppose that the impact of the SE "marketing messages send" begins and ends today $i-1$, which is Monday. The next day $i$ is Tuesday with no additional special events. Then the scale

$$\begin{eqnarray} s_i=\frac {Z_{i-1)} }{ Z_{p(i-1)}}. \end{eqnarray}$$

The index $p(i-1)$ denotes the day, which precedes the present one. Here, we regard a multiple SE as two single SE. First single SE is Monday with marketing. The second single SE is Tuesday without marketing.

Another way is to regard the impacts of those two SE by defining a scale $s_{i-1}$ as a product of two scales.

$$\begin{eqnarray} s_{i-1}=s^1_{i-1}\ s^2_{i-1}. \end{eqnarray}$$

Here

$$\begin{eqnarray} s^1_{i-1}=\frac {z_{p^1(i-1)} }{ z_{p(p^1(i-1))}} . \end{eqnarray}$$

$p^1(i-1)$ is the index of the day preceding Monday without marketing. $p(p^1(i-1))$ is the index of the day preceding $p^1(i-1)$.

$$\begin{eqnarray} s^2_i=\frac {z_{p^2(i-1)} }{ z_{p(p^2(i-1))}} . \end{eqnarray}$$

$p^2(i-1)$ is the index of the day preceding the day with marketing. $(p(p^2(i-1))$ is the index of the day preceding $p^2(i-1)$.

Expression (16.77) is an example of empirical approach when multiple SE are regarded as different special events. Expression (16.78) needs less data but is based on the assumption that scales are multiplicative.

2 Examples of PSE

Here, the impact of "marketing" starts at the second working hour of the present Monday $i-1$. Regarding PSE as another SE the scale is

$$\begin{eqnarray} s_{i-1}=\frac {z_{i-1} }{ z_{p(i-1)}} . \end{eqnarray}$$

Here $p(i-1)$ is the index of the day preceding the present one.

Another way is to regard the PSE as a reduced SE. Then the scale $s_{i-1}$ is proportional to the impact duration

$$\begin{eqnarray} s^{2}_{i-1}=1/4 \ s_{i-1}. \end{eqnarray}$$

Here

$$\begin{eqnarray} s_{i-1}=\frac {z_{p^3(i-1)} }{ z_{p(p^3(i-1))}}. \end{eqnarray}$$

$p^3(i-1)$ is the index of the preceding marketing Monday. $p(p^3(i-1))$ is the index of the day preceding $p^3(i-1)$. Applying expression (16.82) one requires less data. However, this expression is based on the assumption that scales are proportional to the impact duration.

Expression (16.82) involves expert opinion indirectly, by assuming proportional scales. Expression (16.81) is an example of empirical approach when a PSE is regarded as another SE.

3 Method of Least Squares

Denote by $\psi_{ij}(s,d,\tau)$ the call rate prediction by an expert model with fixed parameters $s,d,\tau$. One predicts call rates of the part $j$ of the next period $i$ using data up to the period $i-1$. The parameters of the expert model include scales $s$ and delay and duration times $d,\tau$. All special events, including DSE and PSE are regarded. We minimize the sum

$$\begin{eqnarray} \min_{s,d,\tau} \sum_{i=T_0}^T \sum_{j=1}^J (z_{ij}- \psi_{ij}(s,d,\tau))^2. \end{eqnarray}$$

Here $T$ is the end and $T_0$ is the beginning of the "learning" set. One minimizes (16.64) by corresponding global and local optimization methods.

4 Remarks

As usual, expert models predict the weekly and the seasonal graphs. The weekly graphs are given in a form of average calls for each of seven days. In the seasonal graphs, weeks are represented by weekly averages. Often "observed" values of the call rates $z_i$ cannot be defined directly. Then their estimates are used (see Section 5). Missing or uncertain data is replaced by expert estimates. This is important for the first year of operation.

5 Time Series Model (ARMA)

1 Scalar Prediction

Here we consider how to apply the ARMA model to predict the errors $\nu_i$ of the expert scale model $\psi_i(s,d,\tau)$, where $\nu_i=\psi_i(s,d,\tau) -z_i$ [Mockus et al., 1997]. The formal description of the ARMA model is in Chapter 15. The software is on the web-site (see Section 4). If predicting $\nu_i$ we do not know some values of $\nu_s,\ s=i-1,i-2,...$ then we replace the missing data by the expected values of the unknown $\nu_s$ defined recurrently (see Section 10).

ARMA models predict only the difference between the expert models (EM) and the data. Therefore, in the first step, EM models are adapted to the data by defining the right scales. Only then, the parameters $a,b$ of ARMA models are optimized. This means that one minimizes a squared difference between the adapted EM model and the observed data. Predicting the data, the results of both EM and AH models are summed up.

2 Vector Prediction

The vector prediction of EM errors $\nu_i$ is performed using multi-step ARMA predictions described in Section 10 "Confidence cones" could be useful for preliminary evaluation of prediction errors. They estimate probable deviations from expected values with some "confidence level," say ninety percent or ninety nine percent. These cones are defined empirically as upper and lower limits of the results obtained by repetition of Monte Carlo procedures (see Section 10).

6 Application Examples

In the call center example, both the Expert and ARMA models were tested separately. Figure 15.18 shows daily call rates of this center. Table 15.1 illustrates some results of the ARMA model. Results of the expert model are not available, for publication.

11 Call Center Scheduling

There are "Of-the-Shelf" tools for the scheduling of single-skill agents. The scheduling of multi-skill agents is theoretically possible using Monte Carlo simulation (see Section 7). However, the simulation time is too large for on-line scheduling. A convenient approximation to multi-skill scheduling is the reduction of multi-skill problems to the single-skill ones. One can do that by representing each multi-skill agent as a "weighted" single-skill one

$$\begin{eqnarray} N_{single}=v*N_{multi}. \end{eqnarray}$$

Then, one estimates the unknown weight $v$ by minimizing the sum of squared deviations. The deviations are differences between results of approximate single-skill model and the multi-skill one

$$\begin{eqnarray} \min_{v} 1/K \sum_{k=1}^K (Q_k (\lambda, N_{multi})-q_k(\lambda, N_{single}))^2. \end{eqnarray}$$

Here $Q_k$ denotes the results of the $k$ iteration of a Monte Carlo simulation of the multi-skill system. $\lambda$ is the call rate. $N_{multi}$ is the number of multi-skill agents. $q_k$ is the result of the $k$ iteration of a Monte Carlo simulation of a single-skill system. $N_{single}$ is the number of single-skill agents, which replace the multi-skill ones.

Another way is to extend the single-skill simulation (see Section 11) to the multi-skill case. Here one needs more computing time to estimate the "optimal" weights $v$ for different call rates. For example, the method $Exkor$ is used to obtain the best values of $v$ (see Section 4).

12 Software Example

1 ARMA Model

There are two immlementations of the ARMA model, both are on the web-sites n the "Software Systems". The Java version is marked as $ARMA-GMJ$, the C++ version is denoted by $ARMA-C$.

2 Scales Model

To start the program [Jankauskas, 2000], one goes to one of web-sites (see Section 4) and selects the problem "Call Center-SCALES-GMJ" in the section "Global Optimization ". Figure 16.3 shows the operating page.

Figure 16.3: Operating page.

Here is an illustrative example of data file.