{'fraction_non_alphanumeric': 0.0400650738460557, 'fraction_numerical': 0.011132731506587131, 'mean_word_length': 4.325293018515373, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On spherical Monte Carlo simulations for multivariate normal probabilities

The calculation of multivariate normal probabilities is of great importance in many statistical and economic applications. This paper proposes a spherical Monte Carlo method with both theoretical analysis and numerical simulation. First, the multivariate normal probability is rewritten via an inner radial integral and an outer spherical integral by the spherical transformation. For the outer spherical integral, we apply an integration rule by randomly rotating a predetermined set of well-located points. To find the desired set, we derive an upper bound for the variance of the Monte Carlo estimator and propose a set which is related to the kissing number problem in sphere packings. For the inner radial integral, we employ the idea of antithetic variates and identify certain conditions so that variance reduction is guaranteed. Extensive Monte Carlo experiments on some probabilities calculation confirm these claims.


Introduction
Efficient and precise calculation for the multivariate normal probability is of critical importance in many disciplines. To name a few, the multinomial probit model used in econometrics and biometrics has cell probabilities that are negative orthant probabilities. In financial industry implementation of the CreditMetrics model, the joint migration probabilities in credit migration model are rectangle probabilities for bivariate normal distributions, cf. Gupton et al. (1997). The estimation of Value-at-Risk for risk management considered in Glasserman et al. (2000) requires calculation of multivariate probability of an ellipsoid. In multiple comparisons, multivariate normal probabilities are also considered in Hsu (1996). For more examples and applications, the reader is referred to Genz and Bretz (2009) for a recent summary of multivariate normal distribution and multivariate t distribution.
Motivated by these applications, we investigate efficient calculation for the multivariate normal probability. That is, for an indicator function I A (x) with a support set A in R d , we seek more efficient computation of the following probability where ′ denotes the transpose, X = (X 1 , · · · , X d ) ′ is a d-dimensional nonsingular multivariate normal distribution with mean vector µ and covariance matrix Σ, denoted by X ∼ N d (µ, Σ), and |Σ| denotes the determinant of Σ.
Standard approaches of calculating (1) include classical analytic approximation, numerical integration and Monte Carlo method. Instead of using analytic approximation and numerical method, cf. Miwa et al. (2003) and Craig (2008), which are usually more suitable for low dimensional cases, in this paper, we study Monte Carlo method. Although Monte Carlo method is easy to implement and can overcome the curse of dimensionality, its convergence rate is rather slow (proportional to 1/ √ d). Therefore, additional variance reduction methods are required. Typical methods for variance reduction include antithetic variates, Latin hypercube sampling, and primitive Monte Carlo method, cf. Genz (1992), Genz (1993), Hajivassiliou et al. (1996), Vijverberg (1997), Genz and Bretz (2002), and among others. A comparison study of alternative sampling methods can be found in Sándor and András (2004). A survey on existing methods is in Genz and Bretz (2009).
Monte Carlo methods using spherical transformation have been studied in the literature. For instance, Deák (1980) and Deák (2000) used this transformation as the basis for calculating multivariate normal probabilities, and Fang and Wang (1994) proposed a transformation on the unit hypercube to generate points uniformly distributed on the sphere. Monahan and Genz (1997) proposed a Monte Carlo simulation method for Bayesian computation, to which the authors used randomized extended simplex design for the spherical integral and Simpson weights for the radial integral.
In this paper, we propose a spherical Monte Carlo method with both theoretical analysis and numerical simulations. There are two aspects in this study. First, because the spherical integral requires generating unit vectors uniformly on the sphere, one way to improve the Monte Carlo efficiency is to use a randomly rotated predetermined set of unit vectors. For this purpose, we give a criterion to select such a set on the unit sphere, which involves the minimal distance of any two points in the set and the cardinality of the set. Among all sets with the same minimal distance, the desired optimal set is the one with maximal cardinality. Especially, when the minimal distance equals one, finding the desired optimal set is linked to the kissing number problem in sphere packings.
Next, for the radial integral, we offer a variance reduction technique employing the idea of antithetic variates. For this purpose, we introduce the idea of central symmetry for the set of basis points and central antisymmetry for the event of simulation. To the best of our knowledge, this seems to be a first step to provide sufficient conditions for antithetic variates on spheres.
To illustrate the proposed method, simulation studies are given for orthants, rectangles, and ellipsoids probabilities for multivariate normal distributions.
The simulation results confirm these claims.
The rest of this paper is organized as follows. In §2 we propose the spherical Monte Carlo method with antithetic variates. §3 links the proposed set of basis points to a sphere packing problem and related spherical t-designs, and discusses practical implementation for high dimensional cases. Simulation results are given in §4. §5 concludes. The proofs are deferred to Appendix.

The proposed spherical Monte Carlo method
For easy presentation, we split this section as three subsections. §2.1 formulates the problem, §2.2 presents the spherical integral, and §2.3 studies the radial integral.

Problem formulation
Let I A : R d → R be an indicator function with a support set A. Denote the probability density function (pdf) of a multivariate normal random where µ is the mean vector and Σ is the variance-covariance matrix. The desired multivariate normal probability of the region A is an integral of the form in (1).
Note that a d-variate normal random variable X ∼ N d (µ, Σ) can be expressed by X = µ + ΓZ, where Z is a d-variate standard normal random variable, and Γ is the lower triangle matrix such that Σ = ΓΓ ′ , the so-called Cholesky decomposition of Σ. By change of variables, x = µ + Γz, we rewrite (1) as where φ(·) is the standard normal density andÃ = Γ −1 (A − µ). Abuse the notation a little bit, we denoteÃ by A in the rest of the paper.
By using the spherical transformation, a point z ∈ R d can be written as (r, u), where r is the radius and u is the unit vector of z. Then, I A (z) = I A (r, u) and (2) equals where dA denotes the differential area on the unit sphere S d−1 , and the inner radial integral is Here, Area is the surface area of S d−1 , and k d (·) is the pdf of a χ-distribution with degrees of freedom d, denoted by χ(d). A random unit vector u in R d is a vector uniformly distributed on S d−1 , denoted by u ∼ U(S d−1 ). Therefore, generating a sample from a standard normal distribution z is equivalent to independently generating a radius r ∼ χ(d) and a unit vector u ∼ U(S d−1 ), and setting z = ru.
To rewrite (2) using the spherical transformation, an alternative approach is to consider the spherical integral as the innermost integral as in Monahan and Genz (1997). Here, we take the radial integral as the innermost integral because the radial integral is of one dimension and its calculation is simple. For a simple region A like rectangles, orthants and ellipsoids, when the unit vector u is fixed, the inner radial integral has a closed-form formula in terms of cumulative distribution function (cdf) of a χ-distribution, cf. Deák (2000). For general regions, the inner radial integral can be calculated via Monte Carlo methods or numerical methods, cf. Davis and Rabinowitz (1984). Numerical quadrature methods produce biased estimators in general and are out of the scope of this paper. Here, we only focus on Monte Carlo methods.
The outer spherical integral using Monte Carlo method requires generating a random unit vector u. A straightforward method to generate u is to generate d independent standard normal random variables to have a vector in R d and then normalize the vector by its length. A more efficient algorithm can be found in Fang and Wang (1994) by generating just (d − 1) random numbers to get u.
The crude spherical Monte Carlo estimator iŝ where r ∼ χ(d), u ∼ U(S d−1 ), and r and u are independent. The crude spherical Monte Carlo estimator with antithetic variates iŝ

Variance reduction on the spherical integral
One way to obtain an efficient spherical Monte Carlo estimator for the spherical integral is taking the average value of a randomly rotated predeter- The standard algorithm to generate a random orthogonal matrix can be described as follows, cf. Heiberger (1978). First, generate a random d × d matrix whose entries are independently standard normal random variables.
Then, apply Gram-Schmidt method to the column vectors and obtain the desired random orthogonal matrix.
One important feature of random orthogonal matrix is that if we fix a unit vector and let a random orthogonal matrix act on it, the resulting vector is uniformly distributed on S d−1 . In other words, for a continuous function h(u) on S d−1 , and any unit vector v ∈ S d−1 , we have Here dT is the unique left-invariant probability measure on O(d) mentioned above.
A spherical Monte Carlo estimator using V iŝ where r v are independent random variables with χ(d) distribution, and T ∼ U(O(d)) which is independent of r v . When the innermost radial integral can be calculated explicitly, i.e., f A (u) can be expressed in terms of the cdf of χ(d), an estimator using V iŝ where T ∼ U(O(d)). By using the fact of conditioning, indeed, f A (T v) equals , and thus enjoys smaller variance compared withp V .
To have efficient simulation, the crucial step involves the selection of V to minimize the variance ofp V . For this purpose, we propose a criterion of determining the desired set V in the following, and provide a feasible solution that related to the maximal kissing number problem in sphere packings in §3.
To begin with, for a finite subset V of the unit sphere S d−1 , to characterize the variance of the estimator for the spherical integral, we need to define some notations as follows. First, denote d(·, ·) as the usual Euclidean distance function on R d , and let be the minimal distance of any two points in V . For a region A in S d−1 , define the diameter of A as It is easy to see that when Diam(A) < d min (V ), the intersection of V and A is either the empty set or a single point.
Let π be the normalized measure on S d−1 induced from the Lebesgue measure on R d so that π(S d−1 ) = 1. Note that π is a probability measure where Since orthogonal matrices preserve the Euclidean distance function, we have d min (T V ) = d min (V ) and |A ∩ T V | = 1 or 0. On the other hand, the estimator g V A (T ) is unbiased. Therefore, we have If Diam(A) ≥ d min (V ), we can decompose A and obtain the following The proof of Theorem 2.2.1 will be given in Appendix A.
When N is the minimal number such that the decomposition condition in Theorem 2.2.1 holds, denote the upper bound in Theorem 2.2.1 by C(V, A) for short. Note that C(V, A) reduces to (11), the exact variance, when if one of the following two conditions holds: Although there is no general method to construct such set V in either case, useful results are available. Other than the construction of V based on Case 1 1 , here we provide a solution based on Case 2. Note that in the case of d min (V ) = 2, V must have exactly two points, since the distance between any two unit vectors is less than or equal to two, and the equality holds only when they form an antipodal pair. This shows that the method of antithetic variates is optimal under the above criterion for d min (V ) = 2. For the case of d min (V ) = 1, constructing a set with maximal cardinality is related to the kissing number problem in sphere packings. Details of this linkage will be discussed in §3.1.

Variance reduction on the radial integral
The efficiency of simulators may be improved, for a given number of where , and r v 's and T are independent.
We remark that the proposed set V generated by the shortest nonzero vectors of a lattice as will be described in §3 is centrally symmetric so thatp V AT in (12) is well defined.
It is known that certain conditions are required in Cartesian coordinate to insure antithetic variates enjoy lower variance, cf. Ross (2006). Likewise, an additional condition is needed in the spherical case. A region A is called centrally antisymmetric if A and (−A) are disjoint. It is easy to see that the property of centrally antisymmetric is also preserved under the action of O(d). As a result, we have the following The proof of Theorem 2.3.1 will be given in Appendix B.
Theorem 2.3.1 shows that the variance ofp V AT can be reduced when A is centrally antisymmetric. As will be shown in §4.3, for certain central antisymmetric sets, although Var(p V AT ) = Var(p V ) as a mathematical fact,p V AT is preferred because it requires less number of independently drawn radii and thus enjoys lower computational cost compared withp V . For demonstration purposes, §4.3 provides an explicit central antisymmetric set, for whichp V AT enjoys smaller variance thanp V .

Practical implementation
This section first explores the link of an optimal V based on Case 2 in Corollary 2.2.2 to the kissing number problem in sphere packings, and then makes a connection to spherical t-designs. §3.2 discusses possible approaches for high-dimensional cases.

The linkage to the kissing number problem in sphere packings
For a sphere packing, a kissing number is defined as the number of nonoverlapping spheres that can be arranged such that they each touch another given sphere. Here all spheres must have the same radius. The kissing number problem seeks the maximal possible kissing number in a sphere packing.
Figures 1 and 2 demonstrate two sphere arrangements in R 2 with a kissing number equal to four and six, respectively. Moreover, the later one has the maximal kissing number. Now for a finite subset V of S d−1 with d min (V ) = 1, we can construct an associated sphere packing from V as follows. Consider all spheres centred at v ∈ V with radius 1/2 and one sphere centred at the original with radius 1/2.
It is easy to see that the number |V | equals the kissing number of this sphere packing, and, as a result, to find a V with maximal |V | with d min (V ) = 1 is exactly the same as finding a sphere packing producing the maximal kissing number.
A next question is regarding how to obtain a sphere arrangement with maximal kissing number. This question is quite difficult in geometry, but has been solved partially so that we can directly utilize existing results. For example, the integer lattice Z d is generated by the standard basis {e 1 , · · · , e d } and consists of all integer points. In this case, the kissing number of Z d equals 2d, and the set of shortest vectors Another example considers the lattice A 2 , which is generated by the basis v 1 = ( √ 3 2 , 1 2 ) and v 2 = (0, 1). The lattice and the associated sphere arrangement are depicted in Figure 2. Note that the kissing number of A 2 is six, which is indeed the maximal kissing number for d = 2. Moreover, we obtain It is difficult to prove if a kissing number is maximal for an arbitrary d, cf. Conway and Sloane (1999). For d from two to eight, lattices that produce the known maximal kissing number are the lattices A 2 , A 3 , D 4 , D 5 , E 6 , E 7 , and E 8 . These lattices are called root lattices which come from semisimple Lie algebras and these symbols stands for the name of the corresponding Lie algebra. Explicit construction of these lattices can be found in Conway and Sloane (1999). We can thus form the desired V from these lattices, by taking the normalized shortest vectors in these lattices. Table 1 lists the cardinality of each V L from the corresponding lattices for d from two to eight, 16 and 24.
It is worth mentioning that the set V L above also forms a spherical tdesign for some t > 0. That is for all real-valued continuous functions g(u) on the sphere which are restrictions of polynomial functions of degree less than or equal to t in R d . This |V L | g(T v) produces zero variance for estimating a polynomial function of degree less than or equal to t in R d , under the multivariate normal probability density and hence enjoys a small variance if g(·) can be well approximated by a polynomial of degree less than or equal to t in R d . Table 1 lists the associated t for the proposed V L as a spherical t-design.

High-dimensional cases
In principle, our proposed method is feasible for higher dimensions once a suitable V is selected. Theorem 2.2.1 suggests the use of V with the maximal cardinality given d min (V ) = 1, which is related to the maximal kissing number problem in sphere packings, and can be constructed via a lattice in some cases. The lattice with the maximal presently known kissing number, denoted by L max d , for dimension up to 40 can be found on the website http://www.math.rwth-aachen.de/∼Gabriele.Nebe/LATTICES/. For demonstration purposes, we report results for d from two to eight in §4.3.
Nevertheless, some difficulties arise to construct V from the lattice L max d : not all of the optimal lattices are well-studied lattices so that their shortest vectors can not be constructed generically. Moreover, the size of V may be too large for practical implementation in high-dimensional cases. For example, the maximal kissing number is 196,560 for d = 24, which seems be formidable for practical implementation.
For high-dimensional cases, a simple (but not optimal based on our criterion in Theorem 2.2.1) approach is to use a certain family of lattices. For example, we consider three lattices Z d , D d , and A d in this paper, which are defined as follows.
• Z d is generated by the standard basis e 1 , · · · , e d of R d and V Z d = {±e 1 , · · · , ±e d }.
Therefore, spherical packings generating by lattices Z d , A d and D d have kissing numbers of 2d, d(d + 1) and 2d(d − 1), respectively. An additional advantage of using these lattices is that the associated V of these three families of lattices can be constructed generically. However, as will be compared in §4.4, spherical estimators employing V based on these lattices produce larger variances. As a remark, the spherical estimatorp V D d * in our framework equals to the orthonomalized-2 estimator in Deák (2000), although this connection was not revealed in the original paper.

Simulation studies
For easy presentation, §4.1 presents algorithms for each estimator, and discusses the computational cost and penalized variance ratios. A simulation design is given in §4.2, and the associated numerical results are summarized in §4.3. Finally, §4.4 compares the performance of V arising from the proposed lattice L max d with various lattices for d = 16 and 24.

Algorithms and the computational cost
To fairly compare the efficiency of the proposed estimators, we report variance ratios and penalized variance ratios as measures in our simulation studies. The variance ratio is defined as the variance of a crude Monte Carlo estimator divided by that of an estimator of interest. An estimator with variance ratio larger than one is hence more efficient than the crude Monte Carlo estimator. Since the quantity of interest in this paper is the multivariate normal probability and the integrand consists of an indicator function, we restrict the computational cost to be the possible smallest number of independent random numbers required to generate one realization for the estimator. By incorporating the computational cost, we define the penalized variance of an estimator as the product of its variance and associated computational cost, and define the penalized variance ratio of an estimator as the penalized variance of the crude Monte Carlo estimator divided by that of the estimator of interest.
To estimate p A given in (1), we present explicit algorithms for each estimator described in this paper. Let M be the Monte Carlo sample size.
Because the crude spherical Monte Carlo estimator is indeed the same as the usual crude Monte Carlo estimator, a procedure for generatingp is given as follows.
Similarly, a procedure for generatingp AT in (6) is given as follows.

Setp
It is clear that the computational cost ofp andp AT are both equal to d.

For a given centrally symmetric subset
An algorithm for calculatingp V is implemented as follows.

Setp
Because a realization ofp V requires to generate one random orthogonal matrix and |V | independent radii, the computational cost ofp V is (d + 2)(d − 1)/2 + |V |. Recall in §2.2 that the computational cost to generate a d × d Note that the the cardinality of the proposed V is even since it is centrally symmetric. Let V + consists of all vectors in V whose first non-zero coordinate is positive, and V is decomposed as the disjoint union of V + and −V + . Without loss of generality, we assume that V + = {v 1 , · · · , v |V |/2 }. An algorithm for calculatingp V AT is outlined as follows.
2. Generate a radius, r
Because generating a realization ofp V AT requires a generation of one random orthogonal matrix and |V |/2 radii, the computational cost ofp V AT is (d + 2)(d − 1)/2 + |V |/2.

For calculatingp
. When explicit formulas for fÃ(u) exists (in terms of cdf of chi distributions),p V * can be calculated as follows.

Setp
Note that the calculation of fÃ(u) involves the cdf of a χ-distribution, and such additional computational effort can not be reflected in the definition of the penalized variance ratios. In our simulation studies in §4, we do not report the penalized variance ratio forp V * . Table 2 summarizes the computational cost for four estimators:p,p AT ,p V , andp V AT .

Simulation design
In this subsection, we would like to compare the efficiency of various estimators in calculating P {X ∈ A} for some set A and X ∼ N d (0, Σ).
Similar to the simulation settings in Vijverberg (1997), we consider three types of covariance models: the identity covariance matrix, the one-factor model, and the AR (1)  We consider three types of regions: ellipsoids, orthants and rectangles.
For each type of the region, we denote three sets as follows, • ellipsoid regions: Note that the interiors of R 2 , O 1 , O 3 , and E 2 are centrally antisymmetric.
It is clear that the property of centrally antisymmetric is preserved under changing coordinate via a linear transformation. Therefore when the region A is centrally antisymmetric,Ã = Γ −1 (A) remains centrally antisymmetric.
Also, we run the simulation for d from two to eight. Overall, there are seven dimensions, five estimators, nine covariance structures, and nine sets of interest, and a total of 2,835 combinations. The Monte Carlo sample size is 10, 000 for each case. Table 3 summarizes averaged (penalized) variance ratios among all types of the covariance structures as described in §4.2. Overall, the efficiency ofp V * is the highest, followed byp V andp V AT , andp AT . On average,p AT produces variance ratio about two for all dimensions, whereasp V ,p V AT ,p V * produce higher variance ratios for higher dimensions. In addition, the variance ratio ofp V * is dramatically larger in all cases. This numerical evidence suggests that if the innermost radial integral f A (u) can be calculated explicitly,p V * are preferred.

Numerical results
On the other hand, although calculatingp V * does not need generating the radius, its calculation requires the calculation of the cdf of the χ-distribution, which may be computationally demanding. Thus, we omit listing averaged penalized variance ratios forp V * . Again, the averaged penalized variance ratios forp AT are about two for all dimensions, but increase slightly forp V AT andp AT for higher dimensions. Although averaged variance ratios ofp V AT are slightly smaller than those ofp V in most cases, averaged penalized variance ratios ofp V AT are larger. In other words, these numerical results suggest that althoughp V AT do not yield higher variance ratio thanp V , accounting for the expense of drawing independent random variables,p V AT are preferred. Our simulation studies show thatp V * work particularly well for elliptical regions. In fact, when A is a circle centred at the original, f A (u) is about the same for all u, and thusp V * produces a dramatically high variance ratio. (When the covariance matrix is the identity matrix, f A (u) is identical for all u, and is a zero-variance estimator.) This exemplifies a feature that a spherical Monte Carlo estimator,p V * , is very efficient for certain ellipsoid regions (regardless of the selection of V ).
To investigate how the property of centrally antisymmetric affect the results, we average the (penalized) variance ratios among all centrally antisymmetric sets (R 2 , O 1 , O 3 , and E 2 ). See Table 4 for details. In this case, the variance ratios ofp V AT are quite close to those ofp V . This provides a numerical evidence that the property of centrally antisymmetric does matter when incorporating the idea of antithetic variates.
Note that due to an interesting feature in these centrally antisymmetric sets R 2 , O 1 , O 3 , and E 2 , the averaged (penalized) variance ratios in Table   4 forp V andp V AT are indeed the same. To explain why, let A currently be one of these four sets. Given a v ∈ V , if rv belongs in A, −r ′ v does not belong in A for all r ′ . As a result, a constituent of the spherical estimator and a constituent of the spherical estimatorp V , I A (r v1 , T v ) + I A (r v2 , −T v), both become zero or one. This leads to that both estimators have the same variance. Therefore, the inequality in Theorem 2.3.1 is indeed an equality for these four sets.
To provide a reasonable set showing thatp V AT enjoys a lower variance, con- Clearly, this set is centrally antisymmetric, and, in addition, it allows that when using antithetic variate, a constituent of estimatorp V AT , would take values of zero and one. On the contrary, when no antithetic variate is used, a constituent of the estimatorp V , I A (r v1 , T v )+I A (r v2 , −T v), would take values of zero, one, and two. Because these two estimators are unbiased, the   and the Leech lattice for d = 24. The explicit construction of the shortest vectors of these lattices can be also founded in Conway and Sloane (1999). It is clear that the cardinality of the proposed V L max d is dramatically larger than the other three counterparts. The variance reduction by using V max L d is substantial. As for the variance ratio penalized by the cardinality, the spherical Monte Carlo estimator employing the propsoed V L max d remains the most competitive one at a factor of about three to four.

Conclusion
Clearly, the proposed spherical methods exemplify the omnibus method for calculating multivariate normal probabilities that we seek. To explain the superior performance of the proposed integration methods, we look at the differences in integrating the radial and spherical parts. To reduce the variance for the spherical integral,, instead of sampling one point randomly, we randomly rotate a set of points on the sphere which is constructed from a sphere packing with the maximal kissing number. Moreover, we employ Although the proposed spherical methods give superior efficiency, this is only part of the story in the sense that we only study the case of multivariate normal probabilities. For more general cases such as the calculation of multivariate t distribution, Dirichlet distributions, elliptical copulas, and even the high-dimensional integrals with general function, we need to utilize the proposed spherical method with certain features of the distributions, and the best implementation in this paper can be adapted to various circumstances.
After some exploration, the choice of spherical integral rule, the interval size of the radial integral, and perhaps the radial integral rule itself can be adjusted. Of course, the choice of V depends on the development of kissing numbers in sphere packings, an alternative is to employ basis points chosen based on Case 1. Other possibility is to consider unequal weights for the basis points, which is related to weighted spherical t-designs. Using a general method such as sampling from a χ-distribution for the radial part allows for any reasonable sort of tail behaviour, while it remains competitive with other approaches when the tail of the distribution is not normal.
In conclusion, being able to evaluate multivariate normal probabilities of higher dimensions at low cost opens up new avenues of research on limited dependent variables in time series analysis, panel data, and spatial models.
This one-step estimator becomes more efficient when more adjacent periods are included in the construction of the step, at the computational cost of evaluating higher-order integral. In this paper, we assume that the underlying random variables are independent over time. A more challenging project is to model the time dependence used in, for example, financial industry implementation of the CreditMetrics model and the computation of Basel III incremental risk charge (IRC). Since our method highlights the two aforementioned key features of spherical and radial integrations in general settings, it is also expected that this method can be applied to option pricing, Greeks letters calculation, and correlated default probabilities evaluation among others.
Appendix A. Proof of Theorem 2.2.1 The proof of Theorem 2.2.1 depends on the calculation of variance in (10). To this end, we first consider the covariance of any two random variables appeared in (10). Let A and A ′ be two disjoint regions of S d−1 satisfying that Diam(A) and Diam(A ′ ) are both less than d min (V ). Because the covariances of g V A (T ) and g V A ′ (T ) is difficult to calculate, we would like to find an upper bound of the covariance instead. Since g V A (T ) and g V A ′ (T ) are simple step functions, we have Since all points v in V are predetermined and T v is simply a rotation on the sphere, therefore Recall that when v is fixed, T v is uniformly distributed on S d−1 which implies that P {T v ∈ A} = π(A) and P {T v ′ ∈ A ′ } = π(A ′ ) and E[g V A (T )g V A ′ (T )] ≤ 1 |V | 2 min(|V |π(A), |V |π(A ′ )) ≤ 1 2|V | (π(A) + π(A ′ )). Hence, Then, we have Appendix B. Proof of Theorem 2.3.1 To have the proof, we first note that If A is centrally antisymmetric, then for any χ(d)-distributed independent random variables r 1 , r 2 , and χ(d)-distributed random variable r, we have Cov (I A (r 1 , u), I A (r 2 , −u)) − Cov (I A (r, u), I A (r, u)) Followed by (7), we can rewrite (B.1) as v∈V + Cov(I A (r v , u), I A (r v , −u)) − Cov(I A (r v , u), I A (r −v , −u)) = |V | 2 (Cov (I A (r, u), I A (r, −u)) − Cov (I A (r 1 , u), I A (r 2 , u))) which is less than or equal to zero by (B.2).