Frédérique Bec
email: bec@ensae.fr.

Heino Bohn Nielsen
email: heino.bohn.nielsen@econ.ku.dk.

Sarra Saïdi
email: sarra.saidi@u-cergy.fr.

Christian Francq

Alain Guay

Alain Hecq

Søndergaard Rasmus

Anders Pedersen

Jean-Michel Rahbek

Zakoian

Mixed Causal-Noncausal Autoregressions: Bimodality Issues in Estimation and Unit Root Testing *

Keywords: JEL Classification: C12, C22, Q41 Mixed autoregression, non-causal autoregression, maximum likelihood estimation, unit root test, Brent crude oil price

This paper stresses the bimodality of the widely used Student's t likelihood function applied in modelling Mixed causal-noncausal AutoRegressions (MAR). It first shows that a local maximum is very often to be found in addition to the global Maximum Likelihood Estimator (MLE), and that standard estimation algorithms could end up in this local maximum. It then shows that the issue becomes more salient as the causal root of the process approaches unity from below. The consequences are important as the local maximum estimated roots are typically interchanged, attributing the noncausal one to the causal component and vice-versa, which severely changes the interpretation of the results. The properties of unit root tests based on this Student's t MLE of the backward root are obviously affected as well. To circumvent this issues, this paper proposes an estimation strategy which i) increases noticeably the probability to end up in the global MLE and ii) retains the maximum relevant for the unit root test against a MAR stationary alternative. An application to Brent crude oil price illustrates the relevance of the proposed approach.

Introduction

As emphasized by [START_REF] Hanfelt | Eliminating Multiple Root Problems in Estimation: Comment[END_REF] in his comment on [START_REF] Small | Eliminating Multiple Root Problems in Estimation[END_REF]'s paper, most statisticians doubt the proposition that multiple roots pose a serious problem in data analyses. Also, very few published data analyses in scientific journals mention the existence of multiple roots or describe what methods are used to select among them. Our paper is devoted to fill this gap for the case of mixed causal-noncausal autoregressive (MAR) models, and as a by-product we propose an improvement to the test for a unit root against MAR stationary alternatives. Basically, a maximum likelihood estimation strategy is proposed, which has the desirable feature that it selects the relevant maximum for the topic at hand.

Introduced decades ago in the literature of statistics (see for instance [START_REF] Weiss | Time-reversibility of linear stochastic processes[END_REF], [START_REF] Findley | On bootstrap estimates of forecast mean square errors for autoregressive processes[END_REF], [START_REF] Lawrance | Directionality and reversibility in time series[END_REF], [START_REF] Breidt | Time-reversibility, identifiability and independence of innovations for stationary time series[END_REF], Breidt, Davis, [START_REF] Lh | Maximum likelihood estimation for noncausal autoregressive processes[END_REF], Breidt, Davis and Dunsmuir [1992], [START_REF] Rosenblatt | A note on prediction and an autoregressive sequence[END_REF], [START_REF] Cambanis | On prediction of heavy-tailed autoregressive sequences: forward versus reversed time[END_REF], Cambanis and Fakhre-Zakeri [1996] or [START_REF] Rosenblatt | A note on prediction and an autoregressive sequence[END_REF]), mixed causal-noncausal autoregressions have recently known a revival of interest amongst researchers in economics and econometrics (see e.g. Lanne and Saikkonen [2011], [START_REF] Lanne | Bayesian Model Selection And Forecasting In Noncausal Autoregressive Models[END_REF], Lanne and Saikkonen [2013], [START_REF] Hencic | Noncausal Autoregressive Model in Application to Bitcoin/USD Exchange Rates[END_REF], Gouriéroux and Zakoian [2015], [START_REF] Gouriéroux | Filtering, Prediction and Simulation Methods for Noncausal Processes[END_REF], Gouriéroux and Zakoian [2017] or [START_REF] Fries | Mixed Causal-Noncausal AR Processes and the Modelling of Explosive Bubbles[END_REF]).

To fix ideas, let us introduce the MAR(r,s) model as formulated by Lanne and Saikkonen [2011],

φ(B)ϕ(B -1 )y t = t , (1)
where B is the backward shift operator (B k y t = y t-k for k = 0, ±1, ...) and φ(B) = 1 -φ 1 B -... -φ r B r , ϕ(B -1 ) = 1 -ϕ 1 B -1 -... -ϕ s B -s . Finally, t is a sequence of non-Gaussian independent, identically distributed random variables with mean zero, density function f ( t | λ) where λ is a set of parameters to be specialized later, and E( 2 t ) < ∞ unless otherwise mentioned. Indeed, if the error terms were Gaussian distributed, the model could be written indifferently as a backward or a forward autoregression, as these two representations are observationally equivalent asymptotically in this case. 1 Under the assumption that the polynomials φ(z) and ϕ(z) (z ∈ C) have roots outside the unit circle, it is well-known that y t has a stationary solution in terms of the two-sided moving average representation:

y t = ∞ j=-∞ ψ j t-j .
(2)

As can be seen from equation (2), MAR models allow for dependence on both the past and the future, by contrast with the well-known backward-looking autoregression which rules out dependence on future observations. If r > 0 and s = 0, the process defined in equation (1) becomes a purely causal process, namely the MAR(r,0) or AR(r):

φ(B)y t = t . (3)
If s > 0 and r = 0 in equation (1), one obtains the following purely noncausal process MAR(0,s):

ϕ(B -1 )y t = t . (4)
Papers by [START_REF] Hencic | Noncausal Autoregressive Model in Application to Bitcoin/USD Exchange Rates[END_REF], Gouriéroux and Zakoian [2015], [START_REF] Gouriéroux | Filtering, Prediction and Simulation Methods for Noncausal Processes[END_REF], Gouriéroux and Zakoian [2017] or [START_REF] Fries | Mixed Causal-Noncausal AR Processes and the Modelling of Explosive Bubbles[END_REF], assume Cauchy distributed disturbances in (1), i.e. very fat tailed distributions needed to capture bubble-like dynamics. For other macroeconomic variables such as the inflation or interest rates, a popular choice is the Student's t-distribution with scale parameter σ > 0 and ν > 2 degrees of freedom, see inter alia Lanne and Saikkonen [2011], Nyberg, [START_REF] Nyberg | Does Noncausality Help in Forecasting Economic Time Series?[END_REF], Lanne and Saikkonen [2013], [START_REF] Lof | Noncausality and asset pricing[END_REF] or Lof and Nyberg [2017]. With these Student's t-disturbances, the approximate maximum likelihood estimation (MLE) approach has been advocated by [START_REF] Breidt | Time-reversibility, identifiability and independence of innovations for stationary time series[END_REF], Andrews, [START_REF] Andrews | Maximum likelihood estimation for all-pass time series models[END_REF] and further promoted by Lanne and Saikkonen [2011]. The last two papers have derived conditions under which the MLE is consistent and asymptotically normal. The least absolute deviation (LAD) estimator has been proposed and its properties studied by [START_REF] Huang | Quasi-likelihood Estimation of Noninvertible Moving Average Processes[END_REF] and [START_REF] Wu | Least Absolute Deviation Estimation for General Autoregressive Moving Average Time-Series Models[END_REF] for the cases of unknown density function of the t 's in (1): they have characterized the properties of the Quasi-MLE when the objective function is the Laplace likelihood, which gives the LAD estimator. Alternative possible distributions include the Gaussian mixture or the generalized error distribution.

More recently, [START_REF] Hecq | Identification of Mixed Causal-Noncausal Models in Finite Samples[END_REF] have compared the MLE and LAD estimators in finite sample. There are signs of bimodality in their simulated illustration of finite sample likelihood estimation, see Figures 2 and3 2 . For sample sizes of 200 and 800, two peaks can be seen in the distribution of both LAD and ML estimates when the underlying disturbances are Student's t-distributed with 10 degrees of freedom. To our knowledge, this kind of issue has hardly been discussed so far in the literature. Yet, as noticed by [START_REF] Small | Eliminating Multiple Root Problems in Estimation[END_REF], estimating functions such as the first order conditions of the likelihood function or likelihood equations can have more than one root. Unfortunately, in practice it is often not obvious to determine which root is appropriate as a parameter estimate.

In what follows, we will specialize to the MAR(1,1) model to unveil the bimodality issue for ML estimators. Section 2 presents results to evaluate the magnitude of the bimodality phenomenon. Section 3 will evaluate two consequences of the issue, both stemming from the fact that the global and local maxima have approximately the same roots, but their location is interchanged between the causal and noncausal components: First, the interpretation of the processes generated from the global and local maxima are quite different. Second, the size and power of [START_REF] Saikkonen | Testing for a Unit Root in Noncausal Autoregressive Models[END_REF]'s unit root test against a stationary MAR alternative can be strongly impacted by the confusion in the location of the estimated roots. Section 4 proposes an estimation strategy which circumvents this bimodality issue and a modification of [START_REF] Saikkonen | Testing for a Unit Root in Noncausal Autoregressive Models[END_REF]'s unit root test so that its size and power are restored. Section 5 illustrates the benefits of the proposed approach using the Brent crude oil price while Section 6 concludes.

The bimodality issue

Let us first present the likelihood function to maximize in the remainder of the paper. We assume that t , t = 1, 2, ..., T , is a sequence of non Gaussian i.i.d. random variables with density function f ( t | λ). Lanne and Saikkonen [2011] propose to estimate the MAR model in (1) from the approximate log-likelihood function given by:

log L T (θ) = T -s t=r+1 log f ( t | λ), (5)
where t = φ(B)ϕ(B -1 )y t and θ = (φ, ϕ, λ). For the Student's t-case, considered in most of the paper, λ consists of the scale and the degree of freedom parameter, and the density is given by

f ( t | λ) = Γ( ν+1 2 ) Γ( ν 2 ) πνσ 2 -1 2 1 + 2 t σ 2 ν -ν+1 2 , (6)
with θ = (φ, ϕ, σ 2 , ν). The maximisation of log L T (θ) over permissible values of θ, θ ∈ Θ, gives the approximate MLE of θ.

Throughout this paper, we focus on a characterization of the shape of the likelihood function and let the parameter space for estimation, Θ, be unrestricted apart from σ 2 > 0 and ν > 0. An alternative would be to impose the assumptions for the analysis of the stationary case as in Lanne and Saikkonen [2011], |φ| < 1, |ϕ| < 1 and ν > 2, or the assumptions for the unit-root case as in [START_REF] Saikkonen | Testing for a Unit Root in Noncausal Autoregressive Models[END_REF], allowing φ = 1. This may be preferable for some empirical analyses, but to illustrate that bimodality often occurs with interchanged roots, in particular when the true backward root is unity, φ 0 = 1, we prefer to maximize the likelihood function for (φ, ϕ) unrestricted.

Limiting behaviour

First, Figure 1 reports the concentrated log-likelihood function as a function of ϕ, the noncausal parameter of the MAR(1,1) model, for different cases3

The simulated MAR processes are of length T = 50, 000 or larger, so as to consider the limiting behaviour of the likelihood.

Figure (A)-(D) report results for the Student's t-likelihood for a set of true values (ϕ 0 , φ 0 , ν 0 ). Figure (A) is characterized by a unit root in the backward lag polynomial, φ 0 = 1, and it shows very clearly that two maxima can be reached by maximizing the likelihood function-no matter the degrees of freedom of the Student's t-distribution. Moreover, the global maximum is reached for the true value of ϕ, namely ϕ 0 = 0.6, while the second local maximum is reached for ϕ ≈ 1, so that the unit root is wrongly located as the forward component of the MAR model. Finally, it is worth noticing that when ν 0 goes to infinity, corresponding to the Gaussian case, the two maxima reach exactly the same likelihood value. This is the well-known non-identification result, see [START_REF] Cambanis | On prediction of heavy-tailed autoregressive sequences: forward versus reversed time[END_REF] and Remark 1 below. By contrast, for smaller ν 0 values, the second maximum, interchanging the roots, is a local one only.

Figure (B) repeats the experiments in a stationary case which can still be considered as close to the unit root case: φ 0 = 0.9. The same conclusions hold, but it can be seen

ν 0 =4 ν 0 =10 ν 0 =∞ 0.4 0.6 0.8 1.0 -0.15 -0.10 -0.05 0.00 (A): log L c (ϕ), student: φ 0 =1, ϕ 0 =0.6 ν 0 =4 ν 0 =10 ν 0 =∞ ν 0 =4 ν 0 =10 ν 0 =∞
0.4 0.6 0.8 1.0 -0.10 -0.05 0.00 (B): log L c (ϕ), student: φ 0 =0.9, ϕ 0 =0.4

ν 0 =4 ν 0 =10 ν 0 =∞ ν 0 =4 ν 0 =10 ν 0 =∞
ν 0 =4 ν 0 =10 ν 0 =∞ ν 0 =4 ν 0 =10 ν 0 =∞ 0.4 0.6 0.8 1.0 -0.050 -0.025 0.000 (D): log L c (ϕ), student: φ 0 =0.4, ϕ 0 =0.9 ν 0 =4 ν 0 =10 ν 0 =∞ σ 0 2 =11 σ 0 2 =4 σ 0 2 =1
0.4 0.6 0.8 1.0 -0.10 -0.05 0.00 (E): log L c (ϕ), mixture: φ 0 =0.9, ϕ 0 =0.4

σ 0 2 =11 σ 0 2 =4 σ 0 2 =1 ν 0 =4 ν 0 =10 ν 0 =∞ 0.4 0.6 0.8 1.0 -0.10 -0.05 0.00 (F): log L c (ϕ), LAD: φ 0 =0.9, ϕ 0 =0.4 ν 0 =4 ν 0 =10 ν 0 =∞ Figure 1: Profile likelihood function, log L c (ϕ)
, for the Student's t, Gaussian mixture, and Laplace likelihood for different values of φ 0 and ϕ 0 . Calculated for T=50,000 or larger. The Gaussian mixture in (E) is given by an N(0,1) with provability 0.9 and N (0, σ 2 0 ) with probability 0.1, for different values of σ 2 0 . To make comparison easy, the plotted likelihood values are the deviations from the maximum close to the true values of the parameters.

that the width of the set of starting values which will lead to the wrong optimum is much larger than in the unit root case4 .

The simulation experiments reported in Figure (C) reveal that the bimodality issue may disappear for particular parameter combinations. Actually, with an even lower backward root, here set to φ 0 = 0.7, fat-tailed Student's t-distributions make the profile likelihood unimodal. Nevertheless, the bimodality reappears as ν 0 increases.

Figure (D) illustrate the likelihood surface for the opposite case, ϕ > φ. The profile

ϕ 0 × φ 0 0.0 0.2 0.4 0.6 0.8 1.0 0.5 1.0 (A): Student's t, ν 0 =4 ϕ 0 × φ 0 ϕ 0 × φ 0 0.0 0.2 0.4 0.6 0.8 1.0 0.5 1.0 (B): Student's t, ν 0 =10 ϕ 0 × φ 0 ϕ 0 × φ 0 0.0 0.2 0.4 0.6 0.8 1.0 0.5 1.0 (C): Student's t, ν 0 =20
ϕ 0 × φ 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.5 1.0 (D): Gaussian To emphasize that the bimodality issue is related to the MAR model per se, and not to the particular choice of a Student's t-distribution, Figure (E) shows the profile likelihood for the case where f ( t | λ) is chosen as a mixture of Gaussian distributions (which is also guaranteed to be leptokurtic). The three cases have been chosen to mimic the shape of the tails of the Student's t-distribution with ν 0 ∈ {4, 10, ∞}, and the results are very similar to the results for the Student's t-likelihood.

ϕ 0 × φ 0 ϕ 0 × φ 0
Finally, Figure (F) shows the profile likelihood when f ( t | λ) is chosen as a Laplace distribution, leading to the LAD estimator. We generate the data using the Student's t-distribution to emphasize the role of LAD estimation as a robust QML estimator when the true distribution of the errors is unknown. Again we observe a pronounced bimodality.

The two following remarks sum up these simulation experiments and give conditions on the parameters for bimodality to occur. We first establish analytically the bimodality of the likelihood function in the Gaussian case. Here the bimodality is a reflection of φ and ϕ being unidentified as T → ∞, which has already been stressed from the simulations experiments reported in Figure 1 above. We next emphasize some conditions on the parameters of the MAR(1,1) DGP which increase the probability of bimodality occurence.

Let us begin with the following remark, whose derivations are given in the Appendix:

Remark 1 Consider the MAR(1,1) model in (1) with Gaussian errors and

|ϕ 0 | < 1, |φ 0 | < 1. Let log L ∞ (φ, ϕ) denote the limit of the Gaussian likelihood function as T → ∞. If (ϕ, φ) = (ϕ 0 , φ 0 ) is a local maximum of log L ∞ (φ, ϕ), then the interchanged param- eter point, (ϕ, φ) = (φ 0 , ϕ 0 ), is also a maximum, and log L ∞ (ϕ 0 , φ 0 ) = log L ∞ (φ 0 , ϕ 0 ).
Let us now figure out what are the conditions of the MAR(1,1) model parameters which increase the bimodality case probability. Figure 2 plots in red, for ν 0 ∈ {4, 10, 50, ∞}, all the couples (φ 0 ,ϕ 0 ) ∈ {0.0, 0.01, • • • , 0.99, 1.00} 2 for which bimodality with interchanged roots occurs. It amounts to 10,000 gridpoints (φ 0 ,ϕ 0 ) for each ν 0 . This is done with T=5,000 in order to approximate the limiting likelihood function.

Remark 2 As can be seen from Figure 2, the frequency of bimodality occurence increases if: i) the fatness of the tails decreases, ii) the causal and/or noncausal roots approach one from below, iii) the difference between the causal and non causal roots increases.

Finite sample behaviour

Let us now have a look at the finite-sample properties of the likelihood function and the MLE. First, Figure 3 shows profile likelihoods for smaller samples, T ∈ {100, 250, 1000}, corresponding to some of the limiting cases in Figure 1. Unlike the limiting case, here approximated by T=50,000, the profile likelihood for finite T depends on the specific realization of the simulated MAR time series. As such these are only illustrative examples.

It is worth noting in Figure 3, that in small samples the global maximum might correspond to the case of interchanged roots. This can be seen from the examples in panels (A), (B), and (C), where the likelihood for T = 100 is larger for the interchanged parameter points than for points close to the true values.

Also observe that even if the limiting likelihood function is bimodal, the finite sample likelihood may be unimodal. This happens for one realization for T = 100 in panel (B). Panel (C) illustrates the reverse case, where the limiting likelihood is unimodal, while the small sample likelihoods for T = 100 and T = 250 are bimodal.

To illustrate the occurrence of these issues, Table 1 reports the results of simulation experiments for ν 0 = 4 or ν = 10. The same DGPs as the ones considered in Figure 1 have been simulated 5,000 times for smaller sample sizes: T ∈ {100, 250, 500, 1000, 10000}. The global maximum is found using a grid search over initial values, which will be described in Section 4 below. The columns labelled "% Bimodal" report the percentage of the corresponding model estimates leading to two or more maxima. The columns labelled "% Roots interchanged" report the percentage of global maxima which yields interchanged roots, φ > φ.

We observe that the finite sample likelihood function is frequently bimodal, also in cases where the limiting likelihood is unimodel. Even fat tailed disturbances associated with clearly stationary roots (φ 0 = 0.7 and ϕ 0 = 0.4, or φ 0 = 0.4 and ϕ 0 = 0.7 ) raise bimodality and interchanged roots issues in small samples, but the issue disappears as T → ∞. For instance, when T = 100 and ν 0 = 4, 20.4% of the DGPs display two maxima, and the roots are interchanged in around 15% of the global MLE. With ν 0 = 10, these percentages are much higher, even for T = 500 where bimodality is found in about 42% of the cases and still around 15% of the global maxima interchange the roots location. As φ 0 ↑ 1, the bimodality occurrences strongly increase -it reaches 100% or so as T increases -while the percentage of the interchanged roots is not so much affected. Besides, the latter is a decreasing function of the sample size.

Its consequences

As discussed earlier, the roots obtained from global and local maxima are often interchanged between backward and forward components. Interpretation is strongly affected and so is [START_REF] Saikkonen | Testing for a Unit Root in Noncausal Autoregressive Models[END_REF]'s unit root test.

Impact on interpretation

As already mentioned, an obvious consequence of the roots interchange concerns the interpretation of the series at hand. To illustrate this very simply, let us consider two MAR(1,1) processes with T=200, the scale parameter σ 2 0 = 1, and ν 0 = 4 degrees of freedom. In Model 1 we set (φ 0 , ϕ 0 ) = (0.9, 0.4) while for Model 2 we interchange the roots, (φ 0 , ϕ 0 ) = (0.4, 0.9). These simulated processes, submitted to the exact same sequence of t-distributed shocks, are plotted in Figure 4.

Figure (A) represents the process with a large (respectively small) backward (resp. forward) root and as such, one may recognize the usual behaviour of AR processes. Following the strong negative shock which occurs a little bit before the 70th realization, it plunges deeply from a slightly positive value to around -30, and then progressively absorbs the shock over the thirty next periods.

For Model 2, in Figure 4 (B), where the roots are interchanged, the process anticipates the shock by decreasing steadily almost from 20 periods before on, and then recovers immediately after the adverse shock occurs. For the positive large shock which occurs around the 130th realization, Model 2 again anticipates it and build up a bubble like dynamics since the 110th realization approximately, and then burst right after the shock happens. Again, a mirroring dynamics is generated by model one which is "surprised" by the occurrence of the favourable shock and increases abruptly instantaneously, before the shock effect slowly vanishes with the time passing by.

This difference has already been widely commented on when comparing pure causal and non causal dynamics. What we want to emphasize here, is that the same can be Model 1 has (φ 0 , ϕ 0 ) = (0.9, 0.4) while model 2 has (φ 0 , ϕ 0 ) = (0.4, 0.9).

observed in mixed causal-noncausal processes.

The uncertainty in the location of the forward and backward roots implied by the bimodality issue does not only affect the direct interpretation of the dynamics, e.g. the likelihood of bubble-like dynamics occuring, but may also have severe consequences for forecasts based on the estimated MAR models.

Impact on unit root testing

Another severe consequence is related to the size and power of [START_REF] Saikkonen | Testing for a Unit Root in Noncausal Autoregressive Models[END_REF]'s unit root test. Their test statistic is given by:

τ T = T ( φ -1) G φ T ( θ)
, where G φ T ( θ) is the element of the inverse standardized Hessian matrix corresponding to φ. Their Proposition 3, p.107, derives the limiting distribution of this statistics under the unit root null φ 0 = 1 and shows that it is not nuisance parameters free. For the Student's t-case, the limiting distribution depends on

J = ν 0 (ν 0 + 1) (ν 0 -2)(ν 0 + 3) ,
where ν 0 > 2 is the degrees of freedom parameter. With this in mind, the issue raised by bimodality and interchanged roots can be easily explained. Intuitively, every time the largest root is wrongly located in the forward component, such that φ ≈ 1 and φ 1, the unit root test rejects the null φ 0 = 1. Hence, the test's size and power will be spuriously increased. This is illustrated in Table 2, which reports in the empirical rejection rates of the null by the τ T statistics for [START_REF] Saikkonen | Testing for a Unit Root in Noncausal Autoregressive Models[END_REF]'s unit root test. The authors consider three cases depending on the deterministic specification, and we focus here on case 2, for which the series is demeaned before estimation and unit root testing is performed. The critical values used in this Table have been computed using Table 1, p.110, in [START_REF] Saikkonen | Testing for a Unit Root in Noncausal Autoregressive Models[END_REF]. Notes: ERF stands for empirical rejection frequency. All figures are computed from 5,000 simulated MAR models. The critical values for the unit root test have been calculated from Table 1 in [START_REF] Saikkonen | Testing for a Unit Root in Noncausal Autoregressive Models[END_REF] using the estimated degrees of freedom in each replication, ν.

Top and bottom panels of Table 2 show the results obtained when ν 0 = 4 and ν 0 = 10, respectively. The left panel reports the empirical rejection frequency (ERF) when the global maximum estimates are used to perform the τ T unit root test.

The right panel labelled "Preferred maximum" performs the unit root test after the maximum with φ close to one is eliminated, when there is more than one maximum. Indeed, this case doesn't make sense from an economic point of view, as will be further discussed in Section 4.2. below. Again, the MLE is obtained using a grid search over initial values which will be described in the next section. The impact of interchanged roots on this unit root test size is of course bigger when the percentage of reversed roots is higher, which can be seen from the fourth columns under the "Global maximum" label: In the case where ν 0 = 4, with T=100, 14.6% of the global maximum interchange the roots which yields ERFs well above their nominal levels. For instance, considering the τ T unit root test at the 5% level gives an ERF of 16.3% instead of 5%. The issue becomes worse in the bottom panel when ν 0 = 10, since at the 5% level and for T=100, one obtains an ERF of 33.5%. With this ν 0 , the issue persists for larger samples too. For instance, when T=500, the ERF corresponding to the nominal 5% level of the test is still 18.5%.

Regarding power, Table 3 compares the empirical rejection rates of the null obtained from our estimation method -column labeled "Preferred maximum ERF" -to the ones given in [START_REF] Saikkonen | Testing for a Unit Root in Noncausal Autoregressive Models[END_REF], Fig. 2,p.114 therein. As expected, it can be seen from the last two columns that for small sample sizes, the [START_REF] Saikkonen | Testing for a Unit Root in Noncausal Autoregressive Models[END_REF]'s ERFs are larger than the ones of our modified version of the test by 10 to 15 percentage points. Nevertheless, this power issue becomes smaller and vanishes Notes: ERF stands for empirical rejection frequency. All figures are computed from 5,000 simulated MAR models. The ERFs for [START_REF] Saikkonen | Testing for a Unit Root in Noncausal Autoregressive Models[END_REF] unit root test are the ones reported in Fig. 2, p.114 therein.

when the sample size increases, as does the bimodality issue. In Section 5, we illustrate the benefits from implementing our approach by applying it to the Brent crude oil price series.

4 The proposed approach

Estimation strategy

The issue of parameter estimation based on multimodal criteria functions has a long history in the statistical literature. [START_REF] Barnett | Evaluation of the maximum likelihood estimator where the likelihood equation has multiple roots[END_REF] reviews different approaches, and although the computational possibilities have changed enormously over the last fifty years, the main principles of the possible approaches are still the same.

It is well-known that gradient ascent methods, such as the Newton estimation algorithm, will converge to a local maximum of the criteria function -highly depending on the initial values for the estimation, θ (0) say. One approach in the literature is therefore to use estimation procedures that seek the global maximum by using algorithms allowing also downhill movements. Examples of such algorithms include simulated annealing (e.g. [START_REF] Kirkpatrick | Optimization by Simulated Annealing[END_REF]), genetic algorithms (e.g. [START_REF] Goldberg | Genetic Algorithms in Search, Optimization and Machine Learning[END_REF]), or particle swarm methods (e.g. [START_REF] Zhang | A comprehensive survey on particle swarm optimization algorithm and its applications[END_REF]).

In some cases, however, it is of interest to characterize the entire likelihood surface, and to list the different local maxima. An example of this case is the MAR model with a root close to unity, where a forward root in the vicinity of unity can be excluded by economic reasoning, see e.g. [START_REF] Saikkonen | Testing for a Unit Root in Noncausal Autoregressive Models[END_REF]. Selection between local maxima is also sometimes done in Markov-Switching modelling where a preferred "reasonable" model is chosen among the set of local maxima.

To obtain all local maxima for a given criteria function the standard approaches are either i) the use of randomized initial values, e.g. with θ (0) drawn uniformly over a subset of the parameter space or ii) the use of a deterministically specified grid on the parameter space. In this paper we opt for the grid-search, and based on the insight from the characterization of the likelihood function for finite T as well as the limit T → ∞ above, we choose a coarse but wide grid. In particular, we choose5 (ϕ (0) , φ (0) , σ 2(0) , ν (0) ) ∈ {0.05, ..., 0.95} 2 × {0.5q, ..., 2q} × {3, ..., 50},

where q is a robust measure of the scale, here taken to be the inter-quartile range of the residuals obtained from a causal AR(2). For the simulations, we use four points for the grid on ϕ and φ and three points on σ and ν, as it has proven to be a good compromise between a quickly increasing number of starting values and the need to have enough of them to detect the bimodality. Because of the striking tendency of the likelihood estimation to produce interchanged roots, we consider -in addition -a sequence of initial values given by

(ϕ (0) , φ (0) , σ 2(0) , ν (0) ) = ( φ, φ, σ2 , sν), with s ∈ {0.1, 0.5, 1, 2, 5}, (8)
where ( φ, φ, σ2 , ν) are the (possibly local) estimates obtained from our standard initial values: Note that the backward and forward estimates are interchanged on purpose to initialize the new starting values. In total, this design of starting values covers B = 4 2 3 2 + 6 = 150 points in the parameter space6 . For each starting point we apply the quasi-Newton optimization algorithm BFGS, see e.g. [START_REF] Nocedal | Numerical Optimization[END_REF], and record the unique maxima.

As documented in Table 1, the grid-search algorithm often finds two maxima corresponding to cases with interchanged roots. In addition, it sometimes identifies local maxima with very large ν, corresponding to the maxima of the Gaussian likelihood function. The Gaussian local maxima have much lower likelihood in all cases, and will not be discussed further in this paper7 .

A modified unit root test

For estimation of the MAR model in Sections 2 and 3, we used this grid-search over initial values in order to characterize the entire likelihood surface and identify the global maximum. The grid-search procedure yields a list of unique maxima, and the MLE is taken to be the candidate with the highest maximum. As T → ∞, this global maximization would obviously produce a consistent estimator.

For very persistent processes, however, where one root is in the vicinity of unity, the situation is slightly different. According to Table 1, there is a large probability of observing interchanged roots in finite samples, meaning that even if the true parameter set induces a causal unit root, φ 0 = 1 and ϕ 0 1, both case (a) : φ ≈ 1 and φ 1 and case (b) : φ 1 and φ ≈ 1 are likely to occur. Case (b), however, is not reasonable from an economic point of view, because as φ ↑ 1 it implies that the behavior of the estimated model would be given by

y t = φy t-1 + v t , with v t = ∞ i=0 t+i , (9)
with the interpretation that agents would have to look into the infinite future without discounting.

In light of this lack of a reasonable economic interpretation when φ ≈ 1 and in light of the tendency towards interchanged roots, we therefore suggest to select amongst the local maxima whenever max( φ, φ) ≈ 1. In this case we prefer models with φ > φ, even if this maximum is only local and there exists a maximum with φ < φ ≈ 1 and a higher likelihood. This is similar to the approach sometimes retained in Markov switching models, where a "reasonable" solution is chosen within the set of local maxima.

To be precise, we suggest the following procedure:

Estimation procedure: Do a full search of local maxima for the MAR model, see Section 4.1 above. i) If there is only one maximum of the likelihood function, it is selected as the preferred. ii) If there are multiple maxima, and max( φ, φ) ≈ 1, the maxima with φ > φ is the preferred one.

Based on this preferred maximum, we implement the unit root test of [START_REF] Saikkonen | Testing for a Unit Root in Noncausal Autoregressive Models[END_REF] -thereby avoiding the spurious rejections whenever the roots are interchanged by random finite-sample estimation uncertainty. Remark that this procedure is feasible in practice and that the simulated size and power properties for the [START_REF] Saikkonen | Testing for a Unit Root in Noncausal Autoregressive Models[END_REF] unit root test reported above are the ones to be expected in empirical applications. This is unlike the reported results from the simulations in [START_REF] Saikkonen | Testing for a Unit Root in Noncausal Autoregressive Models[END_REF], where the correct maximum is obtained in the majority of replications because the estimation has been initialized at the true values of the DGP, such that the issue of bimodality is unsubstantial. This approach is obviously infeasible in practice.

Right part of Table 2 shows the empirical rejection frequencies for the unit root test based on our chosen maximum with φ > φ. It can be seen that our proposed estimation and test strategy clearly outperforms the results obtained from the global maximum estimation -left panel of the Table. This is especially the case for small sample sizes T ∈ {100, 250} and ν 0 = 10. It is also worth noticing that the size gets close to the nominal as T → ∞.

Application to Brent crude oil prices

In what follows, we provide an example which illustrates the bimodality-induced weaknesses of [START_REF] Saikkonen | Testing for a Unit Root in Noncausal Autoregressive Models[END_REF]'s unit root test and the benefits of our proposed approach.

The data and unit root test

The data we analyse are monthly Brent crude oil prices, expressed in US dollars per barrel, from January 2004 to December 2018. The series8 is plotted in Figure 5.

First, we find out that an AR(2) model is enough to eliminate any serial correlation in the residuals up to order 24. Then, the normality of the residuals is strongly rejected according to the Jarque and Bera statistics, with a test statistics of 8.44, which rejects the null of normality at the 1% level.9

Let us now turn to the unit root tests. From an ADF unit root test with two lags in levels and including an intercept, we find a statistics value of -2.63, corresponding to a p-value of 9%. Hence, the unit root is not rejected at the 5% level, but the conclusion is not so clear cut either. Let us now turn to [START_REF] Saikkonen | Testing for a Unit Root in Noncausal Autoregressive Models[END_REF] unit root test. The results reported in Table 4 illustrate the bimodality of the Student's t ML estimation and the misleading conclusions it can induce.

Actually, as can be seen from this Table, two maxima co-exist for this series from the ML estimation. For this application, we have not used the method described earlier first, but have purposely started from the most intuitive (or neutral) values for the two roots, i.e. 0.5 for each of them. As can be seen from the first line of the results, this leads to a local maximum, with a lower log-likelihood than the one of the global maximum, and interchanged roots. As a result, [START_REF] Saikkonen | Testing for a Unit Root in Noncausal Autoregressive Models[END_REF]'s unit root test wrongly concludes to a strong rejection of the null. As a matter of fact, the global maximum is reached when the values of φ and ϕ are initialized to 0.05 instead. In this case, reported [START_REF] Saikkonen | Testing for a Unit Root in Noncausal Autoregressive Models[END_REF].

in the second line of the results, the largest root is located in the backward component which in turn leads the unit root test not to reject the null. This is also the maximum which is found by our modified version of [START_REF] Saikkonen | Testing for a Unit Root in Noncausal Autoregressive Models[END_REF]'s unit root test.

MAR models estimates

Table 5 reports the estimation results for all Student's t distributed MAR models where r + s = 2, namely the purely backward AR(2,0), the purely forward AR(0,2), and both the global and local ML estimates, denoted respectively MAR G (1,1) and MAR L (1,1). [START_REF] Fries | Mixed Causal-Noncausal AR Processes and the Modelling of Explosive Bubbles[END_REF] suggest to select between candidate models based on the independence of residuals from estimated causal and noncausal models, see also [START_REF] Cavaliere | Bootstrapping Noncausal Autoregressions: With Applications to Explosive Bublle Modeling[END_REF]. As an alternative to the maximum likelihood criterion, we have also tried several independence tests applied to estimated residuals in order to select between models with interchanged roots. For larger samples, some independence tests, e.g. the BDS test, see [START_REF] Brock | A Test for Independence Based on the Correlation Dimension[END_REF] were found to do rather well, but the overall conclusion was that for small and moderate sample sizes (with less than 500 observations, say) the model selection based on maximizing the likelihood outperforms any independence test. Accordingly, all candidate models have been estimated from the exactly same period so that their estimated likelihoods are perfectly comparable, i.e. from 2004m3 to 2018m10. These results confirm the mixed nature of this autoregressive process, with an estimated backward root smaller than the estimated forward root. Even though smaller, the latter is still significantly different from zero. This noncausal component might be interpreted as capturing the epochs of bubble build-up and burst, as well as nonfundamentalness of shocks. The latter can in turn be seen as evidence that the econometrician uses less information than economic agents do.

Conclusion

This paper aimed at emphasizing and circumventing the undesirable consequences of bimodality of the widely used Student's t likelihood function applied in modelling mixed causal-noncausal autoregressions. The main consequence is that the MLE can end up in a local maximum where the backward and forward roots are interchanged, hence impacting severely the interpretation of the series dynamics. Another important consequence of this bimodality concerns the unit root test procedure developed by [START_REF] Saikkonen | Testing for a Unit Root in Noncausal Autoregressive Models[END_REF], as interchanged roots will mislead the test: any time the largest root is wrongly located in the forward component of the MAR model's estimate, the test on the backward root will wrongly reject the unit root null. This affects significantly both the size and power of the test. For this reason, an estimation strategy relying on grid search is proposed so as to characterize the entire likelihood surface and list all local maxima. In the particular context of unit backward root testing, this can in turn be used to eliminate the meaningless maxima where φ ≈ 1. This proposed strategy is shown to correct both size and power of [START_REF] Saikkonen | Testing for a Unit Root in Noncausal Autoregressive Models[END_REF] unit root test very efficiently. The relevance of our proposed approach is illustrated from the Brent crude oil price for which a naive choice of starting values for the roots would have led to a local minimum with interchanged roots, hence leading to reject the unit root null wrongly. Of course, this paper is only a first insight into this bi-or multi-modality issue in the MAR estimation and testing literature. The extension of our results to the more general MAR(r,s) model is on our research agenda.

  : log L c (ϕ), student: φ 0 =0.7, ϕ 0 =0.4

Figure 2 :

 2 Figure 2: Bimodality occurence as a function of the parameters values, T=5,000.

Figure 3 :

 3 Figure 3: Profile likelihood function, log L c (ϕ), for the Student's t likelihood for finite samples, with ν 0 = 10. To make comparison easy, the plotted likelihood values are the deviations from the maximum close to the true values of the parameters.

Figure 4 :

 4 Figure4: Simulated MAR(1,1) processes. Model 1 has (φ 0 , ϕ 0 ) = (0.9, 0.4) while model 2 has (φ 0 , ϕ 0 ) = (0.4, 0.9).

Figure 5 :

 5 Figure 5: Brent crude oil prices, US dollars per barrel.

Table 1 :

 1 The extent of the bimodality issue in finite sample

				ν 0 = 4	ν 0 = 10
			%	% Roots	%	% Roots
	T	φ 0 ϕ 0	Bimodal interchanged	Bimodal interchanged
	100	1 0.6	60.6	14.6	74.3	37.7
	250	1 0.6	85.1	2.6	97.8	27.5
	500	1 0.6	97.4	0.1	99.9	15.8
	1000	1 0.6	99.9	0.0	100.0	4.7
	10000	1 0.6	100.0	0.0	100.0	0.0
	100 0.9 0.4	58.3	12.1	71.8	37.4
	250 0.9 0.4	61.2	1.6	90.0	24.4
	500 0.9 0.4	63.0	0.1	97.2	12.1
	1000 0.9 0.4	65.9	0.0	99.7	3.5
	10000 0.9 0.4	82.2	0.0	100.0	0.0
	100 0.7 0.4	20.4	15.0	42.9	39.1
	250 0.7 0.4	6.4	3.2	45.6	27.7
	500 0.7 0.4	1.0	0.3	41.7	15.0
	1000 0.7 0.4	0.0	0.0	32.7	6.0
	10000 0.7 0.4	0.0	0.0	0.6	0.0
	100 0.4 0.7	20.4	14.6	42.6	38.0
	250 0.4 0.7	6.7	2.8	45.8	26.9
	500 0.4 0.7	0.9	0.4	41.6	15.4
	1000 0.4 0.7	0.0	0.0	32.4	6.0
	10000 0.4 0.7	0.0	0.0	0.7	0.0

Notes: "% Bimodal" refers to the percentage of the model estimates leading to at least two maxima. "% Roots interchanged" reports the percentage of global maxima which yields interchanged roots, φ > φ. All figures are computed from 5,000 drawings of the Student's t MAR DGPs.

Table 2 :

 2 Impact of bimodality on Saikkonen-Sandberg's unit root test size

	T	Global maximum		Preferred maximum
		% Roots	τ T ERF in %	% Roots	τ T ERF in %
		interchanged	1% 5% 10%	interchanged	1% 5% 10%
			(φ 0 , ϕ 0 , ν 0 ) = (1,0.6,4)
	100	14.6	11.8 16.3 20.5	3.8	2.5 6.0	9.9
	250	2.6	3.5 6.2 11.1	0.0	0.9 3.6	8.5
	500	0.1	1.1 4.5	9.0	0.0	1.0 4.4	9.0
	1000	0.0	0.8 4.6	9.6	0.0	0.8 4.6	9.6
			(φ 0 , ϕ 0 , ν 0 ) = (1,0.6,10)
	100	37.7	28.0 33.5 37.2	11.3	4.8 8.5 11.4
	250	27.5	27.5 28.8 31.5	0.4	0.5 1.8	4.5
	500	15.8	16.0 18.5 22.9	0.1	0.3 2.8	7.3
	1000	4.7	5.2 8.2 12.8	0.0	0.5 3.6	8.4

Table 3 :

 3 Impact of bimodality on Saikkonen-Sandberg's unit root test power

	T	Preferred maximum, ERF in %	SS, ERF in %
		(φ 0 , ϕ 0 , ν 0 ) = (0.85,0.5,4)
	100	58	73
	250	98	100
	500	100	100
		(φ 0 , ϕ 0 , ν 0 ) = (0.9,0.5,4)
	100	39	54
	250	96	99
	500	100	100
		(φ 0 , ϕ 0 , ν 0 ) = (0.95,0.5,4)
	100	18	27
	250	68	80
	500	99	100

Table 4 :

 4 Modified τ T unit root test for demeaned series Notes: All estimations use the 176 observations between 2004m3 and 2018m10. Starting values for the scale and ν parameters are 2.5 and 10 respectively. The critical values reported here have been computed using case 2 of Table 1, p.110, in

	Starting values	Estimation results		Unit root test
	φ (0)	ϕ (0)	Log-lik	φ	φ	τ T	5% c.v. 10% c.v.
	0.5	0.5	-555.70 0.46 0.92	-5.26	-2.83	-2.54
	0.05	0.05	-551.44 0.94 0.40	-2.26	-2.81	-2.51

Table 5 :

 5 MAR estimation results for demeaned Brent crude oil price series.Notes: Estimations use 176 effective observations. MAR G and MAR L refer to global and local MLE respectively. Standard errors are given in () and J-B p-val. gives the p-value of the Jarque and Bera normality test.

				Specifications
			AR(2,0)	AR(0,2)	MAR G (1,1)	MAR L (1,1)
	φ 1	1.36		0.94	0.46
			(0.07)		(0.03)	(0.10)
	φ 2	-0.39
			(0.07)
	ϕ 1		1.35	0.40	0.92
				(0.07)	(0.08)	(0.04)
	ϕ 2		-0.38
				(0.07)
	σ	27.74	25.17	22.46	25.83
	√	ν	(4.54) 4.26	(4.64) 2.88	(4.59) 2.56	(4.49) 3.02
			(2.39)	(0.89)	(0.71)	(0.94)
	Log-lik	-551.98	-555.39	-551.44	-555.70
	J-B p-val.	0.03	<0.01	<0.01	<0.01

See e.g. the discussion inCambanis and Fakhre-Zakeri [1996].

See pages

317-318 therein.

Throughout this paper we consider non-zero mean processes and follow[START_REF] Saikkonen | Testing for a Unit Root in Noncausal Autoregressive Models[END_REF] by demeaning the series before estimation. All results below do not depend on the demeaning or detrending of the data.

This stems from the convergence of the estimator at rate T if φ 0 = 1 and at rate √ T if |φ 0 | < 1.

None of the models in this paper have negative roots, but if negative roots are likely, the grid could be easily extended.

We have also considered alternative and less computationally expensive grid-designs, e.g. the interchanged roots in (8) alone, but it did not consistently locate all local maxima.

Anyway, remind that in this case φ and ϕ are unidentified as T → ∞.

The time series is taken from FRED Economic Data website and is named MCOILBRENTEU.

All results regarding normality tests have been double-checked using the robust version of the Jarque and Bera statistic proposed by[START_REF] Gel | A robust modification of the Jarque-Bera test of normality[END_REF] and the classical Anderson-Darling statistic: all conclusions are confirmed.

Bec and S. Saïdi acknowledge financial support from the Labex MME-DII. H. B. Nielsen acknowledges financial support from the Institute for Advanced Studies at the University of Cergy-Pontoise and the Danish Council for Independent Research, DSF Grant 015-00028B.

A Derivations for Remark 1

The density of the Gaussian distribution is given by

and the approximate average likelihood function for the Gaussian MAR(1,1) model is

For simplicity we focus on (φ, ϕ) and fix

Under the stationarity condition, |ϕ 0 | < 1 and

and the limit of the likelihood per observation is

(A.4)

The limiting likelihood function in (A.4) is symmetric in φ and ϕ, which gives the results i Remark 1.


Latifa Greche
email: latifa.greche@usmba.ac.ma

Mohamed Akil
email: mohamed.akil@esiee.fr

Rostom Kachouri

Najia Es-Sbai

A new pipeline for the recognition of universal expressions of multiple faces in a video sequence

Keywords: Facial expression recognition, Histogram of Oriented Gradient, Linear Discriminant Analysis

Facial Expression Recognition (FER) is a crucial issue in human-machine interaction. It allows machines to act according to facial expression changes. However, acting in real time requires recognizing the expressions at video speed. Usually, the video speed differs from one device to another. However, one of the standard settings for shooting videos is 24 fps. This speed is considered as the low-end of what our brain can perceive as fluid video. From this perspective, to achieve a real-time FER the image analysis must be completed, strictly, in less than 0.042 second no matter how the background complexity is or how many faces exists in the scene. In this paper, a new pipeline has been proposed in order to recognize the fundamental facial expressions for more than one person in real world sequence videos. First, the pipeline takes as input a video and performs a face detection and tracking. Regions of Interest (ROI) are extracted from the detected face in order to extract the shape information when applying the Histogram of Oriented Gradient (HOG) descriptor. The number of features yield by HOG descriptor is reduced by means of a Linear Discriminant Analysis (LDA). Then, a deep data analysis was carried out, exploiting the pipeline, for the objective of setting up the LDA classifier. The analysis aimed at proving the suitability of the decision rule selected to separate the facial expression clusters in the LDA training phase. To conduct our analysis, we used ChonKanade(CK+)

Introduction

The early studies and analysis of facial expressions were studied only in the psychological field where Charles Darwin [START_REF] Darwin | The expression of the emotions in man and animals[END_REF] demonstrated for the first time that despite the differences between individuals and cultures, six facial expressions are innate and common among human beings, which are joy, fear, sadness, surprise, anger and disgust. With the emergence of new technologies and robots, the psychologist provided the machine learning and computer science community with labeled databases [START_REF] Georghiades | Yale face database. Center for computational Vision and Control at[END_REF][START_REF] Pantic | Web-based database for facial expression analysis[END_REF][START_REF] Lucey | The extended cohn-kanade dataset (ck+): A complete dataset for action unit and emotion-specified expression[END_REF] which contain a set of face images of universal expressions or a set of micro-expressions [START_REF] Ekman | Measuring facial movement[END_REF] information. Facial expression recognition has many applications in security, human-machine interaction, and medicine. It is an active area of research where different solutions have been proposed over the years. Generally, proposed solutions can be grouped in two categories. The first category uses a video sequence to build dynamic models for facial expressions. This category is based on approaches that extracts spatio-temporal features [START_REF] Zhao | Dynamic texture recognition using local binary patterns with an application to facial expressions[END_REF][START_REF] Ji | Automatic facial expression recognition based on spatiotemporal descriptors[END_REF][START_REF] Hameed Siddiqi | Facial expression recognition using active contour-based face detection, facial movement-based feature extraction, and non-linear feature selection[END_REF][START_REF] Fan | A spatial-temporal framework based on histogram of gradients and optical flow for facial expression recognition in video sequences[END_REF] from the face over time. Although, this first category provides more information that is useful for expression analysis, the dynamic aspect of those models requires considering the neutral expression as reference to predict the exhibited expressions [START_REF] Dornaika | Improving dynamic facial expression recognition with feature subset selection[END_REF]. Therefore, the dependency to the neutral expression limits the use of dynamic models in real world scenes where the expression is unpredictable. The second category of solutions extracts particular features from single images. Such solutions get rid of the dependency to the neutral expression. In Shan et al. [START_REF] Shan | Facial expression recognition based on local binary patterns: A comprehensive study[END_REF] authors investigate the Local Binary pattern (LBP [START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF]) descriptor for low-resolution FER. Some modified versions of LBP were proposed for expression analysis like the Local Tenary Patterns (LTP) [START_REF] Tan | Enhanced local texture feature sets for face recognition under difficult lighting conditions[END_REF] and Compound Local Binary pattern (CLBP) [START_REF] Ahmed | Personindependent facial expression recognition based on compound local binary pattern (clbp)[END_REF]. Uçar et al. [START_REF] Uçar | A new facial expression recognition based on curvelet transform and online sequential extreme learning machine initialized with spherical clustering[END_REF] proposed curvelet transform and online sequential extreme learning machine (OSELM) with radial basis function (RBF) hidden node having an optimal network structure. In the state-of-the-art, a shape descriptor has been used in recent FER studies which is HOG descriptor. It was first proposed by Dalal et al. [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF] for human detection and then studied as part of FER. D. McDuff et al. [START_REF] Mcduff | Affdex sdk: a cross-platform real-time multi-face expression recognition toolkit[END_REF] present a real-time FER toolkit that automatically recognizes the expressions of multiple people simultaneously. Authors extract in that work HOG features from regions of interest defined by a set of landmark points delimiting the face components like the eyebrows, the eyes and mouth. Then, the Support Vector Machine (SVM) classifier is trained on 10000 images, collected from around the world over the Internet [START_REF] Senechal | Facial action unit detection using active learning and an efficient non-linear kernel approximation[END_REF] after filming volunteers while they are watching videos that stimulate their emotions. However, the landmark points have been manually localized in the images. For the real time recognition the toolkit user must be connected to a server in order to exchange and preprocess the frames instantly. The HOG descriptor has been studied and compared to many descriptors in some recent studies. P. Carcagni et al. [START_REF] Carcagnì | Facial expression recognition and histograms of oriented gradients: a comprehensive study[END_REF] presented with details a comprehensive study of the HOG extractor in FER problem. Where, HOG descriptor has been compared to LBP, CLBP, and Spatial Weber's Local descriptor (SWLD). The descriptors' comparison has been realized with an SVM classifier that measures the recognition rate related to each descriptor. The authors showed that HOG features are more relevant than the LBP, CLBP and SWLD features. K. Lekdioui et al. [START_REF] Lekdioui | Facial decomposition for expression recognition using texture/shape descriptors and svm classifier[END_REF] compared the HOG descriptor to other descriptors like the LBP, LTP, and the pair combinations of the descriptors which are HOG-LBP, HOG-LTP, and LBP-LTP. The performance of recognition has been measured using the SVM classifier. The descriptors' evaluation demonstrated that the association of LTP (a texture extractor) and HOG (a shape descriptor) is the appropriate method that describes facial expressions.

Through these studies and comparisons, the authors showed that HOG is one of the useful shape and form descriptor of facial expressions if compared to the texture descriptors. Therfore, our interest goes towards using the HOG descriptor because of its capability to discriminate the shape and the changes of facial expressions, by computing the distribution of edge directions related to the facial traits, on the one hand. On the other hand its capability to get rid of problems related to the variation in shadowing and illumination [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF]. Automatic facial expression recognition is a composite task that involves the expression classification stage in addition to the feature extraction stage. Some studies [START_REF] Carcagnì | Facial expression recognition and histograms of oriented gradients: a comprehensive study[END_REF][START_REF] Lekdioui | Facial decomposition for expression recognition using texture/shape descriptors and svm classifier[END_REF][START_REF] Zhang | Comparison between geometry-based and gabor-wavelets-based facial expression recognition using multi-layer perceptron[END_REF] put emphasis on searching the relevant features to extract from facial expression data. This paper will focus on analysing the hidden structure of facial expression data, after being encoded by HOG features, for the objective of gaining high FER rates and reducing the processing time until reaching a real-time FER that fits the video speed of 24 fps. In fact, it is difficult to decide which classifier is relevant for the facial expression classes before the hidden structure of data classes is analysed and understood. Visualizing the hidden structure of data classes can be helpful to investigate the suitability of different classification methods. However, it is challenging to visualize the data structure when the subjects within the data are represented by more than three features. To solve this problem some dimensionality reduction methods are proposed like the Principal Component Analysis (PCA) [START_REF] Abdi | Principal component analysis[END_REF] and the LDA [START_REF] Peter N Belhumeur | Eigenfaces vs. fisherfaces: Recognition using class specific linear projection[END_REF]. Though the both methods can reduce thousands of features extracted from each face image into only two or three features, there are differences in the manner they transform the dataset when reducing the number of features. As the PCA ignores the class labels during the search of a low-dimension set of axes that summarize data, reducing our data using this reduction method will certainly cause an information lose related to the clusters' separability within the data. The LDA, however, focuses more on finding a low-dimensional set of axes in which the data classes separability increases. It has been widely used. However, in facial expression recognition issue and related to HOG features the LDA has seldom been addressed. M.H. Sidiqi et al. [START_REF] Hameed Siddiqi | Human facial expression recognition using curvelet feature extraction and normalized mutual information feature selection[END_REF] proposed curvelet transform [START_REF] Tang | Facial expression recognition and its application based on curvelet transform and psosvm[END_REF] to extract features from images. The features are then reduced by employing linear discriminant analysis (LDA). Finally, the hidden Markov model (HMM) is used to recognize the expressions. M.H Sidiqi et al. [START_REF] Hameed Siddiqi | Human facial expression recognition using stepwise linear discriminant analysis and hidden conditional random fields[END_REF] used a Stepwise Linear Discriminant Analysis (SWLDA) [START_REF] Dean | Toward enhanced p300 speller performance[END_REF] that focuses in selecting the features from the images. For recognition, the hidden conditional random field model is applied. J. Wang et al. [START_REF] Wang | 3d facial expression recognition based on primitive surface feature distribution[END_REF] investigated the usefulness of 3D fa- Usually, the LDA is used for dimensionality reduction rather than for classification. In fact, the choice of whether to use the LDA only for reduction or for classification too depends on the relationships between the classes when the original data is transformed into a new set of low-dimensional data. In other words, the presence of significant overlaps between the classes, in the scatter diagram of the transformed data, excludes the LDA from the option of being used as a classifier. However, when the data are separable, the LDA makes the class separation easy with less overlap and can be exploited as a supervised classification method if combined with the suitable decision rule allowing the correct classification of the test set. As the main aim is to reach an effective recognition, we are interested in reducing the number of HOG features, related to each subject, with maximizing the between-class variation and minimizing the within-class variation in the reduced data. For the above mentioned purpose the LDA has been used for dimensionality reduction. Furthermore, our preliminary analysis [START_REF] Latifa | A novel tool for automatic exploration of feature extraction and classication methods: A case of study in facial expression recognition[END_REF] for exploring facial expression data, by means of multiple feature descriptors and classifiers, showed that the LDA can be used also as a powerful classifier for HOG data. This conclusion was reached after automatically tuning the parameters of a preselected set of image descriptors and classifiers and comparing the recognition rates of their pair combinations (extractor-classifier). B. Hariharan et al [START_REF] Hariharan | Discriminative decorrelation for clustering and classification[END_REF] showed as well that the use of the whitened HOG descriptor and the LDA classifier instead of the SVM classifier allows building a competitive object detection system that is significantly faster not only at the level of the training time but also the detection time. The work in this study consists then of adapting the processing chain, that was essentially used to explore images and to work on real world videos.

Making FER systems reliable is challenging for practical real-time application. To address this challenge, this paper proposes the following contribution, which is a real-time pipeline based on a processing chain that extracts the shape of the face from images and makes an extreme reduction of the extracted features with preserving the essential needed information for recognition. The pipeline takes as input a video and returns the expression of individuals whom frontal face is correctly detected in the frames. In addition to the previously mentioned contribution, an automatic evaluation technique has been carried out in order to evaluate the pipeline using labelled videos instead of labelled images. Evaluating the pipeline using real world videos allows enhancing the pipeline functioning, Particularly when the frames where the false alarms are detected are registered and added to the training set. This would help the pipeline to make a self-learning while reviewing new images. The remainder of this paper is divided as follows, in Section 2, a new pipeline to automatically recognise the expressions of multiple faces has been proposed. Then, in Section 3 a set of analysis has been carried out for the objective of setting up the LDA classifier and investigating the pipeline suitability in real world application. Finally, Before concluding and proposing some future works in Section 5, a discussion and comparison of the pipeline results to the state-of-the-art has been presented in Section 4.

The methodology overview

The recognition of facial expression requires a pipeline that involves different processing steps. The pipeline in figure 1 has been used in this work: the first step detects the faces within the frame and then the registered one passes through the ROIs extraction step in order to register the most expressive parts of the face. This preliminary operation allows extracting the top region of the face, including the eyes, the eyebrows, and the forehead, the nose region, and the mouth region. Therefore, relevant information like the shape and the wrinkles can be described by applying the HOG extractor. Finally ,for the recognition of the appeared emotions, the vector of features extracted by HOG is classified by the LDA classifier. The pipeline operating steps are detailed in the following subsections.

Face detection and tracking

As people spend time looking straight at their device screen while working, texting, reading or watching videos, the interest goes towards detecting and processing frontal faces in the frames provided by the front-camera of those devices. In this study, static images and video frames are scanned by the fast Viola and Jones face detector [START_REF] Viola | Detecting pedestrians using patterns of motion and appearance[END_REF] in order to find out faces at multiple scales and positions. Detected faces are then tracked over time using an accurate method which is the Kanade-Lucas-Tomasi (KLT [START_REF] Tomasi | Detection and tracking of point features[END_REF]) method. Usually, the tracker is automatically executed once the face detector finds the face in any frame. It extracts the feature points, which are the Harris corners illustrated in figure 2, from the face and tracks them from one frame to another till the face disappears from the scene. The most important role of the tracker is to reduce the image processing time by reducing the image area where the face detector will search for faces. The time is further reduced by skipping some video frames when recalling the face detector. For this, after measuring the pipeline speed while skipping certain number of frames every time the pipeline re-runs, we decided to make the face detector works after every 10 frames. The interesting point here is even if the face detector skips 10 frames, the face tracker operates to mantain the face tracking from one frame to another without skipping any frame, which keeps the pipeline processing works at the real time frame rate. However, The challenge lies in finding new faces. In other words, if a new face appears in the scene, just after calling the face detector, the pipeline will detect and process the new face only after skipping 10 frames. Analysing this point in terms of time, the skip of 10 frames in the case of devices having a shooting video of 30 fps means waiting 0.33 seconds before updating the face detection stage to be able to search new faces, knowing that the tracking is maintained for the already detected ones during the 0.33 seconds. When running the pipeline, the effect of the delay of 0.33 seconds is not visually observable as it is related only to the first face detection and once the face is detected the expression will be analysed from one frame to another till the face disappears.

ROIs extraction

Once the face tracking is ensured in the video, the ROIs containing transient features like the wrinkles, the forehead lines, the frown lines, and the laugher lines and permanent features like the forehead with the eyes, the nose, and the mouth are extracted from each face window. Thus, the irrelevant face parts are eliminated while facial regions containing transient and permanent features that are prone to changes are selected. As long as the face is detected in its frontal position, the ROIs can be localized easily and rapidly by using geometric method [START_REF] Jeng | Facial feature detection using geometrical face model: an efficient approach[END_REF] rather than using detectors for facial components like the eyes, the nose, and the mouth detectors. The geometric method is based on the fact that human faces have similar geometric configuration and there is a proportionality between the face components dimension, such as the proportionality between the face length and the position of the eyes or the mouth. Thus, the three ROIs have been localized using this method as follows: -The face window having L x L dimensions is divided into 6 x 6 sub-regions as shown in figure 3. -For each ROI, specific sub-regions are selected. For example, the region of the eyes is bounded by subregions (2-5) horizontally and (2-3) vertically. Emotions stimulate initially the face muscles. Therefore, the movement of the muscles appears through two types of face shape changes involved when a new emotion is occurring. The first change is related to the movements of permanent features. The second change concerns the transient features that appear only when Fig. 3 The face partitioning and ROIs extraction the emotion is occurring. The best way to describe facial expression is to extract the shape of the face components. By shape, we mean the relative positional information that can be formed while computing particular information like the local distribution of the edge directions [START_REF] Belongie | Shape context: A new descriptor for shape matching and object recognition[END_REF]. For this reason, the powerful shape descriptor which is HOG extractor has been used to construct the proposed pipeline.

HOG descriptor

Allows extracting the local shape from the regions of interest as illustrated in figure 4. It calculates the distribution of local edge directions in small cells of size 8 × 8 pixels. At each pixel, the horizontal g x and the vertical g y gradients are calculated using 1-D centred derivative masks [-1 0 1] and [-1 0 1] T . A gradient orientation θ (x,y) = arctan gx gy is used to vote into eight corresponding orientation histogram bins equally spaced between 0 -360 o . Then, the votes are weighed by a gradient magnitude M (x,y) = g 2

x + g 2 y and accumulated over pixels of each cell to construct one resulting histogram h i . Histograms calculated over blocks are normalized with an overlap of 50% to provide better invariance to the change of illumination. Finally, all histograms h i are concatenated in order to construct the final vector

H = [h 1 , h 2 . . . , h n ].
Once the HOG descriptor extracts the vector of features, the pipeline final step classifies that vector. However, before going through the classification step, a training phase is required in order to construct the LDA classifier.

LDA training

The training and the classification can only be achievable in two distinct phases. The scheme of the training phase is illustrated in figure 5. The scheme consists of the same above mentioned pipeline preprocessing steps in addition to the LDA reduction step. However, no face tracking was needed for the face detection step since labelled images were used instead of video sequences. The labelled images were collected from two different databases in order to diversify the subjects within the training set. The ChonKanade Image Database (CK+) [START_REF] Lucey | The extended cohn-kanade dataset (ck+): A complete dataset for action unit and emotion-specified expression[END_REF] is specifically acquired for facial expression analysis. This choice has been made according to the following reasons: the diversity of samples within the same classes in terms of gender, as well as ethnicity and age; the necessity to compare the performance of our proposed pipeline to the ones of the state-of-the-art. The CK+ contains image sequences of people with seven facial

data m,n =      h 1,1 h 1,2 • • • h 1,n h 2,1 h 2,2 • • • h 2,n . . . . . . . . . . . . h m,1 h m,2 • • • h m,n      (1)
Where data m,n is the resulting labelled data after extracting HOG features from the images, the rows refer the vector related to each image within the database, n is the number of features extracted from the ROIs, and m is the number of images.

LDA reduction [START_REF] Peter N Belhumeur | Eigenfaces vs. fisherfaces: Recognition using class specific linear projection[END_REF]: when the resulting matrix data m,n is given, the number of variables n can be reduced.

The idea is to transform the matrix data m,n having The calculation of S W and S b are the following [START_REF] Marsland | Machine learning: an algorithmic perspective[END_REF] :

S w = classesC M iεC P C .(h i -µ C )(h i -µ C ) T (2)
S b = classesC (µ C -µ)(µ C -µ) T (3)
Where The produced diagram shows that the classes can cluster separately. However, there's an overlap of samples at the boundaries of the seven classes like the class of happy faces with scared faces. Therefore, we must construct our classifier by separating the data clusters with the adequate decision rule. In this context, the pipeline has been exploited to test multiple decision rules over the produced clusters. For that different classification techniques have been associated with the LDA in which the LDA decision rule has been replaced by the decision rule related to other classifiers like the k-NN, the Naïve Bayes, the SVM, and the Binary Tree. Then, the recognition rate has been measured for each combination test, so that the LDA decision rule can be chosen according to the combination having the maximum recognition rate. After finishing the training and fixing the decision rules of the LDA classifier, the last pipeline step uses this classifier to recognize expressions in real world videos.

µ C = 1 M M i=1 h (i,j) and µ = 1

Classification

At this step when a new face image with the vector of features X is projected in a particular location in the Fig. 6 The scatter diagram of the CK+ database after extracting HOG features and reducing dimension by the LDA scatter diagram, the LDA classifier relies on a multivariate Gaussian to estimate the probability that the facial expression X belongs to a particular class (called also the posterior probability P (C = c i /X = x)), where c i refers the ith class and i {1, 2, 3, 4, 5, 6, 7} refers the class number. The posterior probability is computed using the Bayes rule witch is equal to the probability of X knowing the class P (X = x/C = c i ) multiplied by the probability P (C = c i ) of the class c i and divided by the sum of the probabilities of the facial expression X under all the possible classes. So, the posterior probability is defined as:

P (C = c i /X = x) = P (X = x/C = c i )P (C = c i ) 7 j=1 P (X = x/C = c j )P (C = c j ) (5)
As we assume that our data come from a multivariate distribution, the P (X = x/C = c i ) can be computed as follows:

P (X = x/C = c i ) = exp(-1 2 (x -µ i ) T Co -1 i (x -µ i )) 2π | Co i | (6)
Where, µ i is the vector mean and Co i is the covariance of the classe i.

The requested facial expression X is classified as a class c i according to the class that gives the maximum posterior probability, given by:

P max = argmax ci {P (C = c i /X = x)} (7)
After the LDA recognizes the face expression, the pipeline outputs the location of the faces (if there are so many within the frame) as bounding boxes and exhibits the face expression on the top of those bounding boxes.

In addition to that, the pipeline records the result of recognition and the face location for the objective of evaluating the pipeline. On the whole, this section describes the pipeline steps and details the training phase in which further analysis have been realized to set up the LDA classifier.

The pipeline evaluation

In this section many experimentations and tests have been carried out to evaluate our proposed pipeline functioning when using real world videos. Before the evaluation, first, the pipeline has been exploited to conduct deep analysis on the hidden structure of data and set up the decision boundary of the LDA classifier. 1. Looking at the SVM decision boundary, it is obvious that the top right region of the scatter diagram is identified as a part of the sad class while there is no point that belongs to this class in that region. Thus, in this case, if a requested face image is projected within this region the expression of that face will be misclassified. The analysis of table 1 demonstrates that the LDA with the multivariate Gaussian seems to be suitable for the classification of the seven facial expression classes.

That's because it can make good decision boundaries by separating, linearly, the seven classes (as shown in figure 7) and can classify, accurately, requested images with a maximal average recognition rate of 96.44%.

Once the LDA decision boundary is fixed, the pipeline uses the LDA classifier to recognize the face expressions in videos collected from the MMI database and videos As it was mentioned in the subsection 2.4 that the training is made using the Yale Face and the CK+ databases, it is worthwile to make a closer look regarding the generalization ability of the LDA classifier. According to the obtained results the LDA is general enough to be representative of those databases and some real world situations as we can see in the Subsection 3.2. The CK+ database is widely used for algorithm development and evaluation of facial expression recognition systems [START_REF] Carcagnì | Facial expression recognition and histograms of oriented gradients: a comprehensive study[END_REF][START_REF] Lekdioui | Facial decomposition for expression recognition using texture/shape descriptors and svm classifier[END_REF][START_REF] Sl | Automatic facial expression recognition using features of salient facial patches[END_REF][START_REF] Bartlett | Computer expression recognition toolbox[END_REF] The reason of this is in addition to the diversity of samples within this data, it includes also video sequences that vary in duration (i.e. 10 to 60 frames) and incorporate the onset (which is the neutral frame) to peak formation of the facial expressions. On the other hand, the benefit behind using the Yale Face together with the CK+ database is just to increment the diversity of samples and to add images with partially occluded faces to the experimental dataset. However, we still can't generalize in real world data in complicated situations like occluded and rotated faces, because the training data contain only few samples considering those situations. To make the LDA classifier general enough to be representative of the real world it should be trained using an important quantity of real world data. To make this possible, the pipeline can register the results of recognition including the face image and the expression associated with, so that the false recognition ones can be selected and used to retrain the LDA classifier.

Fig. 8 The automatic evaluation technique of our pipeline using real labelled videos

The evaluation of the proposed pipeline

The challenge of evaluating the pipeline in real world videos leads us to build the scheme presented in figure 8. The scheme consists of the same pipeline processing elements in which two additional operating blocks have been added. The first block aims at evaluating the face detector and the second one aims at evaluating the LDA classifier. For the automatic evaluation, first the scheme takes as inputs a set of labelled videos that contains frames with faces displaying different transitions between facial expressions. Those videos were collected from the MMI database [START_REF] Pantic | Web-based database for facial expression analysis[END_REF]. The reason for choosing this data is diversity of subjects and the existence of partially occluded faces by hair and glasses in addition to the various lighting conditions. The face detector scans the video frames in order to find out the faces. At the face detector output, the block 1 has been added to ensure that the face window found by the face detector is a true positive. Block 1: To decide whether the outputs of the face detector are faces, the block 1 compares the face positions that have been manually extracted from the video frames (the reference vector) to those given by the face detector when the same video is being entered then processed by the pipeline. As the face of individuals participating in the MMI videos are in the same position over the frames, it was easy for us to set the face position for all the videos and thereby to construct the reference vectors for the block 1. Therefore, if the face positions that are given by the face detector don't match with the face position in the reference vector then the recognition of the facial expression fails. Then, the output of the Block 1 could be a false positive F pi (the face detector can't find the face while it exists in the image) or a False Negative F ni (the face detector estimates are found the face while it doesn't exist). In this case, the outputs of this block F pi and F ni will be stored in the false recognition metric. Otherwise, if the face position given by the face detector matches the one in the reference vector, the detected face pass throughout the remaining pipeline steps till arriving at the second block that evaluate the recognition results outputted by the LDA classifier. Block 2: allows comparing the pipeline outputs with the reference vectors stored in this block.To do this, the reference vector related to each video sequence used in the experimentation has been manually constructed as follows. As almost videos in the within the MMI database begin with a neutral expression that changes over time to achieve a high-intensity expression and then returns to the neutral expression, the reference vector can be constructed by assigning the value 1 to the frames that contain the high-intensity expression and 0 otherwise, as illustrated in figure 9. To avoid assigning manually the label for each frame in video sequences that lasts on an average of 6 seconds with a speed of 30 fps, the solution is to just localize the transition frames as follows: T 1 where the expression changes from the neutral state to the high-intensity expression and T 2 where the expression changes to go back to the neutral state. So, The frames between the transitions T 1 and T 2 will have the same label that is equal to 1 and the ones at the left of T 1 and the right of T 2 are labelled with the value 0. However, for few videos the expression begins and ends with a minimal intensity like the sequence presented in figure 10. In this case, all the frame are assigned into 1. For us a sequence is considered correctly recognized when the frames containing the neutral expression are correctly recognized and the frames containing the basic expression of the video sequence are also correctly recognized by the classifier. Finally, for each video, the reference vector is constructed and stored in the block 2, so that the block can compare the pipeline outputs y i with the elements of the reference vector y i as follows:

S i = y i if y i = y i y i if y i = y i (8)
Where y i is the face expression in the reference vector at the i-th frame, y i is the face expression classified by the LDA at the i-th frame The block can decide whether the expression is correct or not using the following condition. If S i = y i then the resulting S i is a positive outputs, which means the facial expression is correctly classified (see scheme 8). Otherwise, if S i = y i then S i is a negative outputs, which means the facial expression is misclassified. On the whole, the positive and negative outputs are used to construct the confusion matrix and compute the f-measure [START_REF] Sokolova | A systematic analysis of performance measures for classification tasks[END_REF] which gives a precise calculation about the average recognition rate. Table 2 reports the resulting confusion matrix computed after the automatic pipeline evaluation. Knowing that the LDA has been trained by the CK+ and Yale Face databases, the tests using MMI videos showed that the pipeline is accurate in recognizing the expression of joy, surprise, and anger. Those three expressions have recognition rates, superior than 90%, which is normal for the joy and surprise expressions as they are the easiest expressions to discriminant [START_REF] Michel | Real time facial expression recognition in video using support vector machines[END_REF][START_REF] Pardàs | Facial animation parameters extraction and expression recognition using hidden markov models[END_REF]. However the recognition rates of fear, disgust,and sadness are limited between 90% and 80%. There can be many reasons why certain expressions are not correctly classified. the first reason is the low-intensity of the expressions when it changes from expression to another which makes the recognition particularly difficult at the frames of transitions. Another reason is the similarities in the appearance of some facial traits, like strict lips and when the eyebrows are down, makes confusions between the expressions of disgust, anger, and sadness. The scheme has been adapted to continuously visualize the change of the expression over time. Figures 11(a) and 10 illustrate the FER in two different videos where the visualisation of the expression transitions has been automatically returned after executing the proposed pipeline. In figure 11 (a), a correct prediction is obtained over the frames containing the basic expression, whereas the prediction is unstable due to the pipeline uncertainty during the permutation of the face expression. The proposed pipeline has been designed to recognize the expression of multiple faces in the case where several people appear in the video scene. An example of expression recognition over four frames extracted from a real video has been presented in Figure 12, where three individuals are filmed using the front camera of a cell phone. After testing the pipeline on some real world videos, it should be noted that the proposed pipeline achieves an average processing time of 0.018 second for one face. This means that pipeline is applicable at least for devices having a shooting video of 24 fps. However, the more the faces are multiplied in the scene the more the time is consumed by the pipeline. Indeed, the processing time varies according to many variables which are the face size and the scene background complexity.

Meaning that in the case of scenes containing faces having large dimension or containing many moving objects in the background, the pipeline consumes a lot of time in analysing the scene frames.

Discussion

For fair comparison with the state-of-the-art, many protocols must be respected which limits some choices regarding the database, the manner we split it, and the metrics we use during the pipeline evaluation. For this reason, the pipeline results (reported in section 3) have been compared with the results of the existing works that use in common: the CK+ database, the 10-fold or LOO cross-validation techniques, and the recognition Fig. 12 The resulting recognition after applying the proposed pipeline on a real world video filmed by the front camera of a cell phone metric F-measure. Based on those protocols, the comparison with the state-of-the-art have been done using two studies [START_REF] Lekdioui | Facial decomposition for expression recognition using texture/shape descriptors and svm classifier[END_REF][START_REF] Carcagnì | Facial expression recognition and histograms of oriented gradients: a comprehensive study[END_REF] in which the HOG descriptor has been used to extract facial expression features. Table 3 reports the comparison result demonstrating that our proposed pipeline provides a high average recognition rate of 96.44 and processing time of 18 ms. The pipeline computational cost, in terms of time, has been measured considering an average of 100 frames using Intel core i7 processors, CPU @ 2.20GHz 2.2GHz, and RAM of 16.0 Go. Lekdioui et al. [START_REF] Lekdioui | Facial decomposition for expression recognition using texture/shape descriptors and svm classifier[END_REF] propose a FER system that employs an automatic technique for finding the proper cropping of facial components. Their experimental analysis showed that the performance of HOG descriptor associated with a texture descriptor, like the LTP, and an SVM classifier provide a recognition rate of 96.06%. However, the recognition time demonstrated that this FER system have an expensive computational cost in terms of time and thereby can be applied only to static images. Carcagni et al. [START_REF] Carcagnì | Facial expression recognition and histograms of oriented gradients: a comprehensive study[END_REF] carried out a comprehensive study of the HOG descriptor as a part of the FER issue, showing that a proper tuning of HOG parameters can lead to a powerful representation of facial expressions and thus an effective recognition with an average of 94.1%. The simplicity of the model, which is based on HOG descriptor and the SVM classifier, have allowed to apply the pipeline in real-world operating conditions with a preprocessing time of 43.38 milliseconds. However, with the analysis carried out in subsection 3.1 and the pipeline evaluation in subsection 3.2 we emphasize that a potential reduction of HOG features by the LDA can speed up the FER with maintaining high recognition rates.

Conclusion

The need for building practical facial expression recognition system, lead us to propose a pipeline after making several data analysis. The pipeline first extracts three ROIs isolating some face areas like the forehead with the eyes, the nose, and the mouth. The ROIs are then processed in order to extract the shape information using a HOG feature extractor. An extreme reduction of HOG feature has been realized by the LDA. For classification, the following classification techniques including k-NN, the Naïve Bayes, the SVM, the Binary Tree, and the pair combination of these classifiers with the LDA have been tested. An automatic evaluation technique has been proposed to compute the pipeline performance of recognition when many expression transitions are present in the videos. Already trained by CK+ video sequences and Yale Face database, the pipeline evaluation using MMI video sequences showed that the proposed pipeline can be applied in real-time (i.e. at the frame rate of 24 fps) with an average processing time of 0.18 millisecond. However, the implementation of the pipeline has highlighted several points that should be improved, namely: the adaptation of the pipeline in order to recognize further expressions such as recognizing the driver fatigue, to survey the state of patients and to recognize the micro-expressions. It is necessary to run this pipeline, constantly, in real scenes either if the facial expression is a standard expression or not. Parallelization is also needed in order to ensure fast facial expression recognition in the presence of dozens of faces in the video scene.

Fig. 1

 1 Fig.1The overall view of the proposed pipeline

Fig. 2

 2 Fig.2The feature points extracted from the face and used for the tracking

Fig. 4

 4 Fig. 4 HOG feature extraction from the ROIs

Fig. 5

 5 Fig. 5 The training phase of the LDA classifier

  i,j) are the mean vector in each class and µ is the overall mean, respectively. M and C are the size of samples in each class and the class, respectively. Making S -1 W S b as large as possible is the key to reduce the thousands of data variables into two discriminant variables that contain almost all the classification information embedded in the original images. Those variables are computed by resolving the generalised eigenvectors of the matrix S -1 W S b in the equation 4 and selecting the eigenvectors υ that have the maximal eigenvalues λ: Aυ = λυ (4) Where A = S -1 W S b is the between-class scatter divided by the within-class scatter assuming that S -1 W exists. The two selected eigenvectors, having maximal eigenvalues, constitute the columns of the two variables of the reduced data. Those two variables are then plotted in the scatter diagram of figure 6. This figure represents the spread of 248 face expressions clustered into seven classes after being transformed and reduced by the LDA, where the axes LDA 1 represents the first eigenvector and the axes LDA 2 represents the second eigenvector.

Fig. 7

 7 Fig. 7 The scatter diagram showing the linear decision boundaries of the LDA classifier

Fig. 9 A

 9 Fig.9A simplified scheme showing the construction process of the reference vector

Fig. 10

 10 Fig. 10 The recognition results of the video "S036-005.avi" taken from the MMI database

Fig. 11

 11 Fig. 11 The graphical representation of the facial expression variation over time in a video. (a) The recognition results of the video "S001-100.avi" taken from the MMI database. (b) The reference vector of the same video "S001-100.avi".

  3.1 The LDA evaluationAn effective use of the HOG descriptor in FER requires the use of robust and fast classification method. However, choosing the appropriate classification tools is difficult if one does not have a prior knowledge about the structure of the data classes. By knowing the hidden structure we can understand how the data classes are organized and clustered. The structure means the geometric relationships among the data classes in twodimensional space. An example of a structure is linear and non-linear relationships among the vectors of the classes. Therefore, understanding the geometric relationships among the data classes will help us figure out the suitable decision boundary that separate the classes. We need to choose the appropriate decision boundary for the LDA classifier. For this purpose, table 1 summarizes the analysis after testing different decision boundaries to separate the clusters within the LDA classifier during the training. The generalization ability of the hybrid classifiers, presented in the second column of the table, reveals that associating the LDA with the Multivariate Gaussian or Naïve Bayes makes the recognition better with an average recognition rate of 96.44% and 95%, respectively. However, combining the LDA with the SVM, the k-NN, and the Binary Tree doesn't enhance the recognition of the facial expressions at all.

	The problem with the Quadratic Discriminant Analy-
	sis (QDA) and the other classifiers is that even if the
	decision boundaries between the classes seem to be cor-
	rect in the scatter diagrams at the third column of the
	table 1, those methods fail in classifying new requested
	images and that's because they make wrong decision
	boundaries in regions where no data points are avail-
	able. This problem is noticeable in the case of the scat-
	ter diagram that represents the SVM decision bound-
	ary, in table

Table 1

 1 The recognition results and the visualization of the hybrid classifiers

	Methods	Recognition rate	Decision boudaries visu-alisation
	The LDA
	with
	the Mul-	0.96
	tivariate
	Gaussian
	The LDA
	with	0.43
	the SVM
	The
	Quadratic
	Discrim-inant	0.80
	Analysis
	(QDA)
	the LDA
	with	0.77
	the k-NN
	The LDA
	with the Naive	0.95
	Bayes
	The LDA
	with the Binary	0.43
	Tree

filmed in real word in the presence of more than one

Table 2

 2 The confusion matrix after the pipline evaluation over video sequences. su, fe, di, an, sa, and ne refer to joy, surprise, fear, disgust, anger, sadness, and neutrality respectively

		jo	su	fe	di	an	sa	ne
	Jo	95.7 0	3.1	0	0	0	1.2
	Su	1.3 96	1.9	0	0	0	0.8
	fe	5.12 7.01	85.09 0	0.08 0.1	2.6
	di	0	0	0	89.92 3.58 3.1	3.4
	an	0	0	0.84 1.6	91.1 4.02 2.44
	sa	0	0	1.3	2.9	8	83.49 4.31
	ne	0.3 0.3	5.5	0.9	8.94 10.23 73.83
	with jo,

Table 3

 3 Performance comparison of our approach vs recent State-of-the-Art approaches (with using CK+ data, 10-fold, and F-measure) Av: is the average recognition rate (the F-measure) time: is the whole preprocessing time of the systems ms: refers the millisecond unit

		Carcagni et	Lekdioui et	Proposed
		al. [19]	al. [20]	pipeline
	Av	94.1	96.06	96.44
	time	43.38 ms	236.18 ms	18 ms