Bertrand Candelon
email: b.candelon@maastrichtuniversity.nl

Elena-Ivona Dumitrescu
email: elena.dumitrescu@univ-orleans.fr

Christophe Hurlin
email: christophe.hurlin@univ-orleans.fr

Franz C Palm
email: f.palm@maastrichtuniversity.nl

Multivariate Dynamic Probit Models: An Application to Financial Crises Mutation

Keywords: C35, F37 Non-linear VAR, Multivariate dynamic probit models, Exact maximum likelihood, Impulse-response function, Financial crises

In this paper we propose a multivariate dynamic probit model. Our model can be considered as a non-linear VAR model for the latent variables associated with correlated binary time-series data. To estimate it, we implement an exact maximum-likelihood approach, hence providing a solution to the problem generally encountered in the formulation of multivariate probit models. Our framework allows us to apprehend dynamics and causality in several ways. Furthermore, we propose an impulse-response analysis for such models. An empirical application on three financial crises is finally proposed.

Introduction

Since the pioneer paper of [START_REF] Sims | Macroeconomics and Reality[END_REF], Vector Autoregressive (VAR) models have been the main tool used in the analysis of continuous macroeconomic time-series. This has paved the way for numerous papers proposing improvements to the standard linear VAR model. A research direction focuses on cointegrated time-series (e.g. Engle and Granger, 1987 ;Johansen, 1988Johansen, ,1991)), leading to Vector Error Correction Models, while another line of research concentrates on adapting VAR specifications to the case of non-linear continuous time-series [START_REF] Gallant | Nonlinear Dynamic Structures[END_REF]Potter, 2000, among others).

Apart from continuous indicators, discrete-choice variables are also ubiquitous in macroeconomics and their forecast has become essential for different economic parties. There are two main ways to proceed with their analysis. First, the standard methodology used to forecast the probability of occurrence of such events consists in simple logit / probit models. They are used, for example, by Berg and Patillo, 1999;[START_REF] Estrella | The Yield Curve as a Leading Indicator: Some Practical Issues[END_REF][START_REF] Fuertes | Optimal Design of Early Warning Systems for Sovereign Debt Crises[END_REF] to forecast recessions as well as financial (currency, banking, sovereign debt, stock market,...) crises. A second, more recent approach, consists in a dynamic univariate framework based on an exact maximum-likelihood estimation method (see Kauppi and Saikonnen, 2008;[START_REF] Candelon | Currency Crises Early Warning Systems: Why They Should Be Dynamic[END_REF] or [START_REF] Nyberg | Dynamic Probit Models and Financial Variables in Recession Forecasting[END_REF][START_REF] Nyberg | Forecasting the Direction of the US Stock Market with Dynamic Binary Probit Models[END_REF], inter alii ). These dynamic models have been shown to be parsimonious, as the dynamic forecasts obtained outperform the forecasts issued from standard static models. In a similar vein, [START_REF] Dueker | Dynamic forecasts of qualitative variables : A Qual VAR model of U.S. recessions[END_REF] estimates a 'mixed' Qual-VAR model including as dependent variables the latent variable lying behind the binary business cycle indicator along with standard continuous indicators by using simulation based methods. Modelling the dynamics of qualitative variables hence appears as important as that of continuous variables (for the improvement of Early Warning Systems, for example). Nevertheless, no multivariate extension of these dynamic binary models exists as far as we know. This is intriguing, since it would definitely improve forecasts in view of the existing correlation among the different macroeconomic binary indicators.

To our knowledge, the only multivariate binary models that exist are not dynamic. Indeed, although numerous macroeconomic studies involve discrete-choice dependent variables, this particular type of non-linearity has received little attention in the literature on multivariate time-series. The main reason behind this relies in the difficulty to formulate a multivariate binary model, in particular to evaluate the underlying likelihood function. In this context, Carey, Zeger and Diggle (1993) and [START_REF] Glonek | Multivariate Logistic Models[END_REF] have proposed generalizations of the binary logistic model to multivariate outcomes by selecting a particular parameterization for the correlations. Besides, [START_REF] Ashford | Multivariate Probit Analysis[END_REF] and [START_REF] Amemiya | Bivariate Probit Analysis: Minimum Chi-square Methods[END_REF] have focused on generalizing the binary probit model. Some attempts have subsequently been made to solve the computational difficulty of evaluating multivariate normal integrals. Note the works of Chib and Greenberg (1998) who developed a simulation-based Bayesian and non-Bayesian approach and those of [START_REF] Song | A Multivariate Probit Latent Variable Model for Analyzing Dichotomous Responses[END_REF], who rely on a Monte Carlo EM algorithm to evaluate the likelihood function for a multivariate probit model. More recently, [START_REF] Huguenin | Estimation of Multivariate Probit Models by Exact Maximum Likelihood[END_REF] have shown that a multivariate probit model cannot be accurately estimated using simulation methods, as generally done so far in the literature. Its estimation requires hence to derive an exact maximum-likelihood function.

The objective of this paper is hence to extend the estimation methodology proposed by [START_REF] Huguenin | Estimation of Multivariate Probit Models by Exact Maximum Likelihood[END_REF] for non-dynamic multivariate probit models to the case of dynamic multivariate probit models. We hence introduce a multivariate dynamic probit model which relies on an exact maximum-likelihood estimation approach to produce dynamic forecasts of binary time-series correlated variables. Our model takes the form of a non-linear VAR for the latent variables associated with the binary indicators under analysis. It is easy to implement and provides a solution to the problem encountered in the estimation of multivariate probit models. For this, in a first step we extend the decomposition of high-order integrals into integrals of lower order proposed by [START_REF] Huguenin | Estimation of Multivariate Probit Models by Exact Maximum Likelihood[END_REF] to the case of our multivariate dynamic model. In the second step, we evaluate the lower-order finite-range multiple integrals by using quadrature-rules over bounded intervals.

Our framework allows us to apprehend dynamics and causality in several ways. First, it can be included as a set of lagged binary variables. Notice the existence of threshold effects in this case, as the event under analysis occurs only if the latent variable goes beyond a certain threshold. Second, it can be introduced via the past latent variables associated with the binary indicators. This comes down to an autoregressive (AR) model, where the lagged latent variable summarizes all the past information of the system. Finally, both types of dynamics can be simultaneously considered. We thus generalize the univariate dynamic probit model developed by [START_REF] Kauppi | Predicting U.S. Recessions with Dynamic Binary Response Models[END_REF] to a multivariate level and derive its exact likelihood, hence obtaining consistent and efficient parameter estimates. Furthermore, we are the first ones to propose an impulse-response function (IRF) analysis for a multivariate probit model.

In an empirical application on three types of financial crises, namely currency, banking and sovereign debt crises, we investigate the potential spill-over from one crisis to another within a number of emerging countries. It appears that in the bivariate case mutations of a banking crisis into a currency crisis (and vice-versa) have been quite common, confirming hence other results in the financial crises literature (e.g. Glick and Hutchinson, 1999). More importantly, for the two countries (Ecuador and South Africa) which suffered from the 3 types of crises, the trivariate model turns out to be more parsimonious, thus supporting its implementation anytime when it is feasible. The conditional probability and IRF analyses confirm these findings.

The rest of the paper is organized as follows. Section 2 presents the multivariate dynamic probit model. In section 3 we describe the Exact Maximum Likelihood method as well as some numerical procedures to estimate the multivariate dynamic probit model. Section 4 describes the steps to follow to construct the IRFs and the associated bootstrapped confidence intervals, while in section 5 the multivariate dynamic probit model is estimated for 17 emerging countries in its bivariate or trivariate form.

A Multivariate Dynamic Probit Model

In this section we describe the multivariate dynamic probit model and insist on its similarities with traditional VAR models. Let us denote by y m,t , m = {1, 2, ..., M }, M binary variables taking the value one if the event under analysis occurs at time t and zero otherwise. Let y * m,t be the normal latent continuous variable associated with y m,t , and define

F t-1 = σ[(y s , x v ) |s < t, v ≤ t]
as the information set available at time t.

The corresponding multivariate probit model takes the well known form:

y * m,t = π m,t + m,t , ∀m = {1, 2, ..., M } y m,t =    1, if y * m,t > 0 0, otherwise , (1)
where π m,t denotes the index, and the innovations t = { 1,t , ..., m,t } verify

t |F t-1 ∼ IIN (0, Ω), (2)
such that V ( ) = I T ⊗ Ω, where I T is the identity matrix of order T . Ω stands for the covariance matrix of t , with Ω = (σ m,m ) and σ m,m = ρ m,m σ m σ m , where ρ m,m represents the correlation coefficient between the m th and m th binary processes and σ m and σ m are the associated standard deviations. The m events under analysis can hence be related at time t through the associated innovations, but there is no dependence in time among these innovations (V is block-diagonal). In this model, the probability of occurrence of the m th event is equal to

p m,t = Pr(-m,t ≤ π m,t |F t-1 ) = Φ(π m,t ), ∀m = {1, 2, ..., M }. (3)
The objective of this paper is to propose a dynamic multivariate modelling of these M qualitative variables. We specify the dynamics of each stochastic process through its associated index π m,t . Traditionally, the index only depends on exogenous explanatory variables x t = {x 1,t , ..., x K,t }, where K is the number of exogenous variables in the model. But in a dynamic model, it can also depend on the past information on the dependent variable. For-mally, for a given event m, the dynamics can be introduced in two ways: either through the lagged binary variables y m,t , or through the lagged latent variable y * m,t . The first equation of our multivariate dynamic probit model in ( 1) is hence given by:

y * m,t = α m + β m x t + M m =1 ∆ m,m y m ,t-s + M m =1 Γ m,m y * m ,t-s + ε m,t , (4)
where 4), with θ = (θ 1 ; θ 2 ; ...; θ M ) , that will be used in the estimation of the model (see section 3).

α m ∈ R is the intercept, β m = {β 1,m , ..., β K,m }, β m ∈ R K is
It is clear that our model has the usual VAR-X structure, as it assumes a linear relation between the latent variables y * m,t and their past. Its dynamics is then enriched by the nonlinear relation between the latent variables y * m,t and the observed binary ones y m,t , which in turn depend upon y * m,t-s (see the second eq. in ( 1)). In other words, our dynamic probit model differs from a standard VAR in two ways. First, it introduces both a linear and a non-linear dynamics. Indeed, the dynamics of the m th process / event can be modeled by considering that the latent variable y * m,t depends either on its lagged value via the Γ m,m coefficient, or on the past regime (0/1) through ∆ m,m .

Second, the analysis of the interdependence, i.e. causality between the M qualitative variables is more sensitive, as it passes through several channels.

1. Unobserved common factors can be taken into account through the contemporaneous dependence of the innovation terms (E( m m ) = σ mm = 0 for m = m ).

2. For an event m, the unobservable latent variable y * m,t depends on past values of other processes y * m ,t-s (where m = m ), themselves unobservable, which can be interpreted as a causality phenomenon.

3. The latent variable may depend on past realizations of the other events, i.e. y m ,t-s = 1.

Formally, y * m,t depends on past values of the observable variable y m ,t-s , s > 0, where m = m . 4. It is possible to combine the two previous cases, assuming that y * m,t , depends on both the latent variable y * m ,t-s , and past values of the observable variable, y m ,t-s , for other binary processes (m = m ).

From this perspective, our model in eq. 4 is a multivariate extension of the univariate dynamic probit model recently proposed by [START_REF] Kauppi | Predicting U.S. Recessions with Dynamic Binary Response Models[END_REF].

Another novelty introduced by this non-linear VAR framework relies in the fact that our new specification enables us to compute not only marginal but also joint and conditional probabilities. The traditional marginal probabilities are associated with each binary response, Pr(y m = 1|y * m ) = Φ(y * m ), and rely on univariate discrete-choice models. In contrast, joint and conditional probabilities, i.e. Pr(y 1 = 1, y 2 = 1, ..., y M = 1|y * ) = Φ M (y * ), and Pr(y m = 1|y * m ) = Φ M (y * )/Φ M -1 (y * m ), for m, m ∈ {1, 2, ..., M }, where Φ and Φ M represent the univariate and M-variate normal cumulative distribution functions respectively can also be obtained here.

The first-order dynamics, i.e. s = 1, is the most common situation in empirical applications with binary-dependent indicators that stand for regime-switches. It is also the case of the empirical application we present in section 5. In vector notation, the process in (4) hence becomes:

y * t = α + Bx t + ∆y t-1 + Γy * t-1 + t , (5)
where α = (α 1 , ..., α m ) , B is a M ×K matrix, and ∆ and Γ are M ×M matrices of parameters. Note that the matrices ∆ and Γ summarize useful information about the dynamics of the binary processes, in particular about their persistence and causality.

One the one hand, as previously discussed, the diagonal terms of Γ specify the persistence of each process. These parameters correspond to a first order autoregressive representation of each latent variable. An increase in the latent variable during a certain period is always transmitted to the next period, hence always increasing the probability of realization of the event (observing a value of 1). The closer these parameters are to 1, the more persistent the processes are. Notice that the diagonal elements of this matrix will be constrained to be strictly inferior to 1 to exclude the case where the latent variable y * m,t follows a random walk, which is not an interesting case in this context.

At the same time, the diagonal terms of ∆ also deliver information about persistence but that is somewhat different from that inferred from Γ. Indeed, they indicate to what extent the probability of occurrence of an event (defined by the m th binary indicator) depends on the regime prevailing the period before, i.e. the occurrence of the same event in the previous period. In this situation we observe the existence of threshold effects, as the regime defined by a value of one for the binary process lasts more than one spell only if the latent variable soars sufficiently to exceed a threshold which initiates this regime in the previous period.

Altogether, we can distinguish between a linear persistence of the phenomenon, captured through the diagonal terms of Γ, and a non-linear, threshold-based one, apprehended by the diagonal terms of ∆.

On the other hand, causality is taken into account in the off-diagonal elements of the two matrices Γ and ∆. These Granger-causal effects between the M binary stochastic processes play a key role in the analysis of numerous economic events (e.g. financial crises, economic cycles). As in the analysis of the persistence of a binary process, both a linear and a nonlinear, threshold-effect transmission can be identified. A significant off-diagonal Γ element shows that no sooner the latent variable for the m binary indicator soars, than the one for another process, say m , rises. By contrast, a ∆ term reveals the presence of causality only if the corresponding latent variable is high enough to impact the occurrence of another binary event.

Exact Maximum Likelihood Estimation

The exact maximum likelihood estimator for the multivariate dynamic probit model cannot be obtained as a simple extension from the univariate model. For this reason, the simulated maximum likelihood method is generally considered. Nevertheless, Holly, [START_REF] Huguenin | Estimation of Multivariate Probit Models by Exact Maximum Likelihood[END_REF] prove that it leads to a bias in the estimation of the correlation coefficients as well as in their standard deviations. Therefore, they advocate the exact maximum likelihood estimation. Since the correlations between the different binary variables, i.e. the contemporaneous transmission channels from one process to another one, is essential for disentangling the mechanism linking several binary indicators, asymptotic unbiased estimation of the correlations is of importance and it calls for an explicit form of the likelihood. This section deals with this objective.

The Maximum Likelihood

Following [START_REF] Greene | Econometric Analysis[END_REF], the full information maximum-likelihood (FIML) estimates are obtained by maximizing the log-likelihood LogL(Y |Z; θ, Ω), where θ is the vector of identified parameters and Ω is the covariance matrix. 1 Under the usual regularity conditions2 (Lesaffre and Kauffmann, 1992), the likelihood is given by the joint density of observed outcomes:

L(y|z, θ; Ω) = T t=1 L t (y t |z t-1 , θ; Ω), (6)
where y t = (y 1,t , ..., y M,t ) and y = [y 1 , ..., y T ]. The individual likelihood L t (.) is given in Lemma 1.

Lemma 1. The likelihood of observation t is the cumulative density function, evaluated at the vector w t of a M-variate standardized normal vector with a covariance matrix Q t ΩQ t :

L t (y t |z t-1 , θ; Ω) = Pr(y 1 = y 1,t , ..., y M = y M,t ) = Φ M,εt (w t ; Q t ΩQ t ), (7)
where Q t is a diagonal matrix whose main diagonal elements are q m,t = 2y m,t -1 and thus depends on the realization or not of the events (q

m,t = 1 if y m,t = 1 and q m,t = -1 if y m,t = 0, ∀ m = {1, 2, ..., M }).
Besides, the elements of the vector w t = [w 1,t , ..., w M,t ] are given by w m,t = q m,t π m,t , where π m,t is the index associated with the m th binary variable.

For a complete proof of Lemma 1, see Appendix 1. Thus, the FIML estimates are obtained by maximizing the log-likelihood:

LogL(y|z, θ; Ω) = T t LogΦ M,ε (w t ; Q t ΩQ t ) (8)
with respect to θ and Ω.

The main problem with FIML is that it requires the evaluation of high-order multivariate normal integrals while existing results are not sufficient to allow accurate and efficient evaluation for more than two variables (see [START_REF] Greene | Econometric Analysis[END_REF], page 714). Indeed, [START_REF] Greene | Econometric Analysis[END_REF] argues that the existing quadrature methods to approximate trivariate or higher-order integrals are far from being exact. To tackle this problem in the case of a static probit, Huguenin, Pelgrin and Holly (2009) decompose the triple integral into simple and double integrals, leading to an Exact Maximum Likelihood Estimation (EML) that requires computing double integrals. Most importantly, they prove that the EML increases the numerical accuracy of both the slope and covariance parameters estimates, which outperform the maximum simulated likelihood method (McFadden,1989) which is generally used for the estimation of multivariate probit models. Therefore, we extend the decomposition proposed by Huguenin, Pelgrin and Holly, (2009) in the case of our multivariate dynamic model so as to obtain a direct approximation of the trivariate normal cumulative distribution function.

The EML log-likelihood function is given by:

LogL(y|z, θ; Ω) = T t=1 Log M m=1 Φ(w m,t ) + G , (9)
where Φ(w t ) is the univariate normal cumulative distribution function of w t . Indeed, the log-likelihood function depends on the product of the marginal distributions (w t ) and the correction term G which captures the dependence between the m events analyzed. The maximum likelihood estimators { θ; Ω} EM L are the values of θ and Ω which maximize (9):

{ θ; Ω} EM L = Arg max θ;Ω M m=1 LogL(.), (10)
with LogL(.) given in [START_REF] Carey | Modelling Multivariate Binary Data with Alternating Logistic Regressions[END_REF]. Under the regularity conditions of Lesaffre and Kaufman (1992), the EML estimator of a multivariate probit model exists and is unique. Besides, the estimates { θ; Ω} EM L are consistent and efficient estimators of the slope and covariance parameters and are asymptotically normally distributed. It is worth noting that in a correctly specified model for which the error terms are independent across the m equations the EML function corresponds to T t=1 M m=1 Φ(w m,t ), since the probability correction term G in eq. ( 9) tends toward zero.

The Empirical Procedure

Nonetheless, most of the empirical applications involving correlated time-series binary data (e.g. financial crises, economic cycles, etc.) tackle only two or three such events at a time. Accordingly, without any loss of generality, henceforth we restrict our attention to the bivariate and trivariate form of the model. This presentation simplifies the comprehension of the exact maximum-likelihood estimation method and corresponds to the empirical analysis performed in section 5. Further details are provided in Appendices 1-3:

Φ 2 (w t ; Q t ΩQ t ) = Φ(w 1,t )Φ(w 2,t ) 1 2π ρ 12 0 exp - 1 2 w 2 1,t + w 2 2,t -2w 1,t w 2,t 1 -λ 2 12 dλ 12 1 -λ 2 12 (11)
9 for a bivariate model and

Φ 3 (w t ; Q t ΩQ t ) = 3 m=1 Φ(w m,t ) + G = Φ(w 1,t )Φ(w 2,t )Φ(w 3,t ) + Φ(w 3,t ) ρ 12 0 φ 2 (w 1,t ,
for a trivariate model, where ρ are the non-diagonal elements of the Q t ΩQ t matrix and λ are the non-diagonal elements of a theoretical 2 × 2 matrix and respectively a 3 × 3 matrix in which one of the correlation coefficients is null. Moreover, ẇt is a vector of indices obtained by changing the order of the elements to (w 2,t , w 3,t , w 1,t ). Similarly ẅt corresponds to a vector of indices of the form (w 3,t , w 1,t , w 2,t ). Finally, ẇt corresponds to w t , ẇt or ẅt respectively, depending on the way the last integral is decomposed. The computation of the last term is not trivial. However, this integral can be decomposed in a non-unique way as follows:

ρ 12 0 ρ 13 0 ρ23 0 ∂ 3 φ
These finite-range multiple integrals are numerically evaluated by using a Gauss-Legendre Quadrature rule3 over bounded intervals. In such a context, two possibilities can be considered: whether the likelihood function is directly maximized, or the first order conditions4 are derived so as to obtain an exact score vector. As stressed by Huguenin, Pelgrin and Holly (2009), the two methods may not lead to the same results if the objective function is not sufficiently smooth. We also tackle the autocorrelation problem induced by some binary time-series variables by considering a Gallant correction for the covariance matrix of the parameters.

Impulse-Response Analysis

Since Sims' (1980) seminal paper, the concept of impulse-response function (IRF) has been playing a major role in the analysis of linear dynamic models. It has later been extended to the general case of non-linear time-series by [START_REF] Gallant | Nonlinear Dynamic Structures[END_REF] and further improved by [START_REF] Potter | Nonlinear impulse response functions[END_REF]. However, as far as we know, no non-linear VAR analysis focusing specifically on the IRFs associated with a dynamic multivariate binary model has been proposed to date. This is intriguing, as given the characteristics of this type of model, namely the presence of a threshold effect, it is particularly relevant to analyze the IRFs associated with this non-linear VAR.

Recall that in the classical linear VAR, the IRFs are immediately obtained from the vector moving average (VMA) representation of the model and their significance can be studied by reporting parametric or non-parametric confidence intervals. By contrast, dynamic multivariate binary models cannot be written in an equivalent VMA representation as it is usually done for linear VAR models by relying on Wold's theorem [START_REF] Wold | Data sample[END_REF]. To overcome this problem we hence propose to express the IRFs in terms of the latent model, i.e. the probability of being in a certain regime (0 or 1) and the observed binary processes by elaborating on eq. 4. Orthogonal impulse-response functions are considered, where the order of the variables (shown to be crucial for such an analysis) is given by economic theory. The exogenous variables are fixed to a certain value, denoted by x, which is usually given by the unconditional mean (x t ) or the last observation in the sample (x T ). The initial value of the vector of lagged latent variables can be set to 0 (which is equivalent to a marginal probability of 0.5 to know each of the events under analysis) or to any other value that is reasonable for a specific application. Additionally, as in any non-linear model, the initial state of the world must be defined for each binary process. Let us denote these M × 1 vectors by ỹ * 0 and ỹ0 , respectively. It follows that the M IRFs are obtained recursively for each period h, h ≥ 1 by iterating on:

ỹ * m,h = αm + βm x + M m =1 ∆m,m ỹm ,h-1 + M m =1 Γm,m ỹ * m ,h-1 + εm , Pr(ỹ m,h = 1) = Φ(ỹ * m,h ), ỹm,h = 1 (ỹ * m,h >0) = 1 (Pr(ỹ m,h =1)>0.5) , (14)
where α, β, ∆ and Γ result from the estimation of the multivariate binary model in [START_REF] Bussière | Balance of Payment Crises in Emerging Markets. How Early Were the "Early" Warning Signals?[END_REF] and the correlated shocks εm are set to 0 for h > 1 and transformed via the Choleski decomposition of the covariance matrix Ω for h = 1, as generally done in linear VAR literature.

To be more precise, if P denotes the lower triangular matrix obtained from the Choleski decomposition, and P m is the m th row of this matrix, then εm = P m * ξ for h = 1, where ξ is the vector of the magnitudes of the orthogonal shocks. Note that this vector includes only one value different than 0, which indicates the equation (the binary process) that is exposed to the shock. It is worth noting here that in the context of dynamic multivariate binary models, that shock should be set to a higher value than one standard deviation (as generally done in VAR models), because it should allow for a potential regime switch in at least one of the binary processes. Indeed, a shock at time h = 1 on the innovations of the binary event m will modify the latent variable y * m,h and it will always be transmitted to the next period y * m ,h+1 for the same event m and not only. But if the shock is strong enough, it also impacts the observable variable y m,h , which switches to the other regime and leads to a non-linear transmission of the shock to y * m ,h+1 . Besides, the length of the IRF path should be chosen so that the impact of the shock on each of the processes completely disappears or converges to its log-term value up to the last period considered.

To analyze the significance of the IRFs we propose to construct bootstrapped confidence intervals. For this, the following four steps are followed:

Step 1: Draw a vector of parameters θ = ( θ 1 ; θ 2 ; ...; θ M ) , with θm = (α m ; βm ; ∆m,m ; Γm,m ) from a M -variate normal distribution with mean θ, i.e. the estimated vector of parameters, and variance Σ, i.e. the estimated covariance matrix associated with these parameters.

Step 2: Proceed to the Choleski decomposition of the covariance matrix Ω to write the correlated shocks as a function of the orthogonal ones, as for the computation of the IRFs.

Step 3: For the same set of initial values as in the case of the IRFs, i.e. xm , ỹ * b 0 = ỹ * 0 and ỹb 0 = ỹ0 , construct a simulated path of the impulse-response functions:

ỹ * b m,h = αm + βm x + M m =1 ∆m,m ỹb m ,h-1 + M m =1 Γm,m ỹ * b m ,h-1 + εm , Pr(ỹ b m,h = 1) = Φ(ỹ * b m,h ), ỹb m,h = 1 (ỹ * b m,h >0) = 1 (Pr(ỹ b m,h =1)>0.5) . (15)
Step 4: Repeat this procedure a large number of times (usually 10000 repetitions are considered) and take the 2.5% and 97.5% quantiles of the IRF distribution to define the 5% confidence interval for each qualitative variable.

In this context, it is important to distinguish between a significant IRF and a significant shift from a regime to the other. First, to simplify interpretation, IRFs are demeaned, i.e. the unconditional mean αm + βm x is subtracted for each of the M processes. It follows that they are significant if the corresponding confidence intervals do not include the value of 0. Second, the probability to shift from one regime to the other or the probability of remaining in the same regime is significantly different from zero at time h if the IRF (m, h), i..e ỹ * m,h , is positive, i.e. the centered IRF (m, h) is significantly lying above the opposite of its unconditional mean (-α -βm x). While the first analysis is common to all vector autoregressive (VAR) models, the second one is specific to non-linear (threshold) time series models.

Empirical Application

This section aims at implementing the multivariate dynamic probit methodology presented above to a system composed by three types of financial crises, i.e. currency, banking and sovereign debt crises. As historical events have proven,most of the time crises do not remain restricted to a single market, but tend to spill-over into another one. Several perspectives are hence needed to get an exhaustive picture of a turmoil. We thus evaluate the probability of mutation of one type of crisis into another one. After a short data description, we estimate bivariate models by excluding sovereign debt crises. This constitutes a benchmark for the second part where the sovereign debt crises are included in the system.

The Database

Monthly macroeconomic indicators expressed in US dollars covering the period from January 1985 to June 2010 have been extracted for 17 emerging countries 5 from the IMF-IFS database as well as from the national banks of the countries under analysis via Datastream. 6 The three types of crises are identified by relying on popular measures, generally considered in the literature. 7 To be more precise, as in [START_REF] Lestano | A Comparison of Currency Crisis Dating Methods: East Asia 1970-2002[END_REF] 2004) is considered for the monthly identification of banking crises, while the non-parametric method based on sovereign debt spread introduced by Pescatori and Sy (2007) is used to detect debt-servicing difficulties. To this aim, the government bond returns are obtained via the JPMorgan EMDB database.

We have selected the main leading indicators used in the literature for the three types of crises that we analyze (see [START_REF] Candelon | How to Evaluate an Early Warning System? Towards a Unified Statistical Framework for Assessing Financial Crises Forecasting Methods[END_REF][START_REF] Jacobs | Indicators of Financial Crises do work! An Early Warning System for six Asian Countries[END_REF][START_REF] Glick | Banking and Currency Crises: How Common Are The Twins?[END_REF][START_REF] Hagen | Money Market Pressure and the Determinants of Banking Crises[END_REF][START_REF] Pescatori | Are Debt Crises Adequately Defined?[END_REF], namely, the one-year growth rate of international reserves, the growth rate of M2 to reserves ratio, one-year growth of domestic credit over GDP ratio, one-year growth of domestic credit, one-year growth of GDP, government deficit, debt service ratio and external debt ratio. Remarks 1. As in [START_REF] Kumar | Predicting Emerging Market Currency Crashes[END_REF], we dampen the magnitude of every variable using the formula : f (x t ) = sign(x t )log(1 + |x t |), so as to reduce the impact of extreme values. 8 2. It should also be noted that the entire sample is used for the identification of currency and banking crises, while the identification of debt crises is realized by using data from December 1997 (See Table 1) since the CDS spread series used for the identification of sovereign debt crises are not available before 1997 in the JPMorgan EMDB database.

Consequently our empirical analysis will consist of two parts, the first one analyzing the case of twin crises (currency and banking) for which the entire database can be used, while the second part focuses on the interactions between the three types of crises 5 Argentina, Brazil, Chile, Colombia, Ecuador, Egypt, El Salvador, Indonesia, Lebanon, Malaysia, Mexico, Panama, Peru, Philippines, South Africa, Turkey and Venezuela. 6 We choose not to include any European country, as i) only few of them have suffered from the three types of crises and ii) if this is the case it corresponds to a single episode: the recent turmoil. 7 For a more detailed description of the three dating methods see the Candelon et al. 2011. 8 Missing values of the series are replaced by cubic spline interpolation.

and is thus based on data from 1997 onwards. The data sample actually used for each of the 17 countries and the two types of analyses is available in Table 1.

3. We only retain the countries for which the percentage of crisis periods is superior to 5% (See Table 2). 9

4. As mentioned in section 2, there are three dynamic multivariate specifications that can be used. However, as shown by Candelon et al. ( 2010), the dynamic model including the lagged binary variable seems to be the best choice according to model selection using the Akaike information criterion. However, since we cannot expect a crisis to have a certain impact on the probability of emergence of another type of crisis from one month to another, which would justify the notation y m,t-1 from the theoretical part, in the empirical application we consider a response lag l of 3, 6 and respectively 12 months for the bivariate models and one of 3 or 6 months for the trivariate models 10 . Therefore, for each type of crisis we build a lagged variable y m,t-l which takes the value of one if there was crisis in the past l periods or at time t, and the value of 0 otherwise:

y m,t-l =      1, if l j=0 y m,t-j > 0 0, otherwise. (16)
5. The significance of the parameters of each model is tested by using simple t-statistics based on robust estimates of standard-errors (which rely on a Gallant kernel, as in [START_REF] Kauppi | Predicting U.S. Recessions with Dynamic Binary Response Models[END_REF]. A special attention is given to the interpretation of cross-effects which stand for the transmission channels of the shocks/crisis. Besides, the joint nullity of the contemporaneous correlations between shocks is tested using a log-likelihood ratio test for the trivariate models.

Bivariate Analysis

Along the lines of Kaminsky et al. (1998) it is possible to find a large number of explanatory variables that may signal the occurrence of a crisis. Nevertheless, Candelon et al. (2010) showed that a univariate dynamic probit model presents the advantage of yielding plausible 9 Argentina, Chile, Ecuador, Egypt, Indonesia, Lebanon, Mexico, South Africa and Venezuela are included in the bivariate analysis, whereas a trivariate model is specified for Ecuador and South Africa. Since the threshold has been arbitrarily set to 5%, we have also checked the borderline countries, like Colombia or Turkey in the bivariate analysis and Egypt in the trivariate analysis respectively, and similar results have been obtained. 10 A 12 months lag is not used in the case of trivariate models since it would significantly reduce the already small number of observations we have at our disposal.

results while being fairly parsimoniously parametrized. Indeed, a large part of the information is integrated either in the past state variable or in the lagged latent variable and thus, only a few explanatory variables turn out to be significant. In this context, we expect their multivariate (bivariate or trivariate) extension to be even more parsimonious. Therefore, we consider the four explanatory variables which are significant in Candelon et al. (2010), i.e. one-year growth of international reserves, one-year growth of M2 to reserves for currency crises as well as one-year growth of domestic credit over GDP and one-year growth of domestic credit for banking crises, resulting in four different specifications including one explanatory variable for each type of crisis. Three different lags (3 months, 6 months and 12 months) are considered for the lagged binary variable y m,t-l . The dynamic probit model is estimated country-by-country using the Exact maximum likelihood. 11 It is indeed a simplification as contagion (or spill-overs) from one country to another are not taken account. A panel version of the model would lead to several problems. First, as shown by [START_REF] Berg | A Cautious Note on the Use of Panel Models to Predict Financial Crises[END_REF] heterogeneity due to country specificities would have to be accounted for. Second, the estimation of a fixed effect panel would be biased without a correction on the score vector. 12Third, in a country by country analysis contagion has to be ignored. For all these reasons, we consider this extension to be beyond the scope of this paper and leave it for future research.

Each model is estimated via maximum-likelihood, the bivariate normal cumulative distribution function being approximated using the Gauss-Legendre quadrature, as proposed by Huguenin, Pelgrin and Holly, (2009). However, the quadrature specified in Matlab by default, i.e. the adaptive Simpson quadrature, has been considered as a benchmark.

Information criteria, namely AIC and SBC, are used to identify the best model for each country; the specification with the lagged binary variable turns out to be preferred. Optimal lag lengths are determined similarly. It is nevertheless worth stressing that the results are generally robust to the choice of explanatory variables and even to the choice of lags.

A summary of the results for the selected models is given in Table 3.

insert Table 3 First of all, it seems that most of the models exhibit dynamics, whatever the lag used to construct the 'past crisis' variable is. This result confirms the findings of [START_REF] Candelon | Currency Crises Early Warning Systems: Why They Should Be Dynamic[END_REF] and [START_REF] Bussière | Balance of Payment Crises in Emerging Markets. How Early Were the "Early" Warning Signals?[END_REF], showing that crises exhibit a regime dependence: if the country is proven to be more vulnerable than investors had initially thought, investors will start withdrawing their investments, thus increasing the probability of a new crisis. More precisely, most of the countries are found to have experienced banking and currency crises, with a significant autoregressive coefficient, i.e. the crisis variable depends on its own past, e.g. Argentina, Egypt, Lebanon, Mexico, South Africa, Venezuela. Besides, only for a small number of cases, only one of the two types of crises is best reproduced by a dynamic model (currency crises in Chile (3 and 12 months), Mexico (6 and 12 months); banking crises in Argentina (6 and 12 months), Ecuador, Lebanon (6 months), South Africa (12 months) and Venezuela (12 months)). Actually, in Chile a past currency crisis has only a short term positive impact on the emergence of another currency crisis, whereas a banking crisis has just a long term effect on the probability of occurrence of another banking crisis. Mexico, however, seems to be more prone to recurring currency crises than banking crises as the former type of crisis has a long-term impact on the probability of experiencing a new crisis, whereas the latter has a positive effect only in the short run. On the contrary, for Argentina, South Africa and Venezuela the impact of past banking crises on currency crises is longer (up to one year) as opposed to that of past currency crises on banking ones (up to three and six months, respectively).

Second, for the majority of the countries (Argentina, Chile, Lebanon, Mexico and Venezuela), currency and banking crises are interconnected. This link between crises can take two forms. On the one hand, a certain type of crisis increases (or diminishes) the probability of occurrence of the other type of crisis. This strong link between banking to currency crisis was emphasized by Glick and Hutchinson (1999) within a panel framework. Nevertheless, there is no reason for the transmission of shocks to be symmetric. Indeed, our country per country analysis reveals that for some countries like Argentina (3 and 6 months) a banking crisis in the past months increased the probability of a currency crisis at time t. At the same time, a banking crisis in Chile in the last 12 months reduced the probability of experiencing a currency crisis. Conversely, a currency crisis in Egypt and in Lebanon (3 months) diminished the probability of a banking crisis.

On the other hand, crisis shocks can be contemporaneously positively correlated. This feature seems to be very stable across models (independent of the lag used). The only exceptions are Egypt and Lebanon, for which there is no instantaneous correlation in the model with 3-months lagged binary variables and Mexico, for which such a correlation appears only for the 12-months lag. To sum up, but for Egypt, all countries are characterized by a positive instantaneous correlation between shocks of currency and banking crises variables, corroborating the previous findings of Glick and Hutchinson (1999).

Third, the macroeconomic variables are rarely significant. 13 These results corroborate our previous findings (see [START_REF] Candelon | Currency Crises Early Warning Systems: Why They Should Be Dynamic[END_REF]) that the dynamics of crises captures most of the information explaining the emergence of such phenomena. Furthermore, when these coefficients are significant, they have the expected sign (an increase in the growth of international reserves diminishes the probability of a crisis, while a surprise in the rest of indicators 13 These results are available upon request.

soars the probability of a crisis).

To summarize, these results confirm the presence of interaction between the banking and currency crises. The twin crisis phenomenon is thus confirmed empirically. Our findings are also robust to the quadrature choice and the lags considered when constructing the dynamic binary variables.

Trivariate Analysis

But is it really enough to look at two crises only? This subsection extends the previous analysis to the trivariate case by modeling simultaneously the occurrence of currency, banking and debt crises. However, only two countries experienced these three events during a sufficiently long period. Ecuador presents for our sample an ex-post probability larger than 5% for each type of crisis. Such a result is not surprising if one remembers that Ecuador faced a strong financial turmoil in the late 1990, affecting first the banking sector, 14 then the Sucre 15 , and the government budget. [START_REF] Jacone | The Late 1990s Financial Crisis in Ecuador: Institutional Weaknesses, Fiscal Rigidities, and Financial Dollarization at Work[END_REF] showed that institutional weaknesses, rigidities in public finances, and high financial dollarization have amplified this crisis. South Africa constitutes a borderline case as the sovereign debt crisis probability is slightly below 5%.

Each of the models is estimated for these countries using both the methodology proposed by [START_REF] Huguenin | Estimation of Multivariate Probit Models by Exact Maximum Likelihood[END_REF] based on the Gauss-Legendre quadrature and the direct approximation of a triple integral based on the adaptive Simpson quadrature that Matlab uses by default. Similar results are obtained for the two methods. 16 However, the latter implies a significant gain in time without any loss in accuracy proving that recently developed quadrature methods are good approximations of the normal cumulative distribution function. Besides, 6 and 12 month-lags of the dynamic crisis variable are considered.

insert Table 4 In the case of Ecuador, the results corroborate our bivariate findings: the banking crises are persistent, while currency crises are not. Nevertheless, it is clear that the bivariate model is misspecified, since it cannot capture the impact of a banking crisis on the occurrence of a currency crisis when using the 6-months lagged binary variables to account for the dynamics of these phenomena (see Table 4).

Moreover, the trivariate model turns out to be more parsimonious in terms of parameters to be estimated since the latent variable of past debt crisis has a positive effect on the probability of occurrence of both currency and debt crises. Therefore it supports the implementation of a trivariate crisis model whenever when it is feasible. We also observe that the contemporaneous correlation matrix is diagonal, ruling out common shocks. Crises in Ecuador turn out to be exclusively driven by transmission channels, as in the late 1990, when the banking distress was diffused to the currency and the government budget.

In the case of South-Africa, both currency and debt crises are dynamic. There is no evidence of causality between the different types of crises, but significant contemporaneous correlation. It highlights the fact that contrary to Ecuador, South African crises did not mutate but they originated from a common shock. It is worth noting that the results are found to be robust in the sensitivity analyses performed, namely to the choice of macroeconomic variables and the use of different lags for the past crisis variables.

Further results

To grasp better the properties of the models estimated and selected, a conditional probability as well as an IRF analysis are provided. For sake of space, we only report the results obtained for Ecuador. 17First, Figure 1 reports the conditional probabilities for each type of crisis obtained from both the bi-and trivariate models considering a forecast horizon of 3 and 6 months. To allow a fair comparison, both models are estimated from the same sample, i.e. from 1997 onwards. It goes without saying that the bivariate model does not provide any conditional probabilities for sovereign debt crisis.

It turns out that the trivariate model outperforms the bivariate one whatever the forecast horizon is, i.e. the conditional probabilities issued from the trivariate model are higher than those obtained from the bivariate model during observed crisis periods, while they appear to be similar for calm periods. Such results corroborate hence our previous findings, stressing that a crisis model should take into account the whole sequence of crises to be accurate. Besides, the conditional probabilities obtained from the trivariate model do not immediately collapse after the occurrence of the crisis, which is the case for the bivariate model. It stresses hence the vulnerability of the economy after the exit from a turmoil in particular if it affects the foreign exchange market.

Second, to evaluate the effect of a crisis, considered here as a shock, an IRF analysis is performed for the trivariate model. As the order of the variables has been shown to be crucial, we consider the historical sequence of crises observed in Ecuador, i.e. banking crises (the most exogenous ones), debt crises and currency crises (the most endogenous ones).

Orthogonal impulse response functions are considered on the latent variable for a 3 monthhorizon. The exogenous variables are fixed to the unconditional mean (x m,t ). Additionally, as in any non-linear model, the IRFs are calculated for two initial states: a tranquil one, ỹt-1 = 0, i.e. "no type of crisis is observed at time t = 0 or in the previous 3 months" and a turmoil regime, ỹt-1 = 1, i.e. "all types of crisis are observed in t = 0". Confidence intervals are built taking the 2.5% and 97.5% percentiles of IRF's distribution obtained from 10,000 simulations of the model. The magnitude of the shock is fixed to 5,18 allowing for a potential mutation of the crisis.

Figures 2 to 4 report the diffusion through the system of banking, currency and debt crises respectively. First, it appears in figure 2 that a banking crisis shock has almost no persistence in a calm initial state, as the IRF function reverses to mean after a single period. On the contrary the persistence jumps to 5 months for an initial crisis state. Similarly, the diffusion of a banking crisis shock to another type of turmoil is exclusively observed in a crisis initial state. Besides, the shift probability from calm to crisis period is significant only for the banking crisis and up to the second period (see the left part of figure 2), whereas the probability of remaining in a crisis period is significant for all three types of crises until t = 2 (see the right part of figure 2). This underlines the uncertainty surrounding the duration of a crisis beyond one month after the shock. Overall, these first results clearly correspond to the path exhibited by the crisis sequence faced by Ecuador in the late 90's. Figure 3 reports the response of the three latent variables to a debt crisis shock. In such a case, the impact of the shock on the banking and currency crises vanishes almost instantaneously in the case of a calm initial state, while it disappears after 4 or 5 months, if the economy is facing initially a joint crisis. As for the debt crisis, the impact of the shock lasts at least 5 months even though we are certain of being in a crisis period during the first two periods (the confidence interval is in the grey area at that time). Finally, Figure 4 presents the IRF after currency crisis shock. As in the previous cases, the impact on the banking crisis is not important if we depart from a calm situation, while it becomes significant during 4 periods for an initial crisis period. At the same time, the response of the debt crisis is slowly dampened towards the baseline for a calm initial state, whereas it is significant during the first 4 periods if the shock occurs while being in a crisis state. It seems that the persistence of the effect of this shock is around two months for a calm initial period while it dies away only after 5 months in the alternative situation.

Overall, the conditional probability and the IRF analyses stress the superiority of the trivariate model to scrutinize the diffusion mechanisms that occurred in Ecuador after the banking crisis in 1998. Strong interactions between the three types of crises are clearly present in particular between banking and other crises.

Conclusion

This paper proposes a multivariate dynamic probit model to produce dynamic forecasts of binary time-series correlated variables. It is easy to implement and relies on an exact maximum-likelihood estimation approach, hence providing a solution to the problem generally encountered in the estimation of multivariate probit models. For this, higher-order integrals are decomposed into lower-order finite-range multiple integrals, that are subsequently evaluated using quadrature-rules over bounded intervals. Our framework allows us to apprehend dynamics and causality in several ways, namely through the lagged binary indicators or the lagged latent variables associated with the qualitative variables. Besides, we propose an impulse-response function (IRF) analysis adapted to this type of models.

To illustrate this methodology we consider an application on three types of financial crises for a sample of emerging countries. To be more precise, we investigate the potential mutations between currency, banking and sovereign debt crises within each country and find that in the bivariate case causality from banking crises to currency crises (and vice-versa) are quite common. More importantly, for the two countries, Ecuador and South Africa, which suffered from the three types of crises, the trivariate model turns out to be the best performing in term of conditional probabilities and comprehension of the reasons why a specific crisis mutates to another one: this can be due to either common shocks (as in South Africa) or to a strong causal structure (as in Ecuador). The conditional probability and IRF analyses reinforce these findings. More generally, this application advocates the use of trivariate probit crisis models whenever it is possible, so as to have a better insight on the financial turmoils.

Appendix 1: Proof of lemma 1

By definition, the likelihood of observation t is given by: L t (y t |z t-1 , θ; Ω) = Pr((-q 1,t y * 1,t ≤ 0), ..., (-q M,t y * M,t ≤ 0)) = Pr(-q 1,t ε 1,t ≤ q 1,t π 1,t , ..., -q M,t π M,t ≤ q M,t π M,t )

= Φ M,-Qtεt (w t |0 M ; Ω) = w M,t -∞ ... w 1,t -∞ φ M,-Qtεt (Q t ε t , Ω) M m=1 dε m,t .
Since each q m,t takes only the values {-1, 1}, it is straightforward to show that Q t = Q -1 t and |Q t ΩQ t | = |Ω|. Moreover, the density of an M-variate standardized normal vector -Q t ε t with covariance matrix Ω may be re-written as the density of an M-variate standardized normal vector ε t with variance-covariance matrix Q t ΩQ t :

φ M,-Qtεt (Q t ε t ; Ω) = |2πΩ| -1 2 exp -1 2 (-Q t ε t ) Ω -1 (-Q t ε t ) = |2π(Q t ΩQ t )| -1 2 exp -1 2 ε t (Q t ΩQ t ) -1 ε t = φ M,εt (ε t ; Q t ΩQ t ).
Therefore, the likelihood of observation t is given by:

L t (y t |Z t-1 , θ; Ω) = q M,t π M,t -∞ ... q 1,t π 1,t -∞ φ M,εt (ε t ; Q t ΩQ t ) M m=1 dε m,t = Φ M,εt (Q t π t ; Q t ΩQ t ).
Appendix 2: The Gauss-Legendre Quadrature rule

The goal of the Gauss-Legendre Quadrature rule is to provide an approximation of the following integral:

b a f (x)dx. ( 17
)
In a first step, the bounds of the integral must be changed from [a, b] to [-1,1] before applying the Gaussian Quadrature rule:

b a f (x)dx = b -a 2 1 -1 f (z)dz, (18)
where z i = b-a 2 abs i + b+a 2 and the nodes abs i , i ∈ {1, 2, ..., p} are zeros of the Legendre polynomial P p (abs). Definition 1. Then, the standard p-point Gauss-Legendre quadrature rule over a bounded arbitrary interval [a,b] is given by the following approximation:

b a f (x)dx ≈ b -a 2 p i=1 v i f (z i ) + R p , (19)
where v i are the corresponding weights, v i = Appendix 3: The EML score vector for a trivariate dynamic probit model

For ease of notation, let us denote by ρ i,j , i, j = {1, 2, 3}, i = j the correlation coefficients associated to the Ω matrix. The likelihood of observation t may be written as:

P t = Φ 3 (q 1 π 1,t , q 2 π 2,t , q 3 π 3,t , q 1 q 2 ρ 12 , q 1 q 3 ρ 13 , q 2 q 3 ρ 23 )

= Φ(q 1 π 1,t )Φ(q 2 π 2,t )Φ(q 3 π 3,t )

+ q 1 q 2 Φ(q 3 π 3,t )Ψ 2 (π 1,t , π 2,t , ρ 12 )

+ q 1 q 3 Φ(q 2 π 2,t )Ψ 2 (π 1,t , π 3,t , ρ 13 ) + q 2 q 3 Φ(q 1 π 1,t )Ψ 2 (π 2,t , π 3,t , ρ 23 ) + q 1 q 2 q 3 Ψ 3 (π 3,t , π 1,t , π 2,t , ρ 13 , ρ 23 , 0) + q 1 q 2 q 3 Ψ 3 (π 2,t , π 3,t , π 1,t , ρ 23 , ρ 12 , 0)

+ q 1 q 2 q 3 Ψ 3 (π 1,t , π 2,t , π 3,t , ρ 12 , ρ 13 , ρ 23 ), (20)
where Ψ 2 (π Therefore, the first order partial derivatives can be obtained as follows:

∂ ∂π 1 P t = q 1 φ(π 1,t )Φ(q 2 π 2,t )Φ(q 3 π 3,t ) + q 1 q 2 Φ(q 3 π 3,t ) ∂ ∂π 1 Ψ 2 (π 1,t , π 2,t , ρ 12 ) + q 1 q 3 Φ(q 2 π 2,t ) ∂ ∂π 1 Ψ 2 (π 1,t , π 3,t , ρ 13 ) + q 1 q 2 q 3 φ(π 1,t )Ψ 2 (π 2,t , π 3,t , ρ 23 ) + q 1 q 2 q 3 ∂ ∂π 1 Ψ 3 (π 3,t , π 1,t , π 2,t , ρ 13 , ρ 23 , 0) + q 1 q 2 q 3 ∂ ∂π 1 Ψ 3 (π 2,t , π 3,t , π 1,t , ρ 23 , ρ 12 , 0) + q 1 q 2 q 3 ∂ ∂π 1 Ψ 3 (π 1,t , π 2,t , π 3,t , ρ 12 , ρ 13 , ρ 23 ), (21)
∂ ∂π 2 P t = q 2 φ(π 2,t )Φ(q 1 π 1,t )Φ(q 3 π 3,t ) + q 1 q 2 Φ(q 3 π 3,t ) ∂ ∂π 2 Ψ 2 (π 1,t , π 2,t , ρ 12 ) + q 1 q 2 q 3 φ(π 2,t )Ψ 2 (π 1,t , π 3,t , ρ 13 ) + q 2 q 3 Φ(q 1 π 1,t ) ∂ ∂π 2 Ψ 2 (π 2,t , π 3,t , ρ 23 ) + q 1 q 2 q 3 ∂ ∂π 2 Ψ 3 (π 3,t , π 1,t , π 2,t , ρ 13 , ρ 23 , 0) + q 1 q 2 q 3 ∂ ∂π 2 Ψ 3 (π 2,t , π 3,t , π 1,t , ρ 23 , ρ 12 , 0) + q 1 q 2 q 3 ∂ ∂π 2 Ψ 3 (π 1,t , π 2,t , π 3,t , ρ 12 , ρ 13 , ρ 23 ), (22)
∂ ∂π 3 P t = q 1 φ(π 3,t )Φ(q 1 π 1,t )Φ(q 2 π 2,t ) + q 1 q 2 q 3 φ(π 3,t )Ψ 2 (π 1,t , π 2,t , ρ 12 )

+ q 1 q 3 Φ(q 2 π 2,t ) ∂ ∂π 3 Ψ 2 (π 1,t , π 3,t , ρ 13 ) + q 2 q 3 Φ(q 1 π 1,t ) ∂ ∂π 3 Ψ 2 (π 2,t , π 3,t , ρ 23 ) + q 1 q 2 q 3 ∂ ∂π 3 Ψ 3 (π 3,t , π 1,t , π 2,t , ρ 13 , ρ 23 , 0) + q 1 q 2 q 3 ∂ ∂π 3 Ψ 3 (π 2,t , π 3,t , π 1,t , ρ 23 , ρ 12 , 0) + q 1 q 2 q 3 ∂ ∂π 3 Ψ 3 (π 1,t , π 2,t , π 3,t , ρ 12 , ρ 13 , ρ 23 ), (23)
∂ ∂ρ 12 P t = q 1 q 2 Φ(q 3 π 3,t ) ∂ ∂ρ 12 Ψ 2 (π 1,t , π 2,t , ρ 12 ) + q 1 q 2 q 3 ∂ ∂ρ 12 Ψ 3 (π 2,t , π 3,t , π 1,t , ρ 23 , ρ 12 , 0) + q 1 q 2 q 3 ∂ ∂ρ 12 Ψ 3 (π 1,t , π 2,t , π 3,t , ρ 12 , ρ 13 , ρ 23 ), (24)
∂ ∂ρ 13 P t = q 1 q 3 Φ(q 2 π 2,t ) ∂ ∂ρ [START_REF] Fuertes | Optimal Design of Early Warning Systems for Sovereign Debt Crises[END_REF] Ψ 2 (π 1,t , π 3,t , ρ 13 ) + q 1 q 2 q 3 ∂ ∂ρ [START_REF] Fuertes | Optimal Design of Early Warning Systems for Sovereign Debt Crises[END_REF] Ψ 3 (π 3,t , π 1,t , π 2,t , ρ 13 , ρ 23 , 0)

+ q 1 q 2 q 3 ∂ ∂ρ 13 Ψ 3 (π 1,t , π 2,t , π 3,t , ρ 12 , ρ 13 , ρ 23 ), (25)
∂ ∂ρ 23 P t = q 2 q 3 Φ(q 1 π 1,t ) ∂ ∂ρ 23 Ψ 2 (π 2,t , π 3,t , ρ 23 ) + q 1 q 2 q 3 ∂ ∂ρ 23 Ψ 3 (π 3,t , π 1,t , π 2,t , ρ 13 , ρ 23 , 0) + q 1 q 2 q 3 ∂ ∂ρ 23 Ψ 3 (π 2,t , π 3,t , π 1,t , ρ 23 , ρ 12 , 0) + q 1 q 2 q 3 ∂ ∂ρ 23 Ψ 3 (π 1,t , π 2,t , π 3,t , ρ 12 , ρ 13 , ρ 23 ), (26)
where Note: Three different lags of the dependent variable are used, namely 3, 6 and 12 months. '∆' stands for the parameters of the lagged crisis variables, while Ω represents the covariance matrix. A'+'/'-' sign means that the coefficient is significant and positive/ negative, while a '.' indicates its non-significance. For example, in the case of Argentina, 3 months, all the parameters are positive and significative except for the impact of a currency crisis on the probability of occurrence of banking crises. Similarly, the correlation coefficient between currency and banking crises is significative.

∂ ∂π 1 Ψ 3 (π 3,t , π

  

Note: Two different lags of the dependent variable are used, namely 3 and 6 months. '∆' stands for the parameters of the lagged crisis variables, while Ω represents the variance-covariance matrix. A'+'/'-' sign means that the coefficient is significant and positive/ negative, while a '.' indicates its non-significance. For example, in the case of Ecuador, 3 months, sovereign debt crises have a positive and significative impact on the probability of occurrence of currency crises.

  the vector of parameters associated with the exogenous variables and ∆ m,m ∈ R and Γ m,m ∈ R are the parameters of the predetermined variables giving the dynamics of the m th equation of the model. s > 0 is the lag associated with the predetermined right-hand-side variables. Denote by θ m = (α m ; β m ; ∆ m,m ; Γ m,m ) the vector of parameters for equation / event m in (

  and Candelon et al. (2012), we use a modified version of the pressure index proposed by Kaminski et al.(1998) to date currency crises. Besides, the money market pressure index proposed by Hagen and Ho (

i=1 v i = 2 ,

 2 and R p is theerror term, R p = Q p f (2p) (ξ) = (b-a) 2p+1 (p!) 4 (2p+1)(2p!) 3 f 2p (ξ), with ξ ∈ (a, b).

Figure 1 :

 1 Figure 1: Conditional crisis probabilities -EcuadorNote: Probabilities at time t are calculated including observed information prior 3 or 6 months.

32

 32

Figure 2 :

 2 Figure 2: IRF after a banking crisis shock -Ecuador 3 months 33

Figure 3 :

 3 Figure 3: IRF after a debt crisis shock -Ecuador 3 months 34

Figure 4 :

 4 Figure 4: IRF after a currency crisis shock -Ecuador 3 months 35

  3 ( ẇt ; λ 12 , λ 13 , λ 23 ) ∂w 1,t ∂w 2,t ∂w 3,t dλ 12 dλ 13 dλ 23 (w t ; λ 12 , λ 13 , λ 23 ) ∂w 1,t dλ 12 dλ 13 -

		ρ 13	ρ 23		ρ 13		ρ 23
	=		∂φ 3 ( ẅt ; λ 12 , λ 13 , λ 23 ) ∂w 3,t	dλ 13 dλ 23 -		∂φ 3 ( ẅt ; 0, λ 13 , λ 23 ) ∂w 3,t	dλ 13 dλ 23
		0	0		0	0
		ρ 12	ρ 23		ρ 12		ρ 23
	=		∂φ 3 ( ẇt ; λ 12 , λ 13 , λ 23 ) ∂w 2,t	dλ 12 dρ 23 -		∂φ 3 ( ẇt ; λ 12 , 0, λ 23 ) ∂w 2,t	dλ 12 dλ 23
		0	0		0	0
	=	λ 12	λ 13	∂φ 3 ρ 12	ρ 13	∂φ 3 (w t ; λ 12 , λ 13 , 0) ∂w 1,t	dλ 12 dλ 13 .
		0	0		0		0

  1,t , π 2,t , ρ 12 ) = (π 2,t , π 3,t , λ 23 )dλ 23 , and Ψ 3 (π 3,t , π 1,t , π 2,t , ρ 13 , ρ 23 , 0) = (π 3,t , π 1,t , π 2,t , λ 13 , λ 23 , 0)dλ 13 dλ 23 Ψ 3 (π 2,t , π 3,t , π 1,t , ρ 23 , ρ 12 , 0) = (π 2,t , π 3,t , π 1,t , λ 23 , λ 12 , 0)dλ 23 dλ 12 2 23 + 2λ 12 λ 13 ρ 23 × φ 3 (π 1,t , π 2,t , π 3,t , λ 12 , λ 13 , ρ 23 )dλ 12 dλ 23 .

		ρ 12
		ρ 23
	0 0	0 φ 2 ρ 13 ρ 23 0 -π 3,t + λ 13 π 1,t + λ 23 π 2,t 1 -λ 2 13 -λ 2 23 φ 3 ρ 23 ρ 12 0 -π 2,t + λ 23 π 3,t + λ 12 π 1,t 1 -λ 2 23 -λ 2 12 φ 3 ρ 12 ρ 13
		0

0 φ 2 (π 1,t , π 2,t , λ 12 )dλ 12 Ψ 2 (π 1,t , π 3,t , ρ 13 ) = ρ 13 0 φ 2 (π 1,t , π 3,t , λ 13 )dλ 13 Ψ 2 (π 2,t , π 3,t , ρ 23 ) = Ψ 3 (π 1,t , π 2,t , π 3,t , ρ 12 , ρ 13 , ρ 23 ) = 0 -(1 -ρ 2 23 )π

1,t + (λ 12 -λ 13 ρ 23 )π 2,t + (λ 13 -λ 12 ρ 23 )π 3,t 1 -λ 2 12 -λ 2 13 -ρ

  1,t , π 2,t , ρ 13 , ρ 23 , 0) = (π 3,t , π 1,t , π 2,t , λ 13 , λ 23 , 0)dλ 13 λ 23 (π 3,t , π 1,t , π 2,t , ρ 13 , λ 23 , 0)dλ 23 , ∂ ∂π 2 Ψ 3 (π 3,t , π 1,t , π 2,t , ρ 13 , ρ 23 , 0) = (π 3,t , π 1,t , π 2,t , λ 13 , λ 23 , 0)dλ 23 λ 13 (π 3,t , π 1,t , π 2,t , λ 13 , ρ 23 , 0)dλ 13 , ∂ ∂π 3 Ψ 3 (π 3,t , π 1,t , π 2,t , ρ 13 , ρ 23 , 0) = (π 3,t , π 1,t , π 2,t , λ 13 , λ 23 , 0)dλ 13 dλ 23 , ∂ ∂ρ 13 Ψ 3 (π 3,t , π 1,t , π 2,t , ρ 13 , ρ 23 , 0) = (π 2,t , π 3,t , π 1,t , λ 23 , λ 12 , 0)dλ 12 λ 23 (π 2,t , π 3,t , π 1,t, , λ 23 , ρ 12 , 0)dλ 23 , (π 2,t , π 3,t , π 1,t , λ 23 , λ 12 , 0)dλ 23 dλ 12 , (π 2,t , π 3,t , π 1,t , ρ 23 , ρ 12 , 0) = -π 2,t + λ 23 π 3,t + ρ 12 π 1,t 1 -λ 2 23 -ρ 2 12 φ 3 (π 2,t , π 3,t , π 1,t , λ 23 , ρ 12 , 0)dλ 23 , ∂ ∂ρ 23 Ψ 3 (π 2,t, , π 3,t , π 1,t , ρ 23 , ρ 12 , 0) = (π 2,t , π 3,t , π 1,t , ρ 23 , λ 12 , 0)dλ 12 , ∂ ∂π 1 Ψ 3 (π 1,t , π 2,t , π 3,t , ρ 12 , ρ 13 , ρ 23 ) = -(λ 12 -λ 13 ρ 23 )π 2,t -(λ 13 -λ 12 λ 23 )π 3,t ] 2 + 2λ 12 λ 13 ρ 23 ) 2 φ 3 (π 1,t , π 2,t , π 3,t , λ 12 , λ 13 , ρ 23 )dλ 12 dλ 13 , ∂ ∂π 2 Ψ 3 (π 1,t , π 2,t , π 3,t , ρ 12 , ρ 13 , ρ 23 ) = (π 1,t , π 2,t , π 3,t, , λ 12 , λ 13 , ρ 23 )dλ 12 λ 13 (π 1,t , π 2,t , π 3,t , ρ 12 , λ 13 , ρ 23 )dλ 13 , ∂ ∂π 3 Ψ 3 (π 1,t , π 2,t , π 3,t , ρ 12 , ρ 13 , ρ 23 ) = (π 1,t , π 2,t , π 3,t , λ 12 , λ 13 , ρ 23 )dλ 13 λ 12 (π 1,t , π 2,t , π 3,t , λ 12 , ρ 13 , ρ 23 )dλ 12 , ∂ ∂ρ 12 Ψ 3 (π 1,t , π 2,t , π 3,t , ρ 12 , ρ 13 , ρ 23 ) = -ρ 2 23 )π 1,t + (ρ 12 -λ 13 ρ 23 )π 2,t + (λ 13 -ρ 12 ρ 23 )π 3,t 1 -ρ 2 12 -λ 2 13 -ρ 2 23 + 2ρ 12 λ 13 ρ 23 × φ 3 (π 1,t , π 2,t , π 3,t , ρ 12 , λ 13 , ρ 23 dλ 13 ), ∂ ∂ρ 13 Ψ 3 (π 1,t , π 2,t , π 3,t , ρ 12 , ρ 13 , ρ 23 ) = + (λ 12 -ρ 13 ρ 23 )π 2,t + (ρ 13 -λ 12 ρ 23 )π 3,t 1 -λ 2 12 -ρ 2 13 -ρ 2 23 + 2λ 12 ρ 13 ρ 23 × φ 3 (π 1,t , π 2,t , π 3,t , λ 12 , ρ 13 , ρ 23 dλ 12 ), ∂ ∂ρ 23 Ψ 3 (π 1,t , π 2,t , π 3,t , ρ 12 , ρ 13 , ρ 23 ) = ∂ 2 ∂π 2,t ∂λ 13 φ 3 (π 1,t , π 2,t , π 3,t , λ 12 , λ 13 , ρ 23 )dλ 12 dλ 13 + (λ 12 -ρ 13 ρ 23 )π 1,t + (ρ 23 -λ 12 ρ 13 )π 3,t 1 -λ 2 12 -ρ 2 13 -ρ 2 23 + 2λ 12 ρ 13 ρ 23 × φ 3 (π 1,t , π 2,t , π 3,t , λ 12 , ρ 13 , ρ 23 )dλ 12 ∂ 2 ∂π 3,t ∂λ 12 φ 3 (π 1,t , π 2,t , π 3,t , λ 12 , λ 13 , ρ 23 )dλ 12 dλ 13 + (λ 13 -ρ 12 ρ 23 )π 1,t + (ρ 23 -ρ 12 λ 13 )π 2,t 1 -ρ 2 12 -λ 2 13 -ρ 2 23 + 2ρ 12 λ 13 ρ 23 × φ 3 (π 1,t , π 2,t , π 3,t , ρ 12 , λ 13 , ρ 23 )dλ 13 .

	ρ 23 φ 3 = 0 ρ 13 0 ∂ ∂λ 13 ρ 23 0 φ 3 ρ 13 0 ρ 23 0 ∂ ∂λ 23 φ 3 = ρ 13 0 φ 3 ρ 13 0 ρ 23 × 1 (1 -λ 2 13 -λ 2 23 ) 2 φ 3 ρ 23 ρ 13 0 -π 3,t + λ 13 π 3,t + ρ 23 π 2,t 1 -λ 2 13 -ρ 2 23 Ψ ρ 23 0 ρ 12 0 ∂ ∂λ 12 φ 3 = ρ 23 0 φ 3 ∂ ∂ ∂π 1 ∂π 2 Ψ 3 (π 2,t , π 3,t, , π 1,t , ρ 23 , ρ 12 , 0) = ρ 23 0 ρ 12 0 [(π 2,t -λ 23 π 3,t -λ 12 π 1,t ) 2 -(1 -λ 2 2 13 -λ 2 23 )] φ 3 (π 3 , π 1 , π 2 , λ 13 , ρ 23 , 0)dλ 13 , 23 -λ 2 12 )] × 1 (1 -λ 2 23 -λ 2 12 ) 2 φ 3 ∂ ∂π 3 Ψ ρ 12 0 ρ 23 0 ∂ ∂λ 23 0 φ ∂ ∂ρ 12 Ψ 3 ρ 23 0 ρ 12 0 -π 2,t + ρ 23 π 3,t + λ 12 π 1,t 1 -ρ 2 23 -λ 2 12 φ 3 ρ 12 0 ρ 13 0 {[(1 -ρ 2 23 )π 1,t -(1 -ρ 2 23 )(1 -λ 2 12 -λ13 3 -ρ 2 23 + 2λ 12 λ 13 ρ 23 )}× 1 (1 -λ 2 12 -λ 2 13 -ρ 2 23 ρ 13 0 ρ 12 0 ∂ ∂λ 12 φ 3 = ρ 13 0 φ 3 ρ 12 0 ρ 13 0 ∂ ∂λ 13 φ 3 = ρ 12 0 φ 3 ρ 13 0 (1 ρ 12 (1 -ρ 2 0 ρ 13 0 = ρ 12 0 -(1 -ρ 2 13 )π 2,t = ρ 12 0 ρ 13 0 0 23 )π 1,t ρ 12 = ρ 13 -(1 -ρ 2 12 )π 3,t

0 [(π 3,t -λ 13 π 1,t -λ 23 π 2,t ) 2 -(1 -λ 0 -π 3,t + ρ 13 π 3,t + λ 23 π 2,t 1 -ρ 2 13 -λ 2 23 φ 3 (π 3 , π 1 , π 2 , ρ 13 , λ 23 , 0)dλ 23 , ∂ ∂ρ 23 Ψ 3 (π 3,t , π

1,t , π 2,t , ρ 13 , ρ 23 , 0) = 3 (π 2,t , π 3,t , π 1,t , ρ 23 , ρ 12 , 0) = 3 (π 2,t , π 3,t , π 1,t , ρ 23 , ρ 12 , 0) = 3 (π 2,t , π 3,t, , π 1,t , λ 23 , λ 12 , 0)dλ 23 λ 12 = ρ 12 0 φ 3 (π 2,t , π 3,t , π 1,t , ρ 23 , λ 12 , 0)dλ 12 , 29

Table 2 :

 2 Percentage of crisis periodsThe entries represent the proportion of crises period over the whole sample. It is indicated in bold as it exceeds 5%.

		Bivariate model		Trivariate model
		Currency crisis Banking crisis	Currency crisis Banking crisis Debt crisis
	Argentina	5.13	8.90	4.00	6.67	10.0
	Brazil	3.77	7.19	0.00	3.33	2.67
	Chile	6.07	10.0	5.79	5.79	3.31
	Colombia	4.95	9.90	9.22	12.8	0.00
	Ecuador	5.73	9.93	6.67	10.8	6.67
	Egypt	6.76	9.96	4.17	7.30	7.30
	El Salvador	3.65	9.85	0.00	0.00	2.50
	Indonesia	5.30	9.90	0.00	14.0	6.25
	Lebanon	9.62	9.96	1.38	8.97	2.76
	Malaysia	3.10	10.0	4.05	6.08	4.73
	Mexico	6.50	9.93	0.00	9.33	0.00
	Panama	0.00	9.89	0.00	6.38	0.00
	Peru	4.45	8.22	0.00	10.7	0.00
	Phillipines	4.90	9.80	5.69	6.50	3.25
	South Africa	6.71	9.89	7.09	7.80	4.26
	Turkey	4.80	8.56	4.00	6.67	0.00
	Venezuela	7.33	10.1	4.17	7.64	2.78

Note:

Table 3 :

 3 Bivariate Analysis

			3 months	6 months	12 months
	Country		∆	Ω	∆	Ω	∆	Ω
	Argentina	currency banking	+ + . +	1 + + 1	. + . +	1 + + 1	. . . +	1 . . 1
	Chile	currency banking	+ . . .	1 + + 1	. . . .	1 + + 1	. -. +	1 + + 1
	Ecuador	currency banking	. . . +	1 . . 1	. . . +	1 . . 1	. . . +	1 . . 1
	Egypt	currency banking	+ . -+	1 . . 1	+ . -+	1 --1	+ . . +	1 --1
	Lebanon	currency banking	+ . -+	1 . . 1	. . . +	1 + + 1	+ . . +	1 + + 1
	Mexico	currency banking	+ . . +	1 . . 1	+ . . .	1 . . 1	+ . . .	1 + + 1
	South Africa	currency banking	+ . . +	1 . . 1	+ . . +	1 . . 1	. . . +	1 . . 1
	Venezuela	currency banking	+ . . +	1 + + 1	+ . . +	1 + + 1	. . . +	1 + + 1

Table 4 :

 4 Trivariate Analysis

Note that to identify the slope and covariance parameters, we impose that the diagonal elements of Ω to be standardized, i.e. equal to one.

If the parameters θ are estimated while the correlation coefficients are assumed constant, the loglikelihood function is concave. In this case the MLE exists and it is unique. Nevertheless, when θ and ρ are jointly estimated (as in our model), the likelihood function is not (strictly) log-concave as a function of ρ. Thus, the MLE exists only if the log-likelihood is not identically -∞ and E(z T z|ρ) is upper semi-continuous finite and not identically 0. Furthermore, if no θ = 0 fulfills the first order conditions for a maximum, the MLE of (θ, ρ) for the multivariate probit model exists and for each covariance matrix not on the boundary of the definition interval, the MLE is unique.

Details about this quadrature are available in Appendix 2.

The score vector of the trivariate probit model is presented in Appendix 3.

Initial conditions are introduced as given by the univariate static probit.

See Candelon et al., (2010) for a discussion about this point.

[START_REF] Glonek | Multivariate Logistic Models[END_REF] out of the 40 banks existing in 1997 faced liquidity problems.

The Ecuadorian currency has been replaced by the U.S. dollar on March,13, 2000.

The results for Ecuador when considering a 6-months lag have been obtained with Matlab's quadrature since the model based on the Gauss-Legendre Quadrature did not converge.

For South Africa, crisis mutation is exclusively driven by the contemporaneous correlation matrix as indicated in Table4. Otherwise we can see that currency and sovereign debt crises are more persistent than banking ones. All figures are available from the authors upon request.

Results for shocks of magnitude 10 are available upon request.