Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
57,699 characters · 19 sections · 40 citation commands
Practical volume computation of structured convex bodies, and an application to modeling portfolio dependencies and financial crises
Modern finance has been pioneered by Markowitz who set a framework to study choice in portfolio allocation under uncertainty, see M52.\footnote{for which he was awarded the Nobel Prize in economics in 1990.} Within this framework, Markowitz characterized portfolios by their return and their risk which is defined as the variance of the portfolios' returns. And an investor would build a portfolio that will maximize its expected return for a chosen level of risk.\footnote{In the same way, by choosing a level of expected return, an investor can construct a portfolio which minimizes the risk.} It has since be common for asset managers to optimize their portfolio within this framework. And it has led a large part of the empirical finance research to focus on the so-called efficient frontier which is defined as the set of portfolios presenting the lowest risk for a given expected return. Figure (ref) (left panel) presents such an efficient frontier. The region on the left of the efficient frontier represent the portfolios domain.
Interestingly, despite the fact that this framework considers the whole set of portfolios, no attention has been given to the distribution of portfolios. Figure (ref) (middle panel) presents such distribution\footnote{10.000.000 portfolios have been sampled as presented later in Section (ref).}. When comparing the contour of the empirical portfolios distribution\footnote{Region over which at least 1 random portfolio lies.} and the portfolio domain bounded by the efficient frontier in Figure (ref) (right panel), we observe that the density of portfolios along the efficient frontier is dim and that most of the portfolios are located in a small region of the portfolios domain.
We also know from the financial literature that financial markets exhibit 3 types of behavior. In normal times, stocks are characterized by slightly positive returns and a moderate volatility, in up-market times (typically bubbles) by high returns and low volatility, and during financial crises by strongly negative returns and high volatility, see e.g. BGP12 for details. So, following Markowitz' framework, in normal and up-market times, the stocks and portfolios with the lowest volatility should present the lowest returns, whereas during crises those with the lowest volatility should present the highest returns. These features\footnote{also called “stylized facts" in the financial literature} motivate us to describe the time-varying dependency between portfolios' returns and volatility.
However this dependency is difficult to capture from the usual mean-variance representation, as in Figure (ref) (middle panel), so we will rely on the copula representation of the portfolios distribution. A copula is a bivariate probability distribution for which the marginal probability distribution of each variable is uniform. As we following Markowitz' framework, the variables considered are the portfolios' return and variance. Figure (ref) illustrates such a copula and shows a positive dependency between portfolios returns and variances. Each line and column sum to 1% of the portfolios.
The methods introduced here can be used to study other dependencies such as the momentum effect JT93 which is implied by the dependencies of asset returns with their past returns.
The dependencies mentioned here are important because
Interestingly, the copulas can be computed over a single period of time making the information available as early as the sample allows. The copula for the momentum dependency can be computed over very short periods (even intra-daily). The copula for the return/volatility dependency requires the estimation of the stock returns variance-covariance matrix which has to be estimated over a sufficiently large period of time to be reliable thus delaying the detection of crisis.\footnote{Methods exist to estimate the stock returns variance-covariance matrix over short periods, see e.g. the range-based estimation method BGP12. However they usually requires high-frequency data and are not widely used. These methods are beyond the scope of this paper.}
In the general case, the framework to describe the dependencies is as follows. First, as the set of portfolios, we consider the canonical $d$-dimensional simplex $\Delta^d\subset{\mathbb R}^{d+1}$ where each point represents a portfolio and $d+1$ is the number of assets. The vertices represent portfolios composed entirely of a single asset. The portfolio weights, i.e.\ fraction of investment to a specific asset, are non-negative and sum to 1. This is the most common investment set in practice today, as portfolio managers are typically forbidden from short-selling or leveraging. Second, considering some asset characteristic $ac$ quantified by $C\in{\mathbb R}^{d+1}$, we define a corresponding quantity $f_{ac}(\omega, C)$ for any portfolio $\omega \in \Delta^d$. For instance, considering the vector of asset returns $R\in{\mathbb R}^{d+1}$, $\omega$ has the return $f_{ret}(\omega, R) = R^T \omega$. Then, we define the cross-sectional score of a given portfolio $\omega^*$ as $$ \rho_{ac} = \frac{ \mbox{vol}(\Delta^*)}{\mbox{vol}(\Delta^d)} ,\; \mbox{ where } \Delta^*=\{ \omega \in \Delta^d \, :\, f(\omega, C) \leq f(\omega^*, C)\} , $$ which corresponds to the share of portfolios with a return lower or equal to $R^*=R^T \omega^*$. This score corresponds to the cumulative distribution function of $f_{ac}(\omega,C)$ where the portfolios are uniformly distributed over the simplex. In the following, we consider the cases where $f_{ac}$ is a linear combination or a quadratic form of $C$. Finally, the relationship between two asset characteristics $ac_1$ and $ac_2$ is presented in the form of a {\em copula} whose marginals are $\rho_{ac_1}$, $\rho_{ac_2}$. In our applications, the asset characteristics considered are the assets' returns and variances, and their values correspond to a linear combination of the returns and a quadratic form of the returns, respectively.
A copula is computed by slicing a simplex, i.e. the set of portfolios, along the asset characteristics. Thus, these questions are formulated in terms of convex bodies defined by intersecting simplices on one hand by a family of parallel hyperplanes and, on the other hand, by another family of parallel hyperplanes in the linear case or a family of concentric ellipsoids in the quadratic case. Furthermore, the latter case yields non-convex bodies between two ellipsoids.
The cross-sectional score of portfolio returns has been introduced in P05 and it is estimated by means of a quasi-Monte Carlo method. The applications have been limited in terms of dimensions: the 30 DAX components and the 24 MSCI Netherlands components in P05, the 35 components of the IBEX in PST04. This score has also been proposed in BCG11, for the set of long/short equally weighted zero-dollar portfolios and whose estimation relying on combinatorics and statistics is computationally limited to around 20 dimensions, and in BH11 where the focus was not on a precise score.
Given that volume computation of polytopes is $\#$P-hard for both V- and H-representations DyerFr88 and no poly-time algorithm can achieve better than exponential error Elekes86, the problem is not expected to admit of an efficient deterministic algorithm in general dimension. Developing algorithms for volume computation has received a lot of attention in the exact setting Bueler00. In the approximate setting, following the breakthrough polynomial-time algorithm by random walks DyerFrKa91, several algorithmic improvements ensued. The current best theoretical bounds are in Lee17b and for polytope sampling in Lee17a. Interestingly, only two pieces of software offer practical algorithms in high dimension: {\tt VolEsti}, a public-domain C++ implementation that scales to a few hundred dimensions VolEsti, based on the Hit-and-Run paradigm Lovasz99, and the Matlab implementation of CousinsV13, which treats hyperplanes as an ellipsoid, and seems competitive to {\tt VolEsti}\ in very high dimensions. Sampling from non-convex bodies appears in experimental works, with very few methods offering theoretical guarantees, e.g.\ in star shaped bodies Chandrasekaran09 or, more recently, in Abbasi17.
We design and implement the following different approaches for volume computation: Efficient sampling from the simplex and using rejection to approximate the target volume, which is fast but inaccurate for small volumes. Exact formulae of integrals of appropriate probability distribution functions, which are implemented for the case of a single hyperplane. Optimizing the use of Lawrence's sign decomposition method, since the polytopes at hand are shown to be simple with extra structure; a major issue here is numerical instability. Extending state-of-the-art random walks, based on the hit-and-run paradigm, to convex bodies defined as the intersection of linear halfspaces and ellipsoids. The latter is experimentally generalized to non-convex bodies defined by two ellipsoids with same quadratic form, and accurate approximations are obtained under certain mild conditions.
Our randomized algorithms for volume approximation extend {\tt VolEsti} , where the main problem to address is to compute the maximum inscribed ball of the convex body $P$ a.k.a.\ Chebychev ball. This reduces to a linear program when $P$ is a polytope. For a convex body defined by intersecting a polytope with $k$ balls, the question becomes a second-order cone program (SOCP) with $k$ cones. When interchanging input balls with ellipsoids, the SOCP yields a sufficiently good approximation of the Chebychev ball.
Our implementations are in C++, lie in the public domain ({\tt github}), are based on CGAL, rely on Eigen for linear algebra, on {\tt Boost} for random number generators, and experiment with two SOCP solvers for initializing random walks. Our software tools are general and of independent interest. They are applied to allow us to extend the computation of a portfolio score to up to 100 dimensions, thus doubling the size of assets studied in financial research. We thus provide a new description of asset characteristics dependencies. Our methods allow us to propose and to effectively compute a new indicator of financial crises, which is shown to correctly identify all past crises with which we experimented. More importantly, it allows us to establish that periods of momentum nearly never overlap with the crisis events, which is a new result in finance.
The rest of the paper is organized as follows. In Section (ref) we describe the convex bodies that intersect the canonical simplex and arise from the financial modeling we work on. Moreover we overview methods for representing and uniform sampling from the simplex and give theoretical guarantees for the sampling - rejection method accuracy. Section (ref) considers volumes defined as the intersection of a simplex and one hyperplane or more hyperplanes, the latter being organized in at most two families of parallel hyperplanes. Section (ref) studies convex and non-convex bodies defined as the intersection of a simplex and an ellipsoid, for which random walk methods are developed. The implementations are discussed in Section (ref), along with experiments that show the validity of our approach in answering open questions in finance. We conclude with current work and open questions. Figures and tables that do not fit are given in the Appendix.
We analyze real data consisting of regular interval (e.g.\ daily) returns of stocks such as the constituents of the Dow Jones Stoxx 600 Europe\texttrademark (DJ600). These are points in real space of dimension $d=600$, respectively: $r_i=(r_{i,1},\dots,r_{i,d}) \in \mathbb{R}^{d}$, $i\ge 1$.
We apply the methodology to a subset of assets drawn from the DJ 600 constituents\footnote{The data used is from Bloomberg\texttrademark.. It is daily and ranges from 01/01/1990 to 31/11/2017.}. Since not all stocks are tracked for the full period of time, we select the 100 assets with the longest history in the index\footnote{This implies a survivor bias, but we use it to assess the effectiveness of the methodology. One would wish to keep 600 constituents, replacing the exiting stocks with the entering ones along the sample.}, and juxtapose:
In financial applications, one considers compound returns over periods of $k$ observations, where typically $k=20$ or $k=60$; the latter corresponds to roughly 3 months when observations are daily. Compound returns are obtained using $k$ observations starting at the $i$-th one where the $j$-th coordinate corresponds to asset $j$ and the component j of the new vector equals: $$ (1+r_{i,j})(1+r_{i+1,j})\cdots (1+r_{i+k-1,j}) -1, \quad j=1,\dots, d. $$ This defines the normal vector to a family of parallel hyperplanes, whose equations are fully defined by selecting appropriate constants. The second family of parallel hyperplanes is defined similarly by using an adjacent period of $k$ observations.
The covariance matrix of the stock returns is computed using the shrinkage estimator of LW04,\footnote{Matlab code on \url{http://www.econ.uzh.ch/en/people/faculty/wolf/publications.html}.} as it provides a robust estimate even when the sample size is short with respect to the number of assets. A covariance matrix $C$ defines a family of ellipsoids centered at the origin $0\in\mathbb{R}^{d}$ whose equations $x^TCx=c$ are fully specified by selecting appropriate constants $c$.
To compute the copulas, we determine constants defining hyperplanes and ellispoids so that the volume between two consecutive such objects is 1% of the simplex volume. The former are determined by bisection using the Varsi's exact formula. For ellipsoids $E(x)=c_i$, we look for the $c_i$'s by sampling the simplex, then evaluating $E(x)$ at each point. The values are sorted and the $c_i$ selected so as to define intervals containing 1% of the values. Two consecutive ellipsoids intersecting the simplex and the family of parallel hyperplanes define a non-convex body for which we practically extend {\tt VolEsti} \ algorithm.
The volume between two consecutive hyperplanes and two consecutive ellipsoids defines the density of portfolios whose returns and volatilities lie between the specified constants. We thus get a copula representing the distribution of the portfolios with respect to the portfolios returns and volatilities. Fig. (ref) illustrates such copulae, and shows the different relationship between returns and volatility in good (left, dot-com bubble) and bad (right, bubble burst) times.\footnote{We consider 100 components of DJ 600 with longest history, over 60 days ending at the given date.}
The main problem is to compute all the volumes that arise from the intersection of the two families with the unit simplex. We totally have to handle three types of full dimensional bodies and thus we develop or use existing methods for three different problems. The first is to compute the volume of the polytope defined by the intersection of the unit simplex with four hyperplanes which are pairwise parallel. The second arises when an ellispoid intersects with the unit simplex and a family of parallel hyperplanes. The third is to compute the volume of a non-convex body defined by the intersection of two concentric ellipsoids with a simplex and a family of parallel hyperplanes.
We develop and use four methods in total. The first (M1) is an exact formula for the volume defined by the intersection of simplex with a hyperplane. The second (M2 or s/r) is to sample the unit simplex and approximate all the volumes directly. The third method (M3) is the optimized Lawrence formula for simple polytopes and is used for the first problem. The fourth method (M4) is the generalization of the {\tt VolEsti} \ algorithm to non-linear and non-convex bodies.
This subsection sets the notation, surveys methods for uniform sampling of the simplex, and discusses their efficient implementation.
The $d$-dimensional simplex $\Delta^{d}\subset{\mathbb R}^{d+1}$ may be represented by barycentric coordinates $\lambda=(\lambda_0,\dots,\lambda_d)$ s.t.\ $\sum_{i=0}^{d} \lambda_i=1, \, \lambda_i \ge 0$. The points are $\sum_{i=0}^{d} \lambda_i v_i$, where $v_0, \dots, v_d \in \mathbb{R}^{d}$ are affinely independent. It is convenient to use a full-dimensional simplex, by switching to Cartesian coordinates $x =(x_1, \dots, x_{d})$: $$ m_{bc}: \mathbb{R}^{d+1} \mapsto \mathbb{R}^{d} : \lambda \rightarrow x = M (\lambda_1,\dots,\lambda_n)^T + v_0,\; \mbox{ where } M = [ v_1 -v_0\, \cdots \ v_{d} -v_0 ] , $$ is a $d\times d$ invertible matrix. The inverse transform is:
where $0_{d}, 1_{d}$ are $d$-dimensional column vectors of 1's and 0's, respectively, and $I_{d}$ is the $d$-dimensional identity matrix.
A number of algorithms exist for sampling, where some have been rediscovered, while others contain errors; see the survey SmTr. Let us start with a unit simplex in Cartesian coordinates. A $O(d \log d)$ algorithm is the following David81,D86,RbKr07: Generate $d$ distinct integers uniformly in $\{1,\dots ,K-1\}$, where $K$ is the largest representable integer. Sort them as follows: $x_0=0<x_1<\cdots <x_{d+1}=K$. Now ${(x_i-x_{i-1})}/K$, $i=1,\dots ,d$, defines a uniform point. Assuming we possess a perfect hash-function, the choice of distinct integers takes $O(d)$. For $d>60$ we implement a variant of Bloom filter to guarantee distinctness.
A linear-time algorithm is given in RbMel98, which is is generally the algorithm of choice, although it is slower for $d<80$: \if 0 Generate $d+1$ independent unit-exponential random variables $y_i$ and normalize by their sum $s$: $y_i/s,\, i=1,\dots,d$ is a uniform point. For this, we uniformly sample real $x_i\in (0,1)$, and set $y_i=- \log x_i$. In practice, this algorithm is slower in low dimensions. \fi
To sample an arbitrary simplex, we can map sampled points from the unit simplex by transformation ((ref)), which preserves uniformity. Due to applying the transformation, the complexity is $O(d^2)$ to generate a uniform point. The same complexity, though slower in practice, is achieved in grimme.
Sampling could be used in order to approximate all the volumes that arising when two families of parallel hyperplanes intersect with a simplex. One can sample the simplex and count the percentage of points in the region of interest. The complexity is $O(kd)$ to generate $k$ points. In the case of a family of $\ell$ parallel hyperplanes, all sample points are evaluated at the hyperplane linear polynomials in time $O(kd)$. Given the $\ell$ constant terms characterizing the hyperplanes, for each point we perform a binary search so as to decide in which layer it lies. Hence the total complexity is $O(k\log \ell)$, which is dominated since $\ell\le 100$ typically. Given a family of $\ell$ ellipsoids with same quadratic form intersecting a simplex, the method requires $O(kd^2)$ to evaluate all sample points and $O(k\log \ell)$ to assign them to layers.
In this subsection we obtain how we could guarantee a bounded error for the sampling - rejection method with high probability. Let $B$ be a convex or non convex full dimensional body in dimension $d$ and let $S$ be an enclosing simplex such that $B\subseteq S$ and let $p=\dfrac{Vol(B)}{Vol(S)}$. Then if we uniformly sample a point from $S$ it lies in $B$ with probability $p$. So if we sample $N$ points from $S$ the random variable $X$ which gives the number, $k$, of points that lie in $B$ follows the binomial distribution. So, $P(X=k)={{N}\choose{k}} p^k(1-p)^{N-k}$. Moreover if we sample $N$ points and reject,
is the probability that the sampling - rejection method error is at most $e$. From Poison Limit Theorem we know that if $\lim_{N\rightarrow \infty}Np=\lambda$ is a constant independent of $N$ then for any fixed $k$, $lim_{N\rightarrow\infty}P(X=k)=e^{-\lambda}\dfrac{\lambda^k}{k!}$. If we set $N=m_1\cdot 10^x$, where $x=m_2+\lceil -\log_{10}p\rceil$, we notice that $Np\approx m_1\cdot 10^{m_2}$. So we can use Poisson limit theorem to approximate the random variable $X$ as we usually set large enough values for $m_1$ and $m_2$. Then from Poison cumulative distribution function we can approximate probability (ref).
In Table (ref) we give for several values of errors the $m1,m2$ in order to get probabilities higher than 0.94. Notice that we use the order of $p$ at the exponent. In practice we could estimate $p$ while we are sampling and obtain the right order of $p$ with very high probability.
This section considers computing the volume of the intersection of a simplex and one or more linear halfspaces. The most general case is to be given two families of parallel hyperplanes and consider all created polytopes. We assume that the simplex is given in V-representation, i.e.\ as a set of vertices, and the hyperplanes by their equations.
We can always transform the simplex to be a unit full-dimensional simplex with the origin as one vertex by the transformation of Sect. (ref). The same transform applies to the hyperplanes, and volume ratios as preserved. \if 0 by shifting them using simplex vertex $V_0$ and use the linear transformation that is defined by matrix $A=[V_1-V_0|\dots |V_{d+1}-V_0]$. If $H_1: c_1x_1\dots c_dx_d<z$ is a hyperplane intersecting a simplex then $H'_1:c_1(x_1+V_{01})+\dots +c_d(x_d+V_{0d})<z\Leftrightarrow c^Tx<z'$ is the shifted one. Then $H:c^TAx<z'$ is the final hyperplane. So in practice is much faster to transform hyperplanes and sample from the unit simplex than sample to the arbitrary simplex. \fi
Surprisingly, there exist an exact, iterative formula (M1) for the volume defined by intersecting a simplex with a hyperplane. A geometric proof is given in VARSI, by subdividing the polytope into pyramids and, recursively, to simplices. We implement a somewhat simpler formula Ali73, which also requires $O(d^2)$ operations. Let $H = \{(x_1,\dots ,x_d) | \sum_{i=1}^{d} a_ix_i\leq z\}$ be the linear halfspace.
If $\Delta^d\subset{\mathbb R}^d$ is the unit simplex then, for $h =J$, $A_K= {\mbox{vol}(\Delta^d\cap H)} / {\mbox{vol}(\Delta^d)}$.
Recall, from Section (ref), that sampling uniformly over the simplex can be obtained by drawing exponential random variables. Thus, an alternative formula follows from computing the cumulative distribution of a linear combination of exponential random variables. In Mathai, they propose an exact method to compute the distribution $f$ of such linear combination. It consists in representing $f$ as its moment generating function, analogous to a Laplace transform, simplifying it with a generalized partial-fraction technique of integration, before inverting its terms. However, in double precision, the method showed numerical discrepancies above 20 dimensions and was thus abandoned. However, it has the advantage of being generalizable to nonlinear combinations (see Sect. (ref)).
This section considers simple polytopes defined by a constant number of families of parallel hyperplanes; in our application there are two such families. The defined polytopes are simple, i.e., all vertices are defined at the intersection of $d$ hyperplanes, assuming that no hyperplane contains any of the simplex vertices and, moreover, two hyperplanes do not intersect on a simplex edge at the same point.
For a simple polytope $P$, the decomposition by Lawrence Lawrence picks $c \in \mathbb{R}^d$, $q\in\mathbb{R}$ such that $c^T x + q$ is not constant along any edge, i.e.\ $c,-c$ do not lie on the normal fan of any edge. For each vertex $v$, let $A(v)$ be the $d \times d$ matrix whose columns correspond to the equations of hyperplanes through $v$. Then $A(v)$ is invertible and vector $\gamma(v)$ such that $A(v) \gamma(v) = c$ is well defined up to a permutation. The assumption on $c$ assures no entry vanishes, then $$ \mbox{vol}(P) = \frac1{d!} \, \sum_{v} \frac{ (c^T v + q )^ d }{ |\det A(v) | \; \prod_{i=1}^d \gamma(v)_i }. $$ The computational complexity is $O (d^3 n)$, where $n$ is the number of vertices. We set $q=0$ for simplicity in the implementation. An issue is to choose $c$ so as to avoid that $c^T x + q$ be nearly constant on some edge, because this would result in very small entries in the denominator and numerical issues. A theoretical choice is given in Lawrence, but its practical importance is very small. The main drawback of Lawrence's decomposition remains numerical instability when executed with floating point numbers, and high bit complexity, when executed over rational arithmetic. The latter is indispensible for $d>30$ in our applications, because then numerical results become very unstable.
To compute the volume defined by the intersection of a simplex and two arbitrary hyperplanes, we exploit the fact that the simplex is unit in order to compute more effectively the determinants and the solutions of the linear system. The hardest case is when vertex $v$ is defined by the two arbitrary hyperplanes $H_a, H_b$, the supporting hyperplane $H_{0}:\sum_{i=1}^{d}x_i=1$, and $d-3$ hyperplanes of the form $H_i:x_i=0$. Then, up to row permutations,
where the $i_j,\, i=a,b$ are the coefficients of the equation of $H_i$ up to permutation. Then we solve the lowest right $3\times 3$ linear system in $O(1)$ and then the computation of each remaining unknown $\gamma(v)_i$, $i=1,\dots ,d-3$ requires $O(1)$ operations for a total of $O(d).$ The corresponding determinant is computed in $O(1)$.
Lawrence's formula requires both H- and V-representation. In our setting, the H-representation is known, but the previous lemma allows us to obtain vertices as well.
The entire discussion extends to polytopes defined by two families of parallel hyperplanes. The matrices $A(v)$ remain of the same form because each vertex is incident to at most one hyperplane from each family.
This section considers more general convex bodies, defined as a finite, bounded intersection of linear and nonlinear halfspaces. For this, we extend the polynomial-time approximation algorithm in {\tt VolEsti} VolEsti so as to handle nonlinear constraints. Our primary motivation here is computing the volume of the intersection of a simplex with an ellipsoid in general dimension.
The method in VolEsti follows the Hit-and-Run algorithm in Lovasz99, and is based on an approximation algorithm in $O^*(d^5)$. It scales in a few hundred dimensions by integrating certain algorithmic improvements to the original method. We have to generalize the method because the input is not a polytope but a general convex body, while {\tt VolEsti} \ works for $d$-polytopes. It suffices to solve two subproblems: Compute the maximum inscribed ball of the convex body a.k.a.\ Chebychev ball, and compute the intersection points of a line that crosses the interior of the convex body $P$ with the boundary of $P$.
The first problem is treated in the next subsection. For the second one, when the body is the intersection of linear and quadratic halfspaces, it suffices to solve systems of linear or quadratic equations. In our case where $P$ has few input hyperplanes we can optimize that procedure by transforming a base of our polytope to an orthonormal base thus obtaining very simple linear systems. One heuristic is to first compute the intersection of the line with all hyperplanes and test whether the intersection points lie inside the ellipsoid so as to avoid intersecting the line with the ellipsoid. Formally, every ray $\ell$ in Coordinate Direction Hit-and-Run is of the form ${p}+\lambda{e}_k$ and parallel to $d-1$ simplex facets. The roots of $\lambda^2+2\lambda p_k+|p|^2-R^2$ define the intersection of a sphere with radius $R$, centered at the origin, and a coordinate direction ray $\ell$. If $C$ is the matrix of an ellipsoid centered at the origin its intersections with $\ell$ are roots of:
Computing the roots, and keeping the largest negative and smallest positive $\lambda$ is quite fast.
In our application, there are non-convex bodies defined by the intersection of two parallel hyperplanes and two concentric ellipsoids. We thus modify {\tt VolEsti}\ in order to compute the non convex volume. We make two major changes. First, in ray shooting, we have to check whether one quadratic equation has only complex solutions, which implies the ray does not intersect the ellipsoid. For $\lambda$, we take the largest negative and the smallest positive root in every step as well. Second, for the initial interior point, we sample from the unit simplex and when we find a point inside the intersection we stop and use it for initialization. We define an inscribed ball with this center and radius equal to some small $\epsilon>0$. We stop the algorithm when we find the first inscribed ball as described in the next subsection. So we can set $\epsilon$ sufficiently small so it always defines an inscribed ball in practice, but the enclosing ball is enough to run the algorithm and do not stop until we find an inscribed ball.
The method works fine for $d<35$ using the same walk length and number of points as for the convex case, and has time complexity and accuracy competitive to running {\tt VolEsti} \ on the convex set defined by one ellipsoid. For $d> 35$, the method fails to approximate volume for most of the cases. This should be due to inaccurate rounding bodies and the inscribed ball we define. Given these first positive results, various improvements are planned.
This section offers methods for computing a ball inside the given convex region. Ideally, this is the largest inscribed ball, aka Chebychev ball, but a smaller ball may suffice. Computing the Chebychev ball reduces to a linear program when $P$ is a polytope (p. 148 in Boyd_opt). For general convex regions, more general methods are proposed.
At the very least, one point must be obtained inside the convex region. When we do not have the Chebychev ball, an issue is that concentric balls with largest radii will again be entirely contained in the convex region, thus wasting time in the computation. In practice we use the one interior point as center of an enclosing ball, then reduce the radius until the first inscribed ball. To decide whether a given ball is inscribed, with high probability, we check whether all boundary points in Hit-and-Run belong to the sphere instead of any other constraint.
We start with some simple approaches. Let us consider the case of intersecting a simplex with an ellipsoid. If there are $z_1$ simplex vertices inside the ellipsoid and $z_2$ outside, then we have $(z_2+1)z_1$ vertices on the boundary of the convex intersection. Since $z_1+z_2=d+1$, then $(z_2+1)z_1\geq d+1$ and a new inscribed simplex is defined. In this case we take its largest inscribed ball and start hit and run. More generally, we sample from the unit simplex until we have $d+1$ points inside our section and then take the largest inscribed ball of this new simplex that is defined by the $d+1$ points. Another approach is to consider the transformation mapping the ellipsoid to a sphere and apply it both to the simplex and to the ellipsoid. We compute the distance from the sphere's center to the new simplex and compare it with the sphere's radius.
For a convex body that comes from intersecting a polytope with $k$ balls the problem becomes a Second-Order Cone Program (SOCP) with $k$ cones. However in our case we need to consider input ellipsoids. Assume that we transformed the ellipsoid to a ball $B'=\{x'_c+u' \, :\, \|u'\|\leq r'\}$, and applied the same transformation to the simplex to have $a_i x\leq b_i$ for $i\in [d+1], a_i\in \mathbb{R}^d, b_i\in\mathbb{R}$. The following SOCP computes the maximum ball $B=\{x_c+u\ :\ \|u\|\leq r\}$ in the intersection of the simplex and $B'$:
There are several ways to solve SOCP's such as to reformulate it to as a semidefinite program or perform a quadratic program relaxation. Moreover, since in our case we only have a single cone we could utilize special methods as in ErIy. However, for our case it suffices to use the generic SOCP solver from DoCh as it is very efficient; for a random simplex and ball, it takes $0.06$ sec in $d=100$ and $<20$ sec in $d=1000$, on Matlab using {\tt ecos} and {\tt yalmip} packages.
It is possible to apply the inverse transformation and get an inscribed ellipsoid, which is not necessarily largest possible. However we can use the maximum inscribed ball in that ellipsoid as an approximation of the Chebychev ball, by taking the center of that ellipsoid and the minimum eigenvalue of its matrix as the radius.
In our financial application, portfolios are points in the unit $d$-dimensional simplex $\Delta^{d}\subset {\mathbb R}^{d+1}$ defined as the convex hull of $v_0,\dots,v_d\in{\mathbb R}^d$, where $v_i$ lies on the $i$-th axis. The simplex lies in hyperplane $\sum_{i=0}^d \lambda_i=1$. To model levels of volatility, a family of full-dimensional ellipsoids in ${\mathbb R}^{d+1}$, centered at the origin, is defined by the covariance matrix $C$ of asset returns. We wish to compute the volume of intersections of this family with the simplex and, moreover, with a family of hyperplanes on the simplex. Rejection sampling would work in this context, however methods employing random walks require a full-dimensional convex body. Given a full $(d+1)$-dimensional ellipsoid $G:\lambda^TC\lambda-c\leq 0$ centered at the origin, where $C\in{\mathbb R}^{(d+1) \times (d+1)}$ is symmetric positive-definite, we compute the equation of the ellipsoid defined $G\cap \Delta^d\subset {\mathbb R}^d$, by imposing the constraint $\sum_{i=0}^d \lambda_i=1$ by transform $m_{cb}$ in expression ((ref)), thus obtaining: $$ (x-v_0)^T \left( M^{-T} \, [-1\, I_d] \, C \left[
\right] M^{-1} \right) (x-v_0) + A(x-v_0) = c', $$ where the expression in parenthesis is the matrix defining the new $d$-dimensional ellipsoid in Cartesian coordinates, and $A\in{\mathbb R}^{d \times d}, c'\in{\mathbb R}$ are obtained by direct calculation. Similarly the simplex maps to Cartesian coordinates.
\if 0 We shift the simplex and the ellipsoid by the vertex $V_d$ and define the matrix $M=[V_1-V_d|\cdots |V_{d-1}-V_d]\in{\mathbb R}^{d\times d-1}$. If ${x}$ is a Cartesian point in the full dimensional simplex in ${\mathbb R}^{d-1}$ then the linear transformation $M^{-1}x$ maps it to the shifted facet simplex.
Let $G_1:y^TCy\leq c$ be a full $d$-dimensional ellipsoid centered at the origin, where $C\in{\mathbb R}^{d\times d}$ is symmetric positive-definite. Hence, there is a matrix $A\in{\mathbb R}^{d\times d}$, such that $C=(AA^T)^{-1}$ and maps the set $\{y\in{\mathbb R}^d:y^TCy\leq c\}$ to the shifted $G_1$ using $R(y)=Ay-V_d$. So all $x^TM^TMx\leq c,\, x\in{\mathbb R}^{d-1}$ are maped by $R$ to the shifted $G_1$. But the last is an ellipsoid centered at the origin. So if we apply $R^{-1}$, using pseudoinverse $(AM)^{+}$, to the shifted vertices of the facet simplex we get an arbitrary simplex in ${\mathbb R}^{d-1}$ which intersects with the ellipdoid $\{G_2:x^TM^TMx\leq c\}$. Then we can transform both and get a simplex and a sphere in order to solve the conical problem and get an inscribed ball. \fi
Our implementations are in C++, lie in the public domain\footnote{https://github.com/TolisChal/volume_approximation}, and are using CGAL and Eigen. All experiments of the paper have been performed on a personal computer with Intel Pentium G4400 3.30GHz CPU and 16GB RAM. Times are averaged over 100 runs. Some resulting tables and figures are given in the Appendix.
We test the following convex bodies: a $d$-simplex intersected with: (1) two arbitrary halfspaces, (2) two parallel halfspaces, (3) an ellipsoid, (4) two parallel halfspaces and two cocentric ellipsoids (non convex body).
In general, M1 is preferred when available. Method M2 is the fastest and scales easily to 100 dimensions, so it is expected to be useful for larger dimensions. However, for small volumes its accuracy degrades; sampling more points makes it slower than M4. The latter is thus the method of choice for volumes $<1\%$ of the simplex volume, but it is not clear whether it would be fast beyond $d=100$. Method M3 is useful, even for small volumes, but it cannot scale to $d=100$ due to numerical instability; if we opt for exact computing, it becomes too slow.
The formula M1 is used in all Tables where exact computation is needed between two parallel hyperplanes intersecting the simplex.
Table (ref) considers the intersection of an arbitrary simplex with two hyperplanes. The vertices of each simplex are randomly chosen uniformly from the surface of a ball with radius 100, using {\tt CGAL} random point generator. All hyperplanes' coefficients are randomly chosen in $[-10,10]$ with {\tt Boost} (mt19937) random generator. For {\tt VolEsti}\ we do not use the rounding option for the input polytope. This means that skinny polytopes have low accuracy since the random walk mixes slow, cf.\ row 10 of Table (ref). On the other hand, M2 is not affected by polytope shape. Up to $d= 30$ and for large volume ratio, namely $>1\%$, M2 yields very accurate and fastest results. The last two experiments show that {\tt VolEsti}\ achieves the most accurate approximation when the ration of accepted sampled points is small.
In Table (ref) we use same runtime for M2 and M4 (analogous numbers of sampled points) and compare their accuracy. We perform two experiments per dimension. For the first, for each dimension we compute the volume between two parallel hyperplanes which is $1\%$ of the simplex volume. For exact volume computation we use (M1). For the second experiment, for each dimension we compute volumes defined by the intersection of 4 hyperplanes which are pairwise parallel with the simplex, which is close to $0.01\%$ of the simplex volume. For exact computation we used {\tt vinci} default method, {\tt rlass}. M2 is fast but inaccurate for small volumes; M4 is most accurate but should not scale beyond $d=100$.
In Table (ref) we have an arbitrary simplex and two arbitrary hyperplanes that intersect with it. We compare our Lawrence implementation in Sect. (ref), using floating-point and rational computation, with {\tt rlass} and M2. We have two parallel hyperplanes intersect the unit simplex. {\tt vinci} fails to compute the volume for $d>31$. Our exact computation works even in $d=100$ but becomes very slow.
Table (ref) compares M2 (s/r) with two variants of M4 for ellipsoid intersection. The only difference for the latter is the way we construct an inscribed ball: In s/V we implement random sampling until $d+1$ points are found, and in o/V we use SOCP. We see M2 is significantly faster than either variant of M4. All methods yield similar output values.
Table (ref) compares s/r with Hit-and-Run for non-convex bodies, as in Sect. (ref). Very small values of Volume means the method failed to approximate the volume. \if 0
We analyze real data consisting of regular interval (e.g.\ daily) returns corresponding to indices such as the Dow Jones Europe Stoxx 600\texttrademark (DJ600). These are points in real space of dimension $d=600$, respectively: $r_i=(r_{i,1},\dots,r_{i,d}) \in \mathbb{R}^{d}$, $i\ge 1$.
We apply the methodology to assets drawn among the DJ 600 constituents\footnote{The data used is daily and ranges from 01/01/1990 to 31/11/2017. It is from Bloomberg\texttrademark.}. Since not all assets are tracked for the full period of time, we select the 100 assets with the longest history in the index\footnote{This implies a survivor bias, but we use it to assess the effectiveness of the methodology. One would wish to keep 600 constituents, by replacing exiting by incoming stocks along the sample.}, and juxtapose:
In financial applications, one considers compound returns over periods of $k$ observations, where typically $k=30$ or $k=60$; the latter corresponds to roughly 3 months when observations are daily. Compound returns are obtained using $k$ observations starting at the $i$-th one where the $j$-th coordinate corresponds to asset $j$ and the component j of the new vector equals: $$ (1+r_{i,j})(1+r_{i+1,j})\cdots (1+r_{i+k-1,j}) -1, \quad j=1,\dots, d. $$ This defines the normal vector to a family of parallel hyperplanes, whose equations are fully defined by selecting appropriate constants. The second family of parallel hyperplanes is defined similarly by using an adjacent period of $k$ observations.
Volatility is modeled using the shrinkage estimator of the covariance matrix\footnote{Matlab code on \url{http://www.econ.uzh.ch/en/people/faculty/wolf/publications.html}.} LW04, as it provides a robust estimate even when the sample size is short with respect to the number of assets. A covariance matrix $C$ defines a family of ellipsoids centered at the origin $0\in\mathbb{R}^{d}$ whose equations $x^TCx=c$ are fully specified by selecting appropriate constants $c$.
We determine constants defining hyperplanes and ellispoids so that the volume between two consecutive such objects is 1% of the simplex volume. The former are determined by bisection using the formula (M1). For ellipses $E(x)=c_i$, we look for the $c_i$'s by sampling the simplex, then evaluating $E(x)$ at each point. The values are sorted and the $c_i$ selected so as to define intervals containing 1% of the values.
The volume between two consecutive hyperplanes and two consecutive ellipsoids defines the density of portfolios whose returns and volatilities lie between the specified constants. We thus get a copula representing the distribution of the portfolios with respect to the portfolios returns and volatilities. Fig. (ref) illustrates such copulae, and shows the different relationship between returns and volatility in good (left, dot-com bubbled) and bad (right, bubble burst) times.\footnote{We consider 100 components of DJ 600 with longest history, over 60 days ending at the given date.} \fi
When we work with real data in order to build the indicator, we wish to compare the densities of portfolios along the two diagonals. In normal and up-market times, the portfolios with the lowest volatility present the lowest returns and the mass of portfolios should be on the up-diagonal. During crisis the portfolios with the lowest volatility present the highest returns and the mass of portfolios should be on the down-diagonal, see Fig. (ref) as illustration. Thus, setting up- and down-diagonal bands, we define the indicator as the ratio of the down-diagonal band over the up-diagonal band, discarding the intersection of the two. The construction of the indicator is illustrated in Fig. (ref) where the indicator is the ratio of the mass of portfolios in the blue area over the mass of portfolios in the red one.
In the following, the indicator is computed using copulae estimated using the sampling method, drawing 500000 points. Computing the indicator over a rolling window of $k=60$ days and with a band of $\pm 10\%$ with respect to the diagonal, we report in Table (ref) all the periods over which the indicator is greater than 1 for more than 60 days. The periods should be more than 60 days to avoid the detection of isolated events whose persistence is only due to the auto-correlation implied by the rolling window. All these periods offer warnings, but only the longest ones correspond to crises.
We compare these results with the database for financial crises in European countries proposed in ESRB17. The first crisis (from May 1990 to Dec. 1990) corresponds to the early 90's recession, the second one (from May 2000 to May 2001) to the dot-com bubble burst, the third one (from Oct. 2001 to Apr. 2002) to the stock market downturn of 2002, the fourth one (from Nov. 2005 to Apr. 2006) is not listed and it is either a false signal or it might be due to a bias in the companies selected in the sample, and the fifth one (from Dec. 2007 to Aug. 2008) to the sub-prime crisis.
Regarding the momentum effect, i.e.\ the effect of the compound returns of the last 60 days on the following 60-day compound returns, we report the indicator in Fig (ref). We observe that there were only 10 events of lasting momentum effect, mostly around the 1998-2004 period. We remark that they nearly never overlap with the crisis events, with the exception of the end of 2011. To the authors' knowledge, this result is new in finance.
Since runtimes are very reasonable, we plan to extend our study to larger subsets of assets of DJ 600 and eventually the whole index in $d=600$. Another extension is to consider polytopes defined by intersections of both families of parallel hyperplanes and the family of ellipsoids, thus creating 3-D diagrams of dependencies, which have never been studied in finance: one difficulty is to model the outcome since visualization becomes intricate.
Random sampling follows a Monte Carlo (MC) approach by relying on C/C++ functions such as {\tt random} which implement pseudorandom generators. We are experimenting with quasi-MC generators which require fewer points to simulate the uniform distribution. Our preliminary experiments indicate that this may yield a speedup of about 2. An obvious enhancement is to parallelize our algorithms, which seems straightforward. Then results can be obtained for larger classes of assets such as the entire DJ 600.
One challenge is to extend the volume formula to the intersection with a ellipsoid. In MP13, they propose a method to approximate the distribution $f$ of quadratic forms in gamma random variables which is a similar problem to that in Mathai (see Sect. (ref)). It consists in fitting $f$ with a generalized gamma distribution by matching its first 3 moments with those of $f$ and to adjust the distribution with a polynomial in order to fit the higher moments. To get an approximation with a polynomial of degree $k$, the method requires the first $2k$ moments. In the case of a quadratic form in $d$ random variables, the moment of order $m$ is obtained by a sum over all the partitions of $m$ into $d^2$ terms. The number of partitions makes the computation of moments challenging even for $d \geq 5$.