Source-linked AI summary
The AEP algorithm for the fast computation of the distribution of the sum of dependent random variables
Philipp Arbenz, Paul Embrechts, Giovanni Puccetti
TL;DR
Computing the distribution of sums of dependent random variables is difficult even when their joint distribution is known, often requiring problem-specific simulation methods. The paper introduces the AEP algorithm, which numerically calculates the distribution at a fixed threshold using only the joint distribution, with convergence under weak assumptions but exponential memory and runtime growth in iteration count.
Problem
Exact calculation of the distribution of sums of dependent random variables is onerous and often requires expertise and problem-specific Monte Carlo or quasi-Monte Carlo methods.
Method
The AEP algorithm numerically computes the distribution function at a fixed threshold using only the joint distribution, with approximations constructed through simplex decompositions and converging under weak assumptions.
Results
AEP shows good convergence across studied dimensions and thresholds; in two dimensions, extrapolation increases accuracy by two decimal digits without increasing computational time.
Takeaways & Limitations
AEP can support numerical distribution and quantile-range calculations for dependent portfolio losses, although high-level quantile estimation is insufficient in dimensions d = 4,5.
Takeaways & Limitations
Memory required for the nth iteration increases exponentially in n, making estimates for d ≥6 prohibitive and limiting higher-dimensional applications.
Abstract
from arXiv · showhide
We propose a new algorithm to compute numerically the distribution function of the sum of $d$ dependent, non-negative random variables with given joint distribution.
1. Motivations and preliminaries
The paper introduces AEP to numerically compute the distribution function of a sum of dependent random variables from their joint distribution, addressing a task that is otherwise onerous and often requires tailored simulation methods. The problem is relevant to risk-management applications such as capital requirements and stress-testing.
- Exact computation of the sum of dependent random variables is onerous even when their joint distribution is known.Existing approaches may rely on Monte Carlo or quasi-Monte Carlo methods tailored to the specific problem.
- AEP numerically calculates the distribution at a fixed threshold s using only the joint distribution H, without problem-specific adaptation.
- Computing the distribution of portfolio losses supports VaR-based capital requirements under risk-management applications.The sum S_d represents overall portfolio risk, and its quantile determines the associated capital requirement.
- Stress-testing can use AEP when marginal distributions are known but dependence must be varied through a copula.
- The paper focuses on non-negative random variables with known joint distribution H and introduces notation for hypercubes, simplex regions, and their measures.
2. Description of the AEP algorithm for d = 2
For d = 2, AEP approximates the simplex representing the event that the sum is at most s by recursively decomposing it into simplexes and approximating their measures with hypercubes. The construction uses signed additions and subtractions, with α controlling the decomposition.
- The method is illustrated under the paper’s assumption of non-negative random variables and a fixed positive threshold s.
- AEP approximates the simplex S(0,s) by hypercubes whose VH-measures are straightforward to compute.
- For d = 2, the initial decomposition expresses the simplex measure through three component simplexes and a subtraction term.
- Each iteration decomposes the remaining simplexes, producing a refined approximation whose signed terms track whether measures are added or subtracted.At the second iteration, nine simplexes form the error representation and are passed to the next iteration.
- The decomposition holds for every α in [1/d,1), but Equation (2.2) depends on α and the simplexes are generally non-disjoint for d > 2.
3. Description of the AEP algorithm for arbitrary d
In arbitrary dimension d, AEP recursively generates signed simplexes and estimates their joint-distribution measures through hypercubes. Under weak assumptions, the resulting sequence converges to the target simplex measure, with convergence ensured when the volume factor is below one.
- At iteration n, AEP receives N^(n−1) simplexes, decomposes each into N new simplexes, and passes them to the next iteration.
- The sequence P_n(s) sums hypercube VH-measures weighted by signed coefficients to approximate VH[S(0,s)].
- Under weak assumptions on H, P_n(s) converges to VH[S(0,s)], and its terms can be calculated from hypercube measures.
- If VH has a bounded density and the volume factor satisfies f(α) < 1, the AEP sequence converges.
- The approximation error decreases through the residual simplex volumes, which shrink exponentially when f(α) < 1.
- For α in (0,2/d), simplex distance from Γ_s decreases geometrically, so bounded density is needed only near Γ_s.
4. Choice of α
The paper recommends choosing α∗ to minimize the distribution-independent volume factor f(α), yielding computational and convergence benefits. This choice reduces generated simplexes in odd dimensions and supports convergence under stated assumptions, but memory grows exponentially with iteration.
- Selecting α∗: α∗ minimizes f(α), providing a distribution-independent choice of the algorithm parameter.The authors suggest this choice because the algorithm should remain independent of the distribution H.
- Computational efficiency: For odd dimensions, α∗ makes d(d + 1)/2 simplex coefficients zero and therefore reduces computational work.These simplexes can be omitted from the decomposition.
- Computational efficiency: In dimension d = 3, α∗ generates 4 new simplexes per iteration instead of 7 for other feasible α values.The reduction follows from the zero coefficients produced by α∗.
- Convergence: α∗ also minimizes γ; when α∗ ∈ (0,2/d), γ < 1 and γ^n s tends to zero as n increases.Under these conditions, bounded density is required only near Γs for convergence.
- Convergence: Theorem 3.2 establishes convergence of Pn(s) for d ≤5, while financial and actuarial examples generally show good performance.Slow convergence is mainly associated with more pathological cases and may occur only at a few thresholds s.
- Computational limits: Memory requirements increase exponentially with iteration, and numerical complexity becomes substantial for d ≥6.For d ≥6, the computation can quickly exhaust the memory of a standard computer.
5. An improvement of the numerical accuracy of the algorithm via extrapolation
The extrapolation technique improves AEP accuracy by using a linear Taylor approximation of a smooth joint density, yielding faster convergence and applicability in higher dimensions.
- Method: AEP's extrapolation technique uses a linear Taylor expansion of a twice continuously differentiable joint density with bounded derivatives to estimate approximation error.The method provides a local estimator of simplex mass from the corresponding hypercube volume after several AEP iterations.
- Choice of α: α* = 2/(d+1) makes the simplex-to-hypercube mass ratio independent of b, s, a, and the linear-density coefficients c_k.This choice motivates the extrapolation construction and also minimizes the contraction factor γ.
- Convergence: The extrapolated estimator P*_n(s) converges to V_H[S(0,s)] faster, and in higher dimensions, than the original estimator P_n(s).Under the stated smoothness assumptions, convergence is established for d ≤ 8 because the extrapolation error factor is less than 1.
- Convergence: For d ≤ 8, the extrapolation error converges to zero under the theorem's smoothness assumptions, while convergence is not assured for d = 9.The error-ratio and simplex-growth quantities are summarized in Table 2.
- Computational cost: The extrapolation error admits an upper bound as a function of the number of evaluations, with a polynomial convergence rate depending only on dimension d.The paper presents these rates as a basis for comparing AEP efficiency with Monte Carlo methods.
6. Applications
The applications evaluate AEP on dependent Pareto portfolios coupled by Clayton and Gumbel copulas, reporting convergence, extrapolation gains, computational cost, and practical limits.
- Applications: For a two-dimensional portfolio, iteration n = 7 already preserves the first eight decimal digits, with a computation time of 0.01 seconds.Extrapolation increases accuracy by two decimal digits on average without increasing computational time.
- Applications: AEP shows good convergence across the studied dimensions and thresholds, while extrapolation provides additional accuracy in higher dimensions.Memory constraints make estimates for d ≥6 prohibitive.
- Applications: For dimensions 2 ≤d ≤5, average computation time increases exponentially with the number of iterations n.The timing study uses Pareto portfolios coupled by a Clayton copula.
- Applications: The extrapolated estimator P*_n(s) provides accurate estimates within a reasonable computational time for the tested Clayton–Pareto portfolios.Tables 3–6 report convergence across dimensions d = 2 to 5.
- Applications: The extrapolation technique increases accuracy for two- and three-dimensional independent Pareto portfolios under a smooth model for H.The authors therefore prefer P*_n(s) over P_n(s) in this setting.
- Applications: AEP can compute numerical VaR ranges for sums of four and five variables, but its accuracy is insufficient for high-level quantiles in dimensions d = 4,5.The resulting numerical error is described as extremely small relative to statistical inference error.
7. A comparison with Monte Carlo, quasi-Monte Carlo and quadrature methods
The paper compares deterministic AEP with Monte Carlo, quasi-Monte Carlo, and quadrature methods, emphasizing different convergence behavior, assumptions, and computational trade-offs.
- Comparison methods: Traditional MC converges at O(M^-1/2), while QMC can approach O(M^-1) in optimal cases but has dimension-dependent worst-case rates.For smooth densities, QMC convergence is at least O((log M)^d M^-3/2).
- Comparison methods: AEP is deterministic and does not require analytic density formulas, overall smoothness, or adaptation to the specific probabilistic model.This contrasts with the randomness and model-specific adaptation used by MC and randomized QMC methods.
- Comparison methods: AEP’s behavior is not affected by the threshold s or marginal tail properties, whereas MC and QMC relative error increases in the tail under heavy-tailedness.This makes the threshold dependence a central distinction between the approaches.
- Comparison methods: Well-designed quadrature rules or specific quasi-random sequences may outperform AEP in particular examples in both accuracy and computational effort.The paper nevertheless presents AEP as competitive under weak assumptions.
- Comparison methods: AEP provides accurate estimates for sums up to five dimensions in reasonable time without adapting the computational rule to each example.The method can handle joint distributions represented through copulas and marginal distributions.
8. Final remarks
The final remarks state AEP’s convergence conditions and practical scope, while identifying unresolved higher-dimensional geometry and future extensions.
- Final remarks: AEP may fail to converge when VH[Γs] > 0; for point masses at 1/2, the sequence P_n(1) alternates between 0 and 1.The convergence guarantee therefore excludes such discontinuous cases.
- Final remarks: If H has a bounded density near Γs, P_n(s) converges to VH[S(0,s)], although convergence speed depends on probability mass near Γs.The paper motivates efficiency improvements for cases with slower convergence.
- Final remarks: AEP converges for d ≤5 with bounded density and for d ≤8 under extra smoothness assumptions, with conditions needed only near Γs.The authors state that these conditions cover most practically relevant examples.
- Final remarks: Convergence in arbitrary dimensions remains unproved because P_n(s) and P*_n(s) can be non-monotone and generated simplexes can overlap for d > 2.For d > 6, the sum of simplex volumes increases with iterations even though their union remains near Γs.
- Final remarks: Memory constraints limit practical use of AEP to dimension d ≤5, so a general convergence theorem would not itself improve current practice.The authors nevertheless expect better convergence rates in these manageable dimensions than the upper bounds in Table 2.
- Final remarks: Future work will study more general aggregating functions and adaptive decompositions that reduce the number of new simplexes per iteration.The proposed extensions target both scope and efficiency.
Appendix: Proof of (3.1)
The appendix proves the simplex decomposition underlying AEP by establishing set identities for shifted simplexes and hypercubes, then combining them algebraically.
- Appendix: Proof of (3.1): The proof introduces binary index vectors and the Kronecker delta to describe coordinate shifts in d dimensions.The notation includes all vectors in {0,1}^d and the canonical basis.
- Appendix: Proof of (3.1): The resulting identities cover the cases α = 1/d and α > 1/d, with transformed parameters used when the latter requires a complementary representation.The final substitution identifies equation (A.1) with equation (3.1).
- Appendix: Proof of (3.1): The proof decomposes set differences into unions and intersections of shifted simplexes, using disjointness where the set algebra permits it.Equations (A.7)–(A.12) assemble the component identities into the target decomposition.