A Joint Quantile and Expected Shortfall Regression Framework

We introduce a novel regression framework which simultaneously models the quantile and the Expected Shortfall (ES) of a response variable given a set of covariates. This regression is based on a strictly consistent loss function for the pair quantile and ES, which allows for M- and Z-estimation of the joint regression parameters. We show consistency and asymptotic normality for both estimators under weak regularity conditions. The underlying loss function depends on two specification functions, whose choice affects the properties of the resulting estimators. We find that the Z-estimator is numerically unstable and thus, we rely on M-estimation of the model parameters. Extensive simulations verify the asymptotic properties and analyze the small sample behavior of the M-estimator for different specification functions. This joint regression framework allows for various applications including estimating, forecasting, and backtesting ES, which is particularly relevant in light of the recent introduction of ES into the Basel Accords.


Introduction
Measuring and forecasting risks is essential for a variety of academic disciplines. For this purpose, risk measures which are formally defined as a map (with certain properties) from a space of random variables to a real number, are applied to condense the complex nature of the involved risks to a single number (Artzner et al., 1999). In the context of financial risk measurement, to date the most commonly used risk measure is the Value-at-Risk (VaR), which is the α-quantile of the return distribution. Its popularity is mainly due to its simple nature and the fact that up to now, the Basel Accords stipulate its use for the calculation of capital requirements for banks. Besides being not coherent (Artzner et al., 1999), the main drawback of the VaR is its inability to capture tail risks beyond itself. This deficiency is overcome by the risk measure Expected Shortfall (ES) at level α, which is defined as the mean of the returns which are smaller than the α-quantile of the return distribution. The ES has the desired ability to capture information from the whole left tail of the return distribution, which is particularly important for measuring extreme financial risks. Over the past few years, ES has increasingly become the object of interest for practitioners, academics, and regulators, especially since its recent introduction into the Basel Accords (Basel Committee, 2016).
A major drawback of the ES (regarded as a statistical functional) is that it is not elicitable, which means that there exists no loss function (scoring function, scoring rule) which the ES uniquely minimizes in expectation (Gneiting, 2011;Weber, 2006). This result has two main consequences. First, consistent ranking of competing forecasts for the ES based on such a loss function is infeasible. Second, and more substantial for this paper, modeling the conditional ES given a set of covariates through a regression model without specifying the full conditional distribution is infeasible since estimation of the regression parameters through M-estimation requires such a loss function. Consequently, and in contrast to quantile regression (which can be used to model the VaR), to date, there exists no such regression framework which models the ES based on a set of covariates. Nadarajah et al. (2014) provide an overview of estimation methods for the ES. However, the reviewed approaches are only applicable for univariate data and not suitable for estimating the conditional ES based on covariates such as in mean and quantile regression. Nevertheless, there are some approaches for the ES which incorporate explanatory variables through indirect estimation procedures. Taylor (2008b) proposes an implicit approach for forecasting ES using exponentially weighted quantile regression and Taylor (2008a) introduces a procedure based on expectile regression and a relationship between the ES and expectiles. Taylor (2017) suggests a joint modeling technique for the quantile and the ES based on maximum likelihood estimation of the asymmetric Laplace distribution. Barendse (2017) proposes generalized method of moments (GMM) estimation for a regression framework for the interquantile expectation.
Even though the ES is not elicitable stand-alone,  show in their seminal paper that the quantile (the VaR) and the ES are jointly elicitable by introducing a class of joint loss functions, whose expectation is minimized by these two functionals. This joint elicitability result and the associated class of loss functions gives rise to a growing literature in both, joint estimation (Zwingmann and Holzmann, 2016) and in joint forecast evaluation (Acerbi and Szekely, 2014;Nolde and Ziegel, 2017;Ziegel et al., 2017) for the risk measures VaR and ES.
In this paper, we utilize the class of loss functions of  for the introduction of a novel simultaneous regression framework for the quantile and the ES and propose both, an M-and a Z-estimator for the joint regression parameters. These strictly consistent loss functions facilitate the opportunity to introduce M-and Z-estimation of the regression parameters without specifying the full conditional distribution of the model, as opposed to maximum likelihood estimation. We show consistency and asymptotic normality for both estimators under weak regularity conditions which are typical for such a regression framework. To the best of our knowledge, we are the first to propose such a joint regression framework for the quantile and the ES together with the joint M-and Z-estimation and the associated results of consistency and asymptotic normality. Furthermore, we are the first to propose a joint semiparametric regression framework for two different functionals based on joint M-estimation without specifying the full conditional distribution.
The employed joint loss function, the estimating equations (for the Z-estimator) and the resulting parameter estimates depend on two specification functions, which can be chosen from some class of functions. Even though consistency and asymptotic normality hold for all applicable choices of these specification functions, they affect the necessary moment conditions, the resulting asymptotic covariance matrices of the estimators, the numerical stability of the optimization algorithm, and the computation times. We discuss the choice of these functions in a theoretical context with respect to asymptotic efficiency and necessary regularity conditions, and with respect to the numerical properties of the optimization algorithm.
The estimation of the asymptotic covariance matrix imposes some difficulties. The first occurs in the estimation of the density quantile function, analogous to quantile regression (cf. Koenker, 2005) and thus, we utilize estimation procedures stemming from this literature. The second issue is the estimation of the variance of the negative quantile residuals conditional on the covariates, a nuisance quantity which is new to the literature. We introduce several estimators for this quantity which are able to cope with limited sample sizes and which can model the dependency of the negative quantile residuals on the covariates. Furthermore, we estimate the covariance matrix using the bootstrap. For ease of application, we provide an R package (Bayer and Dimitriadis, 2017a) which contains the implementation of the M-and Z-estimator. The user can choose the specification functions, the numerical optimization procedure and the estimation method for the covariance matrix of the parameter estimates.
We conduct a Monte-Carlo simulation study where we consider three data generating processes with different properties. We numerically verify consistency and asymptotic normality of the M-estimator for a range of different choices of the specification functions. Furthermore, we find that the Z-estimator is numerically unstable due to the redescending nature of the utilized estimating equations and consequently, we rely on M-estimation of the regression parameters. Moreover, we find that the performance of the M-estimator strongly depends on the specification functions, where choices resulting in positively homogeneous loss functions (Nolde and Ziegel, 2017;Efron, 1991) lead to a superior performance in terms of asymptotic efficiency, computation times, and mean squared error of the estimator. This joint regression technique for the quantile and ES has a wide range of potential applications as it generalizes quantile regression to the pair consisting of the quantile and the ES. In the context of financial risk management, it opens up the possibility to extend the existing applications of quantile regression on VaR in the financial literature to ES, such as e.g. in Chernozhukov and Umantsev (2001), Engle and Manganelli (2004), Koenker and Xiao (2006), Gaglianone et al. (2011), Halbleib andPohlmeier (2012), Komunjer (2013), Xiao et al. (2015) and Žikeš and Baruník (2016). Such estimation, forecasting, and backtesting methods for the ES are particularly sought-after in light of the recent shift from VaR to ES in the Basel Accords. As an illustration, we present an empirical application where we use our regression framework to jointly forecast VaR and ES based on the realized volatility.
The rest of the paper is organized as follows. In Section 2, we introduce the joint regression framework, the underlying regularity conditions together with the asymptotic properties of our estimators and discuss the choice of the specification functions. Section 3 provides details on the numerical implementation of the estimators and on the estimation of the asymptotic covariance matrix. Section 4 presents an extensive simulation study and Section 5 contains an empirical application. Section 6 provides concluding remarks. The proofs are deferred to Appendices B and C.

The Joint Regression Framework
Following Lambert et al. (2008), Gneiting (2011) and , we introduce the concept of (multivariate) p-elicitability. We consider a random variable Z : Ω → R d , defined on some complete probability space Ω, F, P , a class of distributions P on R d , equipped with the Borel σ-field and a functional T : P → D with its domain of action D ⊆ R p , p ∈ N. We call an integrable loss function ρ : R d × D → R strictly consistent for the functional T relative to the class of distributions P, if T is the unique minimizer of E ρ(Z, ·) for all distributions F ∈ P, where F is the distribution of Z. Furthermore, we call a p-dimensional functional T p-elicitable relative to the class P, if there exists a loss function ρ which is strictly consistent for T relative to P. If the dimension p is clear from the context, we simply call the functional elicitable instead of p-elicitable.
Given the generalized α-quantile Q α (Z) = F −1 (α) = inf z ∈ R : F(z) ≥ α for some α ∈ (0, 1), the ES of the random variable Z at level α is defined as If the distribution function of Z is continuous at its α-quantile, this definition can be simplified to the conditional tail expectation Gneiting (2011) shows that the ES is not 1-elicitable with respect to any class P of probability distributions on intervals I ⊆ R, which contain measures with finite support or finite mixtures of absolutely continuous distributions with compact support (see also Weber, 2006). This result has several consequences for the risk measure ES. First, consistent and meaningful ranking of competing forecasts for the functional ES is infeasible. Second, and more consequential for this work, estimating the parameters of a stand-alone regression model for the functional ES in the sense that ES α (Y |X) = X θ e by means of M-estimation, i.e. by minimizing some strictly consistent loss function, is infeasible. Even though the ES is not 1-elicitable,  show that the pair consisting of the ES and the quantile at common probability level α is 2-elicitable relative to the class of distributions with finite first moments and unique α-quantiles and they characterize the full class of strictly consistent loss functions for this pair subject to some regularity conditions. Since the definition of the ES already depends on the respective quantile, the fact that the ES is only elicitable jointly with the quantile is not surprising.
We utilize this joint elicitability result for the introduction of a new joint regression framework for the quantile and the ES where the aforementioned class of strictly consistent loss functions serves as the basis for the M-estimation of the joint regression parameters. For this, let Y : Ω → R and X : Ω → R k be random variables defined on the same probability space Ω, F, P as above. Henceforth, the transpose of X will be denoted by X , the cumulative distribution function of Y given X by F Y |X and the conditional density function by f Y |X . For a k-times differentiable real-valued function G : R → R, we denote the k-th derivative by G (k) (·).
Assumption 2.1 (The joint regression model). The regression framework which jointly models the conditional quantile and ES of Y given X for some fixed level α ∈ (0, 1) is given by We propose both, an M-estimation and a Z-estimation procedure for the compound regression parameter vector θ 0 . For the M-estimation, we adapt the class of strictly consistent joint loss functions1 for the quantile and ES as given in  such that it can be used in a regression framework, where the function G 1 is twice continuously differentiable, G 2 is three times continuously differentiable, G (1) 2 = G 2 , G 2 and G (1) 2 are strictly positive, G 1 is increasing and a and G 1 are integrable. We discuss the choice of the specification functions G 1 and G 2 in a theoretical context in Section 2.3 and by their numerical performance in Section 4.2. The corresponding (ρ-type) M-estimator is defined by a sequencê θ ρ,n , such thatθ ρ,n = argmin θ ∈Θ Instead of minimizing some objective function ρ(Y, X, θ) such as in (2.2), we can also define the corresponding Z-estimator (or ψ-type M-estimator), which sets a vector of estimating equations (moment conditions), denoted by ψ(Y, X, θ), to zero. More generally, it suffices that these estimating equations converge to zero almost surely. Formally, the Z-estimator is a sequenceθ ψ,n , such that which is obtained by differentiating2 (2.2) and where the functions G 1 and G 2 are given as above. When the loss function ρ(Y, X, θ) is continuously differentiable in θ, it is obvious that the M-and Z-estimation approaches are equivalent. However, in this case the loss function ρ(Y, X, θ) is not differentiable and ψ(Y, X, θ) is discontinuous at the points where Y = X θ q . Thus, we treat these two estimation approaches as different estimators and show their asymptotic behavior separately.
1One can interpret the structure of this loss function as follows : The first summand in (2.2) is a strictly consistent loss function for the quantile (Gneiting, 2011) and hence only depends on the quantile, whereas the second summand cannot be split into a part depending only on the quantile and one depending only on the ES. This illustrates the fact that the ES itself is not 1-elicitable, but 2-elicitable together with the respective quantile.
2 Note that the function ρ(Y, X, θ), given in (2.2) is only differentiable for Y X θ q . However, the points of nondifferentiability, Y = X θ q form a nullset with respect to the absolutely continuous distribution of Y given X.

Asymptotic Properties
In this section, we present the asymptotic properties of the M-and Z-estimator of the regression parameters. Consistency and asymptotic normality hold under the following set of weak regularity conditions, which are natural for this regression framework.

Assumption 2.2 (Regularity Conditions).
(A-1) The data (Y i , X i ) for i = 1, . . . , n is an iid series of random variables, distributed such as (Y, X) given above. Furthermore, the conditional distribution F Y |X has finite second moments and is absolutely continuous with probability density function f Y |X , which is strictly positive, continuous and bounded in a neighbourhood of the true conditional quantile, X θ q 0 . (A-2) The matrix E X X is positive definite.

Remark 2.3 (Finite Moment Conditions).
We further have to assume that certain moments of X are finite. For the sake of space, we specify the Finite Moment Conditions (M-1) -(M-4) in Appendix A. Note that these general moment conditions simplify substantially for sensible choices of the specification functions G 1 and G 2 as further outlined in Section 2.3.
Assumption (A-1) is a combination of typical regularity conditions of mean and quantile regression. Absolute continuity of F Y |X with a strictly positive, bounded and continuous density function in a neighborhood of the true conditional quantile is also imposed for the asymptotic theory of quantile regression. Existence of the conditional moments of Y given X is subject to the conditions of mean regression and is included in our regularity conditions since ES is a truncated mean. The positive definiteness (full rank condition) in (A-2) is common for any regression design with stochastic regressors in order to exclude perfect multicollinearity of the regressors. The conditions for the specification functions G 1 and G 2 in (A-3) mainly originate from the conditions for the joint elicitability of the quantile and ES in . Differentiability of these functions is required in this setup for obtaining the estimating equations and for the differentiations in the computation of the asymptotic covariance in Theorem 2.6 and Theorem 2.7. The existence of certain moments of the explanatory variables as in conditions (M-1) -(M-4) in Appendix A is also standard in any regression design relying on stochastic regressors. Even though compactness of the parameter space Θ in Assumption 2.1 generally simplifies the proofs, in this setup it is crucial for consistency of the Z-estimator as the estimating equations ψ 2 are redescending to zero for many reasonable choices of the G 2 function such as e.g. the choices resulting in positively homogeneous loss functions. For details on this, we refer to Section 3.1.
Theorem 2.4. Assume that Assumption 2.1, Assumption 2.2 and the Moment Conditions (M-1) in Appendix A hold true. Then, for every sequenceθ ψ,n ∈ Θ satisfying 1 Theorem 2.5. Assume that Assumption 2.1, Assumption 2.2 and the Moment Conditions (M-2) in Appendix A hold true. Then, for every sequenceθ ρ,n ∈ Θ such that 1 Theorem 2.6. Assume that Assumption 2.1, Assumption 2.2 and the Moment Conditions (M-3) in Appendix A hold true. Then, for every sequenceθ ψ,n ∈ Θ satisfying 1 (2.10) Theorem 2.7. Assume that Assumption 2.1, Assumption 2.2 and the Moment Conditions (M-4) in Appendix A hold true. Then, for every sequenceθ ρ,n ∈ Θ such that 1 where the matrices Λ and C are given as in Theorem 2.6.

Remark 2.8 (Quantile Regression).
Notice that the asymptotic covariance matrix of the quantile-specific parameter estimatesθ q is given by and (2.12) This simplifies to the covariance matrix of quantile regression parameter estimates by setting G 1 (z) = z and G 2 (z) = 0, which means ignoring the ES-specific part of our loss function and estimating equations. This demonstrates that the quantile regression method is nested in our regression procedure, also in terms of its asymptotic distribution.

Remark 2.9 (Asymptotic Covariance of the ES and the Oracle Estimator).
The ES-specific part of the asymptotic covariance is mainly governed by the term C 22 , which depends on the quantity 14) It is reasonable that the asymptotic covariance of ES regression parameters depends on the truncated variance of Y given X as the asymptomatic covariance of mean regression parameters is driven by the conditional (non-truncated) variance of Y given X. The second term X θ q 0 − X θ e 0 2 in (2.14) is included since the ES represents a truncated mean where the truncation point itself is a statistical functional (the quantile). In comparison, we consider an oracle M-estimator for the ES-specific regression parameters θ e , given by the loss function where we assume that the true quantile regression parameters θ q 0 are known. The resulting asymptotic covariance is given by which shows that the additional term X θ q 0 − X θ e 0 2 is not included for this estimator with fixed truncation point X θ q 0 .
Remark 2.10 (Joint Estimation of the Sample Quantile and ES). We can use this regression framework to jointly estimate the quantile and ES of an identically distributed sample Y 1 , . . . , Y n by regressing on a constant only. The asymptotic covariance matrix given in Theorem 2.6 and Theorem 2.7 then simplifies to Σ with components where θ q 0 and θ e 0 are the true quantile and ES of Y . The same result is obtained by Zwingmann and Holzmann (2016), who further allow for a distribution function for Y which is not differentiable at the quantile with strictly positive derivative. Notice that in this simplified case without covariates, the asymptotic covariance matrix is independent of the specification functions G 1 and G 2 used in the loss function and in the estimating equations. Furthermore, (2.17) implies that quantile estimates stemming from our joint estimation procedure have the same asymptotic efficiency as quantile estimates stemming from minimizing the generalized piecewise linear loss (Gneiting, 2011) and as sample quantiles (cf. Koenker, 2005). The same holds true for the efficiency of the sample ES estimators (based on the sample quantile) of Brazauskas et al. (2008) and Chen (2008).

Remark 2.11 (Pseudo-R 2 and the choice of a(Y )). By choosing a(Y
in (2.2), we can guarantee non-negative losses ρ(Y, X, θ) ≥ 0. This choice enables us to define a pseudo-R 2 for our joint regression framework in the sense of Koenker and Machado (1999), whereθ denotes the parameter estimates of the full regression model andθ denotes the parameter estimates of a regression model restricted to an intercept term only. However, this choice of a(Y ) comes at the cost of more restrictive moment conditions, since we need to impose that

Choice of the Specification Functions
The loss functions and the estimating equations given in (2.2) and (2.3) depend on two specification functions, G 1 and G 2 (with derivative G 2 ), which have to fulfill the regularity conditions (A-3) in Assumption 2.2.  already mention the feasible choices G 1 (z) = 0, G 1 (z) = z, G 2 (z) = exp(z) and G 2 (z) = exp(z)/ 1 + exp(z) in order to show that this class is non-empty. In contrast to the loss functions of mean, quantile and expectile regression, there is no natural choice for these specification functions for the quantile and ES yet (Nolde and Ziegel, 2017). However, as the choice of these functions strongly influences the performance of our regression procedure in terms of its asymptotic efficiency, the necessary moment conditions of the regressors and the numerical performance of the optimization algorithm, we discuss sensible selection criteria in the following. Efron (1991) and Nolde and Ziegel (2017) argue that for M-estimation of regression parameters it is crucial that the utilized loss function is positively homogeneous of some order b ∈ R in the sense that for all c > 0. This is an important property for loss functions since the ordering of the losses should be independent of the unit of measurement, e.g. the currency we measure the prices and risk forecasts with. Loss functions following this property guarantee that we can change the scaling and still obtain the same optima and consequently the same parameter estimates. For the pair consisting of the quantile and the ES, Nolde and Ziegel (2017) characterize the full class of positively homogeneous3 loss functions of order b for the case where we restrict the domain of G 2 , i.e. the conditional ES to the negative real line4, There are no positively homogeneous loss functions for the cases b ≥ 1. Our numerical simulations show that there is no gain in efficiency or numerical accuracy by deviating from the choice G 1 (z) = 0 (see also Nolde and Ziegel, 2017;Ziegel et al., 2017), which is also consistent with the homogeneity result. Consequently, we use G 1 (z) = 0 in the following. A different natural guiding principle for selecting the specification functions is induced by choosing G 2 (and G 1 ) such that the moment conditions (M-1) -(M-4) in Appendix A are as least restrictive and as parsimonious as possible. For instance, choosing G 2 such that G 2 and its first and second derivatives are bounded functions (and G 1 (z) = 0) results in the moment condition This motivates the usage of bounded func-tions5 for G 2 such as e.g. the second example of , G 2 (z) = exp(z)/ 1 + exp(z) , which is the distribution function of the standard logistic distribution. Further examples of bounded G 2 functions include the distribution functions of absolutely continuous distributions on the real line. In the simulation study in Section 4.2, we compare the performance of different specification functions in terms of mean squared error, asymptotic efficiency of the estimator and computation times.

Numerical Estimation of the Model
In this section, we discuss the difficulties one encounters and the solutions we propose for estimating the joint regression model. Section 3.1 illustrates the numerical optimization procedure we employ for estimating the regression parameters and Section 3.2 discusses different estimation methods for the covariance matrix of the estimator.

Optimization
Theorem 2.6 and Theorem 2.7 imply that both, M-estimation and Z-estimation of the regression parameters θ have the same asymptotic efficiency and consequently, we discuss these estimation approaches in terms of their numerical performance in the following. The numerical implementation of the Z-estimator relies 3For b = 0, only the loss differences are positively homogeneous. However, the ordering of the losses is still unaffected under this slightly weaker property.
4Since the conditional ES of financial assets for small probability levels is always negative, this is no critical restriction. However, for the numerical parameter estimation, we have to restrict the parameter space Θ such that X i θ e < 0 for all θ ∈ Θ and for all X i in the underlying sample. For details on this, we refer to Section 3.1.
5 Note that the positively homogeneous loss functions exhibit unbounded G 2 functions. However, as the function G 2 (z) does not grow faster than linear as z tends to infinity, the resulting finite moment conditions are not too restrictive. on root-finding of the estimating equations given in (2.3), which we implement as in GMM-estimation by minimizing the inner product . However, the estimating equations are redescending to zero for many attractive choices of G 2 in the sense that ψ 2 (Y, X, θ) → 0 for X θ e → −∞. Consequently, for θ such that θ q = θ q 0 and X θ e → −∞, we get the same minimal value of the Z-estimation as for the true regression parameters θ 0 . Thus, the Z-estimator is numerically unstable and diverges in many setups.
Consequently, we rely on M-estimation of the regression parameters in the following. As the loss functions given in (2.2) are not differentiable and non-convex for all applicable choices of the specification functions (Fissler, 2017), we apply a derivative-free global optimization technique. More specifically, we use the Iterated Local Search (ILS) meta-heuristic of Lourenço et al. (2003), which successively refines the parameter estimates by repeated optimizations with iteratively perturbed starting values. Our exact implementation consists of the following steps. First, we obtain starting values for θ q and θ e from two quantile regressions of Y on X for the probability levels α andα, where we chooseα such that thẽ α-quantile and the α-ES coincide under normality. Second, using these starting values we minimize the loss function with the derivative-free and robust Nelder-Mead Simplex algorithm (Nelder and Mead, 1965). Third, we perturb the resulting parameter estimates by adding normally distributed noise with zero mean and standard deviation equal to the estimated asymptotic standard errors of the initial quantile regression estimates. Fourth, we re-optimize the model with the perturbed parameter estimates as new starting values. If the loss is further decreased by this re-optimization, we update the estimates and otherwise, we retain the previous ones. Fifth, we iterate over the previous two steps until the loss does not decrease in m = 10 consecutive iterations. Our numerical experiments indicate that this repeated optimization procedure yields estimates very close to the ones stemming from other global optimization techniques such as e.g. simulated annealing, whereas the major advantage of ILS is the considerably lower computation time.
For the choices of the specification functions which result in positively homogeneous loss functions, we have to restrict the domain of G 2 to the negative real line as already discussed in Section 2.3. Thus, we have to restrict Θ such that X i θ e < 0 for all θ ∈ Θ and for all i = 1, . . . , n during the optimization process. Even though in financial risk management the response variable Y is usually given by financial returns where the true (conditional) ES is strictly negative, there might still be some outliers X i such that X i θ e 0 ≥ 0. In such a case, imposing the restriction X i θ e < 0 for all i = 1, . . . , n during the optimization process generates substantially biased estimates for θ e . In order to avoid this, we estimate the regression model for the transformed dependent variables Y − max(Y ) for the positively homogeneous loss functions and add max(Y ) to the estimated intercept parameters to undo the transformation6.
We provide an R package for the estimation of the regression parameters (see Bayer and Dimitriadis, 2017a). This package contains an implementation of both, the M-and the Z-estimator, where different optimization algorithms can be chosen (ILS, simulated annealing). The package allows for choosing the specification functions G 1 and G 2 and it includes an option to estimate the model either with or without the translation of the dependent variable. Furthermore, the covariance matrix of the parameter estimates can be estimated either by using the asymptotic theory and the resulting techniques we discuss in the next section, or by using the nonparametric iid bootstrap (Efron, 1979). We recommend applying the M-estimator with the ILS algorithm as this procedure exhibits the best performance in our numerical experiments with respect to accuracy, stability and computation times.
6 Note that this data transformation changes the average loss function as the applied loss functions are in general not translation invariant. Thus, optimizing the translated loss function can lead to different parameter estimates. However, we do not face the risk of obtaining substantially biased estimates in cases where X i θ e 0 ≥ 0 for some i ∈ {1, . . . n}. Our numerical experiments indicate that the difference between estimating the model for Y and for Y − max(Y ) is small when X i θ e 0 < 0 for all i ∈ {1, . . . n}, but can be quite substantial if there is an outlier for X i such that X i θ e 0 ≥ 0.

Asymptotic Covariance Estimation
While most parts of the asymptotic covariance matrix given in Theorem 2.6 and Theorem 2.7 are straightforward to estimate, two nuisance quantities impose some difficulties. The first is the density quantile function f Y |X (X θ q 0 ), which is already well investigated in the quantile regression literature. In particular, we consider the estimators proposed by Koenker (1994), henceforth denoted by iid and by Hendricks and Koenker (1992), henceforth denoted by nid. The main difference between these is that the first is based on the assumption that the quantile residuals are independent of the covariates, whereas the second allows for a linear dependence structure. Both approaches depend on a bandwidth parameter which we choose according to Hall and Sheather (1988).
The second nuisance quantity is the variance of the quantile residuals, conditional on the covariates and given that these residuals are negative, (3.1) Estimation of this quantity is demanding for two reasons. First, for very small probability levels which are typical in financial risk management such as e.g. α = 2.5%, the truncation u q ≤ 0 cuts off all but very few (about α · n) observations. Second, modeling this truncated variance conditional on the covariates X is challenging, especially considering the very small sample sizes. Under the assumption of homoscedasticity, i.e. that the distribution of u q is independent of the covariates X, we can simply estimate (3.1) by the sample variance of the negative quantile residuals and we refer to this estimator as ind in the following. We propose two further estimators which allow for a dependence of the quantile residuals on the covariates. For this purpose, we assume a location-scale process with linear7 specifications of the conditional mean and standard deviation in order to explicitly model the conditional relationship of u q on X, for some parameter vectors ζ, φ ∈ R k and where ε ∼ G(0, 1) follows a zero mean, unit variance distribution, such that u q |X ∼ G X ζ, (X φ) 2 with distribution function F G and density f G . As we need to estimate the truncated variance of u q given u q ≤ 0, i.e. a truncated variant of (X φ) 2 , one possibility is to estimate (3.2) only for those observations where u q ≤ 0. However, this approach particularly suffers from the very few negative quantile residuals as we need to estimate additional parameters compared to the ind approach. We present a feasible alternative by estimating the parameters ζ and φ using all available observations of u q and X by quasi generalized pseudo maximum likelihood (Gourieroux and Monfort, 1995, Section 8.4.4) and we obtain the truncated conditional variance by the scaling formula Var (u q |u q ≤ 0, is the truncated conditional density of u q given X and u q ≤ 0. We propose one parametric estimator, henceforth denoted by scl-N, where we assume that the distribution G is the normal distribution and apply a closed-form solution to the scaling formula. We further propose a semiparametric estimator, henceforth denoted by scl-sp, where we estimate the distribution G nonparametrically and then apply the scaling formula for this estimated density by numerical integration.

Simulation Study
In this section, we investigate the finite sample behavior of the M-estimator and verify the asymptotic properties derived in Section 2.2 through simulations. Furthermore, we compare the performance of different choices for the specification functions and evaluate the precision of the different covariance matrix estimators described in Section 3.2.

Data Generating Process
In order to assess the numerical properties of estimating the joint regression model, we simulate data from a linear location-scale data generating process (DGP), where v ∼ F(0, 1) has zero mean and unit variance, X = 1, X 2 , . . . , X k and γ, η ∈ R k . For this process, the true conditional quantile and ES are linear functions in X, given by where z α and ξ α are the quantile and ES of the distribution F(0, 1), which implies that θ q 0 = γ + z α η and θ e 0 = γ + ξ α η. Furthermore, the conditional distributions of the quantile-and ES-residuals are given by For the simulation study, we want to assess the performance of our regression procedure in various setups. Thus, we specify γ, η and F in the following such that we get data which is homoscedastic (DGP-(1)) and heteroskedastic (DGP- (2)). Furthermore, we include a regression setup with multiple, correlated regressors and a leptocurtic conditional distribution (DGP-(3)), We simulate all three processes 25,000 times with varying sample sizes of n = 250, 500, 1000, 2000 and 5000 observations. For each replication and for each of the sample sizes we regress the simulated Y 's on the covariates X using our joint regression method for the probability level α = 2.5%.

Comparing the Specification Functions
We start the discussion of the simulation results by investigating the numerical performance of the M-estimator based on different choices of the specification function8 G 2 used in the loss function in (2.2). We use three natural examples resulting in positively homogeneous loss functions of order b = −1, b = 0 and b = 0.5 respectively9, a bounded G 2 function and the (unbounded) exponential function: (4.4) Figure 1 presents the sum (over the 2k regression parameters) of the mean squared errors (MSE) of the regression parameters for the three DGPs described above, different sample sizes and for the five choices of the specification functions given in (4.4). As implied by the asymptotic theory, we obtain consistent parameter estimates for all five choices of the specification functions as the MSEs converge to zero for all three DGPs. However, they differ substantially with respect to their small sample properties. The three positively homogeneous specifications result in the most accurate estimates, whereas the choices G 2 (z) = − √ −z and G 2 (z) = − log(−z) tend to perform slightly better than the choice G 2 (z) = −1/z.
Furthermore, the bounded choice G 2 (z) = log 1 + exp(z) still performs better than the unbounded exponential function.   Table 1 reports the Frobenius norms of the lower triangular parts of the true asymptotic covariance matrices and of the respective (lower triangular) quantile-specific and the ES-specific sub-matrices for the three DGPs and for the five choices of the specification functions given in (4.4). For comparison, we also report the Frobenius norm of the lower triangular part of the asymptotic covariance of the quantile regression estimator. We approximate the true asymptotic covariance matrix through Monte-Carlo integration with a sample size of 10 9 using the formulas in Theorem 2.6 and by using the true density and conditional truncated variance. On average, the specification functions G 2 (z) = − log(−z) and G 2 (z) = − √ −z exhibit the smallest asymptotic covariances, closely followed by the third choice for a positively homogeneous loss function, G 2 (z) = −1/z. The non-homogeneous choices lead to considerably larger asymptotic variances for all considered DGPs and sub-matrices. Furthermore, by comparing the quantile-specific parameters of the joint estimation approach (from the positively homogeneous loss functions) to quantile regression estimates, we roughly obtain the same asymptotic efficiency. Table 1: This table reports the Frobenius norms of the lower triangular parts of the asymptotic covariance matrices and the respective quantile-specific and the ES-specific sub-matrices for the three DGPs and for the five choices of the specification functions given in (4.4). For comparison, we report the same quantity for the asymptotic covariance of the quantile regression estimator.

Comparing the Variance-Covariance Estimators
In this section, we compare the empirical performance of the asymptotic covariance estimators discussed in Section 3.2. For the comparison of their precision, Figure 2 reports the average of the Frobenius norm of the lower triangular part of the differences between the estimated covariances and the empirical covariance of the estimated parameters. We report results for the three homogeneous loss functions and the three  DGPs, where each of the plots presents the average norm differences for the four covariance estimators (iid/nid, nid/scl-N, nid/scl-sp and the iid bootstrap) depending on the sample size.
We find that the iid/nid estimator performs well for the first, homoscedastic DGP whereas for the other two DGPs, it fails to capture the underlying more complicated dynamics of the data. The nid/scl-N estimator outperforms the other estimation approaches in the first two DGPs, where the underlying conditional distribution follows a normal distribution whereas its performance drops for the third DGP, which follows a Student-t distribution. The performance of the flexible nid/scl-sp estimator is the most stable throughout all three DGPs. Eventually, the bootstrap estimator accurately estimates the covariance for all three DGPs, whereas in comparison to the other estimators, it is particularly good in small samples. The provided R package contains all four covariance estimators.

Empirical Application
In this empirical application, we use our joint regression framework for forecasting the VaR and ES of the close-to-close log returns of the IBM stock. For that purpose, we adopt the forecasting framework of Žikeš and Baruník (2016) and jointly forecast the VaR and ES of daily financial returns r t by Q α (r t |RV t−1 ) = θ q 1 + θ q 2 RV t−1 and ES α (r t |RV t−1 ) = θ e 1 + θ e 2 RV t−1 , (5.1) where RV t = ( i r 2 t,i ) 1/2 denotes the realized volatility estimator (Andersen and Bollerslev, 1998) for day t, where r t,i denotes the i-th high-frequency return of day t. Our dataset consists of the five minute returns of the IBM stock from January 3, 2001 to July 18, 2017 with total of 4120 days, which we obtain from the TAQ database. We estimate the model parameters using a rolling window of 1000 days and evaluate the forecasts on the remaining 3120 days.
We compare the predictive power of this model against three standard models from the literature. The first is the historical simulation (HS) approach, which forecasts the VaR and ES for day t as the sample quantile and ES of the daily returns of the past 250 trading days. The second is an AR(1)-GARCH(1,1)-t model (Bollerslev, 1986), and the third is the Heterogeneous Auto-Regressive (HAR) model of Corsi (2009), based on the realized volatility estimates given above. Forecasts of the VaR and ES for the HAR model are obtained from the volatility forecasts and by assuming a Gaussian return distribution. While the first two of these approaches rely on daily data only, the third one incorporates the same high frequency information as our approach.
We evaluate the forecasting power of the VaR and ES of these models by the class of strictly consistent loss (scoring) functions for the VaR and ES of . We  Figure 3 displays the average of the elementary score differences of the joint VaR and ES regression model against the three alternative models together with the respective 95% pointwise confidence bands for the elementary scores provided in Ziegel et al. (2017) for the pair VaR and ES. Using this graphical method, we can see that the elementary score differences for the joint regression forecasting model against the historical simulation and AR(1)-GARCH(1,1)-t model are significantly negative for the vast majority of threshold values. This implies that the joint regression forecasting model significantly dominates these other two forecasting approaches. Even though we also observe strictly negative elementary score differences in comparison against the HAR model, these differences are not significant and consequently, we cannot significantly outperform this model.

Conclusion
In this paper, we introduce a joint regression technique for the quantile (the VaR) and the ES. This regression approach relies on the class of strictly consistent joint loss functions introduced by , which permits the joint elicitation of the quantile and the ES. We introduce an M-and a Z-estimator for the parameters of the joint regression model. Given a set of standard regularity conditions, we show consistency and asymptotic normality for both estimators, which we also verify numerically through extensive simulations. The underlying loss functions, the estimating equations and the asymptotic covariance matrices of the estimators depend on the choice of two specification functions, which we investigate in terms of the resulting moment conditions, asymptotic efficiency, numerical performance and computation times. In our numerical simulations, we find that choices resulting in positively homogeneous loss functions dominate other choices with respect to the aforementioned criteria. Furthermore, we propose several estimation methods for the asymptotic covariance matrix, which are able to cope with different properties of the underlying data. We provide an R package (see Bayer and Dimitriadis, 2017a), which implements the M-and Z-estimation procedures where one can choose the underlying specification functions, the numerical optimization approach and the estimation method for the asymptotic covariance matrix.
Our new joint regression technique allows for a wide range of applications for the risk measures VaR and ES. This regression approach can be used to model the ES (jointly with the VaR) by generalizing existing applications of quantile regression on VaR, such as e.g. in Koenker and Xiao (2006), Engle and Manganelli (2004), Chernozhukov and Umantsev (2001), Žikeš and Baruník (2016), Halbleib and Pohlmeier (2012), Komunjer (2013) and Xiao et al. (2015). As an illustration, we present an empirical application in this paper where we use this regression framework to jointly forecast VaR and ES based on realized volatility estimates. Furthermore, Bayer and Dimitriadis (2017b)
(M-1) For Theorem 2.4, we assume that the following moments are finite for some d 0 > 0: (M-2) For Theorem 2.5, we assume that the following moments are finite: (M-3) For Theorem 2.6, we assume that the following moments are finite for some constant d 0 > 0 and for all θ ∈Ū d 0 (θ 0 ): (M-4) For Theorem 2.7, we assume that the following moments are finite for some constant d 0 > 0:

Appendix B Proofs
Henceforth, ||v|| denotes the maximum norm for a vector v ∈ R k and for a matrix A, || A|| denotes the row-sum matrix norm which is induced by the maximum norm for vectors. For convenience of the supremum notation, for all θ ∈ int(Θ) and for some d > 0, we define the open neighborhood U d (θ) = {τ ∈ Θ : ||τ − θ|| < d} and its closureŪ d (θ) = {τ ∈ Θ : ||τ − θ|| ≤ d}. All references to Appendix C refer to the online supplement Dimitriadis and Bayer (2017).
Proof of Theorem 2.4. We apply Theorem 2 from Huber (1967) and show that the function ψ(Y, X, θ) as given in (2.3) satisfies the respective assumptions of this theorem. Note that the parameter space Θ is assumed to be compact and thus, we do not have to show condition (B-4) in the notation of Huber (1967). As the product of continuous functions and the indicator function 1 {Y ≤X θ q } , the function ψ is measurable and regarded as a stochastic process in θ, ψ is separable in the sense of Doob as it is almost surely continuous in θ (Gikhman and Skorokhod, 2004, p.164). This condition assures measurability of the suprema10 given below and in Lemma C.1. In oder to show that ψ has a unique root at θ 0 , let us first define the sets for all θ ∈ Θ such that Ω = W θ ∪ U θ and W θ ∩ U θ = ∅. We first show that P(U θ ) > 0 for all θ θ 0 . In order to see this, we assume the converse, i.e. let us assume that for a fixed θ θ 0 , it holds that However, since θ q θ q 0 , this contradicts the assumption that the matrix E[X X ] is positive definite and we can conclude that P(U θ ) > 0.
The quantity exists under the moment conditions (M-1) in Appendix A and if θ q = θ q 0 , it holds that λ 1 (θ) = 0. Now, we assume that θ ∈ Θ such that θ q θ q 0 . By splitting the expectation, we get that The first summand is obviously zero since for all ω ∈ W θ , F Y |X (X θ q ) − F Y |X (X θ q 0 ) = 0. Since the distribution of Y given X has strictly positive density in a neighbourhood of X θ q 0 , we get that F Y |X is strictly increasing in a neighbourhood of X θ q 0 and thus for all ω ∈ U θ . Furthermore, since αG (1) 1 (X θ q ) + G 2 (X θ e ) > 0 for all θ ∈ Θ and P(U θ ) > 0, we get that and consequently λ 1 (θ) 0. This implies that λ 1 (θ) = 0 if and only if θ q = θ q 0 . Furthermore, Thus, (B.4) simplifies to E (X X )G (1) 2 (X θ e ) θ e − θ e 0 and by applying Lemma C.2, we get that the matrix E (X X )G (1) 2 (X θ e ) is positive definite for all θ ∈ Θ. Consequently, λ 2 (θ) = 0 if and only if θ e = θ e 0 and together with the arguments for λ 1 , we get that λ(θ) = 0 if and only if θ = θ 0 . Eventually, assumption (B-2)' from Theorem 2 of Huber (1967) follows directly from Lemma C.1, which concludes this proof.
Proof of Theorem 2.5. For this proof, we apply Theorem 5.7 from van der Vaart (1998) and show that the respective assumptions of this theorem hold. As in the proof of Theorem 2.6, we can conclude measurability of the suprema since the process ρ is continuous and consequently separable in the sense of Doob. Thus, we do not have to rely on outer probability measures such as in van der Vaart (1998). We start by showing uniform convergence in probability of the empirical mean of the objective function by the help of Lemma 2.4 of Newey and McFadden (1994). Since we have iid data, a compact parameter space Θ and ρ(Y, X, θ) is continuous for all θ ∈ Θ, it remains to show that there exists a dominating function and it holds that d(Y, X) ≥ ρ(Y, X, θ) for all θ ∈ Θ and consequently, we can conclude uniform convergence in probability. We now show that E ρ(Y, X, θ) has a unique and global minimum at θ = θ 0 . For this, we assume that θ ∈ Θ such that θ θ 0 and we define the sets such that Ω = U θ ∪ W θ and U θ ∩ W θ = ∅. We first show that P(U θ ) > 0 for all θ θ 0 . In order to see this, we assume the converse, i.e. we assume that P(W θ ) = 1, which implies that However, since θ θ 0 and consequently either θ q θ q 0 or θ e θ e 0 , this contradicts the assumption that the matrix E[X X ] is positive definite and it follows that P(U θ ) > 0.
From the joint elicitability property of the quantile and ES of , Corollary 5.5 we get that for all x ∈ R k such that x θ q x θ q 0 or x θ e x θ e 0 , it holds that since the distribution of Y given X has a finite first moment and a unique α-quantile. Thus, for all ω ∈ U θ , We now define the random variable and (B.9) implies that h X, θ, θ 0 (ω) < 0 for all ω ∈ U θ . Since P(U θ ) > 0, this implies that E h(X, θ, θ 0 )1 {ω ∈U θ } < 0. Furthermore, for all ω ∈ W θ , it obviously holds that h(X, θ, θ 0 )(ω) = 0 and consequently E h(X, θ, θ 0 )1 {ω ∈W θ } = 0. Thus, we get that for all θ ∈ Θ such that θ θ 0 , which shows that E ρ(Y, X, θ) has a unique minimum at θ = θ 0 .
Proof of Theorem 2.6. We apply Theorem 3 of Huber (1967) for the ψ-function as given in (2.3) and show the respective assumptions of this theorem. Consistency of the Z-estimator is shown in Theorem 2.4. For the measureability and separability of the ψ function, we refer to the proof of Theorem 2.4. It is already shown in the proof of Theorem 2.4 that there exists a θ 0 ∈ Θ such that λ(θ 0 ) = 0. For the technical conditions (N-3), we apply Lemma C.3, Lemma C.1 and Lemma C.4. It remains to show that E ||ψ(Y, X, θ 0 )|| 2 < ∞, which follows from the subsequent computation of C and the Moment

Conditions (M-3) in Appendix A. The asymptotic covariance matrix is given by
Straightforward calculations yield the matrix C as given in (2.8) -(2.10). For the computation of Λ, we first notice that the function is continuously differentiable for all θ in some neighborhood U d (θ 0 ) around θ 0 , since the distribution F Y |X has a density which is strictly positive, continuous and bounded in this area. Let us choose a valuẽ θ ∈ U d (θ 0 ) such that X θ ≤ X θ. Then, (B.14) We consequently get that for all θ ∈ U d (θ 0 ), In order to conclude that ∂ ∂θ E E ψ(Y, X, θ) X = E ∂ ∂θ E ψ(Y, X, θ) X , we apply a measure-theoretical version of the Leibniz integration rule, which requires that the derivative of the integrand exists and is absolutely bounded by some integrable function d(Y, X), independent of θ. For the first term, this can easily be obtained by defining which has finite expectation by the Moment Conditions (M-3). The other two terms follow the same reasoning. Inserting θ = θ 0 eventually shows (2.6) and (2.7).
Proof of Theorem 2.7. For this proof, we apply Theorem 5.23 from van der Vaart (1998) and show that the respective assumptions of this theorem hold. Theorem 2.5 shows consistency of the M-estimator. The map (Y, X) → ρ(Y, X, θ) is obviously measurable as the sum of measurable functions. Furthermore, the map θ → ρ(Y, X, θ) is almost surely differentiable since the only point of non-differentiability occurs where Y = X θ q , which is a nullset with respect to the joint distribution of Y and X and for all θ ∈ Θ such that Y X θ q , its derivative is given by ψ(Y, X, θ). Local Lipschitz continuity with square-integrable Lipschitz-constant follows from Lemma C.5. We have already seen in the proof of Theorem 2.5 that the function E ρ(Y, X, θ) is uniquely minimized at the point θ 0 and is twice continuously differentiable and consequently admits a second-order Taylor expansion at θ 0 . Thus, we have shown the necessary assumptions of Theorem 5.23 from van der Vaart (1998). For the computation of the covariance matrix, we notice that the distribution of Y given X has a density f Y |X in a neighborhood of X θ 0 , which is strictly positive, continuous and bounded. Therefore, by the same arguments as in (B.14), we get that ∂ ∂θ q E G 1 (Y )1 {Y ≤X θ q } X = XG 1 (X θ q ) f Y |X (X θ q ). Thus, straightforward calculations yield that for all θ ∈ U d (θ 0 ), it holds that ∂ ∂θ E ρ(Y, X, θ) X = E ψ(Y, X, θ) X and by applying the Leibniz integration rule such as in the proof of Theorem 2.6, we finally get that Consequently, the asymptotic covariance matrix equals the one given in Theorem 2.6.

Appendix C Technical Results
and assume that Assumption 2.1, Assumption 2.2 and the Moment Conditions (M-1) in Appendix A hold. Then, there are strictly positive real numbers b and d 0 , such that and for all d ≥ 0.
Proof of Lemma C.1. For measurability of the suprema, we refer to the proof of Theorem 2.4. Let in the following d > 0 and θ ∈ Θ such that ||θ − θ 0 || + d ≤ d 0 . We first notice that for some fixed X ∈ R k and for all τ ∈Ū d (θ), it holds that for all Y ∈ R and for some θ q − , θ q + ∈Ū d (θ). SinceŪ d (θ) is compact, we get that for all Y ∈ R and for some values θ q − , θ q + ∈Ū d (θ). Note that the values θ q − and θ q + depend on X and θ, however they are independent of Y . Consequently, it holds that where we apply the mean value theorem for someθ q on the line between θ q − and θ q + , i.e.θ q ∈Ū d (θ).
For the first component of ψ, we get that (C.6) The first term in (C.6) is O(d) since G (1) 1 (X θ q ) and G 2 (X θ e ) are continuously differentiable functions w.r.t θ and thus, by the mean value theorem we get that and the respective moments are finite by assumption. The same arguments hold for the function G 2 . For the second term in (C.6), we apply (C.5) and thus get that (C.8) Since the density f Y |X is bounded in a neighborhood of X θ q 0 and the respective moments are finite by assumption, we get that this term is also O(d).
For the second component of ψ, we get that The first, third and fifth term are linearly bounded by (C.7) since the functions (X θ e − X θ q )G (1) 2 (X θ e ) and (X θ q )G (1) 2 (X θ e ) and G (1) 2 (X θ e ) are continuously differentiable. For the second term, we use the arguments from (C.5). For the fourth term, we use similar arguments as in (C.5), and get that there exist some θ q − , θ q + ∈Ū d (θ) and a valueθ q on the line between θ q − and θ q + , such that since f Y |X is bounded in a neighborhood of X θ 0 and the respective moments exist by assumption. This concludes the proof of the lemma.
Lemma C.2. Let the random variable X ∈ R k with distribution P be such that its second moments exist and the matrix E[X X ] is positive definite. Furthermore, letΘ ⊂ R k be a compact subspace with nonempty interior and let g : R k ×Θ → R be a strictly positive function. Then, the matrix is also positive definite.
Proof of Lemma C.2. Since E[X X ] is positive definite, we know that for all z ∈ R k with z 0, it holds that 0 < z E[X X ]z = E[z (X X )z] = E[(X z) 2 ] and consequently P X z 0 > 0. Since g(X, θ) is a strictly positive scalar for all θ ∈Θ, it also holds that P (X z) g(X, θ) 0 > 0 and thus, for all z 0, This positivity statement holds since X z g(X, θ) 2 is a non-negative random variable and P (X z) g(X, θ) 0 > 0. This shows that the matrix E (X X )g(X, θ) is positive definite.
Proof of Lemma C.4. Let in the following d > 0 and θ ∈ Θ such that ||θ − θ 0 || + d ≤ d 0 . It holds that and consequently, we show that for both components j = 1, 2 and for some d > 0 small enough.
For the first squared component, we get that where we apply (C.5) for the second summand. The remaining two summands can be bounded linearly by the arguments given in (C.7) since G (1) 1 and G 2 are continuously differentiable functions and the respective moments are finite.

(C.30)
Thus, in order to evaluate E sup τ ∈Ū d (θ) ψ 2 (Y, X, τ) − ψ 2 (Y, X, θ) 2 , we have to consider all the cross products out of the five summands in (C.30). Since the techniques applied are very similar, we only show details for two of the cross products.
Proposition C.6. Let Y be a real-valued random variable with distribution function F, finite first and second moments and a unique α-quantile q α = F −1 (α). Then, Proof. We first notice that for a distribution F with finite second moment und unique α-quantile, it holds that which can be obtained by using the identity and by taking expectations on both sides. By applying (C.41), we get that (C.45) and by rearranging the order of integration for the first term in (C.45), we get that (C.46) Thus, by first using (C.45) and (C.46) and by plugging in (C.41) and (C.44), we obtain (C.47) Eventually, using (C.44) and (C.47), straight-forward calculations yield that {x|ψ(x, θ) ∈ F, ∀θ ∈ G ∩ I} (D.2) differ from each other at most by a subset of N. (2004)). Let Θ and Y be metric spaces, Θ be a separable space. The sets (D.1) and (D.2) coincide for all x ∈ Ω for which the stochastic process ψ(x, θ) is continuous in θ.
Thus, letθ ∈ G \ I. Since I is a dense set in Θ, there exists a sequence (θ n ) n∈N ∈ Θ ∩ I, such that θ n →θ and since G is an open set in Θ andθ ∈ G, we can conclude that for m ∈ N large enough, θ n ∈ G for all n ≥ m. Furthermore, by continuity at θ, it holds that ψ(x, θ n ) → ψ(x,θ) and since θ n ∈ G ∩ I for all n large enough, ψ(x, θ n ) ∈ F by assumption. Eventually, since F is a closed set, ψ(x,θ) ∈ F which proves the proposition.