Rouba Ibrahim
email: rouba.ibrahim@ucl.ac.uk

Han Ye
email: hanye@illinois.edu

Pierre L'ecuyer
email: lecuyer@iro.umontreal.cahaipeng@email.unc.edu

Haipeng Shen

MODELING AND FORECASTING CALL CENTER ARRIVALS: A LITERATURE SURVEY

The effective management of call centers is a challenging task mainly because managers are consistently facing considerable uncertainty. Among important sources of uncertainty are call arrival rates which are typically time-varying, stochastic, dependent across time periods and across call types, and often affected by external events. Accurately modeling and forecasting future call arrival volumes is a complicated issue which is critical for making important operational decisions, such as staffing and scheduling, in the call center. In this paper, we review the existing literature on modeling and forecasting call arrivals. We also discuss the key issues in building good statistical arrival models. Additionally, we evaluate the forecasting accuracy of selected models in an empirical study with real-life call center data. We conclude by summarizing future research directions in this important field.

Introduction

Call center is a large and important service industry with more than 2.7 million agents working in the United States and 2.1 million agents working at Europe, the Middle East, and Africa ( [START_REF] Akşin | The modern call center: A multi-disciplinary perspective on operations management research[END_REF]). Efficiently managing a call center is a challenging task because managers have to make staffing and scheduling decisions to balance between staffing cost and service quality, which always contradict, in the presence of uncertain arrival demands. Most staffing or scheduling plans start with forecasting customer call arrivals, which is highly stochastic.

Accurately forecasting call arrivals is one of the keys to achieve optimal operational efficiency, since under-forecasting leads to under-staffing and then long customer waiting, while on the other hand over-forecasting results in a waste of money on over-staffing.

The process of customer arrivals is nontrivial. This process can be modeled as a Poisson arrival process and has been shown to possess several features ( [START_REF] Akşin | The modern call center: A multi-disciplinary perspective on operations management research[END_REF][START_REF] Cezik | Staffing multiskill call centers via linear programming and simulation[END_REF][START_REF] Gans | Telephone call centers: Tutorial, review, and research prospects[END_REF][START_REF] Garnett | Designing a call center with impatient customers[END_REF][START_REF] Wallace | A staffing algorithm for call centers with skill-based routing[END_REF]). One of the most important features is that the arrival rate is time varying, which adds complexity to the forecasting process. Call arrival rates may exhibit intra-day, weekly, monthly, and yearly seasonalities. Given the time-varying arrival rate, the doubly stochastic arrival process can be modeled as a non-homogeneous overdispersed Poisson process. Call center arrivals also show different types of dependencies including intra-day (within-day), inter-day, and intertype dependence. A reasonable forecasting model needs to appropriately account for some of all types of dependencies that exist in real data.

In the presence of intra-day and inter-day dependence of call arrival rates, standard time series models may be applied to forecast call arrivals, for example Autoregressive Integrated Moving Average Models and Exponential Smoothing ( [START_REF] Hyndman | Forecasting with Exponential Smoothing: the State Space Approach[END_REF]). In addition, some recent papers have proposed Fixed Effect Models and Mixed Effect Models to account for within-day dependence, inter-day dependence, and inter-type dependence of call arrivals. Dimension reduction or Bayesian techniques are also adopted in the existing literature. Detailed review of various existing forecasting approaches is given in Sections 3 and 4.

We then conduct a case study and implement several recently proposed forecasting approaches on a Canadian call center data set, which reveals the practical features of those approaches.

The remainder of the paper is organized as follows. In Section 2, we discuss the key features of call center arrival process. In Section 3, we review some relevant theories that can be applied to forecasting call center arrivals. In Section 4, we review forecasting methods that are proposed in the existing literature. We then conduct a case study to compare several models proposed in the recent literature in Section 5.

Key Properties of Call Center Arrival Processes

A natural model for call arrivals is the Poisson arrival process ( [START_REF] Akşin | The modern call center: A multi-disciplinary perspective on operations management research[END_REF][START_REF] Cezik | Staffing multiskill call centers via linear programming and simulation[END_REF][START_REF] Gans | Telephone call centers: Tutorial, review, and research prospects[END_REF][START_REF] Garnett | Designing a call center with impatient customers[END_REF][START_REF] Wallace | A staffing algorithm for call centers with skill-based routing[END_REF]). This model is theoretically justified by assuming a large population of potential customers where each customer independently makes a call with a very small probability; the total number of calls made is then approximately Poisson. As mentioned in [START_REF] Kim | Are call center and hospital arrivals well modeled by nonhomogeneous poisson processes[END_REF], the so-called Poisson superposition theorem is a supporting limit theorem, e.g., see [START_REF] Barbour | Poisson Approximation[END_REF].

Recent empirical studies have shown multiple important properties of the call arrival process, many of which are not consistent with the Poisson modeling assumption. In this section, we describe those properties in detail; for a more abridged description, see §2 in [START_REF] Ibrahim | On the modeling and forecasting of call center arrivals[END_REF].

Time dependence of call arrival rates. One of the most important properties of call arrival rates is that they vary with time. In particular, call arrival rates typically exhibit intraday (within-day), daily, weekly, monthly, and yearly seasonalities. We illustrate this time-dependence property in Figures 1,2, and 3 (taken from [START_REF] Ibrahim | Forecasting call center arrivals: Fixed-effects, mixedeffects, and bivariate models[END_REF]), which show arrival patterns that are commonly observed in call centers.

In Figure 1, we plot the number of calls per day arriving to the call center of a Canadian company between October 19, 2009 and September 30, 2010. Figure 1 shows that there exist monthly fluctuations in the data. For example, the moving average line in the plot, which is computed for each day as the average of the past 10 days, suggests that there is an increase in call volume during the months of January and February, i.e., days 54 to 93 in the plot.

In Figure 2, we illustrate weekly seasonality by plotting daily arrival counts, of the same call type as in Figure 1, over two consecutive weeks in the call center. The call center is closed on weekends, so we have a total of 10 workdays in the plot. Figure 2 clearly shows that there is a strong weekly seasonality in the data. Such weekly patterns are very commonly observed in practice, e.g., see Figure 1 in [START_REF] Taylor | A comparison of univariate time series methods for forecasting intraday arrivals at a call center[END_REF] and Figure 2 in [START_REF] Taylor | Density forecasting of intraday call center arrivals using models based on exponential smoothing[END_REF].

For a more microscopic view of arrivals, we plot half-hourly average arrival counts, per weekday, in Figure 3. These intraday averages constitute the daily profile of call arrivals. Figure 3 shows that call volumes are higher, on average, on Mondays than on the remaining weekdays. Figure 3 also shows that all weekdays have a similar daily profile: there are two major daily peaks for call arrivals. The first peak occurs in the morning, shortly before 11:00 AM, and the second peak occurs in the early afternoon, around 1:30 PM. (There is also a third "peak", smaller in magnitude, which occurs shortly before 4:00 PM on Mondays, Tuesdays, and Wednesdays.) Such intraday arrival patterns are also characteristic of call center arrivals; e.g., see [START_REF] Aldor-Noiman | Workload forecasting for a call center: Methodology and a case study[END_REF][START_REF] Avramidis | Modeling daily arrivals to a telephone call center[END_REF][START_REF] Channouf | The application of forecasting techniques to modeling emergency medical system calls in Calgary, Alberta[END_REF][START_REF] Gans | Telephone call centers: Tutorial, review, and research prospects[END_REF][START_REF] Tanir | Call center simulation in Bell Canada[END_REF]. Given that arrival rates are time-varying, which is not accounted for in a Poisson arrival process, a natural extension is to consider a nonhomogeneous Poisson process with a deterministic and time-varying arrival-rate function. For simplicity, it is commonly assumed that call arrival rates are constant in consecutive 15 or 30 minute intervals during a given day; e.g., see [START_REF] Brown | Statistical analysis of a telephone call center: A queueing-science perspective[END_REF][START_REF] Green | Coping with time-varying demand when setting staffing requirements for service system[END_REF][START_REF] Liao | Staffing a call center with uncertain non-stationary arrival rate and flexibility[END_REF].

Nevertheless, it is important to perform statistical tests to confirm that it is appropriate to model call center data as a nonhomogeneous Poisson process. [START_REF] Brown | Statistical analysis of a telephone call center: A queueing-science perspective[END_REF] proposed a specific test procedure based on Kolmogorov-Smirnov test and did not reject the null hypothesis that arrivals of calls are from a nonhomogeneous Poisson process with piecewise constant rates. [START_REF] Kim | Choosing arrival process models for service systems: Tests of a nonhomogeneous poisson process[END_REF] examined several alternative test procedures which have greater power compared to the one suggested in [START_REF] Brown | Statistical analysis of a telephone call center: A queueing-science perspective[END_REF]. [START_REF] Kim | Are call center and hospital arrivals well modeled by nonhomogeneous poisson processes[END_REF] applied Kolmogorov-Smirnov tests to banking call center and hospital emergency department arrival data and showed that they are consistent with the nonhomogeneous Poisson property, but only if certain common features of data have been accounted for including: data rounding, interval partition and overdispersion caused by combining data.

Overdispersion of arrival counts. A consequence of the Poisson modeling assumption is that the variance of the arrival count in each time period is equal to its expectation during that period. However, there is empirical evidence which invalidates this assumption. Indeed, it has been observed that the variance of an arrival count per time period is usually much larger than its expected value; see [START_REF] Aldor-Noiman | Workload forecasting for a call center: Methodology and a case study[END_REF][START_REF] Avramidis | Modeling daily arrivals to a telephone call center[END_REF][START_REF] Jongbloed | Managing uncertainty in call centers using Poisson mixtures[END_REF][START_REF] Steckley | Performance measures for service systems with a random arrival rate[END_REF]. One way of dealing with this overdispersion of count data is to assume that the Poisson arrival process is doubly stochastic, i.e., that the arrival rate itself is a stochastic process; e.g., see [START_REF] Aldor-Noiman | Workload forecasting for a call center: Methodology and a case study[END_REF][START_REF] Avramidis | Modeling daily arrivals to a telephone call center[END_REF][START_REF] Ibrahim | Forecasting call center arrivals: Fixed-effects, mixedeffects, and bivariate models[END_REF][START_REF] Jongbloed | Managing uncertainty in call centers using Poisson mixtures[END_REF][START_REF] Shen | Statistical analysis of call-center operational data: Forecasting call arrivals, and analyzing customer patience and agent service[END_REF][START_REF] Soyer | Modeling and analysis of call center arrival data: A Bayesian approach[END_REF][START_REF] Steckley | Service system planning in the presence of a random arrival rate[END_REF][START_REF] Weinberg | Bayesian forecasting of an inhomogeneous Poisson process with applications to call center data[END_REF].

A doubly stochastic Poisson process can be viewed as a two-step randomization: A stochastic process (for the arrival rate) is used to generate another stochastic process (for the call arrival process) by acting as its intensity. We now illustrate why a doubly stochastic

Poisson process is a way to deal with a higher variance in the arrival count data. Denote by X j the number of arrivals in a given period j, and let Λ j denote the cumulative arrival rate (its integral) over period j. Then, assume that conditional on Λ j , X j has a Poisson distribution with mean Λ j . To simplify notation, we assume in this paper that all periods have the same length and also that the time unit is equal to one period. Then, when the arrival rate is constant over each period, this rate is the same as the cumulative rate Λ j , and we denote both by Λ j . By conditioning on Λ j , the variance of X j is given by:

Var[X j ] = E[Var[X j |Λ j ]] + Var[E[X j |Λ j ]] = E[Λ j ] + Var[Λ j ] . (2.1)
With a random arrival rate function, we have that Var[Λ j ] > 0 on the right-hand side of (2.1), which accounts for the additional variance in Var[X j ].

To model the doubly stochastic Poisson process, [START_REF] Maman | Queues with random arrival rates: Inference, modelling and asymptotics (c-staffing)[END_REF] proposed a Poisson mixture model with a parametric form of the random Poisson arrival rate. They then incorporated the Poisson mixture model into the M/M/n + G queue and derived asymptotically optimal staffing levels (c-staffing).

Interday and Intraday dependencies. In real-life call centers, there is typically evidence for dependencies between the arrival counts, or arrival rates, in different time periods within a single day, or across several days; e.g., see [START_REF] Aldor-Noiman | Workload forecasting for a call center: Methodology and a case study[END_REF][START_REF] Avramidis | Modeling daily arrivals to a telephone call center[END_REF][START_REF] Channouf | A normal copula model for the arrival process in a call center[END_REF][START_REF] Shen | Interday forecasting and intraday updating of call center arrivals[END_REF][START_REF] Tanir | Call center simulation in Bell Canada[END_REF][START_REF] Whitt | Improving service by informing customers about anticipated delays[END_REF]. Those interday (day-to-day) and intraday dependencies typically remain strong even after correcting for detectable seasonalities. Indeed, it is important to do such a correction to avoid erroneously overestimating dependencies in the data.

In Tables 1 and2 we illustrate interday and intraday correlations in the same call center as in Figures 123. Tables 1 and2 remaining weekdays; (iv) correlations are strong and positive between successive half-hourly periods within a day; and (v) intraday correlations are slightly smaller with longer lags.

There are different measures that could be used to capture interday and intraday dependencies in call arrival data. The most commonly used measure is Pearson's correlation coefficient which captures linear dependence in the data; e.g., see [START_REF] Aldor-Noiman | Workload forecasting for a call center: Methodology and a case study[END_REF][START_REF] Avramidis | Modeling daily arrivals to a telephone call center[END_REF][START_REF] Ibrahim | Forecasting call center arrivals: Fixed-effects, mixedeffects, and bivariate models[END_REF][START_REF] Shen | Interday forecasting and intraday updating of call center arrivals[END_REF]. However, since dependencies may not be linear, it is also useful to consider alternative measures such as rank correlation coefficients; see [START_REF] Channouf | A normal copula model for the arrival process in a call center[END_REF] and references therein. For example, Spearman's rank correlation coefficient measures how well the relationship between two variables can be described using a monotonic, but not necessarily linear, function.

Mixed-effects models ( [START_REF] Aldor-Noiman | Workload forecasting for a call center: Methodology and a case study[END_REF][START_REF] Ibrahim | Forecasting call center arrivals: Fixed-effects, mixedeffects, and bivariate models[END_REF]) and, more generally, copulas ( [START_REF] Channouf | A normal copula model for the arrival process in a call center[END_REF][START_REF] Jaoua | Call type dependence in multiskill call centers[END_REF]) are ideally suited to easily capture interday and intraday dependencies in call center arrival data. Models that fail to account for positive interday and intraday dependence in call arrivals may give an overoptimistic view of call center performance measures, and the resulting errors can be very significant; see [START_REF] Avramidis | Modeling daily arrivals to a telephone call center[END_REF][START_REF] Avramidis | Modeling and simulation of call centers[END_REF][START_REF] Steckley | Performance measures for service systems with a random arrival rate[END_REF][START_REF] Steckley | Forecast errors in service systems[END_REF].

X X X X X X X X X X X X Type B Type A (
Intertype dependencies. In multi-skill call centers, there may be positive dependencies between the arrival counts, or arrival rates, corresponding to different call types. For one example, this could occur in multilingual call centers where the same service request is handled in two or several languages. For another example, this may be due to promotions or advertisements which affect several services offered by the same call center. Neglecting dependencies between different call types may lead to overloads, particularly when the same agent handles several correlated call types.

In Table 3 (taken from [START_REF] Ibrahim | Forecasting call center arrivals: Fixed-effects, mixedeffects, and bivariate models[END_REF]), we present estimates of correlations between half-hourly arrival counts for two different call types, Type A and Type B. In Table 3, we focus on the same consecutive half-hour periods as in Table 2. Table 3 illustrates that intertype correlations can be strong and positive. Here, call arrivals to the Type A queue originate in the province of Ontario, and are mainly handled in English, whereas arrivals to the Type B queue originate in the province of Quebec, and are mainly handled in French. Otherwise, arrivals to both queues have similar service requests. Thus, it is reasonable that there exist correlations between their respective arrival processes. There has been some recent effort to model intertype dependencies in the data; see [START_REF] Ibrahim | Forecasting call center arrivals: Fixed-effects, mixedeffects, and bivariate models[END_REF] and [START_REF] Jaoua | Call type dependence in multiskill call centers[END_REF].

Using auxiliary information. Auxiliary information is often available in call centers to improve point or distributional forecasts considerably. For example, when a company sends notification letters to customers, or makes advertisements, this may trigger a large volume of calls; see [START_REF] Landon | Modeling latent sources in call center arrival data[END_REF]. Also, large sporting events or festivals can bring a significant increase of calls to emergency systems; see [START_REF] Channouf | The application of forecasting techniques to modeling emergency medical system calls in Calgary, Alberta[END_REF].

The past service level in the call center may also be a valuable source of information for predicting future arrivals. For example, long previous delays may lead to a high call abandonment rate, which in turn may lead to more redials in the future. Moreover, when the quality of service is poor, callers may not have their problems resolved during the first call that they make, and they may need to reconnect later. Ignoring such redials and reconnects may lead to considerably underestimating call arrival counts; see [START_REF] Ding | On the estimation of the true demand in call centers with redials and reconnects[END_REF].

Finally, in certain types of call centers, for example where people may call to report power outages or those designated to emergency services, bursts of high arrival rates over short periods of time do occur. In this context, an important accident may trigger several dozen different calls within a few minutes, all related to the same event, resulting in a much larger than expected number of calls during that time frame; e.g., see [START_REF] Kim | Forecasting hourly peak call volume for a rural electric cooperative call center[END_REF] for the modeling of peak periods in a rural electric cooperative call center.

In recent years, there have been a few studies on forecasting call arrivals. In §3, we review some relevant theoretical background, and in §4 we review the relevant literature.

Theoretical Background

In this section, we review some relevant theory. Let X t denote observations taken at equallyspaced intervals. Usually, X t is the number of call arrivals per time period, such as a half-hour interval or a day. Alternatively, X t may also be the rate of call arrivals per time period.

Autoregressive Integrated Moving Average (ARIMA) Models

ARMA Models. Since an ARIMA model is a generalization of an autoregressive moving average (ARMA) model, we begin by describing the latter.

We say that the process {X t , t ≥ 0} follows an ARMA model of orders p and q, denoted by ARMA(p, q), if it can be written in the following form:

Φ p (B)(X t -µ) = Θ q (B) t , (3.1)
where: p, q are non-negative integers,

B is the backshift operator defined by BX t ≡ X t-1 , Φ p (B) = 1 -φ 1 B -... -φ p B p , Θ q (B) = 1 -θ 1 B -... -θ q B q , t
are independently and identically distributed (i.i.d.) as N (0, σ 2 ).

We call φ i , 1 ≤ i ≤ p, the autoregressive parameters and θ j , 1 ≤ j ≤ q, the moving-average parameters. We call µ the location parameter. We assume that all roots of the polynomials Φ p (•) and Θ q (•) lie outside the unit circle. This will guarantee that the process X t is both stationary and invertible; see [START_REF] Box | Time Series Analysis: Forecasting and Control[END_REF] for additional details. In what follows, we consider special cases of (3.1) by assigning specific values to the parameters p and q.

MA Models. Letting p = 0 and q = 1 in (3.1) yields:

X t -µ = (1 -θ 1 B) t = t -θ 1 t-1 . (3.2)
This is called a moving-average model of order 1, and is denoted by MA(1). Intuitively, t and t-1 in (3.2) can be interpreted as shocks that disturb X t and move it away from the level µ. The model in (3.2) means that t and t-1 affect X t in the proportions 1 and -θ 1 , respectively. However, for a fixed s, each s will have no effect on the process beyond time s + 1. That is, the effect of each shock in the model does not persist with time.

In order to obtain MA models of higher orders, simply let q > 1 and p = 0 in (3.1). For example, the expression for an MA(q) model is given by:

X t = t -θ 1 t-1 -θ 2 t-2 -... -θ q t-q . (3.3)
AR Models. Letting p = 1 and q = 0 in (3.1), along with |φ 1 | < 1, yields:

(1 -φ 1 B)(X t -µ) = X t -µ -φ 1 X t-1 = t . (3.4)
This is called an autoregressive model of order 1, and is denoted by AR(1). Expanding (3.4) by exploiting the relation between X t and X t-1 yields the following expression:

X t -µ = t + φ t-1 + φ 2 t-2 + ... , (3.5)
which is a special case of an MA model where the effects of all previous shocks, s for s ≤ t, persist at time t; however, their respective influences decay exponentially since |φ 1 | < 1.

In order to obtain AR models of higher orders, simply let p > 1 and q = 0 in (3.1). For example, the expression for an AR(p) model is given by:

X t = µ + θ 1 X t-1 + θ 2 X t-2 + ... + θ p X t-p + t . (3.6)
Differencing and ARIMA Models. The models in (3.1), (3.2), and (3.4) are all stationary models. This means that, for every k, the distribution of the vector (X t , X t+1 , ..., X t+k ) is independent of t and depends only on k. For example, the assumption that |φ| < 1 in (3.4) ensures for the stationarity of the process. For example, to see why it is a necessary condition, consider the same AR(1) model but let φ 1 = 1. Then, based on (3.5), we obtain that:

X t -µ = X t-1 + t = t + t-1 + t-2 + .. , (3.7)
where the effect of each s , s ≤ t, is permanent. It is not difficult to see that (3.7) implies that the distribution of X t changes with t (e.g., the variance of X t increases with t), which means that the process is nonstationary.

Although the process defined in (3.7) is nonstationary, the process Y t defined as

Y t ≡ X t -X t-1 = t ,
clearly is. This subtraction operation is called differencing and, more generally, it is common to difference a nonstationary time series to transform it into a stationary one.

We are now ready to give the definition of an ARIMA(p, d, q) model, which is given by:

Φ p (B)(1 -B) d (X t -µ) = Θ q (B) t , (3.8)
using the same notation as in (3.1) and letting d be the differencing degree parameter; for example, with d = 2, we first difference the original series X t and then difference the resulting differenced series Y t . In general, d is chosen so that the differenced time series is stationary.

Time series that can be made stationary by differencing are called integrated processes.

It is also possible to accommodate for multiple seasonalities (e.g., daily, weekly, monthly, etc.) in an ARIMA model. The resulting model is a seasonal ARIMA model. For example, the general form of an ARIMA(p, q, d)×(p 1 , q 1 , d 1 ) model with one seasonality period, s, is:

Φ p (B)Γ s (B s )(1 -B) d (1 -B s ) d 1 (X t -µ) = Θ q (B)Ω s (B s ) t , (3.9)
where Γ s (•) is a polynomial of some order p 1 , Ω s (•) is a polynomial of some order q 1 , and d 1 is a differencing parameter.

Exponential Smoothing

Another important special case of ARIMA modelling emerges from letting p = 0 and d = q = 1 in (3.8). In this case, we obtain the following series:

Y t ≡ X t -X t-1 = t + θY t-1 + θ 2 Y t-2 + θ 3 Y t-3 + ... ;
in other words, Y t is equal to an exponentially weighted sum of its past observations. This corresponds to the well-known and widely used practice of exponential smoothing.

Holt-Winters Smoothing. The Holt-Winters method is an extension of exponential smoothing which accommodates both a trend and a seasonal pattern; see [START_REF] Winters | Forecasting sales by exponentially weighted moving averages[END_REF]. The Holt-Winters method has two versions, additive and multiplicative, the use of which depends on the characteristics of the particular time series at hand. To illustrate, here are the smoothing equations for the additive Holt-Winters method:

M t = α 0 (X t -S t-s ) + (1 -α 0 )(M t-1 + B t-1 ) , B t = α 1 (M t -M t-1 ) + (1 -α 1 )B t-1 , S t = α 2 (X t -M t ) + (1 -α 2 )S t-s , (3.10)
where B t is the slope component, M t is the level component, S t is the seasonal component, and s is the period of seasonality. The constants α 0 , α 1 , and α 2 are smoothing parameters, whose values are between 0 and 1. In [START_REF] Taylor | A comparison of univariate time series methods for forecasting intraday arrivals at a call center[END_REF], Taylor extended the Holt-Winters method to accommodate multiple seasonalities in the data.

The early work on call forecasting relied mostly on standard techniques such as ARIMA modeling and exponential smoothing; we review these papers in §4.1. More recently, such traditional forecasting methods are often used as benchmarks against more advanced mod-elling approaches. For a newly proposed call arrival model to be worth serious consideration, it should at least outperform those standard forecasting techniques.

Fixed-Effects and Mixed-Effects Models

Gaussian linear fixed-effects and mixed-effects models are useful models which usually build on ARMA models. Indeed, the residuals in fixed-effects and mixed-effects models are often modeled using some ARMA(p, q) process with appropriately chosen p and q. Fixed-effects and mixed-effects models are used in several recent papers on modeling call arrivals; see §4.2.

In general, a linearly additive fixed-effects (FE) model for X t is given by:

X t = k u=1 α u Z u + t , (3.11)
where t are i.i.d. N (0, σ 2 ) variables and Z u are explanatory variables assumed to be known with certainty. The parameters α u , 1 ≤ u ≤ k, are to be estimated from data.

For a linearly multiplicative FE model, simply replace the sums in (3.11) by products.

Multiplicative models are appropriate when the variance of the variable X t increases with its mean. In order to model a series using a multiplicative model, it is common to proceed indirectly by modeling its logarithm using an additive model.

A linearly additive mixed-effects (ME) model is an FE model which includes additional random effects: These are normally distributed random variables which quantify random deviations. The general expression for such an ME model is given by:

X t = k u=1 β u Z u + l v=1 β v γ v + t , (3.12)
where γ v are normally distributed variables which are assumed to be independent from , and the parameters β u , 1 ≤ u ≤ l, are to be estimated from data. For example, γ v may represent the daily deviation for X t , whereas t may represent the normal intradaily noise associated with X t . Then, imposing specific covariance structures on γ and would allow modeling interday and intraday dependencies in the data, respectively. For additional background on linear mixed models, see [START_REF] Muller | Linear Model Theory: Univariate, Multivariate, and Mixed Models[END_REF]. In §4.2, we illustrate how linear fixed-effects and mixed-effects models are used in the context of modelling call arrivals.

Call Forecasting Approaches

Ideally, we want arrival models that seek to reconcile several objectives. For an arrival model to be realistic, it needs to reproduce the properties that we described in section 2.

Simultaneously, for an arrival model to be practically useful, it needs to be computationally tractable. That is, it needs to rely on a relatively small number of parameters so as to avoid overfitting. Moreover, these parameters need to be easy to estimate from historical data.

Finally, parameter estimates should not be hard to update (e.g., via Bayesian methods)

based on newly available information, e.g., throughout the course of a day. These updated estimates would then be used to update operational decisions in the call center.

In this section, we review alternative models proposed in the literature which aim to reconcile those objectives. We first review early papers which rely mostly on standard forecasting methods ( §4.1). Then, we focus on more recent models for arrivals over several days or months ( §4.2). Finally, we move to models for arrivals over a single day ( §4.3).

Standard Forecasting Techniques

The early work on forecasting call arrivals usually focused on modeling daily or even monthly total call volumes. Part of the reason for this is due to the lack of relevant data. In addition, only point forecasts of future arrival rates or counts were produced.

One of the earliest papers on forecasting call arrivals is [49], where the authors modeled monthly call arrivals for two different call types. Interestingly, they noted that there may be an interdependence between the arrival streams of these two call types, but they did not explore this issue further. They used seasonal ARIMA models to forecast future call volumes, and relied solely on the past history of call arrivals in their models.

In [START_REF] Mabert | Short interval forecasting of emergency phone call (911) work loads[END_REF], Mabert relied on multiplicative and additive regression models, including covariates for special events and different seasonalities, to forecast daily call arrivals to an emergency call center. He also considered model adjustments which exploit previous forecasting errors to yield more accurate forecasts. He found that such models yield the most accurate forecasts, and are superior to ARIMA models.

Other early papers also relied on standard time series models. For example, [START_REF] Andrews | Bean improves call-center forecasting[END_REF] modeled daily call arrivals to the call center of a retailer. The authors considered ARIMA models with transfer functions and incorporated covariates for advertising and special-day effects. They

showed that using such information can dramatically improve the accuracy of their forecasts, and may have a significant impact on the operational decision-making in the call center.

Similarly, [START_REF] Bianchi | Improving forecasting for telemarketing centers by ARIMA modeling with intervention[END_REF] used ARIMA models with intervention analysis to forecast telemarketing call arrivals. In this paper, the authors found that such models are superior to additive and multiplicative Holt-Winters exponentially weighted moving average models.

More recently, [START_REF] Antipov | Forecasting call frequency at a financial services call centre[END_REF] modeled the daily number of applications for loans at a financial services telephone call center. The authors also went beyond standard ARIMA models by including advertising response and special calender effects; they did so by adding exogenous variables in a multiplicative model. In [START_REF] Channouf | The application of forecasting techniques to modeling emergency medical system calls in Calgary, Alberta[END_REF], the authors developed simple additive models for the (small) number of ambulance calls during each hour, in the city of Calgary. Their models capture daily, weekly, and yearly seasonalities, selected second-order interaction effects (e.g., between the time-of-day and day-of-week), special-day effects (such as the Calgary Stampede which leads to increased call volumes), and autocorrelation of the residuals between successive hours. Their best model outperforms a doubly-seasonal ARIMA model for the residuals of a model which captures only special-day effects.

Models Over Several Days

To describe more recent arrival modelling approaches, we need some additional notation.

Let X i,j denote the number of call arrivals during period j, 1 ≤ j ≤ P , of day i,

1 ≤ i ≤ D.
The standard assumption is that call arrivals follow a Poisson process with a (potentially) random arrival rate Λ i,j , which is taken to be constant over each period j. The cumulative arrival rate over period j is also Λ i,j if a unit time period is assumed. Thus, conditional on the event Λ i,j = λ i,j , X i,j is Poisson distributed with rate λ i,j .

Several papers ([3, 11, 23, 51]) exploit the following "root-unroot" variance-stabilizing data transformation:

Y i,j ≡ (X i,j + 1/4) 1/2 . (4.1)
Conditional on the event Λ i,j = λ i,j , and for large values of λ i,j , Y i,j is approximately normally distributed with mean λ i,j and variance 1/4; see [START_REF] Brown | The root-unroot algorithm for density estimation as implemented via wavelet block thresholding[END_REF]. The unconditional distribution, with random Λ i,j , is then a mixture of such normal distributions; therefore, it has a larger variance.

Nevertheless, it can be assumed (as an approximation) that the square-root transformed counts Y i,j are normally distributed, particularly if Var[Λ i,j ] is not too large. The resulting normality is very useful because it allows fitting linear Gaussian fixed-effects and mixedeffects models to the square-root transformed data; for background, see §3.

A better alternative than modelling the arrival counts X i,j would be to model the rates Λ i,j directly. The reason is that it is considerably easier to simulate the system with a distributional forecast for the rates rather than one for the counts. Indeed, to simulate arrivals based on a distributional forecast for counts, one has to generate the number of arrivals in each period, and then generate the arrival times by splitting the counts uniformly and independently over the given time period. (This is consistent with the Poisson assumption.)

In contrast, given a distributional forecast for the rates, one can generate the arrival times directly. Nevertheless, most arrival models in the literature are for the counts X i,j , rather than the rates Λ i,j . The reason being that, in practice, we do not observe the arrival rates themselves but only the counts which give only partial information on the rates. This makes estimating arrival rates a more complicated task.

Multiple papers, such as [START_REF] Ibrahim | Forecasting call center arrivals: Fixed-effects, mixedeffects, and bivariate models[END_REF][START_REF] Shen | Interday forecasting and intraday updating of call center arrivals[END_REF][START_REF] Taylor | A comparison of univariate time series methods for forecasting intraday arrivals at a call center[END_REF][START_REF] Weinberg | Bayesian forecasting of an inhomogeneous Poisson process with applications to call center data[END_REF], consider a linear fixed-effects model as a benchmark for comparison. To illustrate, let d i be the day-of-week of day i, where i = 1, 2, . . . , D.

That is, d i ∈ {1, 2, 3, 4, 5} where d i = 1 denotes a Monday, d i = 2 denotes a Tuesday, . . . , and d i = 5 denotes a Friday. In [START_REF] Ibrahim | Forecasting call center arrivals: Fixed-effects, mixedeffects, and bivariate models[END_REF], the authors considered the following fixed-effects model for the square-root transformed arrival counts:

Y i,j = 5 k=1 α k I k d i + 22 l=1 β l I l j + 5 k=1 22 l=1 θ k,l I k d i I l j + µ i,j , (4.2)
where I k d i and I l j are the indicators for day d i and period j, respectively. That is,

I k d i (I l j ) equals 1 if d i = k (j = l)
, and 0 otherwise. The products I k d i I l j are indicators for the cross terms between the day-of-week and period-of-day effects. The coefficients α k , β l , and θ k,l are real-valued constants that need be estimated from data, and µ i,j are independent and identically distributed (i.i.d.) normal random variables with mean 0. Equation (4.2) simplifies to

y i,j = α d i + β j + θ d i ,j + µ i,j . (4.3)
Fixed-effect models seem hard to beat in terms of accuracy of long-term point forecasting (e.g., 2 weeks or more); see [START_REF] Ibrahim | Forecasting call center arrivals: Fixed-effects, mixedeffects, and bivariate models[END_REF] and [START_REF] Taylor | A comparison of univariate time series methods for forecasting intraday arrivals at a call center[END_REF]. Nevertheless, with short forecasting times, one can exploit interday and intraday dependencies in the data to obtain more accurate forecasts.

As an improvement, and based on real call center data analysis, [START_REF] Aldor-Noiman | Workload forecasting for a call center: Methodology and a case study[END_REF] proposed the following linear mixed-effects model:

Y i,j = α d i + β j + θ d i ,j + γ i + i,j ,
where γ i denotes the daily volume deviation from the fixed weekday effect on day i. Then, γ i is the random effect on day i. Let G denote the D × D covariance matrix for the sequence of random effects. The random effects, γ i , are identically normally distributed with expected value E[γ] = 0 and variance Var[γ] = σ 2 G . The authors assume that these random effects follow an AR(1) process. Considering an AR(1) covariance structure for G is both useful and computationally effective, because it requires the estimation of only two parameters, σ G and ρ G . The residuals i,j are also assumed to have an AR(1) structure within each day. As such, this model captures both interday and intraday dependencies in the data. [START_REF] Brown | Statistical analysis of a telephone call center: A queueing-science perspective[END_REF] proposed an earlier version of this model, also based on call-center data, without intraday correlations and without special-day effects.

In [START_REF] Ibrahim | Forecasting call center arrivals: Fixed-effects, mixedeffects, and bivariate models[END_REF], the authors extended this ME model to two bivariate ME models for the joint distribution of the arrival counts to two separate queues, which exploit correlations between different call types. These models account for the dependence between the two call types by assuming that the vectors of random effects or the vectors of residuals across call types are correlated multinormal. This corresponds to using a normal copula; see [START_REF] Kim | Forecasting hourly peak call volume for a rural electric cooperative call center[END_REF]. The choice of copula can have a significant impact on performance measures in call centers, because of the strong effect of tail dependence on the quality of service [START_REF] Jaoua | Call type dependence in multiskill call centers[END_REF]. A strong upper tail dependence for certain call types, for example, means that very large call volumes tend to arrive together for these call types. When this happens, this produces very large overloads.

To reduce the dimensionality of the vectors (Y i,1 , . . . , Y i,P ), [START_REF] Shen | Analysis of call centre arrival data using singular value decomposition[END_REF] proposed the use of singular-value decomposition to define a small number of vectors whose linear transformations capture most of the information relevant for prediction. Based on this, [START_REF] Shen | Interday forecasting and intraday updating of call center arrivals[END_REF] then developed a dynamic updating method for the distributional forecasts of arrival rates. [START_REF] Shen | Forecasting time series of inhomogeneous Poisson processes with application to call center workforce management[END_REF] proposed a method to forecast the latent rate profiles of a time series of inhomogeneous Poisson processes to enable forecasting future arrival rates based on a series of observed arrival counts. [START_REF] Weinberg | Bayesian forecasting of an inhomogeneous Poisson process with applications to call center data[END_REF] also used Bayesian techniques in their forecasts. They exploited the (normal)

square-root transformed counts to include conjugate multivariate normal priors, with specific covariance structures. They used Gibbs sampling and the Metropolis Hastings algorithm to sample from the forecast distributions, which requires long computational times. Moreover, it is unclear how to incorporate exogenous covariates in such a model. [START_REF] Aktekin | Call center arrival modeling: A Bayesian state space approach[END_REF] recently proposed a model based on a Poisson-Gamma process, where Λ i,j = W i,j λ i,j

for fixed λ i,j 's, and where the multiplicative factors W i,j have a gamma distribution and obey a gamma process. [START_REF] Soyer | Modeling and analysis of call center arrival data: A Bayesian approach[END_REF] analyzed the effect of advertisement campaigns on call arrivals.

Theirs is a Bayesian analysis where they model the Poisson rate function using a mixed model approach. This mixed model is shown to be superior to using a fixed-effects model instead. [START_REF] Weinberg | Bayesian forecasting of an inhomogeneous Poisson process with applications to call center data[END_REF] propose an adaptation of [START_REF] Brown | Statistical analysis of a telephone call center: A queueing-science perspective[END_REF] to enable it to update the forecasts of a day defined from the previous days using newly available observations during this day.

[51] also used Bayesian techniques in their forecasts. They exploited the (normal) squareroot transformed counts to include conjugate multivariate normal priors, with specific covariance structures. They used Gibbs sampling and the Metropolis Hastings algorithm to sample from the forecast distributions, which unfortunately requires long computational times. Moreover, it is unclear how to incorporate exogenous covariates in such a model.

In the empirical analysis of [START_REF] Taylor | A comparison of univariate time series methods for forecasting intraday arrivals at a call center[END_REF], several time-series models are compared including ARMA and Holt-Winters exponential smoothing models with multiple seasonal patterns. The latter method was adapted by [START_REF] Taylor | Short-term electricity demand forecasting using double seasonal exponential smoothing[END_REF] for modeling both the intraday and intraweek cycles in intraday data. In [START_REF] Taylor | Density forecasting of intraday call center arrivals using models based on exponential smoothing[END_REF], Taylor extended his model and considered the density forecasting of call arrival rates. To this aim, he developed a new Holt-Winters Poisson count data model with a gamma distributed stochastic arrival rate. He showed that this new model outperformed the basic Holt-Winters smoothing model. [START_REF] Shen | Exponentially weighted methods for forecasting intraday time series with multiple seasonal cycles: Comments[END_REF] comments about Taylor's work, highlighting the difference between modeling arrivals as a single time series, and as a vector time series where each day is modeled as a component of that vector.

Models Over a Single Day

In this section, we focus on modeling arrivals over a single day. The day is divided into p time periods. We denote by X = (X 1 , . . . , X p ) the vector of arrival counts in those periods.

It is commonly assumed that intraday arrivals follow a Poisson process with a random arrival rate. [START_REF] Whitt | Dynamic staffing in a telephone call center aiming to immediately answer all calls[END_REF] proposed to do that by starting with a deterministic arrival rate function {λ(t), t 0 ≤ t ≤ t e }, where t 0 and t e are the opening and closing times of the call center for the considered day, and to multiply this function by a random variable W with mean E[W ] = 1, called the busyness factor for that day. The (random) arrival rate process for that day is then Λ = {Λ(t) = W λ(t), t 0 ≤ t ≤ t e }.

Under this model, the arrival rates at any two given times are perfectly correlated, and

Corr[Λ j , Λ k ] = 1 for all j, k. We also expect the X j 's to be strongly correlated. More specifically, let I j denote the time interval of period j, let λj = I j λ(t)dt, and let X j be the number of arrivals in I j . Using variance and expectation decompositions, one can find that Var[X j ] = λj (1 + λj Var[W ]) and for j = k:

Corr[X j , X k ] = Var[W ] (Var[W ] + 1/ λj )(Var[W ] + 1/ λk ) -1/2 .
This correlation is zero when Var[W ] = 0 (a deterministic rate) and approaches 1 when [26] examined a similar model, but with independent busyness factors, one for each period of the day. Under their model, the Λ j 's are independent, as are the X j 's, which is inconsistent with intraday dependence of call center arrivals. [START_REF] Channouf | Modélisation et optimisation d'un centre d'appels téléphoniques: étude du processus d'arrivée[END_REF] considered a variant of the model where λ(t) is defined by a cubic spline over the day, with a fixed set of knots, and also shows how to estimate model parameters. This can provide a smoother (perhaps more realistic) model of the arrival rate. [START_REF] Channouf | Modélisation et optimisation d'un centre d'appels téléphoniques: étude du processus d'arrivée[END_REF] and [START_REF] Channouf | A normal copula model for the arrival process in a call center[END_REF] proposed models that account for time dependence, overdispersion, and intraday dependencies with much more flexibility to match the correlations between the X j 's, by using a normal copula to specify the dependence structure between these counts. In principle, similar copula models could be developed for the vector of arrival rates, (Λ 1 , . . . , Λ p ), instead of for the vector of counts. [START_REF] Oreshkin | Rate-based daily arrival process models with application to call centers[END_REF] examined the relationship between modeling for the vector of counts and for the vector of rates. In particular, they gave explicit formulas for the relationship between the correlation between rates and that between counts in two given periods, which implied that for a given correlation between rates, the correlation between counts is much smaller in low traffic than in high traffic.

Case Studies

In this section, we present empirical results from a case study using real data collected at a Canadian call center as described in Section 2. We use the data based on two call types and 200 consecutive workdays (excluding weekends). For each call type, we implement the following four methods to forecast the arrival counts based on 6-week historical data. The four methods used in our case study are discussed in [START_REF] Aldor-Noiman | Workload forecasting for a call center: Methodology and a case study[END_REF], [START_REF] Ibrahim | Forecasting call center arrivals: Fixed-effects, mixedeffects, and bivariate models[END_REF] and [START_REF] Gans | Parametric forecasting and stochastic programming models for call-center workforce scheduling[END_REF].

• MU: the multiplicative univariate forecasting model in [START_REF] Gans | Parametric forecasting and stochastic programming models for call-center workforce scheduling[END_REF] • ME: univariate mixed effect model in [START_REF] Aldor-Noiman | Workload forecasting for a call center: Methodology and a case study[END_REF] • BME1: bivariate mixed effect model in [START_REF] Ibrahim | Forecasting call center arrivals: Fixed-effects, mixedeffects, and bivariate models[END_REF] • BME2: bivariate mixed effect model in [START_REF] Ibrahim | Forecasting call center arrivals: Fixed-effects, mixedeffects, and bivariate models[END_REF].

Each method is applied to the data for an out-of-sample rolling forecasting experiment.

To compare different models, we use the Root Mean Squared Error (RMSE) to assess the point forecast accuracy as defined below:

RM SE = 1 K i,j (X i,j -Xi,j ) 2 ,
where Xi,j is the predicted value of X i,j by the model, and K is the total number of predictions. We also report the coverage probability for the 95% prediction interval to evaluate the forecasting distribution, which is defined as:

Cover = 1 K i,j
I(X i,j ∈ ( Li,j , Ûi,j )),

where ( Li,j , Ûi,j ) is the 95% prediction interval for X i,j given by the model.

Tables 4 and5

Conclusions

Forecasting call center arrivals plays a crucial role in call center management such as determining staffing level, scheduling plan and routing policy. Call center arrival process is complex and has to be modeled appropriately to achieve better forecasting accuracy, and as a result, more efficient operational decisions. In this survey paper, we reviewed the existing literature on modeling and forecasting call center arrival process. We also conducted a case study to evaluate several recently proposed forecasting methods with real-life call center data.

An interesting future research direction is to extend the existing forecasting models or develop new models to forecast more than two call types simultaneously. As some stochastic optimization models for staffing and scheduling rely on the joint forecasting distribution of multiple types of arrivals, such multi-type forecasting models with full distributional forecasts have the potential to better meet the quality of service level and improve operational efficiency.

Another research question worth pursuing is to examine the operational impact of improved forecasts, as most existing literature about call center forecasting evaluate forecasting approaches based on only traditional statistical measures such as RMSE and coverage probability without looking at how those improved forecasting models affect call center operations.

By looking at the operational effect of forecasting models, managers can obtain more insight regarding forecasting model selection and system performance evaluation. [START_REF] Gans | Parametric forecasting and stochastic programming models for call-center workforce scheduling[END_REF] has tried to tackle this problem for one call type. More research is needed in this direction.

Figure 1 :Figure 2 :

 12 Figure 1: Daily call arrival counts over successive months in a Canadian call center.

Figure 3 :

 3 Figure 3: Intraday profile of call arrival counts per weekday in a Canadian call center.

  Var[W ] → ∞.[START_REF] Avramidis | Modeling daily arrivals to a telephone call center[END_REF] studied this model in the special situation where W has a gamma distribution with E[W ] = 1 and Var[W ] = 1/γ. Then, each Λ j has a gamma distribution, X has a negative multinomial distribution, the parameters of this distribution are easy to estimate, and the variance of the arrival counts can be made arbitrarily large by decreasing γ toward zero. The model's flexibility is rather limited, because given the λj 's, Var[X j ] and Corr[X j , X k ] for j = k are all determined by a single parameter value, namely Var[W ]. In an attempt to increase the flexibility of the covariance matrix Cov[X], and in particular to enable a reduction of the correlations,[START_REF] Avramidis | Modeling daily arrivals to a telephone call center[END_REF] introduced two different models for X, based on the multivariate Dirichlet distribution.

Table 2 :

 2

	illustrate several properties which are very commonly observed
	in practice: (i) correlations are strong and positive between successive weekdays; (ii) interday
	correlations are slightly smaller with longer lags; (iii) Mondays are less correlated with the

Correlations between arrivals in consecutive half-hour periods on Wednesday morning in a Canadian call center.

Table 3 :

 3 Correlations between Type A and Type B arrivals in consecutive half-hour periods on Wednesday in a Canadian call center.

	(10, 10:30)	0.75	0.72	0.67	0.60	0.59
	(10:30, 11)	0.76	0.73	0.72	0.64	0.62
	(11, 11:30)	0.66	0.65	0.67	0.67	0.63
	(11:30, 12)	0.60	0.56	0.63	0.63	0.63
	(12, 12:30)	0.58	0.54	0.58	0.65	0.62

10, 10:30) (10:30, 11) (11, 11:30) (11:30, 12)

(12, 12:30)

Table 4 :

 4 summarize the comparisons among the four methods. For both call types, ME is most accurate in point forecasts in most scenarios. BME1 and BME2 has better coverage probability when the leading period is one day or one week, and MU has better coverage probability when the leading period is 2 weeks. Forecasting comparison among five methods for call type A.

		Type A
			MU	ME BME1 BME2
	1 day ahead forecast	RMSE 23.51 21.76	22.59	22.69
		Cover	0.89 0.91	0.93	0.93
	1 week ahead forecast RMSE 30.63 29.59	31.37	31.10
		Cover	0.88 0.88	0.88	0.86
	2 weeks ahead forecast RMSE 37.64 37.51	38.04	37.07
		Cover	0.84 0.82	0.80	0.79

Table 5 :

 5 Forecasting comparison among five methods for call type B.

		Type B
			MU	ME BME1 BME2
	1 day ahead forecast	RMSE 16.80 16.31	16.46	16.49
		Cover	0.93 0.92	0.95	0.95
	1 week ahead forecast RMSE 18.55 17.95	17.99	18.12
		Cover	0.95 0.91	0.94	0.94
	2 weeks ahead forecast RMSE 21.05		20.81	20.75
		Cover	0.93 0.86	0.90	0.90


Quentin Heath

Dale Miller

A framework for proof certificates in finite state exploration

Model checkers use automated state exploration in order to prove various properties such as reachability, non-reachability, and bisimulation over state transition systems. While model checkers have proved valuable for locating errors in computer models and specifications, they can also be used to prove properties that might be consumed by other computational logic systems, such as theorem provers. In such a situation, a prover must be able to trust that the model checker is correct. Instead of attempting to prove the correctness of a model checker, we ask that it outputs its "proof evidence" as a formally defined document-a proof certificate-and that this document is checked by a trusted proof checker. We describe a framework for defining and checking proof certificates for a range of model checking problems. The core of this framework is a (focused) proof system that is augmented with premises that involve "clerk and expert" predicates. This framework is designed so that soundness can be guaranteed independently of any concerns for the correctness of the clerk and expert specifications. To illustrate the flexibility of this framework, we define and formally check proof certificates for reachability and non-reachability in graphs, as well as bisimulation and nonbisimulation for labeled transition systems. Finally, we describe briefly a reference checker that we have implemented for this framework.

Introduction

Model checkers are one way in which logic is implemented. While one of the strengths of model checkers is to aid in the discovery of counterexamples and errors in specifications [START_REF] Clarke | Model Checking[END_REF], they can also be used to prove theorems. Furthermore, such theorems might be of interest to other computational logic systems such as more general theorem provers. One then encounters the problem of whether or not such a theorem prover is willing to trust that model checker or at least a particular theorem it proves. Formally verifying a model checker might be both extremely hard to do and undesirable especially if that checker is being revised and improved. A more plausible option might be to have a model checker output its "proof evidence" as a document (a certificate). If that proof certificate can be formally checked by a trusted checker, one might then be willing to use the theorem in a theorem prover.

Of course, model checkers are asked to solve many kinds of problems so their proof evidence might take many different forms, ranging from decision procedures to paths in graphs, bisimulations, traces, and winning strategies. If we need to have trusted checkers for all these different kinds of proof evidence, then maybe we have not really improved the situation of trust. Here, we contribute to the foundational proof certificate (FPC) effort [START_REF] Miller | A proposal for broad spectrum proof certificates[END_REF] by providing a framework for defining the semantics of a range of proof evidence that naturally arises in model checking. Such a formal semantic model for proof evidence allows anyone to build a proof checker of any formally defined evidence. Furthermore, it is possible to have an implementation of the entire framework of FPC so that this one system could check a wide range of proof evidence.

While this paper has a number of parallels with FPCs for first-order logic in [START_REF] Chihani | Foundational proof certificates in first-order logic[END_REF], that work was limited to first-order logic without fixed points and, as a result, that work was not directly applicable to topics of model checking and inductive and co-inductive theorem proving.

Proof theory for fixed points and certificates

Having proof certificates that are foundational here means that we need to find proof theoretic descriptions of model checking. We shall now describe a few recent developments in proof theory that we bring together in this paper. Of course, the topic of model checking is mature and varied. In order to lay down a convincing and direct proof theory for model checking, we eschew many of its more advanced topics-e.g., predicate abstraction and partial order reduction-for later consideration.

Fixed points as defined predicates

One of the earliest applications of sequent calculus to computational logic was to provide an execution model for logic programming [START_REF] Miller | Uniform Proofs as a Foundation for Logic Programming[END_REF]. That analysis, however, supported only the "open world assumption" of logic programming: negation-as-finite-failure was not touched by that work. Schroeder-Heister [START_REF] Schroeder-Heister | Rules of Definitional Reflection[END_REF] and Girard [START_REF] Girard | A Fixpoint Theorem in Linear Logic[END_REF] showed how sequent calculus could be extended with inference rules for fixed points (or defined predicates), thereby embracing important aspects of the closed world assumption and negationas-finite-failure. The key additions to sequent calculus were rules for unfolding fixed point expressions as well as dealing with equality over the Herbrand universe. A series of papers [START_REF] Gacek | Nominal abstraction[END_REF][START_REF] Mcdowell | Cut-elimination for a logic with definitions and induction[END_REF][START_REF] Miller | A proof theory for generic judgments[END_REF] added induction and co-induction to the sequent calculi for intuitionistic and classical logics. Those papers have been used to design the Bedwyr model checker [START_REF] Baelde | The Bedwyr system for model checking over syntactic expressions[END_REF][START_REF] Tiu | Mixing Finite Success and Finite Failure in an Automated Prover[END_REF] and the Abella interactive theorem prover [START_REF] Baelde | Abella: A System for Reasoning about Relational Specifications[END_REF].

Fixed point expressions will be written as µB t or νB t, where B is an expression representing a higher-order abstraction, and t is a list of terms. The unfolding of the fixed point expression µB t is written as B(µB) t. It is important to understand that we shall treat both µ (least fixed point operator) and ν (greatest fixed point operator) as logical connectives since they will have introduction rules: they are also de Morgan duals of each other.

Fixed points in linear logic

Surprisingly, it is linear logic and not intuitionistic or classical logics per se that is most relevant to our exposition on model checking. The logic MALL (multiplicative additive linear logic) is an elegant, small logic that is, in and of itself, not appropriate for formalizing mathematics and computer science since it is not capable of modeling unbounded behaviors (for example, it is decidable). While Girard extended MALL with the "exponentials" (! and ?) [START_REF] Girard | Linear Logic[END_REF], Baelde [START_REF] Baelde | Least and greatest fixed points in linear logic[END_REF] extended it by adding the least (µ) and greatest (ν) fixed points operators as logical connectives. The resulting logic, called µMALL, forms the proof theoretic foundation of this paper.

To make the use of linear logic easier to swallow for those more familiar with model checking, we adopt the following shallow changes to its presentation. First, we use a two-sided sequent calculus instead of the one-sided calculus used for µMALL. While this change will double the size of our proof system, it will make inference rules look more familiar. Second, we replace the linear logic connectives with familiar connectives (although with annotations). In particular, we replace ⊗, &, ⊕ and their units 1, and 0 with ∧ + , ∧ -, ∨, t + , t -and f + , respectively. (Truth functionally, the two versions of these operators are equivalent: their differences only influence the structure of focused proofs.) We also replace the negatively biased false ⊥ with f -, and instead of the multiplicative disjunction A℘B, we use the implication

A ⊥ ⊃ B: the de Morgan dual of A ⊃ B is A ∧ + B ⊥ . Negation is written as • ⊃ f -.
In addition, we consider µ as positive and ν as negative; this arbitrary choice has been shown to give a convenient natural interpretation to the structure of focused proofs [START_REF] Baelde | Least and greatest fixed points in linear logic[END_REF]. We therefore have the negative connectives f -, ⊃, t -, ∧ -, ∀, = and ν, and the positive connectives t + , ∧ + , f + , ∨, ∃, = and µ.

Focused proof systems

In order to have the kind of control we need to support a definable notion of proof certificate, we make use of a focused proof system. Such sequent calculus proof systems are built from alternating phases which allow us to define flexible proof building protocols that can be used to drive proof search. During the asynchronous phase of proof building, simple (invertible) computations build a proof and during the synchronous phase, information needed for the construction of a proof (such as which branch of a disjunction to prove) must be found.

Focusing requires polarizing all formulas as being either negative or positive. A formula is negative or positive according to its top-level connective, and it is purely negative (resp. purely positive) when its connectives are positive if, and only if, they occur under an odd (resp. even) number of implications. Notice that the de Morgan dual of a positive (resp. purely positive) formula is a negative (resp. purely negative) formula. We call a formula bipolar when it is made of purely negative (resp. positive) subformulas occurring under an even (resp. odd) number of implications in a purely negative context.

Focusing also relies on the sequents having additional storage zones on each side of the turnstile, where formulas can be stored and left untouched by logical inference rules. For instance, the usual onesided focused presentation of µMALL [START_REF] Baelde | Least and greatest fixed points in linear logic[END_REF] has one of these zones, similarly to the focused proof system for linear logic given by Andreoli [START_REF] Andreoli | Logic Programming with Focusing Proofs in Linear Logic[END_REF]. A two-sided subsystem of µMALLF, called µF, makes use of two storage zones, noted N and P, which are lists of, respectively, negative and positive formulas. (Appendix B contains an example of a µF proof.) Between the arrows and the turnstile, are the contexts Γ and ∆: these are lists of formulas in (unfocused) ⇑-sequents, and sets of up to one formula in (focused) ⇓-sequents. The sequents of the µF system are therefore:

N ⇑ Γ ∆ ⇑ P unfocused, similar to the µMALLF sequent N ⊥ , P ⇑ Γ ⊥ , ∆ ⇓ A left-focused, similar to ⇓ A ⊥ A ⇓
right-focused, similar to ⇓ A

Foundational proof certificates

If we think of the implementers of computational logic systems (e.g., model checking systems) as our clients, our job in this project is to formally check our client's proof evidence for formal correctness.

Our approach is to have this evidence translated into a sequent calculus proof. Of course, we would not dream of asking our clients to supply a sequent calculus proof in the first place: such proofs are often huge, too messy, and too esoteric. Instead, we want to take from our clients objects with which they are familiar (e.g., paths, simulations, etc.) and find flexible and high-level means to have our framework extract information from those objects in order to trace out a complete formal sequent calculus proof.

To this intent, Figures 1 and2 present µF a , which is a version of µF augmented with a term Ξ (encoding an actual certificate) as well as with clerk and expert predicates (examples of which we provide soon). This augmentation has two components. First, every sequent (either ⇑ or ⇓) is given an extra argument we write as Ξ. Thus, sequents now display as

Ξ : N ⇑ Γ ∆ ⇑ P, Ξ : ⇓ A , and
Ξ : A ⇓ .
Second, every inference rule is given an additional premise. In all cases, this premise is an atomic formula with either a clerk or expert predicate as its head symbol: if the conclusion of the inference rule is a ⇓-sequent, then the premise atom uses an expert predicate (noted e (. . . ) for the rule ); otherwise, the conclusion is an ⇑-sequent and the atom uses a clerk predicate (noted c (. . . )).

In the case of the clerk rules, the premise atom relates the Ξ value of the concluding sequent with the corresponding value of Ξ for all premises: e.g.,

Ξ 1 : N ⇑ A 1 , Γ ∆ ⇑ Ξ 2 : N ⇑ A 2 , Γ ∆ ⇑ ∨ c (Ξ 0 , Ξ 1 , Ξ 2 ) Ξ 0 : N ⇑ A 1 ∨ A 2 , Γ ∆ ⇑ ∨ L
In this way, the certificate Ξ, intended to aid in the proof of the concluding sequent, can be transformed into two certificates that are used to prove the two premise sequents. We refer to the predicates used in the asynchronous phase as clerks since these predicates do not need, in general, to examine the actual information in the proof certificate (except for the induction and co-induction rules, there is no consumption of information during the asynchronous phase). Instead, the clerks are responsible for keeping track of how a proof is unfolding: for example, Ξ 1 might be a copy of Ξ 0 but with the fact that checking has moved to the left premise instead of the right premise.

Experts are responsible for extracting information from a certificate. For example, µF a contains the inference rule

Ξ 1 : C t ⇓ ∃ e (Ξ 0 , Ξ 1 ,t) Ξ 0 : ∃x.C x ⇓ ∃ R
Notice here that the exists-expert ∃ e (•, •, •) not only computes the continuation certificate Ξ 1 but also a term t to be used to witness the existential variable.

The exact nature of both certificate terms Ξ and of the clerk and expert predicates is not important to guarantee soundness of this system. That is, no matter how certificates, clerks, and experts are specified, if there is a proof in µF a then there is a proof in µF of the same sequent, which can be obtained by deleting from the proof in µF a all references to Ξ, including the additional premises. Notice also that experts are not required to act particularly expertly: it is entirely possible for the ∃ e (Ξ 0 , Ξ 1 ,t) premise to functionally determine one t from Ξ 0 , or to relate all terms t to Ξ. In the latter case, the actual value of t used in a successful µF a proof is determined from other aspects of the proof checking process (typically implemented using unification).

A proof system underlying model checking

FPCs were first proposed in [START_REF] Chihani | Foundational proof certificates in first-order logic[END_REF][START_REF] Miller | A proposal for broad spectrum proof certificates[END_REF] in the context of first-order logic and were used successfully to define and check proof evidence in the form of resolution refutations, Herbrand instances (expansion trees), natural deduction (λ -terms), Frege proofs, etc. We shall now adapt this approach to formally define the semantics of a range of proof evidence that can arise in simple but real model checking problems.

We shall later illustrate just how such a formal semantics can be provided for the following four kinds of proof evidence. These particular examples have been selected for their universality: numerous problems in model checking are related to them.

1. The fact that two nodes are related by the transitive closure of a graph's adjacency relation can be witnessed by an explicit path through the graph.

2. The fact that two nodes are not related by transitivity can be witnessed by pointing out that the reachable set of one does not contain the other.

3. Given an LTS (labeled transition system), the fact that two nodes are similar/bisimilar can be witnessed by a set of pairs called simulation/bisimulation.

4.

If two nodes in an LTS are not bisimilar, then there is a Hennessy-Milner logic (HML) formula that is satisfied by one but not by the other.

ASYNCHRONOUS CONNECTIVE INTRODUCTIONS Ξ 1 θ : N θ ⇑ Γθ ∆θ ⇑ = s c (Ξ 0 , Ξ 1 ) Ξ 0 : N ⇑ s = t, Γ ∆ ⇑ = s L † Ξ 1 θ : N θ ⇑ ⇑ = f c (Ξ 0 , Ξ 1 ) Ξ 0 : N ⇑ s = t ⇑ = f R † Ξ 1 : N ⇑ Γ ∆ ⇑ t + c (Ξ 0 , Ξ 1 ) Ξ 0 : N ⇑ t + , Γ ∆ ⇑ t + L Ξ 1 : N ⇑ ⇑ f - c (Ξ 0 , Ξ 1 ) Ξ 0 : N ⇑ f -⇑ f - R Ξ 1 : N ⇑ A 1 , A 2 , Γ ∆ ⇑ ∧ + c (Ξ 0 , Ξ 1 ) Ξ 0 : N ⇑ A 1 ∧ + A 2 , Γ ∆ ⇑ ∧ + L Ξ 1 : N ⇑ A 1 A 2 ⇑ ⊃ c (Ξ 0 , Ξ 1 ) Ξ 0 : N ⇑ A 1 ⊃ A 2 ⇑ ⊃ R = f c (Ξ 0 ) Ξ 0 : N ⇑ s = t, Γ ∆ ⇑ = f L ‡ = s c (Ξ 0 ) Ξ 0 : N ⇑ s = t ⇑ = s R ‡ f + c (Ξ 0 ) Ξ 0 : N ⇑ f + , Γ ∆ ⇑ f + L t - c (Ξ 0 ) Ξ 0 : N ⇑ t -⇑ t - R Ξ 1 : N ⇑ A 1 , Γ ∆ ⇑ Ξ 2 : N ⇑ A 2 , Γ ∆ ⇑ ∨ c (Ξ 0 , Ξ 1 , Ξ 2 ) Ξ 0 : N ⇑ A 1 ∨ A 2 , Γ ∆ ⇑ ∨ L Ξ 1 : N ⇑ A 1 ⇑ Ξ 2 : N ⇑ A 2 ⇑ ∧ - c (Ξ 0 , Ξ 1 , Ξ 2 ) Ξ 0 : N ⇑ A 1 ∧ -A 2 ⇑ ∧ - R Ξ 1 y : N ⇑ C y, Γ ∆ ⇑ ∃ c (Ξ 0 , Ξ 1 ) Ξ 0 : N ⇑ ∃x.C x, Γ ∆ ⇑ ∃ L Ξ 1 y : N ⇑ C y ⇑ ∀ c (Ξ 0 , Ξ 1 ) Ξ 0 : N ⇑ ∀x.C x ⇑ ∀ R SYNCHRONOUS CONNECTIVE INTRODUCTIONS = f e (Ξ 0 ) Ξ 0 : ⇓ t = t = f L = s e (Ξ 0 ) Ξ 0 : t = t ⇓ = s R f - e (Ξ 0 ) Ξ 0 : ⇓ f - f - L t + e (Ξ 0 ) Ξ 0 : t + ⇓ t + R Ξ 1 : A 1 ⇓ Ξ 2 : ⇓ A 2 ⊃ e (Ξ 0 , Ξ 1 , Ξ 2 ) Ξ 0 : ⇓ A 1 ⊃ A 2 ⊃ L Ξ 1 : A 1 ⇓ Ξ 2 : A 2 ⇓ ∧ + e (Ξ 0 , Ξ 1 , Ξ 2 ) Ξ 0 : A 1 ∧ + A 2 ⇓ ∧ + R Ξ 1 : ⇓ A i ∧ - e (Ξ 0 , Ξ 1 , i) Ξ 0 : ⇓ A 1 ∧ -A 2 ∧ - L Ξ 1 : A i ⇓ ∨ e (Ξ 0 , Ξ 1 , i) Ξ 0 : A 1 ∨ A 2 ⇓ ∨ R Ξ 1 : ⇓ C t ∀ e (Ξ 0 , Ξ 1 ,t) Ξ 0 : ⇓ ∀x.C x ∀ L Ξ 1 : C t ⇓ ∃ e (Ξ 0 , Ξ 1 ,t) Ξ 0 : ∃x.C x ⇓ ∃ R STRUCTURAL RULES Ξ 1 : N ⇑ Γ ∆ ⇑ store L (Ξ 0 , Ξ 1 ) Ξ 0 : ⇑ N, Γ ∆ ⇑ S L Ξ 1 : ⇑ ⇑ P store R (Ξ 0 , Ξ 1 ) Ξ 0 : ⇑ P ⇑ S R Ξ 1 : ⇓ N decide L (Ξ 0 , Ξ 1 ) Ξ 0 : N ⇑ ⇑ D L Ξ 1 : P ⇓ decide R (Ξ 0 , Ξ 1 ) Ξ 0 : ⇑ ⇑ P D R Ξ 1 : ⇑ P ⇑ release L (Ξ 0 , Ξ 1 ) Ξ 0 : ⇓ P R L Ξ 1 : ⇑ N ⇑ release R (Ξ 0 , Ξ 1 ) Ξ 0 : N ⇓ R R Figure 1:
The µF a 0 proof system. (This proof system is best viewed using color). y stands for a fresh eigenvariable, s and t for terms, N for a negative formula, P for a positive formula, and C for the abstraction of a formula over a variable. The † proviso requires that θ is the mgu of s and t, and the ‡ proviso requires that s and t are not unifiable.

FIXED-POINT

RULES Ξ 1 ȳ : ⇑ B S ȳ S ȳ ⇑ Ξ 2 : N ⇑ S t, Γ ∆ ⇑ ind(Ξ 0 , Ξ 1 , Ξ 2 , S) Ξ 0 : N ⇑ µB t, Γ ∆ ⇑ µ Ξ 1 : N ⇑ S t ⇑ Ξ 2 ȳ : ⇑ S ȳ B S ȳ ⇑ co-ind(Ξ 0 , Ξ 1 , Ξ 2 , S) Ξ 0 : N ⇑ νB t ⇑ ν Ξ 1 : N ⇑ B(µB) t, Γ ∆ ⇑ µ-unfold L (Ξ 0 , Ξ 1 ) Ξ 0 : N ⇑ µB t, Γ ∆ ⇑ µ L Ξ 1 : N ⇑ B(νB) t ⇑ ν-unfold R (Ξ 0 , Ξ 1 ) Ξ 0 : N ⇑ νB t ⇑ ν R Ξ 1 : ⇓ B(νB) t ν-unfold L (Ξ 0 , Ξ 1 ) Ξ 0 : ⇓ νB t ν L Ξ 1 : B(µB) t ⇓ µ-unfold R (Ξ 0 , Ξ 1 ) Ξ 0 : µB t ⇓ µ R Figure 2:
The µF a proof system results from adding these rules to µF a 0 . ȳ stands for an list of fresh eigenvariables, t for an list of terms, and B for the abstraction of a formula over a predicate and a variable list.

Core proof system

Figures 1 and 2 contain the rules for the augmented focused proof systems µF a 0 and µF a . One could obtain the non-augmented systems µF 0 and µF by ignoring the certificates (annotated Ξ variables) and the clerk and expert premises; the resulting rules would be no more than (slightly restricted) two-sided versions of the µMALLF rules. The various clerk and expert predicates are named and displayed in their corresponding inference rules. Notice that those inference rules that involve the use of eigenvariables (∃ L , ∀ R , µ and ν) require the associated clerk predicates to return abstractions over certificates: in this way, premise certificates can be applied to the eigenvariables.

A key element of our proof theoretic treatment of model checking via µF a is the fact that focused sequents contain only one formula. This fact entails that µF a can only be complete with respect to µMALLF on a fragment where derivations satisfy this constraint. In particular, the N and P zones must never contain more than one formula, and never both at the same time. This can be ensured at least for the µF a 0 subsystem by the following restriction on formulas. Definition 1 (switchable formula, switchable occurrence). A µF a formula is switchable if

• whenever a subformula C ∧ + D occurs negatively (under an odd number of implications), either C or D is purely positive;

• whenever a subformula C ⊃ D occurs positively (under an even number of implications), either C is purely positive or D is purely negative.

An occurrence of a formula B is switchable if it appears on the right-hand side (resp. left-hand side) and B (resp. B ⊃ f -) is switchable.

Notice that both a purely positive formula and its de Morgan dual are switchable. The follow theorem is proved by a simple induction on the structure of µF a 0 proofs. Theorem 1 (switchability). Let Π be a µF a 0 derivation of either ⇑ A ⇑ or ⇑ A ⇑ , where the occurrence of A is switchable. Every sequent in Π that is the conclusion of a rule that switches phases (either a decide or a release rule) contains exactly one occurrence of a formula and that occurrence is switchable.

Theorem 1 states that an invariant of the µF a 0 proof system (for switchable theorems) is that the number of non-purely asynchronous formulas (i.e. non-purely positive from N and Γ, and non-purely negative from P and ∆) is one or less. Keeping sequents mostly asynchronous allows the asynchronous phase to deal with most of the context: that way the synchronous phase is left with a single, meaningful formula. (The structure of focused proofs based on switchable formulas is similar to the structure of simple games in the game-theoretic analysis of focused proofs in [7, Section 4].) While the restriction to switchable formulas provides a match to the model checking problems we develop here, that restriction is not needed for using clerks and experts (the examples in [START_REF] Chihani | Foundational proof certificates in first-order logic[END_REF] involve non-switchable formulas).

Encoding of recursively defined predicates

In order to exploit the properties of µF a 0 in model checking problems, we need them to extend to µF a by adding fixed-point rules. As those rules make use of the higher-order variables S (an invariant which is either a pre-fixed point or a post-fixed point) and B (the body of a predicate definition), they cannot be used freely without violating Theorem 1. We propose the following constraints on µF a proofs of switchable formulas so as to have exactly one formula per sequent when phases are switched:

• "arithmetic" restriction: S and B are purely positive (resp. negative);

• "model checking" restriction: S is purely negative (resp. positive), and the context does not trigger synchronous rules (N is empty, Γ is purely positive and ∆ is purely negative). The former restriction would allow to extend the scope of the framework by handling simple theorems involving inductive definitions (e.g. about natural numbers), but is not treated here. The latter restriction better suits our needs (since an asynchronous context fits the spirit of model checking) and is respected by all our examples. Example 2. Horn clauses (in the sense of Prolog) can be encoded as purely positive fixed point expressions. For example, here is the Horn clause logic program (using the λ Prolog syntax, the sigma Y\ construction encodes the quantifier ∃Y ) for specifying the graph in Figure 3

. i ((p = u i ) ∧ + (a = v i ) ∧ + (q = w i ))
and the simulation and bisimulation relations can be defined as the following greatest fixed point expressions (note: the second contains both ∧ -and ∧ + ). Both of these formulas are switchable. ν λ Sλ pλ q. ∀a∀p . p a -→ p ⊃ ∃q . q a -→ q ∧ + S p q (sim)

ν λ Bλ pλ q. (∀a∀p . p a -→ p ⊃ ∃q . q a -→ q ∧ + B p q ) ∧ -(∀a∀q . q a -→ q ⊃ ∃p . p a -→ p ∧ + B q p ) (bisim)

Common proof certificates

The presentation of an FPC now involves the following three steps.

1. Describe how unpolarized formulas should be polarized.

2. Describe the structure of certificates Ξ. This can be done, for example, by describing the signature of constructors for certificates.

Define the clerk and expert predicates.

To ease steps 2 and 3, we define the following certificate constructors (shown together with their types), which describe generic focused proof behaviors. The associated clauses can be included into any subsequent clerks and experts definitions.

The stop:cert certificate authorizes no search; it is to be used as a continuation certificate for other certificate constructors.

The sync:cert->cert certificate constructor authorizes µF a to conduct an unbounded synchronous search for a proof before handing the search over to a continuation certificate. It has no clerks and its experts run an exhaustive non-deterministic search for ∨ and ∃. The experts for the right rules are:

= s e (sync(Ξ)). ∨ e (sync(Ξ), sync(Ξ), 1).

∧ + e (sync(Ξ), sync(Ξ), sync(Ξ)). ∨ e (sync(Ξ), sync(Ξ), 2).

∀T. ∃ e (sync(Ξ), sync(Ξ), T ). µ-unfold R (sync(Ξ), sync(Ξ)). release R (sync(Ξ), Ξ).

The async:cert->cert certificate constructor is the dual of sync; it handles an asynchronous phase and has no experts apart from the decide rules. The clerks for the left rules are:

= s c (async(Ξ), async(Ξ)).
∨ c (async(Ξ), async(Ξ), async(Ξ)).

∧ + c (async(Ξ), async(Ξ)). ∃ c (async(Ξ), λ x. async(Ξ)).

µ-unfold L (async(Ξ), async(Ξ)). store L (async(Ξ), async(Ξ)).

decide L (async(Ξ), Ξ).

bipole n :cert is actually short-hand for a chain of n async(sync(•)) before a final stop. It is used for bounded-depth search when simple search strategies would otherwise not terminate. We also write bipole:cert for bipole 1 = async(sync(stop)).

The decproc:cert constructor is short-hand for bipole ∞ , the unbounded version of bipole n . It is a general purpose decision procedure used for automated and unguided proving. Its rules are similar to those from sync and async, and can be obtained via the equivalence decproc = async(sync(decproc)).

The two constructors inv and co-inv:(i->i->bool)->cert->cert each take an explicit predicate S as parameter. It is expected to be proved to be an invariant with the help of bipole.

∀S. ind(inv(S, Ξ), λ x. bipole, Ξ, S)

∀S. co-ind(co-inv(S, Ξ), Ξ, λ x. bipole, S)

We now turn our attention to describing how to formally define the four kinds of proof evidence mentioned earlier in Section 3. Some of the constructors defined above will be used in those definitions.

Examples: certificates for graphs

We use the notations from Example 2 to define • -→ • and path.

Lists as reachability certificates

The natural choice for a certificate of the proof of path(x, y) is an explicit path, i.e. a list of nodes starting right after x and ending right before y. In fact, this list L can be used directly as the proof certificate. Aside from the initial store R , no clerks are invoked in the process of checking this particular FPC. The following clauses defining the experts only use the provided information to instantiate the logical variables of the proof.

∀X∀L. ∃ e (X :: L, L, X).

∀L. ∧ + e (L, sync(stop), L). ∀L. decide R (L, L).

∀X∀L. ∨ e (X :: L, X :: L, 2). ∀L. ∨ e (nil, sync(stop), 1). ∀L. µ-unfold R (L, L).

In this setting, the sync(stop) certificate will terminate quickly since it is only searching through the term that defines • -→ •.

Invariants as non-reachability certificates

The non-reachability problem comes in two forms: if there are no loops in the graph, then a simple check of the set of nodes reachable from the first node provides a simple decision procedure; if there are loops, then induction is needed.

In the first case, the decision procedure can directly be translated as an FPC for proving ¬path(x, y).

Example 5. In Figure 3, (a) is not reachable from (d), as witnessed by async(stop).

On the other hand, if the underlying graph has loops, then the rules of Figure 1 only do not allow proof search to terminate. As the body B of the path expression (i.e., the displayed formula without µ) is purely positive, a bipole can prove that a chosen purely negative predicate S containing no fixed point expressions is an induction invariant (bipole : ⇑ B S x y S x y ⇑ ), which means that we can use the certificate constructor inv(S, •). Then we use another bipole as continuation certificate for this constructor to check that the invariant is adequate for the refutation of path(t, u) (bipole :

⇑ St u • ⇑ ).
Here, the invariant can be chosen so as to represent the fact of not belonging to the set T × {u}, where T is the reachable set of {t}. Example 6. In Figure 3,(d) is not reachable from (b), as witnessed by inv(S, bipole), where the invariant S (built from the set {b, c} × {d}) is Bisimilarity and similarity are important relationships in the domain of process calculus and model checking. To illustrate how these can be captured as FPCs, we first restrict our attention to the existence of a simulation between finite labeled transition systems; bisimilarity is then addressed by expanding on this presentation. We define •

λ xλ y. ((x = b ∧ + y = d) ∨ (x = c ∧ + y = d)) ⊃ f -
• -→ • (for the LTS), sim (for simulated-by) and bisim (for bisimulated-by) as seen in Example 3.

Invariants as simulation certificates

We shall consider two cases: one where the underlying transition system is noetherian and one where it is not. An LTS is said to be noetherian if there is no infinite sequence of transitions p 1

a 1 -→ p 2 a 2 -→ • • •
(in the setting of finite LTSs, this is equivalent to the absence of loops).

In the noetherian case, there is a decision procedure to determine whether or not one process is simulated by another: one simply attempts to incrementally check simulation at every point. This systematic search can be described using the clerks and experts of the decproc certificate, which allows a proof to be built from any number of bipoles (one for each unfolding of the simulation predicate, which formula is itself bipolar).

Example 7. In Figure 4, the process (1) is simulated by the process (6), as witnessed by the certificate decproc.

In the more general (possibly non-noetherian) setting, we need to recall the formal definition of the simulation relation as a set. A binary relation S is a simulation if whenever p, q ∈ S and whenever p a -→ p holds, then there exists a q such that q a -→ q holds and p , q ∈ S. We say that process p is simulated by process q if there is a simulation S such that p, q ∈ S.

Let S be a finite set of pairs and let Ŝ be the purely positive expression λ xλ y. p,q ∈S (x = p∧ + y = q). As the body B of the sim expression is a bipolar formula, a bipole can prove the closure condition for (finite) simulations (bipole : ⇑ Ŝ x y B Ŝ x y ⇑ ), so we can use the certificate constructor co-inv( Ŝ, •). Once again, we use another bipole as continuation certificate to complete the proof that p is simulated by q (bipole : ⇑ • Ŝ p q ⇑ ). Providing an entire invariant as part of a proof certificate or restricting to the case when an invariant is finite certainly limits what kinds of simulation relationships can be proved. In general, invariants will not be finite and, even when they are, they are large. It is for reasons such as this that there has been a great deal of work on bisimulation-up-to [START_REF] Milner | Communication and Concurrency[END_REF][START_REF] Pous | Enhancements of the bisimulation proof method[END_REF]: generally, it is possible to discover and check a closure property of a much smaller relationship and then via various meta-theoretic properties, ensure that such closure properties entail the existence of a proper (bi)simulation.

Assertions as non-simulation certificates

Hennessy and Milner [START_REF] Hennessy | Algebraic Laws for Nondeterminism and Concurrency[END_REF] provided a characterization of bisimulation in terms of an assertion language over modal operators [a] and a . The characterization states that two processes are bisimilar if and only if they satisfy the same assertion formulas. Thus, if p and q are not bisimilar, there is some assertion formula A which is true for p and not for q. Formally, we write p |= A and q |= A.

It is possible to use such assertion formulas directly as proof certificates in the simpler and related problem of the absence of simulation, i.e. for theorems of the form ¬sim(p, q). In that case, the assertion language needs only the diamond modality • as well as the conjunction. More formally, let Act be a set of actions. The restricted set of assertions over Act is given by the recurrence A := i∈I a i A i , where I is a finite set and a i ∈ Act; that is, we have a strict alternation of (indexed) conjunctions and the diamond modality. The statement p |= i∈I a i A i means that, for every i ∈ I, there exists a q i such that p a i -→ q i and q i |= A i . We shall choose to write true for empty conjunctions and we can drop i∈I when I is a singleton. Thus, a true stands for i∈{•} a j∈{} b i, j A i, j . Some of the clerks and experts needed for this interpretation of an assertion as a certificate are listed below; the rest of the definition can be taken from the async constructor.

∀(a i ) i ∀(A i ) i ∀ j. decide L ( i a i A i , a j A j ).
∀a∀A. ∀ e ( a A, A, a).

∀A. ⊃ e (A, sync(stop), A). ∀T ∀A. ∀ e (A, A, T ).

∀a∀A. ν-unfold L ( a A, a A). ∀A. release L (A, A).

Example 9. In Figure 4, the process (6) is not simulated by the process (1): if Ξ is the assertion formula a ( b true ∧ c true), then 6 |= Ξ but 1 |= Ξ.

Assertions as non-bisimilarity certificates

It is possible to extend the FPC described in Section 5.2 to account for the absence of bisimulation in addition to the absence of simulation. As bisimilarity is finer than similarity, this will require a richer class of assertion formulas. The fact that it is a symmetric relation suggests that assertions should contain negations.

We could use full Hennessy-Milner logic (i.e. any arbitrary mix of • , [•], ∨ and ∧ or, equivalently, • , ∧ and ¬), but instead we choose the smaller but equivalent set of assertions defined by the following recurrence.

A := i∈I B i B := a i A i | ¬( a i A i )
It can be shown that this set characterizes the same relation as full Hennessy-Milner logic. The statement p |= i∈I B i means that, for every i ∈ I, p |= B i ; the statement p |= a A means that there exists a q such that p a -→ q and q |= A; and the statement p |= ¬( a A) means that p |= a A. Very little more is needed to extend the FPC to handle this certificate. We need to make sure that, in addition to certificates with a top-level • , decide L and ν-unfold L allow (and propagate) certificates with a top-level ¬ • . We also need an expert to consume ¬, and an expert to handle the additional ∧ - connective (see the definition of bisimilarity from Example 3). If we give these two roles to the same new expert, namely ∧ - e , the link between reflexivity and negations in the assertions appears clearly. The resulting set of clerks and experts for theorems of the form ¬bisim(p, q) is the following.

∀A. store L (A, A). ∀(B i ) i ∀ j. decide L ( i B i , B j ).
∀B. ν-unfold L (B, B). ∀a∀A. ∧ - e ( a A, a A, 1). ∀a∀A. ∀ e ( a A, A, a).

∀a∀A. ∧ - e (¬ a A, a A, 2). ∀T ∀A. ∀ e (A, A, T ).

∀A. ⊃ e (A, sync(stop), A).

∀A. release L (A, A).

∀A. ∃ c (A, λ x. A).

∀A. ∧ + c (A, A). ∀A. µ-unfold L (A, A).

∀A. ∨ c (A, A, A).

∀A. = s c (A, A).

This FPC extension is conservative, in that it can still check a certificate for non-simulation.

Example 10. In Figure 4, the processes (6) and (10) are similar but not bisimilar: if Ξ is the generalized assertion formula a ¬ b true, then 10 |= Ξ but 6 |= Ξ.

A reference proof checker

The framework for foundational proof certificates described in [START_REF] Miller | A proposal for broad spectrum proof certificates[END_REF][START_REF] Chihani | Foundational proof certificates in first-order logic[END_REF] was based on proof theory without fixed point definitions. In that setting, a standard logic programming language (in that case, λ Prolog [START_REF] Miller | Programming with Higher-Order Logic[END_REF]) was an ideal prototyping language for implementing and testing FPCs. The FPCs described in this paper are not so easily implemented in standard logic programming languages since the unification of eigenvariables must be done alongside the usual unification of "logic variables" that makes proof reconstruction possible. The implementation of λ Prolog, for example, considers eigenvariables as constants during unification.

We have built a prototype proof checker for testing the FPCs described in this paper using the Bedwyr extension to logic programming [START_REF]The Bedwyr model checker[END_REF][START_REF] Baelde | The Bedwyr system for model checking over syntactic expressions[END_REF]. That system, originally designed to tackle various kinds of model checking problems, provides the necessary unification for logic and eigenvariables along with backtracking search and support for λ -terms, λ -conversion and higher-order pattern unification.

One could have imagined implementing the non-augmented proof system µF 0 directly and in a sense, this is already done by Bedwyr itself. For example, if (µB t) is a purely positive fixed point encoding a Prolog predicate, when the system is given the sequent (µB t) ⇓ , it would emulate the Prolog search. Similarly, if it is given the sequent ⇑ (µB t) ⇑ , it would emulate a finite failing proof search. But, as anyone familiar with Prolog-style depth-first search knows, such proof search is limited in its effectiveness. For example, if one is attempting to prove that there is or is not a path between two points, a cycle in the underlying graph can make the search non-terminating. Bedwyr handles this with a loop-detection mechanism that can be embedded in the rules from Figure 2, making it a partial implementation of µF.

However, our goal with the µF proof system is not to use it by itself, but together with clerks and experts, as the engine (as a "kernel") for checking already existing proof evidence. Since the logic of the existing Bedwyr system has no native support for proof objects, we implemented µF a as an "object logic", without using some native features such as loop-detection. The Bedwyr specification files that we use (available directly at http://slimmer.gforge.inria.fr/bedwyr/pcmc/, or from the authors' homepages) are rather direct translations of the inference rules in Figures 1 and2 as well as of the various FPCs listed in the previous few sections. It has thus been easy for us to experiment and test FPCs.

While we have found the Bedwyr system to be useful for prototyping a proof checker, our proposal for FPC is not tied to any one particular implementation. Instead, the framework is defined using inference rules (such as found in Figures 1 and2). Any system that can implement the logical principles required by such inference rules can be used as a proof checking FPC kernel.

Conclusion

We have taken the basic structure of foundational proof certificates that had been developed elsewhere for first-order logic and described how it could be imposed on a logic based on fixed points. The resulting logic is much richer (think of the difference between first-order logic and first-order arithmetic) and additional logic principles need to be accounted for in the description of proof certificates.

In the areas of model checking that we have discussed, proof evidence is often taken to be, say, a path through a graph, a set of pairs of nodes (satisfying certain closure conditions), or a Hennessy-Milner logic assertion formula. We have illustrated how each of these familiar objects can be easily transformed into hints to guide a proof checker though the construction of a detailed and complete sequent calculus proof. The architecture of focused proof systems and the clerk and expert predicates allow this conceptual gap (between familiar proof evidence and sequent calculus proofs) to be bridged in a flexible and natural fashion.

We have also provided a novel look at the proof theory foundations of model checking systems by basing our entire project on the µMALL variant of linear logic and on the notion of switchable formulas. This latter notion seems to provide an interesting demarcation between the logically simpler notion of model checking and the more general notion of (inductive and co-inductive) deduction.

A On the augmented focusing proof system µF a It should be noted that, although the system presents two left rules for the connective =, one with the clerk = s c for success and one with the clerk = f c for failure, the implementation is usually expected to have one single unification facility that, given an equation, will or will not succeed, and which is tied to one single clerk. If the unification fails, the rule succeeds immediately without generating a premise certificate Ξ 1 , and the constraint on the conclusion certificate Ξ 0 is actually the same as for success. Hence the clerk = f c can be defined as an existential closure of = s c . Likewise, = s c can be defined in terms of = f c .

= f c (•) ≡ (∃Ξ 1 . = s c (•, Ξ 1 )) = s c (•) ≡ (∃Ξ 1 . = f c (•, Ξ 1 ))
It is also possible to remove the truth and falsity connectives, as we expect to have the equivalences ). This has the same effect as setting them to be always true. The system presented in Figures 1 and2 does not have these simplifications, but the Bedwyr-based implementation does.

B A simple example of a µF proof

The following proof can be seen as the justification that {1, 3} ⊆ {1, 2, 3}. In particular, encode these two sets as the predicates (i.e., abstractions over formula):

λ x[x = 1 ∨ x = 3] and λ x[x = 1 ∨ x = 2 ∨ x = 3].
The sequent calculus proof of inclusion can then be written as the following focused proof.

1 = 1 ⇓ 1 = 1 ∨ 1 = 2 ∨ 1 = 3 ⇓ ⇑ ⇑ 1 = 1 ∨ 1 = 2 ∨ 1 = 3 ⇑ 1 = 1 ∨ 1 = 2 ∨ 1 = 3 ⇑ ⇑ x = 1 x = 1 ∨ x = 2 ∨ x = 3 ⇑ 3 = 3 ⇓ 3 = 1 ∨ 3 = 2 ∨ 3 = 3 ⇓ ⇑ ⇑ 3 = 1 ∨ 3 = 2 ∨ 3 = 3 ⇑ 3 = 1 ∨ 3 = 2 ∨ 3 = 3 ⇑ ⇑ x = 3 x = 1 ∨ x = 2 ∨ x = 3 ⇑ ⇑ x = 1 ∨ x = 3 x = 1 ∨ x = 2 ∨ x = 3 ⇑ ⇑ [x = 1 ∨ x = 3] ⊃ [x = 1 ∨ x = 2 ∨ x = 3] ⇑ ⇑ ∀x.[x = 1 ∨ x = 3] ⊃ [x = 1 ∨ x = 2 ∨ x = 3] ⇑

  and its transitive closure: step a b. step b c. step c b. path X Y :-step X Y. path X Z :-sigma Y\ step X Y, path Y Z. We can translate the step relation into the binary predicate • -→ • defined by µ λ Aλ xλ y. ((x = a) ∧ + (y = b)) ∨ ((x = b) ∧ + (y = c)) ∨ ((x = c) ∧ + (y = b)) which only uses positive connectives. Likewise, path can be encoded as path: µ λ Aλ xλ z. x -→ z ∨ (∃y. x -→ y ∧ + A y z) In general, it is sensible to view any purely positive least fixed point expression as a predicate specified by Horn clauses. (For example, SOS rules for CCS are easily seen as Horn clauses.) Example 3. Let the ternary predicate • • -→ • describe a labeled transition system. It can be defined as a purely positive fixed point expression of the form µ λ Aλ pλ aλ q

Figure 3 :

 3 Figure 3: (Un)reachability problem

Example 4 .

 4 In Figure 3, (c) is reachable from (a), as witnessed by certificates like [b], [b; c; b], etc.

Figure 4 :

 4 Figure 4: Non-(bi)similar noetherian labeled transition systems

Example 8 .

 8 According to Figure 5, the set {(21, 23), (22, 24)} is a simulation and, therefore, the process (21) is simulated by the process (23). This corresponds to the following certificate co-inv λ xλ y. (x = 21 ∧ + y = 23) ∨ (x = 22 ∧ + y = 24), bipole .

Figure 5 :

 5 Figure 5: Non-noetherian labeled transition systems

t

  + ≡ (a = a) f + ≡ (a = b) f -≡ (a = a) t -≡ (a = b)if a and b are distinct constants, hence the following: is customary to leave clerks and experts out of rules with no premises (i.e.t + e , = s e , f - e , = f e , f + c , = f c , t - c and = f c

Acknowledgment. We thank the reviewers for their detailed and useful comments on an earlier draft of this paper. This work has been funded by the ERC Advanced Grant ProofCert.