extfieldsofstudy:
Computer Science
Economics
Mathematics

text:
Dynamic modeling of mean-reverting spreads for statistical arbitrage

Statistical arbitrage strategies, such as pairs trading and its generalizations, rely on the construction of mean-reverting spreads enjoying a certain degree of predictability. Gaussian linear state-space processes have recently been proposed as a model for such spreads under the assumption that the observed process is a noisy realization of some hidden states. Real-time estimation of the unobserved spread process can reveal temporary market inefficiencies which can then be exploited to generate excess returns. Building on previous work, we embrace the state-space framework for modeling spread processes and extend this methodology along three different directions. First, we introduce time-dependency in the model parameters, which allows for quick adaptation to changes in the data generating process. Second, we provide an on-line estimation algorithm that can be constantly run in real-time. Being computationally fast, the algorithm is particularly suitable for building aggressive trading strategies based on high-frequency data and may be used as a monitoring device for mean-reversion. Finally, our framework naturally provides informative uncertainty measures of all the estimated parameters. Experimental results based on Monte Carlo simulations and historical equity data are discussed, including a co-integration relationship involving two exchange-traded funds.

Introduction
A time series is known to exhibit mean reversion when, over a certain period of time, is "reverting" to a constant mean. In recent years, the notion of mean reversion has received a considerable amount of attention in the financial literature. For instance, there has been increasing interest in studying the long-run properties of stock prices, with particular attention being paid to investigate whether stock prices can be characterized as random walks or mean reverting processes. If a price time series evolves as a random walk, then any shock is permanent and there is no tendency for the price level to return to a constant mean over time; moreover, in the long run, the volatility of the process is expected to grow without bound, and the time series cannot be predicted based on historical observations. On the other hand, if a time series of stock prices follows a mean reverting process, investors may be able to forecast future returns by using past information. Since the seminal work of Fama and French [1988] and Poterba and Summers [1988], who first documented mean-reversion in stock market returns during a long time horizon, several studies have been carried out to detect mean reversion in several markets (e.g. Chaudhuri and Wu [2003]) and many asset classes (e.g. Deaton and Laroque [1992], Jorion and Sweeney [1996]).
Since future observations of a mean-reverting time series can potentially be forecasted using historical data, a number of studies have also examined the implications of mean reversion on portfolio allocation and asset management; see Barberis [2000] and Carcano et al. [2005] for recent works. Active asset allocation strategies based on mean-reverting portfolios, which generally fall under the umbrella of statistical arbitrage, have been utilized by investment banks and hedge funds, with varying degree of success, for several years. Possibly the simplest of such strategies consists of a portfolio of only two assets, as in pairs trading.
This trading approach consists in going long a certain asset while shorting another asset in such a way that the resulting portfolio has no net exposure to broad market moves. In this sense, the strategy is often described as market neutral. Entire monographs have been written to illustrate how pairs trading works, how it can be implemented in real settings, and how its performance has evolved in recent years (see, for instance, Vidyamurthy [2004] and Pole [2007]). The underlying assumption of pairs trading is that two financial instruments with similar characteristics must be priced more or less the same. Accordingly, the first step consists in finding two financial instruments whose prices, in the long term, are expected to be tied together by some common stochastic trend. What this implies is that, although the two time series of prices may not necessarily move in the same direction at all times, their spread (for instance, the simple price difference) will fluctuate around an equilibrium level.
Since the spread quantifies the degree of mispricing of one asset relative to the other one, these strategies are also refereed to as relative-value. If a common stochastic trend indeed exists between the two chosen assets, any temporary deviation from the assumed mean or equilibrium level is likely to correct itself over time. The predictability of this portfolio can then be exploited to generate excess returns: a trader, or an algorithmic trading system, would open a position every time a substantially large deviation from the equilibrium level is detected and would close the position when the spread has reverted back to the its mean.
This simple concept can be extended in several ways, for instance by replacing one of the two assets with an artificial one (e.g. a linear combination of asset prices), with the purpose of exploiting the same notions of relative-value pricing and mean-reversion, although in different ways; some relevant work along these lines has been documented, among others, by  and Montana and Parrella [2009], who describe statistical arbitrage strategies involving futures contracts and exchange-traded funds (ETFs), respectively. One aspect that has not been fully investigated in the studies above is how to explicitly model the resulting observed spread. A stochastic model describing how the spread evolves over time is highly desirable because it allows the analyst to precisely characterize and monitor some of its salient properties, such as mean-reversion. Moreover, improved trading rules may be built around specific properties of the adopted spread process.
Recently, Elliott et al. [2005] suggested that Gaussian linear state-space processes may be suitable for modeling mean-reverting spreads arising in pairs trading, and described how such models can yield statistical arbitrage strategies. Their main observation is that the observed process should be seen as a noisy realization of an underlying hidden process describing the true spread, which may capture the true market conditions; thus, a comparison of the estimated unobserved spread process with the observed one may lead to the discovery of temporary market inefficiencies. Based on the additional assumption that the model parameters do not vary over short periods of time, Elliott et al. [2005] suggested to use the EM algorithm, an iterative procedure for maximum likelihood estimation, for tracking the hidden process and estimating the other unknown model parameters. To make the exposition self-contained, we briefly review their model in Section 2, and state under what conditions the stochastic process is mean-reverting.
In this paper we build upon the model by Elliott et al. [2005] and extend their methodology in a number of ways. First, in Section 3, we introduce time-dependency in the model parameters. The main advantage of this formulation is a gain in flexibility, as the model is able to adapt quickly to changes in the data generating process; Section 3.1 further motivates our formulation and discusses its potential advantages. In Section 3.2 we derive new conditions that need to be satisfied for a model with time-varying parameters to be mean-reverting.
In Section 3.3 we describe a Bayesian framework for parameter estimation which then leads to a recursive parameter estimation procedure suitable for real-time applications. The final algorithm is detailed in Section 4; an analysis and discussion on the convergence properties of the algorithm as well as practical suggestions on how to specify the initial values and prior distributions are provided. Unlike the EM algorithm, our estimation procedure also produces uncertainty measures without any additional computational costs. With a view on statistical arbitrage, in Section 4.5 we add a note discussing how pairs trading may be implemented using the spread models proposed in this work and enumerate other important issues involved in realistic implementations, together with some pointers to the relevant literature. However, an empirical evalutation of trading strategies is beyond the scope of this work. For further discussions on statistical arbitrage approaches based on mean-reverting spreads and many illustrative numerical examples the reader is referred to Pole [2007].
In Section 5, based on a battery of Monte Carlo simulations, we demonstrate that posterior means estimated on-line by our Bayesian algorithm recovers the true model parameters and can be particularly advantageous when the analysts wishes to track sudden changes in the mean-level of the spread and its mean-reverting behavior. For instance, real-time monitoring may be used to derive stop-loss rules in algorithmic trading. Two examples involving real historical data are given in Section 5.2, where the cointegrating relationship between a pair of stocks and a pairs of ETFs are discussed. Final remarks are found in Section 6 and the proofs of arguments in Sections 3.2 and 4.3 can be found in the appendix.
2 Time-invariant state-space models for mean-reverting spreads Throughout the paper we will assume that the trader has identified two candidate financial instruments whose prices are observed at discrete time points t = 1, 2, . . . and are denoted by p (j) t , with j = 1, 2. At any given time t, let y t denote the price spread, defined as for some parameters α and β which are usually estimated by ordinary least squares (OLS) methods using historical data. It seems common practice to select the order of i and j such that y t yields the largest β and the resulting spread captures as much information as possible about the (linear) co-movement of the two assets. In Section 5.2 we briefly mention how a penalized OLS model may be used for recursive estimation of a time-varying β. More generally, the observed spread y t may also be obtained in different ways or may represent the return process of an initial price spread. One of the two component processes {p (j) t } may even be artificially built using a linear combination of a basket of assets. For our purposes, the only requirement is that the process {y t } is assumed to be mean-reverting. Furthermore, following Elliott et al. [2005], we assume that the observed spread y t is a noisy realization of a true but unobserved spread or state x t . The state process {x t } is defined such that where 0 < b < 2, a is an unrestricted real number and x 1 is the initial state. The restriction 0 < b < 2 is imposed, because otherwise {x t } is non-stationary and thus mean reversion has probability zero to occur. The innovation series {ε t } is taken to be an i.i.d. Gaussian process with zero mean and variance C 2 , and ε t+1 is assumed to be uncorrelated of x t , for t = 1, 2, . . .. Conditions for the state process to be mean-reverting are established using standard arguments, as follows. First, rewrite (1) as Expanding on this, we obtain Then, taking expectations and variances, It is observed that, when |1 − b| < 1, and regardless of a, lim t→∞ (1 − b) t−1 = 0 and therefore lim t→∞ E(x t ) = a/b. Therefore, in the long run, the state process fluctuates around its mean level a/b. Otherwise, when |1 − b| ≥ 1, (1 − b) t−1 is unbounded and hence E(x t ) is unbounded too. Analogously, when |1 − b| < 1 and regardless of a, the variance Var(x t ) with geometric speed. It is concluded that the hidden process {x t } is mean reverting when 1 − b lies inside the unit circle. Adopting the notation of Elliott et al. [2005], we define A = a and B = 1 − b, so that the process x t can be rewritten as Without loss of generality, we postulate that {y t } is a noisy version of {x t } generated as where {ω t } is Gaussian white noise with variance D 2 and ω t is uncorrelated of x t , for t = 1, 2, . . .. From (3), it also follows that {y t } is a mean-reverting process.
Note that, together with an initial distribution of the state x 1 , equations (2) and (3) define a Gaussian linear state-space model with parameters A, B, C, D. State-space models were originally developed by control engineers [Kalman, 1960] and are useful tools for expressing dynamic systems involving unobserved state variables. The reader is also referred to Harvey [1989] and West and Harrison [1997] for book-length expositions.
3 Time-varying dynamic models and on-line estimation

Preliminaries
The linear Gaussian state-space model described by equations (2) and (3) contains the unknown parameters A, B, C and D which need to be estimated using historical data. When the parameters are known, the Kalman filter provides a recursive procedure for estimating the state process x t [Kalman, 1960]. Full derivations of the Kalman filter and lucid explanations in a Bayesian framework can be found in Meinhold and Singpurwalla [1983]. In practice, maximum likelihood estimation (MLE) of the unknown parameters is required in order to fully specify the model. MLE for state-space models can be routinely carried out in a missing-data framework using the EM algorithm, as first proposed in the 1980s by Shumway and Stoffer [1982] and Ghosh [1989]; a detailed derivation can also be found in Ghahramani and Hinton [1996]. In the context of pairs trading, Elliott et al. [2005] reports some simulation and calibration studies demonstrating that the EM algorithm provides a consistent and robust estimation procedure for the model (2)-(3), and suggest that the finite-dimensional recursive filter described in Elliott and Krishnamurthy [1999] may also be used for estimation (although no results were provided). Elliott et al. [2005] suggest to base model estimation on data points belonging to a lookback window of size N . A full iterative calibration procedure is then run till convergence every time a new data point is observed and the window has been shifted one-step ahead.
This approach implicitly requires the analyst to select a value of N (the effective sample size) ensuring that, within each time window, the model parameters do not vary. The selection of N may be difficult without a proper model selection procedure in place to test the assumption that the model is locally appropriate. For instance, although a small value of N may guarantee adequacy of the model, it could also lead to notable biases in the parameter estimates.
When N is too large, a number of factors such as special market events, persisting pricing inefficiencies or structural price changes may invalidate the modeling assumptions. Clearly, the question of how much history to use to calibrate a model and the corresponding trading strategy is a critical one. In order to cope with these limitations, in this section we present and discuss our three main contributions. Firstly, we introduce more flexibility and release some of the modeling assumptions by allowing the model parameters to vary over time; in this way, both smooth and sudden changes in the data generating process (such as those created by structural price changes and unusual persistence of market inefficiencies) will be more easily accommodated.
Secondly, we propose a practical on-line estimation procedure that, being non-iterative, can be run efficiently over time, even at high sampling frequencies, and does not inflict the burden of frequent re-calibration and window size selection. Ideally, a model should be able to adapt to changes in the data generating mechanism with minimal user intervention, and should be amenable to on-line monitoring so that the key parameters characterizing the underlying mean-reverting property can always be under continuous scrutiny. These features enable the trader (or trading system) to take fast dynamic decisions. Thirdly, as a result of the Bayesian framework proposed here for recursive estimation, measures of uncertainty extracted from the full posterior distribution can be routinely computed at no extra cost. These measures can be very informative in quantifying and assessing estimation errors, and can potentially be exploited to derive more robust trading strategies; see Section 5 for some practical examples.

The proposed model
In this section we initially propose a variation of the classic state-space model used by Elliott et al. [2005] in which the parameters are not assumed to be constant over time. This modification will then force us to reconsider under which conditions the spread process is mean-reverting.
First, let us rearrange the model (2) and (3) in an autoregressive (AR) form. From (3), note that x t = y t − ω t . Then, from substitution in (2) for t ≥ 2, we obtain model is an AR model of order 1 with parameters A, B and σ 2 .
We achieve time-dependence in the parameters of (4) by replacing A and B with A t and B t , respectively, and postulating that both parameters evolve over time, according to some weakly stationary process. Here we consider the case of A t and B t changing over time via AR models, but more general time series may be considered. These choices lead to the specification of a time-varying AR model of order 1, or TVAR(1). Accordingly, the observed spread is described by the following law, where φ 1 and φ 2 are the AR coefficients, usually being assumed to lie inside the unit circle so that A t and B t may be weakly stationary processes.
with Φ = diag(φ 1 , φ 2 ) and error structure governed by the observation error ǫ t ∼ N (0, σ 2 ) and the evolution error vector ν t = (ν 1t , ν 2t ) ′ ∼ N 2 (0, σ 2 V t ), where N 2 (·, ·) denotes the bivariate Gaussian distribution. It is assumed that the innovation series {ǫ t } and {ν t } are individually and mutually uncorrelated and they are also uncorrelated of the initial state vector θ 1 , i.e. Theorem 1. If {y t } is generated from model (6)-(7), then, conditionally on a realized sequence B 1 , . . . , B t , {y t } is mean reverting if one of the two conditions apply: (b) φ 1 and φ 2 lie inside the unit circle, V t is bounded and |B t | < 1, for all t ≥ t 0 and for some integer t 0 > 0. Some comments are in order. First we note that if A t = A and B t = B (this is achieved by setting φ 1 = φ 2 = 1, and by forcing the covariance matrix of ν t to be zero for all t), the condition |B 1 | = |B| < 1 of Theorem 1 reduces to the known condition of mean reversion for the static AR model, as in the previous section. In the dynamic case, when A t and/or B t change over time, the condition |B t | < 1 enables us to check mean reversion in an on-line fashion. Following the approach of Elliott et al. [2005] for the AR model, we use model (6) in order to obtain estimatesB t of B t and then we check |B t | < 1 in order to declare whether {y t } is mean reverting or not; in the following sections we detail the computations involved in the estimation of B t . Structural changes in the level of y t are accounted through estimates of A t , but these do not affect the mean reversion of {y t } as A t controls only the level of y t .
For the case of A t = A, this is explained in some detail in Elliott et al. [2005] and for more information on structural changes for cointegrated systems the reader is referred to Johansen [1988] and Lütkepohl [2006, Chapter 6]. The following result is a useful corollary of Theorem 1.
The proof of this corollary follows by combining the proofs of (a) and (b) of Theorem 1 (see the appendix). Corollary 1 gives an important case, in which A t = A, for all t as in Elliott et al. [2005], but B t changes according to a weakly stationary AR model. This can be used when it is expected that A t will be approximately constant and benefit may be gained by reducing the tuning of the four parameters φ 1 , φ 2 , δ 1 , δ 2 to tuning of two parameters φ 2 , δ 2 .
For a further discussion on this topic see Sections 4.4 and 5.
In this paper we propose (5) as a flexible time-varying model for the observed spread.
However, more general time-varying autoregressive models may be used. Consider that y t is generated from a time-varying AR model of order d, i.e.
and the time-varying AR parameters A t and B it follow first order AR models, as where d is the order or lag of the autoregression, the innovations ǫ t and ν it are individually and mutually uncorrelated and they are uncorrelated with the initial states A d and B id . Certain Gaussian distributions may be assumed on ǫ t and ν it and on the states A t and B it . It is readily seen that this model can be casted in state space form (6) It is clear that model (5) is a special case of model (8) with d = 1. When the general model is adopted, the conditions of mean reversion of {y t } of Theorem 1 need to be revised, as follows. For φ i (i = 1, . . . , d + 1) being inside the unit circle, for t > t 0 , all B it x i must lie outside the unit circle. This effectively means that after some t 0 , {y t } is locally stationary [Dahlhaus, 1997]. For the remainder of this paper, we consider the situation of d = 1, i.e.
model (5), as this is a simple and parsimonious model.

A Bayesian framework
We adopt a Bayesian formulation that, within the realm of conjugate analysis, allows us to derive fast recursive estimation procedures and naturally compute measures of uncertainty. The analysis we propose in this section has roots in the work of West et al. [1999], Prado and Huerta [2002], and Triantafyllopoulos [2007a]. Initially, we assume that, given the observational variance σ 2 , the initial state θ 1 follows a bivariate Gaussian distribution with mean vector m 1 and covariance matrix σ 2 P 1 . Also, we place an inverted gamma density prior with parameters n 1 /2 and d 1 /2 on σ 2 . In summary, the prior structure is specified as follows θ 1 |σ 2 ∼ N 2 (m 1 , σ 2 P 1 ) and σ 2 ∼ IG(n 1 /2, d 1 /2), where m 1 , P 1 , n 1 , d 1 are assumed known; we comment on their specification in Section 4.4.
Note that, unconditionally of σ 2 , the initial state θ 1 follows a Student t distribution.
With these priors in place, the posterior distribution of θ t |σ 2 and the predictive distribution of y t |σ 2 are routinely obtained by the Kalman filter. We elaborate more on this as follows.
First, assume that at time t − 1 the posteriors are given by θ t−1 |σ 2 , y t−1 ∼ N 2 (m t−1 , σ 2 P t−1 ) and σ 2 |y t−1 ∼ IG(n t−1 /2, d t−1 /2), for some m t−1 , P t−1 , n t−1 and d t−1 . Here the notation y t means that all data points observed up to time t are included. Then, writing the likelihood function (or evidence) for an observation y t as p(y t |θ t , σ 2 ), an application of the Bayes theorem gives p(θ t |σ 2 , y t ) = p(y t |θ t , σ 2 )p(θ t |σ 2 , y t−1 ) p(y t |σ 2 , y t−1 ) , It follows that the posterior density of θ t |σ 2 is Gaussian, and specifically The recurrence equations for updating m t and P t are provided in Section 4. The probability density p(y t |σ 2 , y t−1 ) refers to the one-step ahead forecast density, which is obtained from the prior p(θ t |σ 2 , y t−1 ) as y t |σ 2 , y t−1 ∼ N (f t , σ 2 Q t ). Again, see Section 4 below for recursive equations needed to update f t and Q t .
which gives an inverted gamma density σ 2 |y t ∼ IG(n t /2, d t /2), depending on parameters n t and d t . Here y t |y t−1 follows a t distribution with n t−1 degrees of freedom y t |y t−1 ∼ From the density p(θ t |σ 2 , y t ), the posterior distribution of θ t , unconditionally of σ 2 , is easily obtained by integrating σ 2 out; it It then follows that θ t |y t ∼ t 2 (n t , m t , P t S t ). From where m t = (m 1t , m 2t ) ′ , P t = (P ij,t ) ij=1,2 and t γ denotes the 100γ% quantile of the standard t distribution with n t degrees of freedom. The (1 − γ)% confidence interval for x t+1 is where the recurrence relationships of R t+1 and Q t+1 are given below.
Some references on related time series models are in order. From a frequentist perspective, time varying AR models have been discussed in Dahlhaus [1997], Francq and Zakoan [2001], Francq and Gautier [2004] and Anderson and Meerschaert [2005]. Among other works, recursive estimation of time varying autoregressive processes in a nonparametric setting is discussed in Moulines et al. [2005] and, for non Gaussian processes, in Djurić et al. [2002], using particle filters. Standard Bayesian AR models have been developed since the early 70's, see e.g. Zellner [1972], Monahan [1983], Kadiyala and Karlsson [1997] and Ni and Sun [2003]. Free software for model estimation is widely available 1 .
4 On-line estimation 4.1 An adaptive and recursive algorithm using discount factors In this section we provide the updating equations needed to compute the posterior densities of θ t |y t and of σ 2 |y t at each time step. Starting at time t = 1 with a quadruple of initial values (m 1 , P 1 , n 1 , d 1 ), the calibration algorithm then proceeds as follows: For any t = 2, . . . , T , the above algorithm estimates the target posterior quantities of interest; for instance, we can extract posterior and predictive mean and variances, as well as relevant quantiles and credible bounds of θ t and σ 2 . From θ t = (A t , B t ) ′ and the posterior distribution of θ t |y t , we can extract the posterior distribution of B t |y t . The condition for mean-reversion established in Theorem 1 can be monitored recursively by extracting the posterior mean of B t |y t , sayB t , and assessing whether |B t | is strictly less than one. Credible bounds can also be associated to the posterior mean in order to better assess the possibility that the process is still mean-reverting -see the examples in Section 5.
The full specification of algorithm (10) requires the selection of a covariance matrix V t , which is responsible for the stochastic evolution of the signal θ t and hence the stochastic change of A t and B t . Following West and Harrison [1997, Chapter 6] we advocate a practical and convenient analytical solution which allows us to learn this variance component directly from the data in a sequential way by means of two discount factors, δ 1 and δ 2 ; this is referred to as component discounting. The idea is that by assuming P 1 and V t to be diagonal matrices we can use the two discount factors to discount the precision of the updating of the mean and the variance of θ t as we move from time t − 1 to t. In other words we use δ 1 and δ 2 to specify the covariance matrix V t as This implies that R t = diag(φ 2 1 p 11,t−1 /δ 1 , φ 2 2 p 22,t−1 /δ 2 ) and thus, as we move from t − 1 to t, the prior variance of A t is increased by a factor of 1/δ 1 and of B t by a factor of 1/δ 2 . Of course if δ 1 = δ 2 = 1, then V t = 0 and in this case θ t carries no stochastic evolution. If we allow δ 1 = 1 and δ 2 < 1, then only B t has stochastic evolution over time.

Model comparison and model assessment
The performance of the estimation procedure of Sections 3.3 and 4 can be formally evaluated using model diagnostic and model comparison tools; see, for instance, Li [2004] for a general exposition of time series diagnostics and Harrison and West [1991] for diagnostics in state space models. In this section we briefly discuss three diagnostic tools, namely the mean of the squared standardized forecast errors (MSSE), the likelihood function, and sequential Bayes factors.
From the Student t distribution of y t |y t−1 , i.e. y t |y t−1 ∼ t(n t−1 , f t , Q t S t−1 ), we can define the standardized one-step forecast errors (or standardized residuals) as u t = Q −1 t S −1 t−1 (y t −f t ), so that u t |y t−1 ∼ t(n t−1 , 0, 1) (the standard t distribution with n t−1 degrees of freedom). We can therefore construct diagnostics and outlier detection tools based on the above t distribution of u t . Writing v t = (1 − 2n −1 t−1 )u t we have E(v 2 t |y t−1 ) = 1 and so the MSSE is defined as (T − 1) −1 T t=2 v 2 t , which if the model fit is good, should be close to 1. From the Student t distribution of y t |y t−1 the log-likelihood function of φ 1 , φ 2 , δ 1 , δ 2 based on data y T = {y 2 , . . . , y T } is where Γ(.) denotes the gamma function. Model camparison can be carried out by using either one of the following criteria: likelihood function, Akaike's information criterion (AIC) and Bayesian information criterion (BIC) . In particular, we can choose optimal values of some or all of the hyperparameters φ 1 , φ 2 , δ 1 , δ 2 by maximizing ℓ(.). A discussion on the specification of the hyperparameters of the model can be found in Section 4.4.
For the application of the above diagnostic criteria, all data y T is needed to be available, or historical data can be used. However, sometimes it is useful to construct sequential diagnostics so that the model can be assessed and updated over time in an adaptive way. Such diagnostics tools include sequential likelihood ratios and sequential Bayes factors. Here we briefly discuss the latter, the foundations of which are discussed in detail in West and Harrison [1997, Chapter 11]. Suppose that, given a sample y T = {y 1 , . . . , y T } we have two candidate models of the form of (6) that is they have the same structural form, but they may differ in the values of φ 1 , φ 2 , δ 1 and δ 2 . Suppose that we denote the two models by M 1 and M 2 and for i = 1, 2 we write φ i1 , φ i2 , δ i1 and δ i2 to indicate the dependence of model M i in these parameters. Then the Bayes factor of M 1 versus M 2 is given by the ratio of their respective one-step forecast densities, i.e.
where we have used that y t |y t−1 , M i ∼ t(n t−1 , f it , Q it S i,t−1 ), with the quantities f it , e it , Q it , S i,t−1 being appropriately indexed by i = 1, 2. Given data y T one can either judge the performance of the two models sequentially (by comparing H t to 1, for 2 ≤ t ≤ T ) and thus arriving to a sequential monitoring of the two models, or use the entire data set y T to compare the models globally, e.g. one can extract the mean or other features of the empirical distribution of {H t }.

Convergence analysis
Algorithm (10) is quite similar to the celebrated Kalman filter; conditional on σ 2 , the algorithm exactly reduces to the Kalman filter, but the full algorithm allows for the estimation of σ 2 that results in the Student t posterior distribution for θ t . On the performance of the However, the performance of the posterior covariance matrix P t when the components of the model are made time-dependent has not been investigated; in our system this is conveyed via the time-varying vector F t = (1, y t−1 ) ′ . This aspect is important as instability or divergence of P t could result in instability of the estimation of A t and B t and hence of x t . The next result states that, in our system, P t converges to stable values and we provide an explicit formula for the computation of the limit of P t .
Theorem 2. Suppose that {y t } is generated from model (6). If {y t } is mean reverting and if, for j = 1, 2, it is δ j < φ 2 j , then as t → ∞ the limit P of the covariance matrix P t = Var(θ t |y t ) exists and it is given by P = diag(p 11 , p 22 ), where with a 1,t = 1 and a 2,t = y 2 t−1 .
Some comments are in order. First we note that if φ 1 = φ 2 = 1 and V t = 0 (we have already seen that this setting reduces the model to the time-invariant AR model considered in Elliott et al. [2005]), then the condition δ j < φ 2 j is satisfied for all values of δ j , since 0 < δ j < 1.
From the mean reversion assumption of {y t }, if we write y t ≈ µ, where µ denotes the equilibrium mean of the spread, then we can write the limit covariance matrix P as In the important special case of φ 1 = δ 1 = 1, for which A t = A is time-invariant, we can easily see that where p 11,1 is the prior variance Var(A).
The convergence rate of the limit of Theorem 2 is geometric, since after some appropriately large t L , we can write y t ≈ µ, for all t > t L and the limit of P depends on a geometric series.
The above convergence results for P t are given conditional on the variance σ 2 . Given data up to time t, σ 2 has a posterior inverted gamma distribution σ 2 |y t ∼ IG(n t /2, d t /2); hence, as the time index gets larger, the variance of σ 2 , which is given by converges to 0. Therefore, as t → ∞, σ 2 concentrates about its mode S t = d t /n t asymptotically degenerating.

Hyperparameter specification
The estimation algorithm (10) relies upon the specification of prior distributions and corresponding starting values (m 1 , P 1 , n 1 , d 1 ) and values of the model components (φ 1 , φ 2 , δ 1 , δ 2 ), which are selected by the user. In this brief section, considering weakly informative priors, we provide some guidance on how to choose these values. Of course, depending on the specific application, other specifications may be preferred; for instance, the analyst may want to include stronger prior beliefs regarding the spread being traded, see e.g. Kadane et al. [1996].
Nevertheless, it is important to note that, given a reasonable amount of data, the sensitivity of the calibration procedure on these initial specifications becomes negligible, especially over streaming data, because the initial information is deflated over time. This phenomenon is discussed in some detail in Ameen and Harrison [1984] and in Triantafyllopoulos [2007a].
Detailed studies on prior specification for the estimation of AR models can be found in Kadiyala and Karlsson [1997], Ni and Sun [2003] and in references therein.
The parameter m 1 is the prior mean of the hidden state, given the observational variance, i.e. the mean of θ 1 |σ 2 . A common choice is to set m 1 equal to our prior expectation of (A 1 , B 1 ) ′ , which may be obtained from the availability of historical data. In all examples of Section 5 we have used m 1 = (0, 0) ′ . This setting together with the vague prior P 1 that follows, communicates a prior assumption of mean reversion, but with a large uncertainty placed a priori on (A 1 , B 1 ) ′ . The convergence results reported in Theorem 2 above, guarantee that the choice of m 1 and P 1 are not crucial for accurate estimation and forecasting. The covariance matrix P 1 is chosen to be proportional to the 2 × 2 identity matrix, i.e. P 1 = p 1 I 2 . Here a large value of p 1 reflects a weakly informative or defuse prior specification, since in this case the precision P −1 1 gets close to zero. Finally, values for n 1 and d 1 need to be provided. It can be noted that, having placed an inverted gamma prior on σ 2 , the expected value of the observational variance is given by E(σ 2 ) = d 1 /(n 1 − 2), for n 1 > 2. Based on this observation, a sensible choice is to set n 1 = 3 and use the prior expectation of σ 2 as a starting value d 1 .
Historical data may be used to specify d 1 , but in the examples of Section 5 we have simply used d 1 = 1.
Proceeding now with the specification of φ 1 , φ 2 , δ 1 , δ 2 we can optimize these parameters by maximizing the log-likelihood function, given in Section 4.2, under the condition that δ i < φ 2 i so that Theorem 2 applies. Alternatively, according to Corollary 1 we can set φ 1 = δ 1 = 1 and optimize only φ 2 and δ 2 . In Section 5.1, where we present simulation studies, we use the latter, while in Section 5.2, where we analyze real data, we use the former (full optimization of four parameters). We note that the likelihood function or Bayes factors can be used to compare and optimize models using single discount factors δ 1 = δ 2 , known as single discounting [West and Harrison, 1997], and models using two different discount factors (component or multiple discounting).

Pairs trading
Under the assumption that the observed spread process involving two tradable assets is meanreverting, and that the model of Eqs. (2)-(3) describes well its evolution at discrete observational times t = t 1 , . . . , t N , with N sufficiently small, a simple pairs trading strategy immediately follows [Elliott et al., 2005]. Let us assume thatx t denotes our best estimate of the hidden state, which is obtained by calibrating the model on data collected in the above data window.
At each time t, if the observed spread y t is strictly greater than the true statex t , then a sensible decision would be to take a long position in this portfolio, with the intention of closing this position at a later time, when the spread has reverted back to its mean.
Conversely, if y t <x t , the trader may decide to take a short position in the portfolio; this bet is expected to be a profitable one as soon as the spread process corrects itself again. Realistic implementations of this popular strategy may ask for additional layers of sophistication which in turn require the trader to face a few practical questions; some examples are: • How can transaction costs be included in this simple model? In other words, when is a trade expected to be profitable, so that an 'entry' signal can be generated? For instance, a long position could be initiated when y t −x t > z t , where z t is a threshold that guarantees a profitable trade, after costs. The question then becomes, how should z t be calibrated? For instance, Vidyamurthy [2004] suggests a re-sampling procedure and provides some general guidance. There may exist several other alternative ways in which one could define entry points, perhaps based on empirical modeling of the extreme values of the y t − x t process. Theoretical results on zero-crossing rates for autoregressive processes, as in Cheng et al. [1995], may also be explored. For aggressive strategies that execute a trade at each single time tick,  forecasts the one-step ahead expected spread using dynamic regression methods, whereas Montana and Parrella [2008] embrace the principle of "learning with experts" to deal with the uncertaincty involved in future movements on the spread.
• Analogously to the previous issue, how should an 'exit' signal be generated? And shall a trade be closed at an exit point, or simply reversed so that a long position becomes short, and viceversa?
• What stop-loss mechanism can be implemented to make sure that the assumptions on which the strategy relies are still satisfied? Surely, if the spread process is no longer believed to be mean-reverting, a stop-loss signal should be quickly generated. As will appear clearer later (see, for instance, the examples of Section 5), our estimation procedure can be used to monitor mean-reversion sequentially and flag deviations from the acceptable behaviour of the spread process as soon as they occur. Related co-integration arguments may also be used, as in Lin et al. [2006].
• How can suitable pairs of assets be chosen in the first place, especially when the universe of assets to search from is extremely large? Since arbitrage profits between two assets depend critically on the presence of a long-term equilibrium between them (see, for instance, Alexander et al. [2002]), data mining methods built around co-integration techniques may be explored, as in d'Aspremont [2008]. See also Vidyamurthy [2004] and Pole [2007] for alternative methods including simple correlation analysis, turning point analysis and latent factor models.
As a final note, we mention a technique that may be deployed in a dynamic modeling setting, such as ours, to obtain the spread y t = p t in a recursive fashion. As noted before, the regression coefficient is usually estimated on historical data, but on-line procedures such as recursive least squares may also be used. The assumption of a time-invariant regression coefficient β could also be released so as to allow β to change slightly over time; such a modification would capture a time-varying co-integration relationship between the two asset prices, where this extension deemed necessary. Assuming T historical observations, a regression model with a time-varying regression coefficient β t minimizes a cost function where µ ≥ 0 is a scalar determining how much penalization to place on temporal changes in the regression coefficient. When µ is very large, changes in the coefficient are penalized more heavily and, in the limit µ = ∞, the usual OLS estimate is recovered. A solution to the optimization problem above was originally proposed by Kalaba and Tesfatsion [1988].
Following their approach, called flexible least squares (FLS), a recursive estimator for each β t can easily be derived asβ where we have defined the quantities The recursions are initially started with some arbitrarily chosen values S 1 and s 1 . Montana et al.
[2009] show a clear algebraic connection between FLS and the Kalman filter and use this estimation method to develop a dynamic statistical arbitrage strategy.

Simulated data
In this section we initially report on a Monte Carlo simulation study demonstrating that the fast recursive algorithm (10) described in Section 3.3 accurately estimates the parameters of the proposed model. We have simulated a large number of time series under model (6) using a range of values for A, B and σ 2 . The true parameters are kept constant in these initial simulations for simplicity, so they can be easily compared with the estimated posterior means. We have found that convergence to the true parameters A, B and σ 2 is quickly achieved and the estimated values of these parameters are not sensitive to the initial parameters µ 1 , P 1 , n 1 , d 1 (results not shown).
We have also explored situations in which the parameters are time-varying. First, we have considered the case of a sudden change in the level of the spread; the time series fluctuates around an equilibrium level till t = 1500, and after that time it jumps to a much higher equilibrium. Clearly for 1 ≤ t ≤ 1499 the process is mean reverting, then at t = 1500 it looses mean reversion, but it retains it in the sub-period 1500 ≤ t ≤ 3000; of course the process is not mean reverted for the entire period 1 ≤ t ≤ 3000. Figure 1 shows how the posterior mean of |B t | is tracked using two different values of the discount factor δ 2 . Our focus is on monitoring B t because, as established in Theorem 1, this parameter is the ultimate object of interest. As shown in Figure 1, the algorithm with δ 2 = 1 (which corresponds to a model with time-invariant parameters) does not manage to capture the loss of mean-reversion observed at time t = 1500; in fact the algorithm gives the misleading result of mean reversion throughout the time range. On the contrary, when using a smaller discount factor (which corresponds to a model with time-varying parameters), the algorithm tracks the jump almost in real-time and communicates the result that after t = 1500 the process has locally regained mean reversion.
Furthermore, we have considered a more hypothetical scenario that may be of practical interest: Figure 2 corresponds to a scenario where B t is piece-wise constant and undergoes a large sudden jump at time t = 1500. Again, the algorithm is able to track well mean reversion locally, although the true parameter B t may not be estimated very accurately.

Equity data
In this section we apply our methods to spreads obtained from historical equity data. Each spread is computed using the flexible least squares (FLS) method with a very large µ parameter; this is almost equivalent to ordinary least squares (OLS) regression but allows for recursive estimation. As a simple validation exercise, we also compare the findings obtained from our model to formal cointegration tests which assume the availability of all data points. The very first procedure for the estimation of cointegrating regressions, based on OLS, was proposed by Engle and Granger [1987]. Since then several other procedures have been developed including the maximum likelihood method of Johansen [1988Johansen [ , 1991 and the fully modified OLS of Phillips and Hansen [1990]. Hargreaves [1994] lists eleven categories of procedures, and several more have been added in more recent years. For our analysis we have considere only three popular tests: Engle-Granger's ADF test [Engle and Granger, 1987], Phillips-Perron's PP test [Perron, 1988] and Phillips-Ouliaris's PO test [Phillips and Ouliaris, 1990].
The first data sets we present consists of daily share prices of two companies: Exxon Mobil (XOM) and Southwest Airlines (LUV). We have used all the available data for this pair of stocks, which spans a period from March 23, 1980 to August 6, 2008. Figure  When using all historical data , all three standard cointegration tests cannot reject the null hypothesis of unit roots (p-values: 0.246, 0.219 and 0.15). This is in agreement with the patterns captured by Figure 3, which reveals that after GLD is an ETF that tries to reflect the performance of the price of gold bullion, whereas GDX tries to replicate as closely as possible, before fees and expenses, the price and yield performance of the AMEX Gold Miners index. This is achieved by investing in all of the securities which comprise the index (in proportions given by their weighting in the index).
This analysis is based upon all the historical data available for the pair, which covers a shorted period compared to the previous example, from May 23, 2006until August 06, 2008. Figure   5 shows the observed spread process jointly with the estimated hidden process and confidence bands, while Figure 6 indicates that a co-integrating relationship between the two ETFs does exist in the period from July 19, 2006till 17 December, 2007. For this data set we have used φ 1 = 0.999, φ 2 = 99, δ 1 = 0.95 and δ 2 = 0.98 that maximize the log-likelihood function (11).
When all the historical data is used, the ADF and the PP tests indicate the presence of cointegration at a 5% significance level (p-values: 0.01 and 0.01, respectively) and only the PO test suggest lack of co-integration, a result that also agrees with the pattern reported in Figure   6. Considering the period July 19, 2006 till 17 December, 2007, for which our results suggest mean reversion, we find that all three tests also suggest co-integration (p-values: 0.0201, 0.013 and 0.012). Further formal comparisons and more detailed studies will be needed in order to characterize some of the discrepancies; however, based on this empirical evidence, our suggested time-varying model seems to generally agree with most formal cointegration tests.

Conclusions
In this paper we have proposed a Bayesian time-varying autoregressive model, expressed in time-space form, and an efficient recursive algorithm based on forgetting or discount factors.
The procedure can be used for real-time estimation and tracking of the underlying spread process and may be seen as a more efficient alternative to standard iterative MLE procedures such as the EM algorithm. Conditions for mean-reversion as well as the convergence properties of the on-line estimation algorithm have been studied analytically and discussed. The model seems particularly useful for monitoring mean-reversion using financial data streams and as a building block for statistical arbitrage strategies such as pairs trading. Related algorithmic trading strategies that exploit co-integration of financial instruments, for instance index arbitrage [Sutcliffe and Board, 2006] and enhanced index tracking [Alexander et al., 2002], may also benefit from the methods proposed here. Moreover, although the focus of this work has been on applications in computational finance, we believe that the methods described here are of broader interest and may appeal to other users, within the management science community, who need to model and monitor mean-reverting time series arising in different application domains There are several aspects of the suggested methodology that we would like to explore further in future work. First, purely from an empirical point of view, we would like to better understand how the methodology relates to more formal statistical procedures for testing the hypothesis of mean reversion based on finite sample sizes. As already mentioned, since mean-reversion is closely linked to second order stationarity, many efforts have been directed to constructing unit root tests. These standard econometric procedures may lack the power to reject the null hypothesis of a random walk, and we feel that our method may at least complement them well. Besides, some of the recently suggested procedures, such as the bootstrap methods described by Li and Xiao [2003], are too computationally expensive to be of any use in the real settings and applications that we have described. Another important aspect that we plan to investigate is the question of how to learn the discounting factors needed to specify the V t matrix in a more adaptive fashion, so that they become self-tuning, rather than being kept constant at all times. A number of techniques have been successfully used for training adaptive artificial neural networks and other time-varying stochastic processes using forgetting factors [Saad, 1999, Niedźwiecki, 2000 and there may be scope for improvement along this direction.
We write A t = (A 1 , . . . , A t ) and B t = (B 1 , . . . , B t ), for t = 1, . . . , T . Since {ǫ t } is white noise, we have This is a convergent series if |B t | < 1, for all t > t 0 , for some positive integer t 0 . To see this first write x (1) t = t i=2 B i , which is a decreasing series as |x t } is convergent. For the variance of y t we have where it is used that and Cov(A t , A t−j−1 ) ≤ σ 2 V 11 1 − φ 2 1 for V 11,t ≤ V 11 , since from the hypothesis V t is bounded, and so there exists some V 11 > 0 so that V 11,t ≤ V 11 . Now we show that the series x (2) 1 j i=0 B t−i are both convergent. For the former series we note that given |B t | < 1, we can find some B so that |B t | < |B| < 1, from which it follows that which is proportional to a geometric series that converges for |B| < 1 and since x is convergent as t−3 j=0 |φ 1 B| j+1 is a geometric series with |φ 1 B| < 1 and x (3) t is a positive series. With these convergence results in place, the convergence of Var(y t |B t ) is obvious. Given, B t , we have shown that the mean and the variance of {y t } are convergent and so {y t } is mean reverting.
Proof of Theorem 2. From the diagonal structure of P t = diag(p 11,t , p 22,t ) and the updating of P t as in the calibration algorithm (10) we have i p ii,t−1 We can clearly see that p ii,t > 0, for all t and so we have Now since δ j i φ −2j i a i,t−j is a positive sequence and since from the mean reversion of {y t } and the definition of a it , the above sequence is bounded above by the geometric sequence δ j i φ −2j where M is an upper bound of {y t }, it follows immediately that ∞ j=0 δ j i φ −2j a i,t−j < ∞ and this proves that p ii,t converges to ∞ j=0 δ j i φ −2j a i,t−j . The proof is completed by inverting (A-2), noting that p ii,t > 0 and ∞ j=0 δ j i φ −2j a i,t−j > 0.


Dependence of defaults and recoveries in structural credit risk models

The current research on credit risk is primarily focused on modeling default probabilities. Recovery rates are often treated as an afterthought; they are modeled independently, in many cases they are even assumed constant. This is despite of their pronounced effect on the tail of the loss distribution. Here, we take a step back, historically, and start again from the Merton model, where defaults and recoveries are both determined by an underlying process. Hence, they are intrinsically connected. For the diffusion process, we can derive the functional relation between expected recovery rate and default probability. This relation depends on a single parameter only. In Monte Carlo simulations we find that the same functional dependence also holds for jump-diffusion and GARCH processes. We discuss how to incorporate this structural recovery rate into reduced form models, in order to restore essential structural information which is usually neglected in the reduced-form approach.

Introduction
An accurate description of portfolio credit risk is of vital interest for any financial institution. It is also a prerequisite for realistic ratings of structured credit derivative products. Furthermore, it is a crucial aspect in banking regulations. We can distinguish two conceptually different approaches to credit risk modelling: structural and reduced-form approaches. The structural models go back to Black and Scholes (1973) and Merton (1974). The Merton model assumes that a company has a certain amount of zero-coupon debt which becomes due at a fixed maturity date. The market value of the company is modelled by a stochastic process. A possible default and the associated recovery rate are determined directly from this market value at maturity. In the reduced-form approach default probabilities and recovery rates are described independently by stochastic models. The aim is to describe the dependence of these quantities on common (macroeconomic) covariates or risk factors. For some well known reduced-form model approaches see, e.g., Jarrow and Turnbull (1995), Jarrow et al. (1997), Duffie and Singleton (1999), Hull and White (2000) and Schönbucher (2003). First Passage Models constitute a third approach -which is usually regarded as structural, but is better described as a mixed or pseudo-structural approach. They were first introduced by Black and Cox (1976). As in the Merton model, the market value of a company is modelled as a stochastic process. Default occurs as soon as the market value falls below a certain threshold. In contrast to the Merton model, default can occur at any time. In this approach, the recovery rate is not determined by the underlying process for the market value. Instead recovery rates are modelled independently, for example, by a reduced-form approach (see e.g., Asvanunt and Staal (2009a,b)). In some cases recovery rates are even assumed constant, for instance, in Giesecke (2004). The independent modelling of default and recovery rates leads to a serious underestimation of large losses.
In this paper we take a step back, historically, and revisit the Merton model, where defaults and recoveries are both determined by an underlying process. Hence, they are intrinsically connected. For a correlated diffusion process the Merton model has been treated analytically, e.g., in Bluhm et al. (2002) and in Giesecke (2004). In a straightforward calculation we can also derive the functional relation between expected recovery rate and default probability. This relation depends on a single parameter only. In Monte Carlo simulations we find that the same functional dependence also holds for other processes like jump-diffusion and GARCH. We discuss how to incorporate this relation into reduced-form models, in order to restore essential structural information which is usually neglected in the reduced-form approach.
The paper is organized as follows. We give a short introduction to the Merton model in Section 1. In Section 2 we treat the diffusion case analytically and compare the results to Monte-Carlo simulations. In Sections 3 and 4 we extend the Merton model to jumpdiffusion and GARCH processes. We discuss the applicability of the structural recovery rate beyond the Merton model in Section 6. We summarize our findings in Section 6.

Merton model
The Merton model assumes that a company k has a certain amount of zero-coupon debt; this debt has the face value F k and will become due at maturity time T . The company defaults if the value of its assets at time T is less than the face value, i.e., if V k (T ) < F k . The recovery rate then reads R k = V k (T )/F k and the loss given default is We denote the loss given default with an asterisk to distinguish it from the loss including non-default events. The individual loss can be expressed as where Θ is the Heaviside function. In the Merton model, defaults and losses -and hence also recoveries -are directly determined by the asset value at maturity. Therefore, the stochastic modelling of the market value V k (t) of a company allows to assess its credit risk. Let p V k (V k (T )) be the probability density function (pdf) of the market value at maturity. Then the default probability is given by and the expected recovery rate can be calculated as Let us now consider a portfolio of K credit contracts, where the market value of each company k is correlated to one or more covariates. Conditioned on the values of the covariates we obtain different values for P D,k and R k . In fact, we find a functional dependence between default probability and recovery rate. This is in stark contrast to modelling approaches which assume an independence of these quantities. In section 2 we derive this functional dependence analytically for the diffusion case. We compare the result with Monte Carlo simulations. In section 3 we extend the model by including jump terms in the random process. In section 4 we consider a GARCH process. In both cases we find the same functional dependence between default and recovery rates.

A. Analytical discussion
In the diffusion case we can easily derive all results analytically. To keep the notation simple we consider a homogeneous portfolio of size K with the same parameters for each asset process, and with the same face value, F k = F , and initial market value, V k (0) = V 0 . We model the time evolution of the market value of a single company k by a stochastic differential equation of the form This is a correlated diffusion process with a deterministic term µdt and a linearly correlated diffusion. The parameters of this process are the drift constant µ, the volatility σ and the correlation coefficient c. The Wiener processes dW k and dW m describe the idiosyncratic and the market fluctuations, respectively. For discrete time increments ∆t = T /N , where the time to maturity T is divided into N steps, we arrive at the discrete formulation of Eq. (5). The market value of company k at maturity can be written as We define the market return X m as the average return of all single companies k over the time horizon up to maturity, For K → ∞ we can express the average over k as the expectation value for ε k,t . Due to the independence of ε k,t for different k and t, we can write with ε k,t = 0 .
Thus, expression (9) simplifies to In the last step of the calculation we expanded the logarithm up to first order in ∆t. The random variables η m,t are standard normal distributed. Therefore the variable is normal distributed with mean µT − 1 2 cσ 2 T and variance cσ 2 T . This implies a shifted log-normal distribution for the market return itself, For a single company k we can write Conditioned on a fixed value for the market return X m , all variables V k (T ) are independent and ln V k (T )/V 0 is normal distributed with mean ln (X m + 1) − 1 2 (1 − c)σ 2 T and variance (1 − c)σ 2 T . Since we consider a homogeneous portfolio, we will omit the index k in the following. This allows for a clearer notation. The probability density function for the market value V (T ) is given by the following log-normal distribution We obtain the individual default probability by integrating this pdf from 0 up to the face value F , where Φ is the cumulative standard normal distrubution. The expectation value for the individual loss given default L * = 1 − V (T )/F can be calculated as The expected recovery rate is simply And the portfolio loss for a homogeneous portfolio, or simply the average loss, is obtained as For the sake of clarity, we introduce the function and the compound parameter Now, the expressions for P D (X m ) and R(X m ) simplify to and The relation between default probability and expected recovery rate depends only on B and is parametrized by A(X m ). Thus, the parameter B can be calibrated to default and recovery rate data. Furthermore, by inverting expression (23), we can express A in terms of P D , This leads us to the functional dependence of recovery rate and default probability, It is worth repeating that this functional relation depends on a single parameter B only.
In Figure 1 we demonstrate this parameter dependence of Equation (26). We plot R(P D ) for a set of parameter values. Larger values for B lead to an overall decrease of recovery rates. In addition the dependence on P D becomes steeper. From Equation (26) we obtain the functional relation of portfolio losses and default probabilities, The parameter dependence of this relation is also shown in Figure 1. In order to calculate Value at Risk and Expected Tail Loss we make use of the substitution Here we work with a simplified notation where L refers to L(X m ) . With the cumulative distribution function we can express the Value at Risk VaR α as where the function L refers to the X m -dependence given in Equation (20) and the value F −1 Xm (1−α) is the (1−α)-quantile of the market return. The Expected Tail Loss is calculated as We will now compare these analytic results with Monte Carlo simulations for a finite portfolio size.

B. Monte Carlo results
In the Monte Carlo simulation we consider the stochastic process in Equation (6) for discrete time increments. The simulation is run with an inner loop and an outer loop. In the inner loop we simulate K = 500 different realizations of ε k,t for a single realization of the market fluctuations η m,t with t = 1, . . . , N . The inner loop can be interpreted as a homogeneous portfolio of size K, or simply as an average over the idiosyncratic part of the process. In each run of the inner loop, we calculate the market return X m , the number of defaults N D (X m ) and the expected recovery rate R(X m ) . The market return X m is defined as the average return at maturity, For sufficiently large K the idiosyncratic terms average out and the market return X m is solely defined by the realization of η m,t . This is the reason why we use the market return as a parameter for the other observables. The number of defaults N D (X m ) simply counts how many times the condition V k (T ) < F is fulfilled. We can estimate the default probability as We obtain the portfolio loss as the average of individual losses in Equation (2), Using the relation we can estimate the expected recovery rate as Here, we assume that the number of defaults is strictly non-zero, which is justified for large portfolio size K. The outer loop runs over 10 6 realizations of the market terms, where we obtain different values for the market return X m and, consequently, the number of defaults N D (X m ), the default probability P D (X m ) and the expected recovery rate R(X m ) .
In this paper we do not discuss the parameter dependence of the models or aim at calibrating them to a given portfolio. Instead we only present the results for a single set of parameters with economically sensible values. As correlation coefficient we choose c = 0.5. The parameters for the diffusion process are µ = 0.05 and σ = 0.15. The initial market value is set to V 0 = 100, the face value of the zero-coupon bonds is F = 75 with maturity time T = 1. Figure 2 shows the dependence of recovery rates R(X m ) on default probabilities P D (X m ) in one plot, and the dependence of portfolio losses L(X m ) on default probabilities P D (X m ) in a second plot. In both cases we observe a very good agreement between the Monte Carlo simulations and the analytical results in Equations (26) and (27), respectively. Deviations from the average values are smaller for higher default probabilities. Low default probabilities correspond to positive or small negative market returns X m , see Figure 3. In this case defaults and recoveries are mostly influenced by the idiosyncratic part of the process and the correlation to the market plays only a minor role. This is why we observe a much broader range of recovery rates for low default probabilities.
When only limited historical data on defaults and recoveries is available, it might appear reasonable to assume a constant recovery rate. This corresponds to a linear dependence of the portfolio loss on default probability. For small default probabilities this is a good approximation of the results shown in Figure 2. However, for larger default probabilities a constant recovery rate model severly underestimates portfolio losses. The combination of high default probabilities and low recovery rates strongy influences the tail of the loss distribution. In Figure 3 we demonstrate that the dependence of default probabilities and portfolio losses on the market return is well described by Equations (17) and (20), respectively. Given the functional relation between portfolio loss L and market return X m we can transform the pdf of market returns into the pdf of portfolio losses, A comparison to the Monte Carlo result for the loss distribution is given in Figure 4. We observe an excellent agreement, even for extremely large portfolio losses. Finally, we present the results for Expected Loss, Value at Risk and Expected Tail Loss in Table B. The latter two are calculated for the 0.99-quantile. The analytical values describe the simulation results very well. Since the diffusion process was used both in the Monte Carlo simulations and in the analytical derivation, the good agreement of the respective results comes as no surprise. In the following sections we will explore two different processes in the Monte Carlo simulations and compare these with the analytical results for the diffusion.

A. Analytical discussion
We extend the diffusion model of the previous section by adding two jump terms to the stochastic process, dJ k for idiosyncratic jumps and dJ m for jumps which affect the entire market. The jump terms are not contained in Merton's original model; they ensure    that the default probability does not vanish as the time to maturity becomes very short. For applications of jump-diffusion processes in credit risk modelling see, e.g., Zhou (2001), Schäfer et al. (2007) and Kiesel and Scherer (2007). The stochastic differential equation for the market value reads This is a correlated jump-diffusion process. We model the jumps by a Poisson process with intensity λ. We recall that in such a process the probability function for the event to occur n times between zero and time t is given by The size Λ of the jump, measured in units of the current market value V k (t), is a random variable with a distribution which we have to specify. Jumps can be positive or negative. The largest possible negative jump is 100% of the current market value. Based on this information, a possible distribution of the jump size Λ is a shifted lognormal distribution, Λ + 1 ∼ LN(µ J , σ J ), with mean µ J and standard deviation σ J .
Without the jump term, the distribution of the market value V k (t) is log-normal. The jumps render the tails of the distribution fatter. The parameters of the jump process can be adjusted in order to match the tail behavior of a given empirical time series of the market value. Here, we use the same parameters for idiosyncratic and market wide jumps.
Let us now consider the stochastic process in Eq. (38) for discrete time increments ∆t = T /N , where the time to maturity T is divided into N steps. Then the market value of company k at maturity is given by Akin to the diffusion case, we find for the market return of the jump-diffusion For K → ∞ we can express the average over k as expectation values for ε k,t and dJ k,t . Since the averages are independent for different k and t, we can write where the averages are given as Unless we choose µ J = −σ 2 J /2, the expectation value for the jump term does not vanish. In the diffusion case we expanded the logarithm in Equation (11). An analogous step is not possible for the jump-diffusion, because the jump term dJ m,t does not scale with ∆t. This complicates a further analytical discussion of the jump-diffusion case.

B. Monte Carlo results
In the Monte Carlo simulation we consider the stochastic process in Equation (40) for discrete time increments. As in the diffusion case we simulate K = 500 different realizations of ε k,t and dJ k,t for a single realization of the market terms η m,t and dJ m,t with t = 1, . . . , N . For each realization of the idiosynchratic terms we calculate the market return X m , the number of defaults N D (X m ), the default probability P D (X m ) and the expected recovery rate R(X m ) . The full Monte Carlo simulation consists of 10 6 realizations of the market terms.
For the diffusion terms and the contract details we use the same parameters as in Section 2. The additional parameters for the jump terms are λ = 0.005, µ j = 0.4 and σ j = 0.3. Figure 5 shows the dependence of recovery rates R(X m ) on default probabilities P D (X m ) in one plot, and the dependence of portfolio losses L(X m ) on default probabilities P D (X m ) in a second plot. In both cases we observe a very good agreement between the Monte Carlo simulations of the jump-diffusion and the analytical results for the diffusion, see Equations (26) and (27), respectively. The average recovery rates deviate only for very low default probabilities from the diffusion result. For the portfolio loss this deviation is suppressed and no longer visible. The parameter B is here the same as in the diffusion case, this means it is not influenced by the additional jump terms of the process. The most striking difference to the plots in Figure 2 is the abundance of simulation results with very high default probabilities. This is due to the market-wide jumps which render the tails of the market return distribution fatter. The idiosyncratic jumps lead to slightly larger deviations of individual results from the mean values.
The analytical results for the diffusion case do not only describe the P D dependence of the jump-diffusion results, but also their dependence on the market return X m . This is shown in Figure 6. Both default probabilities and portfolio losses are well described by Equations (17) and (20), respectively. Given the functional relation between portfolio loss L and market return X m we can again use Equation 37 to transform the pdf of market returns into the pdf of portfolio losses. Here we use the simulation data on X m to determine their distribution numerically. In Figure 7 we compare this to the Monte Carlo result for the loss distribution. As in the diffusion case, we observe an excellent agreement, even for extremely large losses. The market-wide jumps lead to a very slow decay in the tail of the loss distribution. The results for Expected Loss, Value at Risk and Expected Tail Loss are listed in Table B. The analytical values describe the simulation results very well.
Our results show that the Merton model with jump-diffusion can be described by the diffusion results if we exchange the distribution of market returns. The dependencies on default probabilities and market returns are nearly the same as in the diffusion case. This is somewhat reminiscent of a copula approach to credit risk.  (1 − c)σ 2 T is determined from the simulation parameters. simulation analytical result EL 1.03 · 10 −3 9.07 · 10 −4 VaR 0.99 1.50 · 10 −2 1.47 · 10 −2 ETL 0.99 3.70 · 10 −2 3.75 · 10 −2

Correlated GARCH process
Finally, we examine the Merton model with a GARCH(1,1) model for the underlying asset value process. The GARCH model was first introduced by Bollerslev (1986). It is able to reproduce many of the stylized facts found in empirical financial time series. In particular, it exhibits volatility clustering and fat-tailed return distributions. The GARCH(1,1) model is a discrete time process with an autoregressive volatility. The return at time t reads where η t and ε k,t are independent normal distributed random variables. The parameters α 0 , α 1 and β 1 are chosen in order to mimic the behavior of a typical empirical time series of daily returns. The initial values for the volatilities have been set homogeneously as σ k,0 = σ √ ∆t, where σ = 0.15 is the same value used in the diffusion and the jump-diffusion case. However, the volatility does not remain the same for all k as t evolves, since σ k,t also depends on the idiosyncratic random part in r k,t−1 . Thus, the homogeneity of the portfolio is lost to some degree. Instead the process covers a wide range of volatilities both within one realization of the market terms, i.e., within one portfolio, and also over different market realizations. The GARCH process is therefore the most general case we examine here. It can provide a good indication of how broadly our analytical results for the diffusion case can be applied.
For the GARCH process the market value of company k at maturity reads The deterministic drift term can also be included in Equation (45) instead. Then the parameters α 0 , α 1 and β 1 have to be adjusted accordingly. The drift constant is again set to µ = 0.05. As in the diffusion and jump-diffusion case we simulate 10 6 portfolios of size K = 500. For each portfolio we calculate the market return X m , the number of defaults N D (X m ), the default probabilities P D (X m ) and the expected recovery rate R(X m ) .
In Figure 8 we present the Monte Carlo results for the P D dependence of the GARCH process. One plot shows the dependence of recovery rates R(X m ) on default probabilities P D (X m ), a second plot shows the dependence of portfolio losses L(X m ) on default probabilities P D (X m ). First, we observe a much broader range of the individual results than in the diffusion and jump-diffusion case. This is to be expected because the fluctuating volatilities lead to pronounced fat tails in the return distribution of individual companies. While the returns are correlated in the GARCH(1,1) model, volatilities are rather uncorrelated, see Schäfer and Guhr (2010). Thus, for the market return the fluctuations of individual volatilities average out to some degree. This is why we scarcely observe very high default probabilities in the simulations. The situation is different for empirical stock return time series, where volatilities are also strongly correlated.
It is a rather remarkable result that the average behavior of the GARCH simulation is so well described by the analytical results for the diffusion, see Equations (26) and (27), respectively. The parameter B has been fitted to the simulation data.
In Figure 9 we demonstrate the dependence of default probabilities and portfolio losses on the market return. Here the behavior of the GARCH model is completely different from the diffusion case. For the analytical results of the diffusion case, described by Equations (17) and (20), we use the parameter value B which has been fitted to the P D dependence in Figure 8. Given the functional relation between portfolio loss L and default probability P D we can transform the pdf of default probabilities into the pdf of portfolio losses, A comparison to the Monte Carlo result for the portfolio loss distribution is presented in Figure 10. We observe a very good agreement, even for extremely large losses. Finally, we present the results for Expected Loss, Value at Risk and Expected Tail Loss in Table 4. We compare the Monte Carlo result for Value at Risk with where F −1 P D (α) is the α-quantile of default probabilities. The Expected Tail Loss can be expressed as The analytical values describe the simulation results very well.
Our results demonstrate that the functional dependence of recovery rates on default probabilities does not depend on the underlying process. The analytical result for the diffusion process describes also other underlying processes in the Merton model. Hence we call Equation (26)

Applications of the structural recovery rate
In a recent study Altman et al. (2005) find a strong negative correlation between default probabilities and recovery rates in empirical data. The structural recovery rate provides such a negative relation between default and recovery rates. It correctly describes the case of zero-coupon bonds for various underlying processes, but may even be applicable to a more general debt structure. An indication for this is given by Chen and Panjer (2003) who show that the Merton model with jump-diffusion can be calibrated to provide the same yield spreads as reduced-form models.
Reduced-form models are also able to reproduce negative correlations between default and recovery rates, if the default model and the recovery model depend on a single common covariate, see Chava et al. (2008). Compared to the structural recovery rate, however, the reduced-form recovery introduces more parameters and lacks a deeper motivation.   In first passage models as described in Giesecke (2004), default may occur before maturity if the market value of the company falls below a default barrier. If defaults before maturity dominate in this model, recovery rates are independent from default rates. The average recovery rate is then constant and does not resemble the negative relation mentioned above. However, if defaults occur mostly at maturity, the structural recovery rate is recovered.
Finally, we give a first empirical indication that the structural recovery rate is indeed applicable in a realistic setting. Figure 11 shows data from Moody's annual default study. Although there are only very few data points, corresponding to the years 1987 to 2008, we can observe a negative relation between the firm-wide ultimate recovery rates and the annual default rates. We fit the structural recovery rate to these data points and find a reasonable agreement.

Conclusions
The interdependence of default and recovery rates has a crucial influence on large credit losses. Yet default probabilities and recovery rates are often modelled independently in current credit risk models. In this paper we revisited the Merton model for different underlying processes with correlations. While the original Merton model was conceived for the diffusion process, we also consider a jump-diffusion and a GARCH process. For the correlated diffusion we derived a functional dependence between default and recovery rates. This functional dependence is determined by a single parameter. In Monte Carlo simulations we showed that it describes not only the diffusion case, but also the jump-diffusion and GARCH process. Due to its independence on the underlying process, we call this functional relation the structural recovery rate. It is straightforward to use this structural recovery rate in addition to a first passage model or any other model for default probabilities. We believe that this has great potential to improve current credit risk models. In further studies we shall address two important questions: How well does the structural recovery rate describe empirical data on defaults and recoveries? And how does it compare with reduced-form models?