Combining Predictive Distributions

Predictive distributions need to be aggregated when probabilistic forecasts are merged, or when expert opinions expressed in terms of probability distributions are fused. We take a prediction space approach that applies to discrete, mixed discrete-continuous and continuous predictive distributions alike, and study combination formulas for cumulative distribution functions from the perspectives of coherence, probabilistic and conditional calibration, and dispersion. Both linear and non-linear aggregation methods are investigated, including generalized, spread-adjusted and beta-transformed linear pools. The effects and techniques are demonstrated theoretically, in simulation examples, and in case studies on density forecasts for S&P 500 returns and daily maximum temperature at Seattle-Tacoma Airport.


Introduction
Probabilistic forecasts aim to provide calibrated and sharp predictive distributions for future quantities or events of interest. As they admit the assessment of forecast uncertainty and allow for optimal decision making, probabilistic forecasts continue to gain prominence in a wealth of applications, ranging from economics and finance to meteorology and climatology (Gneiting 2008). The general goal is to maximize the sharpness of the predictive distributions subject to calibration (Murphy and Winkler 1987;Gneiting, Balabdaoui and Raftery 2007). For a real-valued outcome, a probabilistic forecast can be represented in the form of a predictive cumulative distribution function, which might be discrete, mixed discrete-continuous or continuous, with the latter case corresponding to density forecasts.
In many situations, complementary or competing probabilistic forecasts from dependent or independent information sources are available. For example, the individual forecasts might stem from distinct experts, organizations or statistical models. The prevalent method for aggregating the individual predictive distributions into a single combined forecast is the linear pool (Stone 1961). While other methods for combining predictive distributions are available (Genest and Zidek 1986;Clemen and Winkler 1999;2007), the linear pool is typically the method of choice, with the pioneering work of Winkler (1968) and Zarnowitz (1969), and 1 recent papers by Mitchell and Hall (2005), Wallis (2005), Hall and Mitchell (2007), Jore, Mitchell and Vahey (2010), Kascha and Ravazzolo (2010) and Garratt et al. (2011) being examples in the case of density forecasts. Similarly, linear pools have been applied successfully to combine discrete predictive distributions; for recent reviews, see Ranjan and Gneiting (2010) and Clements and Harvey (2011). Despite the ubiquitous success of the linear pool in a vast number of applications, fragmented recent work points at potential shortcomings and limitations. Hora (2004) proved that any nontrivial convex combination of two calibrated density forecasts is uncalibrated. In the discrete case, Dawid, DeGroot and Mortera (1995) and Ranjan and Gneiting (2010) showed that linear combination formulas with strictly positive coefficients fail to be coherent and demonstrated potential improvement under nonlinear aggregation.
Our goal here is to unify these results and to extend them in various directions. Towards this end, we develop novel theoretical approaches to studying combination formulas and aggregation methods, where we think of an aggregation method as a class of combination formulas. For example, the traditional linear pools comprises the linear combination formulas with nonnegative weights that sum to one. Technically, we operate in terms of cumulative distribution functions, which allows us, in contrast to earlier work, to provide a unified treatment of all real-valued predictands, including the cases of density forecasts, mixed discrete-continuous predictive distributions, probability mass functions for count data and probability forecasts of a dichotomous event, all of which are important in applications. The extant literature compares combination formulas by examining whether or not they possess certain analytic characteristics, such as the strong setwise function and external Bayes properties (Genest and Zidek 1986;French and Ríos Insua 2000). Again in contrast to much of the earlier work, we assess combination formulas and aggregation methods from the perspectives of coherence, calibration and dispersion.
Section 2 sets the stage by introducing the key tool of a prediction space, which is a probability space tailored to the study of forecasts and combination formulas. We revisit the work of  and Ranjan and Gneiting (2010) in the prediction space setting and show, perhaps surprisingly, that if the outcome is binary, probabilistic calibration and conditional calibration are equivalent. Throughout the technical parts of the paper, information sets are represented by σ-algebras, and a predictive distribution is considered ideal if it agrees with the conditional distribution of the predictand, given the information basis or σ-algebra at hand. Section 3 is devoted to the study of specific, linear and non-linear combination formulas and aggregation methods. A key result is, roughly, that dispersion tends to increase under linear pooling. This helps explain the success of linear combination formulas in aggregating underdispersed component distributions, and leads us to unify and strengthen the aforementioned results of Dawid et al. (1995), Hora (2004) and Ranjan and Gneiting (2010). In particular, any linear combination formula with strictly positive coefficients fails to be coherent, and the traditional linear pool fails to be flexibly dispersive. In view of these limitations, we investigate parsimonious nonlinear alternatives, including generalized linear pools (Dawid et al. 1995), the spread-adjusted linear pool, which has been used successfully in the meteorological literature (Berrocal, Raftery and Gneiting 2007;Glahn et al. 2009), and the beta-transformed linear pool (Ranjan and Gneiting 2010), which we demonstrate to be flexibly dispersive.
Section 4 turns to a simulation study and data examples on density forecasts for daily maximum temperature at Seattle-Tacoma Airport and S&P 500 returns. The paper ends in Section 5, where we summarize our findings from applied and theoretical perspectives, suggest directions for future work and, in addition to discussing probabilistic forecasts, hint at the closely related problem of the fusion of expert judgements that are expressed in terms of probability distributions.

Combination formulas and aggregation methods
In a seminal paper, Murphy and Winkler (1987) proposed a general framework for the evaluation of point forecasts, which is based on the joint distribution of the forecast and the observation. Dawid et al. (1995) developed and used a related framework in studying multiple probability forecasts for a binary event. Here we respond to the call of Dawid et al. (1995, p. 288) for an extension, and we start with an informal sketch of a fully general approach, in which the observations take values in just any space.
The most general setting considers the joint distribution of multiple probabilistic forecasts and the observation on a probability space (Ω, A, Q). More explicitly, we assume that the elements of the sample space Ω can be identified with tupels of the form (P 1 , . . . , P k , Y ), where each of P 1 , . . . , P k is a probability measure on the outcome space of the observation, Y . For i = 1, . . . , k, we require the random probability measure P i to be measurable with respect to the sub-σ-algebra A i ⊆ A that encodes the forecast's information set or information basis, consisting of data, expertise, theories and assumptions at hand. The probability measure Q on (Ω, A) specifies the joint distribution of the probabilistic forecasts and the observation.
In this setting, the probabilistic forecasts P 1 , . . . , P k might stem from distinct experts, organizations or statistical models, as commonly encountered in the practice of forecasting. In aggregating them, the theoretically optimal strategy is to combine information sets, that is, to issue the conditional distribution of the observation Y given the σ-algebra σ(A 1 , . . . , A k ) generated by the information sets A 1 , . . . , A k . However, as Dawid et al. (1995, p. 264) note, "this ideal will almost always be rendered unattainable, by the extent of the data, company confidentiality, or the inability of the experts to identify clearly the empirical basis and background knowledge leading to their intuitive opinions." The best that we can hope for in practice is to find the conditional distribution of the observation Y given the σ-algebra σ(P 1 , . . . , P k ) generated by the random probability measures P 1 , . . . , P k . Of course, it is always true that σ(P 1 , . . . , P k ) ⊆ σ(A 1 , . . . , A k ), and in most cases of practical interest the left-hand side constitutes a substantially lesser information basis than the right-hand side.

Prediction spaces
In what follows, we restrict the discussion to the case of a real-valued observation. A probabilistic forecast then corresponds to a Lebesgue-Stieltjes measure on the real line, R, which we identify with the associated right-continuous cumulative distribution function (CDF). We write σ(A 1 , . . . , A m ) and σ(X 1 , . . . , X n ) to denote the σ-algebra generated by the families A 1 , . . . , A m of subsets of Ω, and the random variables X 1 , . . . , X n , respectively. We use the symbol L generically to denote an unconditional or conditional law or distribution and follow standard conventions in identifying the sub-σ-algebras on which we condition.
We now introduce the key tool of a prediction space, which is a probability space tailored to the study of combination formulas for real-valued outcomes, though we allow the case k = 1 of a single probabilistic forecast.
Definition 2.1. Let k ≥ 1 be an integer. A prediction space is a probability space (Ω, A, Q) together with sub-σ-algebras A 1 , . . . , A k ⊆ A, where the elements of the sample space Ω can be identified with tupels (F 1 , . . . , F k , Y, V ) such that (P1) for i = 1, . . . , k, F i is a CDF-valued random quantity that is measurable with respect to the sub-σ-algebra A i , 1 (P2) Y is a real-valued random variable, (P3) V is a random variable that is uniformly distributed on the unit interval and independent of A 1 , . . . , A k and Y .
All subsequent definitions and results are within the prediction space setting. Phrases such as almost surely or with positive probability refer to the probability measure Q on (Ω, A) that determines the joint distribution of the probabilistic forecasts and the observations. While (P1) and (P2) formalize the predictive distributions and the observation, assumption (P3) is purely technical, allowing us to define a generalized version of the classical probability integral transform. The sub-σ-algebra A i encodes the information set for the CDF-valued random quantity F i which may, but need not, be ideal in the following sense. 2 Definition 2.2. The CDF-valued random quantity F i is ideal relative to the sub-σ-algebra In the subsequent examples we write N (µ, σ 2 ) to denote a univariate normal distribution with mean µ and variance σ 2 , and similarly for the bivariate normal distribution. We use the symbols Φ and φ to denote the standard normal cumulative distribution function and density function, respectively. 1 That is, {F i (x j ) ∈ B j for j = 1, . . . , n} ∈ A for all finite collections x 1 , . . . , x n of real numbers and B 1 , . . . , B n of Borel sets.
2 In a recent comment, Tsyplakov (2011) proposes the same terminology.
Example 2.3 (probability forecasts of a binary event). Here we find a prediction space for the simulation example of Ranjan and Gneiting (2010), which considers a dichotomous outcome, Y . For technical consistency later on, we identify a success or occurence with the outcome Y = 0, and a non-success with Y = 1. The CDF-valued random quantities F i thus are of the form F i (y) = p i I(y ≥ 0) + (1 − p i )I(y ≥ 1) and can be identified with the random probability forecasts p i for a success, for i = 1, 2. Writing ω = (ω 1 , ω 2 , ω 3 , ω 4 ) for an elementary event, we define the probability forecasts where σ 1 > 0 and σ 2 > 0 are fixed constants, the observation Y (ω) = ω 3 , and the auxiliary variable V (ω) = ω 4 . To complete the specification of the prediction space, we let Ω = R × R × {0, 1} × (0, 1) and let A = B 4 be the corresponding Borel-σ-algebra, define Q to be the product of N (0, σ 2 1 ), N (0, σ 2 2 ) and standard uniform measures on the first, second and fourth coordinate projections, respectively, and let for Borel sets B 1 , B 2 ⊆ R, where λ denotes the Lebesgue measure. Then p 1 is measurable with respect to the sub-σ-algebra A 1 = σ(ω 1 ), and p 2 is measurable with respect to A 2 = σ(ω 2 ). Moreover, by the arguments in the appendix of Ranjan and Gneiting (2010), F 1 is ideal relative to A 1 , and F 2 is ideal relative to A 2 .
Typically, it suffices to consider the joint distribution of the tuple (F 1 , . . . , F k , Y ), without any need to explicitly specify other facets of the prediction space.
Example 2.4 (density forecasts). To define a prediction space, let where µ ∼ N (0, 1), and let τ attain the values 1 and −1 with equal probability, independently of µ and Y . Table  1 places the density forecasts in the simulation example of Gneiting, Balabdaoui and Raftery (2007) in this setting. The perfect forecast is ideal relative to the sub-σ-algebra generated by µ. The climatological forecast is ideal relative to the trivial sub-σ-algebra. The unfocused and sign-biased forecasts are not ideal, as we will see in Example 2.10 below.

Combination formulas
As noted, in aggregating predictive cumulative distribution functions, the theoretically optimal strategy is to combine information sets, that is, to issue the conditional distribution of the observation Y given the σ-algebra σ(A 1 , . . . , A k ) generated by the information sets A 1 , . . . , A k . However, information aggregation often is not feasible in practice, when individual sources of expertise reveal predictive distributions, rather than information sets. What we can realistically aim at is to model the conditional distribution of the observation Y given the σ-algebra generated by the predictive cumulative distribution functions, namely where we define with Q being the set of the rational numbers. In practice one resorts to empirical combination formulas. Specifically, let F be a class of fixed, non-random cumulative distribution functions such that F 1 , . . . , F k ∈ F almost surely. For example, if we are concerned with density forecasts on on the real line R, we consider the class D R of the cumulative distribution functions that admit a Lebesgue density. Further 6 classes F of interest are listed in Table 2. A combination formula then is a mapping of the form We allow the case k = 1, where the mapping may provide calibration and dispersion adjustments for a single predictive distribution, as described later on in Sections 3.4 and 5.
Here our interest is in the case k ≥ 2, where Dawid et al. (1995) introduced the notion of a coherent combination formula in the context of probability forecasts of a binary event. We now define a more general notion that applies to probabilistic forecasts of general, real-valued quantities.
Definition 2.5. Suppose that k ≥ 2. The combination formula G : F k → F is coherent relative to the class F if there exists a prediction space (Ω, A, Q) with sub-σ-algebras A 1 , . . . , A k ⊆ A and CDF-valued random quantities F 1 , . . . , F k such that It is important to note that condition (C1) has two requirements, the first being that F i = L(Y |A i ) is ideal, and the second that F i ∈ F almost surely. Condition (C2) excludes trivial cases, while (C3) requests the aggregated cumulative distribution function to be ideal relative to the σ-algebra generated by the full set F 1 , . . . , F k of the components. Note that the smaller the class F , the stronger a statement about coherence. Conversely, the larger the class F , the stronger a statement about lack of coherence.
Conceptually, a coherent combination formula is compatible with fully informed and theoretically literate decision making, whereas an incoherent combination formula is not. In this light, coherence is an appealing property of a combination formula. Nevertheless, the applied relevance of the notion of coherence is limited, as probabilistic forecasts are hardly ever ideal in practice, which may invite, or even necessitate, the use of flexible classes of potentially incoherent combination formulas. Motivated by these applied perspectives, we now turn to a discussion of the important notions of calibration and dispersion.

Calibration and dispersion
If F denotes a fixed, non-random predictive cumulative distribution function for an observation Y , the probability integral transform is the random variable Z F = F (Y ). It is well known that if F is continuous and Y ∼ F then Z F is standard uniform (Rosenblatt 1952). If the more general, randomized version of the probability integral transform proposed by Brockwell (2007) is used, the uniformity result applies to arbitrary, not necessarily continuous, but still fixed, non-random cumulative distribution functions. In the prediction space setting, we need the following, further extension that allows for F to be a CDF-valued random quantity.
Definition 2.6. In the prediction space setting, the random variable is the probability integral transform of the CDF-valued random quantity F .
In a nutshell, the probability integral transform is the value that the predictive cumulative distribution function attains at the observation, with suitable adaptions at points of discontinuity. The probability integral transform takes values in the unit interval, and so the possible values of its variance are constrained to the closed interval [0, 1 4 ]. A variance of 1 12 corresponds to a uniform distribution and continues to be the most desirable, as evidenced by Theorem 2.9 below.
We are now ready to define and study notions of calibration and dispersion. In doing so, we use the terms CDF-valued random quantity and forecast interchangeably.
Definition 2.7. In the prediction space setting, let F and G be CDF-valued random quantities with probability integral transforms Z F and Z G .
(e) The forecast F is regular if the distribution of Z F is supported on the unit interval.
In the defining equality E Q [F (y)] = Q(Y ≤ y) for marginal calibration, the left-hand side depends on the law of the predictive distribution, whereas the right-hand side depends on the law of the observation. In parts (c) and (d) of Definition 2.7, we define dispersion in terms of the variance of the probability integral transform, thus involving the joint law of the predictive distribution and the observation. In contrast, the spread of the predictive distribution itself is a measure of sharpness that does not consider the observation.
Our current setting of prediction spaces differs from, but relates closely to, the approach of , who studied notions of calibration from a prequential perspective. Specifically, if the CDF-valued random quantity F is probabilistically calibrated in the sense of Definition 2.7 and we sample from the joint law of F and Y , the resulting sequence is probabilistically calibrated in the sense of . An analogous statement applies to marginal calibration.
Returning to the prediction space setting, the following result is immediate.
Proposition 2.8. A probabilistically calibrated forecast is neutrally dispersed and regular.
The converse is not necessarily true, in that a forecast which is neutrally dispersed need not be calibrated nor regular. However, an ideal forecast is always calibrated.
Theorem 2.9. A forecast that is ideal relative to a σ-algebra ist both marginally calibrated and probabilistically calibrated.
where I denotes an indicator function, thereby proving the statement about marginal calibration. Turning to probabilistic calibration, let Q 0 denote the marginal law of Y under Q, for z ∈ (0, 1), where the final equality uses the key result of Brockwell (2007).
Example 2.10. We revisit the density forecasts in Example 2.4 and Table 1. The perfect forecast and the climatological forecast are ideal, and so by Theorem 2.9 they are both probabilistically calibrated and marginally calibrated. Arguments nearly identical to those in  show that the unfocused forecaster is probabilistically calibrated but not marginally calibrated, and that the sign-biased forecaster is marginally calibrated but not probabilistically calibrated. Hence, there is no sub-σ-algebra or information set relative to which the unfocused or the sign-biased forecaster is ideal. Dawid (1984), Diebold et al. (1998),  and Czado et al. (2009), among others, have argued powerfully that probabilistic calibration is a critical requirement for probabilistic forecasts that take the form of predictive cumulative distribution functions, with Theorem 2.9 lending further support to this approach. Indeed, checks for the uniformity of the probability integral transform have formed a cornerstone of density forecast evaluation. In practice, one observes a sample {(F 1j , . . . , F kj , y j ) : j = 1, . . . , J} from the joint distribution of the probabilistic forecasts and the observation, and the uniformity of the probability integral transform is assessed empirically. The prevalent way of doing this is by plotting histograms of the probability integral transform values for the various forecasting methods, which show the corresponding frequency distribution over an evaluation or test set. U-shaped histograms correspond to underdispersed predictive distributions with prediction intervals that are too narrow on average, while hump or inverse U-shaped histograms indicate overdispersed predictive distributions. Formal tests of uniformity can also be employed; for a review, see Corradi and Swanson (2006).
Example 2.11. Let Y = X + ǫ, where X and ǫ are independent, standard normal random variables, and consider the Gaussian predictive distribution F σ = N (X, σ 2 ). A stochastic domination argument, the details of which we give in Appendix A, shows that F σ is underdispersed if σ < 1, neutrally dispersed if σ = 1 and overdispersed if σ > 1. If σ = 1 then F σ is ideal and thus both marginally calibrated and probabilistically calibrated. A more detailed calculation, which is also given in Appendix A, shows that the probability integral (2) In Figure 1 we plot var(Z σ ) as a function of the predictive standard deviation, σ. Figure  2 shows probability integral transform histograms for a Monte Carlo sample of size 10, 000 from the joint distribution of the observation Y and the forecasts F σ , where σ = 3 4 , 1 and 5 4 . The histograms are U-shaped, uniform, and inverse U-shaped, reflecting underdispersion, neutral dispersion and calibration, and overdispersion, respectively.
Recall from Example 2.3 that in the case of a binary outcome Y , we identify a CDFvalued random quantity F (y) = pI(y ≥ 0) + (1 − p)I(y ≥ 1) with the probability forecast p for a success, that is, Y = 0. The extant literature, including Schervish (1989) and Ranjan and Gneiting (2010) and the references therein, calls p calibrated if Here we refer to this property as conditional calibration. Perhaps surprisingly, our next result shows that if the outcome is binary, the notions of probabilistic calibration and conditional calibration are equivalent. For an illustrating example, see Appendix C.
Theorem 2.12. Consider a prediction space (Ω, A, Q) with a binary outcome Y , where Y = 0 corresponds to a success and Y = 0 to a failure, and a CDF-valued random quantity F (y) = p I(y ≥ 0) + (1 − p)I(y ≥ 1), which can be identified with the probability forecast p for a success. Then the following statements are equivalent: (i) The forecast F is probabilistically calibrated, that is, its probability integral transform Z F is uniformly distributed on the unit interval.
(ii) The probability forecast p is conditionally calibrated, that is, Q(Y = 0| p) = p almost surely.
(iii) The forecast F is ideal relative to the σ-algebra generated by the probability forecast p.
Proof. It is clear that (ii) and (iii) are equivalent, and by Theorem 2.9 the statement (iii) implies (i). To conclude the proof, we show that statement (i) implies (ii). To this end, suppose that the forecast F is probabilistically calibrated. By standard properties of conditional expectations, there exists a measurable function q : [0, 1] → [0, 1] such that Q(Y = 0| p) = q(p) almost surely. Let H denote the marginal law of p under Q. If H has a point mass at 0 or 1, it is readily seen that q(0) = 0 or q(1) = 1, respectively.   A version of the conditional density u(z|x) of the probability integral transform Z F given for Lebesgue-almost all z ∈ (0, 1) and δ ∈ (0, 1 − y). Let 0 < a < b < 1, and consider the signed measure defined by : q(x) < x}, which implies that q(x) = x almost surely with respect to the restriction of H to [a, b]. To summarize, we have shown that q(x) = x almost surely with respect to H, whence Q(Y = 0| p) = p almost surely with respect to Q, as desired.
Theorem 2.12 draws a connection from the probability integral transform histogram to the reliability diagram or calibration curve, which is the key diagnostic tool for assessing the calibration of probability forecasts for a binary event (Dawid 1986;Murphy and Winkler 1992;Ranjan and Gneiting 2010). A reliability diagram plots conditional event frequencies against binned forecast probabilities, with deviations from the diagonal indicating violations of the conditional calibration condition (3). For non-binary outcomes Y and the natural generalization of the conditional calibration criterion, namely the auto-calibration property introduced by Tsyplakov (2011) in a comment on Mitchell and Wallis (2010), the equivalence to probabilistic calibration fails, as demonstrated in Table 3 for a ternary outcome. In general, auto-calibration is a much stronger requirement than probabilistic calibration, and there is a small but colorful, growing strand of recent literature that addresses its applied aspects for continuous outcomes (Hamill 2001;Mason et al. 2007;Held, Rufibach and Balabdaoui 2010;Bröcker, Siegert and Kantz 2011;Mason et al. 2011).

Aggregation methods
An aggregation method is a family of combination formulas of the form that share a common value of k and a common class F of fixed, non-random cumulative distribution functions. For example, if G is the traditional linear pool, we can take F to be any convex class of cumulative distribution functions, and we may identify the index set Θ with the unit simplex in R k . The extant literature studies individual combination formulas by examining whether or not they possess certain analytic characteristics, such as the strong setwise function and external Bayes properties (McConway 1981;Genest 1984aGenest , 1984bGenest and Zidek 1986;Genest, McConway and Schervish 1986;French and Ríos Insua 2000). In contrast, we share the recent perspective of Hora (2010) and put the focus on calibration and dispersion. In particular, we study aggregation methods in terms of the behavior of the probability integral transform under the corresponding family G = {G θ : θ ∈ Θ} of combination formulas. As noted, the probability integral takes values in the unit interval, and so the possible values of its variance lie between 0 and 1 4 . The value 1 12 corresponds to neutral dispersion and is the most desirable.
Following French and Ríos Insua (2000, p. 113) we say that the combination formula G θ is anonymous if for all F 1 , . . . , F k ∈ F and all permutations π on k elements. With this, we are ready to define notions of flexibility for aggregation methods.
Definition 2.13. Consider a family G = {G θ : θ ∈ Θ} of combination formulas of the form (4) that share a common k ≥ 2 and a common class F of fixed, non-random cumulative distribution functions.
(a) The aggregation method G is flexibly dispersive relative to the class F if for all F 0 ∈ F and F 1 , . . . , The applied relevance of the definitions is appreciated as follows. Suppose that the aggregation method G is flexibly dispersive relative to F . Then, given any marginal law F 0 ∈ F for the observation Y and any collection F 1 , . . . , F k of probabilistic forecasts for Y , 13 we can find a combination formula G θ ∈ G such that the aggregated predictive distribution, namely G θ (F 1 , . . . , F k ), is neutrally dispersed. If G is exchangeably flexibly dispersive, we can do so while treating F 1 , . . . , F k exchangeably, which is a frequent requirement in the practice of the combination of expert judgements (Jouini and Clemen 1986).
Note that a positive statement about flexible dispersivity is the stronger, the larger the class F . Conversely, a statement about the lack of flexible dispersivity is the stronger, the smaller the class F .

Linear and nonlinear aggregation methods
In this section we investigate specific combination formulas and aggregation methods from the perspectives of coherence, calibration and dispersion. First, we consider the traditional linear pool, then we move on to discuss non-linear ramifications, namely generalized linear pools, the spread-adjusted linear pool and the beta transformed linear pool. Note that the linearity of a combination formula for cumulative distribution functions is preserved in the corresponding formula for aggregating densities or event probabilities. In contrast, the functional form of a nonlinear combination formula changes if it is expressed in terms of densities or event probabilities.

Linear pool
We proceed to state and prove a simple but powerful result about linear combination formulas that generalizes earlier findings by Hora (2004) and Ranjan and Gneiting (2010). The gist of the statement is that dispersion tends to increase under linear aggregation.
Theorem 3.1. In the prediction space setting, suppose that k ≥ 2 and consider any linearly combined probabilistic forecast F = k i=1 w i F i with weights w 1 , . . . , w k that are strictly positive and sum to 1. For i = j, suppose that F i = F j with positive probability. Then the following holds: (a) The linearly combined forecast F is at least as dispersed as the least dispersed of the components F 1 , . . . , F k .
(b) If the component forecasts F 1 , . . . , F k are regular, then F is more dispersed than the least dispersed of the components.
(c) If the components are neutrally dispersed and regular, then F is overdispersed.
14 which demonstrates part (a). To prove part (b) suppose, for a contradiction, that F 1 , . . . , F k are regular and var(Z) = max 1≤i≤k var(Z i ). Then Z i and Z j are perfectly correlated for i, j = 1, . . . , k, and we conclude that there exist constants a ij > 0 and b ij ∈ R such that Z i = a ij Z j + b ij almost surely. By the assumption of regularity, Z i and Z j are supported on the unit interval, whence a ij = 1 and b ij = 0. Therefore, Z i = Z j almost surely, contrary to the assumption that F i = F j with positive probability. Part (c) concerns the special case of part (b) in which var(Z i ) = 1 12 for i = 1, . . . , k. As noted, Theorem 3.1 yields various extant results as corollaries. For instance, Hora (2004) applied Fourier analytic tools to show that if two distinct density forecasts are probabilistically calibrated, then any nontrivial linear combination is uncalibrated, which is an immediate consequence of part (c). However, the statement of part (c) is considerably stronger, in that it substitutes the weaker condition of neutral dispersion and regularity for the assumption of probabilistic calibration, allows for any number k ≥ 2 of components, allows for cumulative distribution functions rather than the special case of densities, and exposes the direction of the deviation, in that the linearly combined forecast is overdispersed. Each of the four facets is useful in practice. For instance, there are real data situations where density forecasts are approximately neutrally dispersed and regular, but clearly not calibrated. Discrete and mixed discrete-continuous predictive cumulative distribution functions also occur frequently in practice, such as in quantitative precipitation forecasting (Sloughter et al. 2007) and for count data (Czado et al. 2009). Finally, the tendency to increase dispersion helps explain the success of linear pooling in applications, where the component distributions are frequently underdispersed. For a prominent example, see Table 10 of Hoeting et al. (1999).
The next theorem nests extant results in the case of probability forecasts for a binary event, which were proved by Dawid et al. (1995) and Ranjan and Gneiting (2010).
Theorem 3.2. If k ≥ 2, any linear combination formula with strictly positive weights fails to be coherent relative to the class F R .
Proof. Suppose, for a contradiction, that the combination formula G(F 1 , . . . , F k ) = F = k i=1 w i F i with strictly positive weights that sum to 1 is coherent relative to F R . Then there exists a prediction space (Ω, A, Q) with sub-σ-algebras A 1 , . . . , A k ⊆ A and CDF-valued random quantities F 1 , . . . , F k such that properties (C1), (C2) and (C3) of Definition 2.5 hold. By property (C1) and Theorem 2.9, the components F 1 , . . . , F k are probabilistically calibrated, and by Proposition 2.8 they are neutrally dispersed and regular. Thus, part (c) of Theorem 3.1 applies and the linearly aggregated forecast F is overdispersed. However, by property (C3), Theorem 2.9 and Proposition 2.8 F is neutrally dispersed, for the desired contradiction.
Thus far, we have considered individual linear combination formulas. The following result views the traditional linear pool as an aggregation method G = {G θ : θ ∈ Θ}, where we may identfy the parameter space Θ = ∆ k−1 with the unit simplex in R k . We state the theorem relative to the full class F R , even though it remains valid relative to much smaller classes.  (5). The table states assumptions on the weights, w 1 , . . . , w k , and instances of classes F , such that the combination formula maps F k into F . The conditions depend on the domain and the range of the link function, h, which we assume to be continuous and strictly monotone.
Theorem 3.3. If k ≥ 2, the linear pool fails to be flexibly dispersive relative to the class F R .
Proof. In view of part (a) of Theorem 3.1 it suffices to find an F 0 ∈ F R and distinct F 1 , . . . , F k ∈ F R , each of which is an overdispersed as a probabilistic forecast for an observation Y with L(Y ) = F 0 . For example, we can take F 0 to be standard normal and F i to be normal with mean zero and variance i + 1 for i = 1, . . . , k. Dawid et al. (1995) introduced and studied generalized linear combination formulas for combining probability forecasts of a binary event. Here we apply the approach to cumulative distribution functions, where we obtain combination formulas of the form

Generalized linear pools
where h is a continuous and strictly monotone link function. It is interesting to note the formal resemblance to Archimedean copulas (Genest and MacKay 1986;McNeil and Neslehová 2009), as observed and investigated by Jouini and Clemen (1996). Table 4 shows conditions on the weights, w 1 , . . . , w k , along with instances of classes F , so that the generalized linear combination formula (5) maps F k into F . In the first type, the link function is defined on the closed unit interval, and the combination formula operates on the full class F R , with the traditional linear pool, for which h(x) = x is the identity function, being the most prominent example. In the remaining types, the link function is defined on the open unit interval only, and we need to restrict attention to F + I or B + , with the harmonic pool and the geometric pool being key examples, occuring when h(x) = 1/x and h(x) = log x, respectively. While not being exhaustive, the listing in the table is comprehensive, in that most link functions can be adapted to fit one of the types considered.
The following result applies to generalized linear combination formulas with nonnegative weights that sum to at most 1.
Theorem 3.4. Suppose that k ≥ 2, and let I ⊆ R be an interval. Consider a generalized linear combination formula G of the form (5) with a continuous and strictly monotone link function, h, and weights w 1 , . . . , w k that are strictly positive and sum to at most 1. Proof. We proceed similarly to the proofs of Theorems 3.1 and 3.2, which correspond to the special case where h is the identity function and the weights sum to 1. As regards part (a) suppose, for a contradiction, that the generalized linear combination formula G is coherent relative to the class C R . Consider the prediction space (Ω, A, Q) with sub-σalgebras A 1 , . . . , A k ⊆ A, CDF-valued random quantities F 1 , . . . , F k ∈ C R , and observation Y such that conditions (C1), (C2) and (C3) of Definition 2.5 hold. Let X = h(G(Y )) and X i = h(F i (Y )) for i = 1, . . . , k. By (C1), (C3), the continuity of F 1 , . . . , F k , G and h, and Theorem 2.9, the random variables X and X 1 , . . . X k are identical in distribution and bounded, and so they share a common finite value of the variance. Furthermore, condition (C2) implies that if i = j then X i = X j with positive probability, whence X i and X j cannot be linearly dependent. Consequently, for the desired contradiction. In part (b) we argue identically, noting that even though X and X 1 , . . . , X k may now be unbounded, they still share a common finite value of the variance, by the assumption of square-integrability for the link function. Dawid et al. (1995) proved in various special cases that generalized linear combination formulas with nonnegative weights summing to 1 fail to be coherent relative to the class B + of the nondegenerate Bernoulli measures. Furthermore, Dawid et al. (1985, pp. 282-283) conjectured that such results hold in broad generality. While Theorem 3.4 does not apply to Bernoulli measures, as its proof depends on the continuity of F 1 , . . . , F k and G, we share this belief.
However, if we allow for individual weights that exceed 1, the situation may change, as exemplified now in the case of Bernoulli measures.

Hence, the type D generalized linear combination formula
with a normal quantile link function, Φ −1 , and weights w 1 > 1 and w 2 > 1 is coherent relative to the class B + of the nondegenerate Bernoulli measures, or any larger class.
The defining equation (5) for a generalized linear combination formula implies that if the weights w 1 , . . . , w k are nonnegative and sum to at most 1 then In the previous section we applied this type of argument to show that the traditional linear pool fails to be flexibly dispersive. Similarly, we conjecture that generalized linear pools with nonnegative weights that are bounded above fail to be flexibly dispersive.
It is important to note that the extant literature considers generalized linear combination formulas that operate on event probabilities or densities, rather than cumulative distribution functions. For example, Dawid et al. (1995) study geometric and harmonic linear pools for a binary outcome, while Bjørnland et al. (2011) apply the geometric pool to density forecasts. Whether or not the resulting combination formulas and aggregation methods are coherent and flexibly dispersive, respectively, remains to be investigated.

Spread-adjusted linear pool
The aforementioned limitations of linear and generalized linear pools suggest that we consider more flexible, nonlinear aggregation methods. In this section, we focus on the class D + R , so that we may identify the cumulative distribution functions F 1 , . . . , F k with the corresponding Lebesgue densities f 1 , . . . , f k .
In the context of probabilistic weather forecasts and approximately neutrally dispersed Gaussian components f 1 , . . . , f k , Berrocal et al. (2007), Glahn et al. (2009) and Kleiber et al. (2011) observed empirically that linearly combined predictive distributions are overdispersed, as confirmed by Theorem 3.1. In an ad hoc approach, they proposed a nonlinear aggregation method which we now generalize and refer to as the spread-adjusted linear pool (SLP), as follows.
To describe this technique, it is convenient to write where µ i is the unique median of F i ∈ D + R , for i = 1, . . . , k. The SLP combined predictive distribution then has cumulative distribution function and Lebesgue density respectively, where w 1 , . . . , w k are nonnegative weights that sum to 1, and c is a strictly positive spread adjustment parameter. For neutrally dispersed or overdispersed components values of c < 1 are appropriate; for example, Table 2 of Berrocal et al. (2007) reports estimates ranging from 0.65 to 1.03. Underdispersed components may suggest values of c ≥ 1, and the traditional linear pool arises when c = 1. We do not know whether or not there are coherent combination formulas of this type. However, the SLP method performs well in the aforementioned applications, and the following result serves to quantify its flexibility.
Proposition 3.6. Suppose that L(Y ) = F 0 ∈ D + R and that F 1 , . . . , F k ∈ D + R have medians µ 1 ≤ · · · ≤ µ k . Let Z c = G c (Y ) denote the probability integral transform of the SLP aggregated predictive cumulative distribution function. Let v 0 = 0 and p 0 = F 0 (µ 1 ), let v i = i j=1 w j and p i = F 0 (µ i+1 ) − F 0 (µ i ) for i = 1, . . . , k − 1, and let v k = 1 and p k = 1 − F 0 (µ k ). Then as the spread adjustment parameter c > 0 varies, the variance of Z c attains any positive value less than (7) Proof. As c → 0, the probability integral transform Z c converges weakly to the discrete probability measure with mass p 0 , . . . , p k at v 0 , . . . , v k , which has variance (7). As c → ∞, it converges weakly to the Dirac measure in 1 2 . In view of var(Z c ) being a continuous function of the spread-adjustment parameter c > 0, this proves the claim.
Our next result views the spread-adjusted linear pool as an aggregation method with parameter space Θ = ∆ k−1 × R + . While the SLP approach is sufficiently rich in typical applications, where the individual predictive distributions are neutrally dispersed or underdispersed, its flexibility is limited.
Theorem 3.7. The spread-adjusted linear pool fails to be flexibly dispersive relative to the class D + R .
Proof. Let F 0 be standard normal, and for i = 1, . . . , k let F i be normal with mean m + i m and variance 1. As m → ∞, the probability integral transform of the SLP combined forecast G c attains values less than 1 2 with probability tending to one, irrespectively of the values of the SLP weights w 1 , . . . , w k and the spread adjustment parameter c. Thus, if m is sufficiently large, the variance of the PIT remains below the critical value of 1 12 that corresponds to neutral dispersion.
The SLP combination formula (6) can be generalized to allow for distinct spread adjustment parameters for the individual components. However, such an extension does not admit neutral dispersion either, and tends not to be beneficial in applications, unless the component densities have drastically varying degrees of dispersion. The assumption of a common spread adjustment parameter yields a more parsimonious model and stabilizes the estimation.

Beta-transformed linear pool
The beta transformed linear pool (BLP) composites the traditional linear pool with a beta transform. Introduced by Ranjan and Gneiting (2010) in the context of probability forecasts for a binary event, it generalizes readily to the full class F R of the cumulative distribution functions on R. Specifically, the BLP combination formula maps F 1 , . . . , for y ∈ R. Here, w 1 , . . . , w k are nonnegative weights that sum to 1, and B α,β denotes the cumulative distribution function of the beta density with parameters α > 0 and β > 0.
In contrast to the spread-adjusted linear pool, the value of the BLP aggregated predictive cumulative distribution function G α,β at y ∈ R depends on F 1 , . . . , F k only through the values F 1 (y), . . . , F k (y), in a locality characteristic that resembles the strong setwise function property of McConway (1981). If F i has Lebesgue density f i for i = 1, . . . , k, the aggregated cumulative distribution function G α,β is absolutely continuous with Lebesgue density where b α,β denotes the beta density with parameters α > 0 and β > 0. This nests the traditional linear pool that arises when α = β = 1.
The following result concerns the flexibility of the BLP combination formula (8) when the cumulative distribution functions F 0 ∈ C R and F 1 , . . . , F k ∈ C R are continuous and the weights w 1 , . . . , w k ≥ 0 are fixed, while the transformation parameters vary.
The next result views the beta-transformed linear pool as an aggregation method with parameter space Θ = ∆ k−1 × R 2 + .
Theorem 3.9. The beta transformed linear pool is exchangeably flexibly dispersive relative to the class C + I , for every interval I ⊆ R.
As hinted at in Section 2.4, the spread-adjusted and beta transformed linear pools can be applied in the case k = 1 of a single probabilistic forecast to provide calibration and dispersion adjustments. If the original predictive distribution is symmetric about its center, the symmetry is retained by the spead-adjusted linear pool, and broken by the beta transformed linear pool, except when α = β.
In practice, the BLP weights w 1 , . . . , w k and transformation parameters α, β > 0 are estimated from training data, say {(F 1j , . . . , F kj , y j ) : j = 1, . . . , J}. If the predictive cumulative distribution functions F 1j , . . . , F kj are absolutely continuous with Lebesgue densities where B denotes the classical beta function. The logarithmic score is simply the logarithm of the value that the density forecast attains at the realizing observation. It is positively oriented, that is, the higher the score, the better, and it is proper, in the sense that truth telling is an expectation maximizing strategy. Alternatively, the corresponding estimates can be viewed as maximum likelihood estimates under the assumption of independence between the training cases.
The optimization can be carried out numerically using the method of scoring, for which we give details in Appendix B. Approximate standard errors for the estimates can be obtained in the usual way, by evaluating and inverting the Hessian matrix for the mean logarithmic score or log likelihood function. However, the estimates of the weights w 1 , . . . , w k need to be nonnegative. Thus, if unconstrained optimization results in negative weights, we turn to the active barrier algorithm implemented in the constrained optimization routine constrOptim in R (R Development Core Team 2011). Similarly, linear, generalized linear and spreadadjusted linear combination formulas can be fitted by maximizing the mean logarithmic score over training data. For reasons of simplicity and tradition, we frequently refer to the resulting estimates as maximum likelihood estimates.

Simulation and data examples
We now illustrate and complement our theoretical results in simulation and data examples on density forecasts. This corresponds to the prediction space setting, where the CDFvalued random quantities F 1 , . . . , F k are absolutely continuous almost surely, and thus can be identified with random Lebesgue densities f 1 , . . . , f k . Throughout the section, we fit combination formulas by maximizing the mean logarithmic score over training data, in the ways described above and in Appendix B. To lighten the notation, we use the acronyms PIT, TLP, SLP and BLP to refer to the probability integral transform and the traditional, spread-adjusted and beta transformed linear pool, respectively.
The recent work of Ranjan and Gneiting (2010) and Clements and Harvey (2011) contains a wealth of simulation and data examples on the combination of probability forecasts for a binary event. In the concluding Section 5 we summarize these experiences and relate them to the findings in the case studies hereinafter.

Simulation example
In this simulation example, the data generating process for the observation, Y , is the regression model Y = X 0 + a 1 X 1 + a 2 X 2 + a 3 X 3 + ǫ, where a 1 , a 2 and a 3 are real constants, and X 0 , X 1 , X 2 , X 3 and ǫ are independent, standard normal random variables. The individual predictive distributions rest on partial knowledge of the data generating process, in that density forecast f 1 has access to the covariates X 0     and X 1 , but not to X 2 or X 3 , and similarly for f 2 and f 3 . Thus, we seek to combine the density forecasts f 1 = N (X 0 + a 1 X 1 , 1 + a 2 2 + a 2 3 ), f 2 = N (X 0 + a 2 X 2 , 1 + a 2 1 + a 2 3 ) and f 3 = N (X 0 + a 3 X 3 , 1 + a 2 1 + a 2 2 ), where X 0 stands for shared, public information, while X 1 , X 2 and X 3 represent proprietary information sets. The density forecasts represent the true conditional distributions under the regression model (11), given the corresponding partial information, as represented by the σ-algebras A 1 = σ(X 0 , X 1 ), A 2 = σ(X 0 , X 2 ) and A 3 = σ(X 0 , X 3 ), respectively. Hence, the forecasts are ideal in the sense of Definition 2.2, and by Theorem 2.9 they are both probabilistically calibrated and marginally calibrated.
We estimate the TLP, SLP and BLP combination formulas on a simple random sample {(f 1j , f 2j , f 3j , Y j ) : j = 1, . . . , J} of size J = 500 from the joint distribution of the forecasts and the observation, and evaluate on an independent test sample of the same size. The regression coefficients in the data generating model (11) are taken to be a 1 = a 2 = 1 and a 3 = 1.1, so that f 3 is a more concentrated, sharper density forecast than f 1 and f 2 . Table 5 shows maximum likelihood estimates, along with approximate standard errors, for TLP, SLP and BLP combination formulas. For all three methods, the weight estimate is highest for f 3 , whereas the estimates for f 1 and f 2 are smaller and not significantly different from each other. The SLP spread adjustment parameter c is estimated at 0.78, and the BLP transformation parameters α and β at 1.49 and 1.44, respectively.
The PIT histograms for the various types of density forecasts over the test set are displayed in Figure 3, with complementary results shown in Table 6. In addition to the variance of the PIT, which is our standard measure of dispersion, the table quantifies sharpness in terms of the root mean variance (RMV), that is, the square root of the average of the variance of the predictive density over the evaluation set. The component forecasts f 1 , f 2 and f 3 are probabilistically calibrated and thus show uniform empirical PIT histograms, up to sample fluctuations. As mandated by Theorem 3.1, the linearly combined TLP density forecast is overdispersed and lacks sharpness. The SLP and BLP agrregated density forecasts show nearly uniform PIT histograms; they are approximately neutrally dispersed and much sharper than their competitors. Table 7 shows the mean logarithmic score for the various types of density forecasts. The best individual density forecast is f 3 , because it is sharper than f 1 and f 2 . The linearly combined density forecast outperforms the individual density forecasts, even though it is overdispersed. The nonlinearly aggregated SLP and BLP density forecasts show higher scores than any of the individual or linearly combined forecasts, both for the training data, where this is trivially true, as the nonlinear methods nest the traditional linear pool, and for the test data, where such cannot be guaranteed.

Density forecasts for daily maximum temperature at Seattle-Tacoma Airport
With estimates of some one-third of the economy, as well as much of human activity in general, being weather sensitive (Dutton 2002), there is a critical need for calibrated and sharp probabilistic weather forecasts, to allow for optimal decision making under inherent environmental uncertainty. In practice, probabilistic weather forecasts rely on ensemble prediction systems. An ensemble system comprises multiple runs of a numerical weather prediction model, with the runs differing in the initial conditions and/or the details of the mathematical representation of the atmosphere (Palmer 2002;Gneiting and Raftery 2005). Here we consider two-days ahead forecasts of daily maximum temperature at Seattle-Tacoma Airport, based on the University of Washington Mesoscale Ensemble (Eckel and Mass 2005), which employs a regional numerical weather prediction model over the Pacific Northwest, with initial and lateral boundary conditions supplied by eight distinct weather centers. A brief description of the ensemble members is given in Table 8.
Our training period ranges from January 1, 2006 to August 12, 2007, with a few days missing in the data record, for a total of 500 training cases. The test period extends from August 13, 2007to June 30, 2009, for a total of 559 cases.
Each ensemble member is a point forecast, which can be viewed as the most extreme form of an underdispersed density forecast. To address the underdispersion and obtain approximately neutrally dispersed components, we use the maximum likelihood method on the training data to fit, for each ensemble member i = 1, . . . , 8 individually, a Gaussian predictive density of the form . Here x ij is the point forecast from the ith ensemble member on day j, a i and b i are member specific linear bias correction parameters, and σ i is the member specific predictive standard deviation. From Table 9 we see that the estimates for σ 1 , . . . , σ 8 range from 1.958 to 2.214.
Next we combine the eight individual density forecasts. Table 10 shows maximum likelihood estimates for TLP, SLP and BLP combination formulas. For all three methods, the GFS member, f 1 , obtains the highest weight and the ETA member, f 3 , the lowest weight. This can readily be explained, in that both members have a common institutional origin, and thus are highly correlated, whence the more competitive GFS member subsumes the  (Eckel and Mass 2005), with member acronyms and organizational sources for initial and lateral boundary conditions. The United States National Centers for Environmental Prediction supply two distinct sets of initial and lateral boundary conditions, namely, from its Global Forecast System (GFS) and Limited-Area Mesoscale Model (ETA).    Raftery et al. (2005).  weight of the ETA member. The SLP spread adjustment parameter is estimated at 0.768, and the BLP transformation parameters both at 1.467. Figure 4 illustrates the various density forecasts for June 28, 2008, an unusually hot day at Seattle-Tacoma Airport with a verifying maximum temperature of 32.8 degrees Celsius or 91 degrees Fahrenheit. The member specific individual density forecasts are shown by the dotted lines, and the linearly combined TLP forecast by the dash-dotted line. The nonlinearly aggregated SLP and BLP density forecasts, which are shown by the solid and dashed line, respectively, are sharper than the TLP density.
PIT histograms for the test period are shown in Figure 5, along with summary measures of dispersion and sharpness in Table 11. The individual, member specific density forecasts tend to be a bit overdispersed. The linearly aggregated TLP density forecast is much more severely overdispersed, as reflected by an inverse U-shaped and skewed PIT histogram. Of course, the overdispersion is not surprising, as it is a direct consequence of Theorem 3.1. The SLP and BLP aggregated density forecasts show somewhat rough and skewed, yet more nearly uniform PIT histograms.
These results are corroborated by Table 12, which shows the mean logarithmic score for the various types of density forecasts, both for the training period and the test period.
The linearly combined TLP forecast shows a higher score than any of the individual density forecasts, which attests to the benefits of aggregation. Nevertheless, the linearly combined density forecast is suboptimal, because it is overdispersed and lacks sharpness, and thus it is outperformed by the nonlinearly aggregated SLP and BLP density forecasts.  nique, which is a state of the art approach to generating density forecasts from forecast ensembles. The BMA density forecast for day j is of the form with BMA weights, w 1 , . . . , w 8 , that are nonnegative and sum to 1, member specific bias parameters a i and b i for i = 1, . . . , 8, and a common variance parameter, σ 2 . In view of our individual density forecasts being Gaussian, the TLP and BMA densities are of the same functional form. However, there is a conceptual difference, in that the TLP weights are fitted conditionally on the individual density forecasts. Thus, a two-stage procedure is used, in which the member specific component densities are estimated first, and only then the weights, with the components held fixed. In contrast, the BMA method estimates the weights, w 1 , . . . , w 8 , and the common spread parameter, σ, for the component forecasts in the Gaussian mixture model (12) simultaneously. While the BMA method can be employed with member specific spread parameters, the assumption of a common spread parameter stabilizes the estimation algorithm and does not appreciably deteriorate the predictive performance . Table 10 shows maximum likelihood estimates for the BMA parameters, obtained with the R package ensembleBMA (Fraley et al. 2011). The BMA weights echo the SLP weights. The BMA spread parameter σ is estimated at 1.566 and differs from the predictive standard deviations for the member specific density forecasts in Table 9 by factors ranging from 0.707 to 0.800, much in line with our estimate of 0.768 for the SLP spread adjustment parameter, c. Thus, the SLP and BMA density forecasts are very much alike, which is confirmed by the PIT histograms in Figure 5, the summary measures in Table 11 and the logarithmic scores in Table 12. In Figure 4 the graphs for the SLP and BMA density forecasts are nearly identical and lie essentially on top of each other, and so we refrain from plotting the BMA density.

Density forecasts for S&P 500 returns
In this final data example, we follow Diebold et al. (1998) in considering S&P 500 log returns for the period of July 3, 1962 to December 29, 1995. The data record through December 1978 is used as training set, for a total of 4,133 training cases. All estimates reported are maximum likelihood fits on the training period obtained with the R package fGarch (Wuertz and Rmetrics Core Team 2007). The balance of the record is used as test period, for a total of 4,298 one-day ahead density forecasts.
The first component forecast, f 1 , is based on a generalized autoregressive conditional heteroscedasticity (GARCH; Bollerslev 1986) specification for the variance structure. With r t denoting the log return on day t, our GARCH(1,1) model assumes that r t = σ t ǫ t , where ǫ t is Student-t distributed with ν degrees of freedom and variance 1, while σ t evolves dynamically as σ 2 t = ω + αr 2 t−1 + βσ 2 t−1 . The maximum likelihood estimates for the GARCH parameters are ω = 0.000, α = 0.089, β = 0.903 and ν = 9.25.
The second component forecast, f 2 , is based on a standard moving average (MA) model for the mean dynamics, which assumes that r t = Z t + θZ t−1 , where {Z t } is a Gaussian white noise process with mean zero and variance σ 2 . The maximum likelihood estimates for the MA parameters are θ = 0.252 and σ = 0.00736.
Our goal now is to combine the density forecasts f 1 and f 2 . Table 13 shows maximum likelihood estimates for TLP, SLP and BLP combination formulas. For all three methods, the conditionally heteroscedastic density forecast f 1 obtains a much higher weight than the simplistic density forecast f 2 . The SLP spread adjustment parameter is estimated at 0.940, and the BLP transformation parameters α and β at 1.100 and 1.081. This suggests that the overdispersion of the TLP density forecast is quite mild, which is confirmed by the PIT histogram in Figure 6 and the summary measures in Table 14.  Table 15 shows the mean logarithmic score for the various types of probabilistic forecasts. The TLP density forecast performs slightly better than the individual component f 1 , with a score that is very slightly lower than for the nonlinearly aggregated SLP and BLP density forecasts, both for the training and the test period. As observed by Geweke and Amisano (2011), there is little reward for using more elaborate, less parsimonious aggregation methods for density forecasts of S&P 500 returns. 3 Finally, we consider the predictive performance of a more comprehensive predictive model, which addresses both the first and the second order dynamics, in that r t = µ t + ǫ t where {µ t } and {ǫ t } are MA(1) and Student-t GARCH(1,1) processes, respectively. The maximum likelihood estimates in this mixed specification are θ = 0.269 and σ = 0.00736 for the MA parameters, and ω = 0.000, α = 0.098, β = 0.892 and ν = 8.284 for the GARCH parameters. The resulting density forecast can be thought of as combining information sets  with respect to the first and second order dynamics, as opposed to combining the corresponding component forecasts f 1 and f 2 . It outperforms the other types of density forecasts and achieves a mean logarithmic score of 3.638 for the training period and 3.473 for the test period.

Discussion
We have studied methods for combining predictive distributions. From a theoretical perspective, our approach deviates from previous work in major ways. Technically, we operate in terms of prediction spaces and cumulative distribution functions, which allows for a unified treatment of all real-valued predictands including, for example, density forecasts for continuous variables and probability forecasts of a binary event. Conceptually, our work is motivated by applications in probabilistic forecasting, and thus we assess combination formulas and aggregation methods in terms of coherency, calibration, and dispersion.
While it is not difficult to construct coherent combination formulas, the formulas typically used in practice tend to be incoherent. Our analytical results about the linear pool extend  and unify extant work by Dawid et al. (1995), Hora (2004) and Ranjan and Gneiting (2010), and show that any linear combination formula with strictly positive coefficients fails to be coherent. In this sense, combined evidence is nonlinear. An interesting open problem is whether or not parsimonious ramifications, such as spread-adjusted or beta transformed linear combination formulas, are coherent. That said, the applied relevance of the notion of coherency is limited, as predictive distributions or expert opinions issued in practice are hardly ever ideal. In typical practice, underdispersed or approximately neutrally dispersed predictive distributions are to be aggregated. In the case of underdispersed components, the tendency of linear combination formulas to increase dispersion can be beneficial, and helps explain the success of linear pooling in applications (Madigan and Raftery 1994). However, if the components are neutrally dispersed, the failure of the traditional linear pool to be flexibly dispersive is a serious limitation. Berrocal et al. (2007), Glahn et al. (2009) and Kleiber et al. (2011) observed this empirically in the context of probabilistic weather forecasts, and proposed a special case of the spread-adjusted linear pool as an ad hoc remedy. Our theoretical results document the increased flexibility of the spread-adjusted linear pool, and demonstrate that the beta-transformed linear pool is exchangeably flexibly dispersive.
Not surprisingly, the parsimonity principle and the bias-variance tradeoff apply in the practice of the combination of predictive distributions. Thus, in data poor settings, where training data are scarce, the parsimonious traditional linear pool might be the method of choice, despite its theoretical shortcomings, as demonstrated persuasively in the recent simulation study of Clements and Harvey (2011). In data rich settings, where predictive models can reliably be estimated, linear aggregation tends to be suboptimal. Hence, we have studied parsimonious nonlinear alternatives, including the spread-adjusted linear pool (SLP) and the beta transformed linear pool (BLP). Further options include consensus methods (Winkler 1968) and nonparametric approaches, such as isotonic recursive partitioning (Luss, Rosset and Shahar 2011). As Winkler (1986) noted, "different combining rules are suitable for different situations".
The SLP and BLP approaches can also be used to provide calibration and dispersion corrections to a single predictive distribution, similar to the methods described by Cox (1958), Platt (1999), Zadrozny and Elkan (2002) and Primo et al. (2009) in the context of probability forecasts for a binary event. An interesting question then is whether dispersion adjustments ought to be applied to the individual components prior to the aggregation. In situations in which the components show substantially differing degrees of dispersion, or are uniformly under-or overdispersed, we indeed see potential benefits in doing this, with (here unreported) simulation experiments providing partial support to this view. In our temperature example, the components derive from point forecasts, which is the most extreme form of underdispersion, and prior to aggregating the components we apply a simple Gaussian technique, which obtains approximately neutrally dispersed individual density forecasts.
In addition to their relevance in probabilistic forecasting, our findings bear on the related problem of the fusion of expert opinions that are expressed in the form of probability distributions. Ha-Duong (2008) reviews methods for doing this, and applies them to combine expert opinions about the climate sensitivity constant, which is a key quantity in the study of the greenhouse effect. Our analytic results imply that if each individual expert is neutrally dispersed, linear aggregation results in combined assessments that are underconfident and show an unduly wide range of uncertainty, when in fact a sharper assessment could be made. In situations of this type, the parsimonious SLP technique with its single, easily interpretable spread-adjustment parameter, c, might be preferable to the BLP, particularly when the individual components are symmetric and unimodal. Our theoretical results show that neutrally dispersed components require values of c < 1, with more specific recommendations depending on the subject matter and situation at hand. The case c = 1 corresponds to the traditional linear pool and can be a useful choice in situations in which overconfident experts make underdispersed judgements.