QMF Sydney 2006:
Generating Quantile Functions 
using Computer Algebra:
Distributional Alchemy & 
other tools for 
Copula & Monte Carlo Simulations

Introduction

Given a distribution function (Subscript[F, X])(x), a simple means of simulation is to set

X = (Subscript[F, X]^(-1))( U)

where U is a sample from the uniform distribution on [0, 1]. Making U is compsci probem. The inverse CDF: (Subscript[F, X]^(-1))= (Subscript[Q, (Subscript[F, X])]) is the Quantile Function associated with the distribution. We do not have to do this, witness the use of Box-Muller, Polar-Marsaglia methods for the Normal case, and its extension to Student by Bailey. But it is very useful if we can, particularly if we are working with algorithms based on hypercube-filling quasi-Monte-Carlo methods, or in particular copula methods.

Quantile Functions and Copulas

For the case of  copula-based simulation, we make a sample ((Subscript[X, 1]),  (Subscript[X, 2])) from a given bivariate (in general multivariate) distribution, then form a sample from the associated copula:

{(Subscript[U, 1]), 
        (Subscript[U, 2])} = {(Subscript[F, (Subscript[X, 1])])((Subscript[X, 1])), (Subscript[F, (Subscript[X, 2])])(
          (Subscript[X, 2]))}

Then to get samples with marginals with CDFs (Subscript[G, i]):

{(Subscript[Y, 1]), 
        (Subscript[Y, 2])} = {(Subscript[G, 1]^(-1))[(Subscript[U, 1])], 
          (Subscript[G, 2]^(-1))[(Subscript[U, 2])]}  =  {(Subscript[Q, (Subscript[G, 1])])[(Subscript[U, 1])], 
          (Subscript[Q, (Subscript[G, 2])])[(Subscript[U, 2])]}

It is always helpful to know quantile functions. But note also we do not need the (Subscript[U, i]), at least when we are dealing with samples from a real distribution, for it would suffice to understand the composite mappings in the following:

{(Subscript[Y, 1]), 
        (Subscript[Y, 2])}={(Subscript[G, 1]^(-1))[(Subscript[F, (Subscript[X, 1])])((Subscript[X, 1]))], 
        (Subscript[G, 2]^(-1))[(Subscript[F, (Subscript[X, 2])])((Subscript[X, 2]))]}

This is one of several motivations for considering composite maps of the form y = (G^(-1))[F(x)], where F, G are CDFs.

Topics of the talk - see how far we get

Direct Construction of Quantile Functions themselves

Example 1: Recap Direct Student (Shaw, JCF, Vol 9, 4, pp.37-73,2006)

Distributional Alchemy - from one Quantile Function to another

Example 1: Sample Transmutation: Student from Normal via 
The "Exact Cornish-Fisher Map"

Example 2: Rank Transmutation: Introducing Skewness to Any Distribution without Gram-Charlier issues

Construction of Quantile Functions themselves

Construction example 1:
Student Quantile Functions, iCDFs for even n

We can give special methods for the Student T distribution CDF when the degrees of freedom n is an even integer. Note EP's talk emphasized role of (Subscript[T, 4]) case. First we have an iterative specification of the CDF

a[0,  n_] := Gamma[(n + 1)/2]/Gamma[n/2]/
          Sqrt[n Pi];
a[k_,  
          n_] := a[k, 
            n] = (n - 2 k)/n/(2  k + 1) a[k - 1,  
              n];

    
TCDF[n_, x_] := 
        1/2  + 
 
          x*Sum[a[p,  n] x^(2  p),  {p,  0,  
                  n/2 - 1}]/(1 + 
                    x^2/n)^(1/2 (n - 1));

    
{Simplify[TCDF[4, x]],  Simplify[D[TCDF[4, x], x]],  
      TCDF[4, -300] // N,  TCDF[4, 300] // N}

Solve TCDF[n, x]  =  u , a polynomial problem! Let a  =  4  u(1 - u) and p = n + (x^2)

P[n_,  x_]  :=  
        x*Sum[a[p,  n] x^(2  p),  {p,  0,  
              n/2 - 1}];
R[n_,  p_]  :=  
        Expand[PowerExpand[
4*
                  n^(n - 1) (P[n, x]) ^2   +  
                p^(n - 1) (a - 1) /. 
              x  ->  Sqrt[p - n]]] ;
TraditionalForm[
      Table[{2 k,  R[2  k, p] == 0},  {k,  1, 7}]]

and so on. We get interesting family of sparse polynomials. Almost half the coefficients are missing. Two exact solutions n = 2, 4, and n = 6, 8, 10, … easy by Newton-Raphson! (See JCF 2006, E. Platen's new book, for 1,2,4!)

Distributional Alchemy

Sample Transmutation

Given one "base" distribution function, say Φ(x) - the base might typically be Normal - and another distribution F(x), we define a sample transmutation mapping (Subscript[T, S]) by the identity

(F^(-1))[
        U] = (Subscript[T, S])((Φ^(-1))(U)), 
    i . e . ,   (Subscript[T, S])(
        z) = (F^(-1))(Φ(z)) = (Subscript[Q, F])(Φ(
            z))

where 0 ≤ U ≤ 1 and z is in its appropriate range (the real line in the Normal case). So if we have the  Φ Quantile function we can get the (Subscript[Q, F]) quantile by post-applying (Subscript[T, S]). This function "transmutes" samples from one distribution into samples from another. We might have a decent expression for Φ but (Subscript[Q, F]) may well be hard to nail down. But we will show some cunning computer algebra for getting series, and we can compose this with a series for Φ to get a good series for (Subscript[T, S]). 

This is not a new idea. It has previously found expression in the asymptotic setting via the use of Cornish-Fisher expansions, but we will now see how to use it in an essentially exact setting via the use of computer symbolic algebra.

Rank Transmutation

There is no particularly good reason why a transmutation mapping should be applied after applying a standard quantile rather than before. So we can define a corresponding rank transmutation mapping (Subscript[T, R]) by the following relationship

(F^(-1))[
        U] = (Φ^(-1))((Subscript[T, R])(U)),  
    i . e . ,   (Subscript[T, R])(
        u) = Φ((F^(-1))(z)) = Φ((Subscript[Q, F])(
            u))

This will allow us to introduce modulations into a distribution in an exact way, and potentially avoid the use of asymptotic (Edgeworth/Gram-Charlier [EGC]) methods and their problems. Note that equation (6) only makes sense if the two distributions have the same sample space.

Why think about these maps?

There are several reasons, some of which will represent a mathematical formalization of some things already being done by practitioners:

1. It helps us generate new (hard) Quantile functions from old (easy) ones for easy, QMC or copula applications;
2. We do not have to use the unit interval as a standard domain - we can change variables e.g. to Gaussian real line;
3. We can transmute a given sample to assess distributional risk in pricing/risk calculation, and avoid MC noise in much as the same way as Greeks are calculated in Monte Carlo;
4. Introduce skewness to standard distributions in a natural way;
5. Take the idea of QQ plot and transform it from a distributional assessment tool to a PDF generation tool;
6. Avoid pecularities of GC methods such as negative PDFs;s
7. Our original example based on copula simulation.

Example1:  Sample Transmutation: Student from Normal via the "Exact Cornish-Fisher Map"

We consider the Student T distribution in the notation of Shaw (2006). Given 0 < u < 1, we set

\(TraditionalForm\`v == \(\((u - 1\/2)\)\ \@\(n\ π\)\ \(Γ(n\/2)\)\)\/\(Γ(\(n + 1\)\/2)\)\)

The quantile equation is then

v == (∫_0^x)(((1 + 
                      ((s^2)/n)))^(-(1/2) (n + 
                        1))) ds = x (Subscript[, 2])(Subscript[F, 1])((1/2),((n + 1)/2);(3/2);-((x^2)/n))= F(x, n)

Off[General::spell1];
Off[Series::esss];

F[x_,  n_]  :=  
        x*Hypergeometric2F1[1/2, (n + 1)/2, 
            3/2, -(x^2/n)];
SCDF[x_,  n_]  :=  
      1/2 +  1/Sqrt[n Pi]*Gamma[(n + 1)/2]/Gamma[n/2] 
          F[x, n]

Simplify[D[SCDF[x, n], x]]

Integrate[%,  {x,  -Infinity,  Infinity},  
      Assumptions  ->  n > 0]

The inversion of a series for such a CDF can be carried out step by step in any computer algebra system, following the methods described in on-line information at KCL web site. In Mathematica you can sometimes just ask for the inverse as a series:

Map[Factor,  
      InverseSeries[Series[F[x, n],  {x,  0,  9}],  v]]

With a bit more work a corresponding tail series can be developed. Requires a little intervention. (See Shaw, J. Comp Fin. 2006 for details). So we have the Quantile function for general real n. What we want to discuss here is the sample transmutation mapping.

One example of what sample transmutation then gives us is a relatively painless way of extracting a high order Cornish-Fisher expansion, and even a non-asymptotic form of it. Compare the following with the usual rash of high-order moments and Hermite functions in the case of the Student:

Φ(z) - (1/2) == 
      (1/2) erf((z/(Sqrt[2])))

Phi[z_] := 1/2(1+Erf[z/Sqrt[2]]);
PhiMinusHalf[z_] := Erf[z/Sqrt[2]]/2;
v[pmh_,n_]:=(pmh*Sqrt[n*Pi]*Gamma[n/2])/Gamma[(n+1)/2]

and we recall the form of F(x, n) and set

QuantileF[v_,  n_,  
          truncation_] := 
InverseSeries[
          Series[F[x, n],  {x,  0,  truncation}],  
          v];

    Transmutation[z_,  n_,  trunca_,  
        truncb_] :=  
Module[{QF,  toexp}, 
        QF  =  QuantileF[v,  n,  trunca]; 
toexp  = 
          Normal[QF /. {v -> v[PhiMinusHalf[z], n], 
                m -> n}]; 
Series[
          toexp, {z, 0, truncb}]]

rawform  =  
        Transmutation[z, n, 10, 10];
CornishFisherExpansion  = 
      Map[
Together,  
        Normal[Series[rawform,  {n,  Infinity,  4}]]]

That can be found in Abramowitz and Stegun! But you wont find this one anywhere (prove me wrong!). Even my laptop has to think for a bit.

rawform  =  
        Transmutation[z, n, 20, 20];
CornishFisherExpansion  = 
      Map[
Together,  
        Normal[Series[rawform,  {n,  Infinity,  9}]]]

We get Cornish-Fisher to a high order with little effort! This last formula extends (and corrects) Equation (75) of Shaw (JCF, 2006), and extends the published results (see e.g. A&S). It is an example for the Student but we can do this trick for all kinds of pairs. But note also we did NOT have to expand in powers of n: we have a raw form as:

rawseries  =  Series[rawform,  {z,  0,  19}];

Series[rawseries,  {z,  0,  5}]

Much of the complication and lack of precision in the Cornish-Fisher expansion arises from it having an unnecessary expansion of all these gamma functions in inverse powers of n. We have literally added up this part of the expansion. Let's do it for a decidedly non-asymptotic n = 3! This is a good case to consider as the Quantile Functions are all known in closed form for n = 1,2,4 and the cases of even n are simple polynomial-solution problems. We look at the errors compared to the exact quantile function.

squeezer[z_] = 
        N[Normal[rawseries /.  n -> 3],  
          20];
Squeezer[z_]  :=  squeezer[z];
Squeezer[z]

Note that there is NO approximation in the coefficients here. The only approximation is in the number of terms we take.

ExactQuantileF[y_,  n_]  :=  
      Module[{arg  =  
            If[y < 1/2,  2 y,  2 (1 - y)]}, 
Sign[
            y - 1/2] 
          Sqrt[n*(1/InverseBetaRegularized[arg,  n/2,  1/2] - 
                  1)]]

explot  =  
        Plot[{ExactQuantileF[u, 3], 
            Squeezer[Sqrt[2]*InverseErf[2*u - 1]]},  {u,  0.0001,  
            0.9999},  PlotRange  ->  All];

There really are two plots here, but we cannot see the difference from this overlay. So we plot the difference, and we have to look well into the tail to see it:

explot  =  
        Plot[ExactQuantileF[u, 3] - 
            Squeezer[Sqrt[2]*InverseErf[2*u - 1]],  {u,  0.99,  0.9999},  
          PlotRange  ->  All];

These are absolute errors - the relative error is much smaller. So we have the possibility of using truncated polynomial expressions of Cornish-Fisher type in non-asymptotic domain. The errors are under investigation, as is the possibility of using a better sample transmutation mapping in the tail region.

We also have a variation of a multivariate T with general degrees of freedom in the marginals, based on a Gaussian copula - just squeeze the mixed Gaussians using a high-order Cornish-Fisher expansion.

explot  =  
        Plot[100*(Squeezer[Sqrt[2]*InverseErf[2*u - 1]]/
                  ExactQuantileF[u, 3] - 1),  {u,  0.9999,  
            0.9999999},  PlotRange  ->  All];

Example 2: Rank Transmutation and introduction of skewness into any distribution

To define a rank transmutation mapping in complete generality, suppose that we have two distributions with a common sample space, with CDFs (Subscript[F, 1]) and (Subscript[F, 2]). We can form

(Subscript[G, R12])(
        u) = (Subscript[F, 2])((Subscript[F, 1]^(-1))(u)) ,   (Subscript[G, R21])(
        u) = (Subscript[F, 1])((Subscript[F, 2]^(-1))(u))

and this pair of maps the unit interval I = [0, 1] into itself and under suitable assumptions are mutual inverses and satisfy (Subscript[G, i j])(0) = 
        0,  (Subscript[G, i j])(1) = 1. We shall further assume that these rank transmuation maps are continuously differentiable, as otherwise it is easy to check that a transmuted density may be discontinuous. We could consider all kinds of such maps arising from a particular choice of the (Subscript[F, i]), but here, rather, we shall postulate an interesting form. Consider, for \(TraditionalForm\`\(\(|\)\(λ\)\(|\)\(\(≤\)\(\ \)\(1\)\)\)\),

(Subscript[G, R12])(u) = 
      u + λ u (1 - u)

This has the consequence that the CDFs are related by

(Subscript[F, 2])(
        x) = (1 + λ) (Subscript[F, 1])(
            x) - λ (((Subscript[F, 1])(x))^2)

and the sampling algorithm remains tractable as the Quantile functions are related by

\(TraditionalForm\`\(F\_2\%\(-1\)\)(u) = 
      F\_1\%\(-1\)\ \((\(G\_R21\)(u))\), \ \ \ \(G\_R21\)(
        u)\  = \ \(1 + λ - \@\(\((λ + 1)\)\^2 - 4\ λ\ u\)\)\/\(2\ λ\)\)

There are two important extremal cases. First, if λ = -1, then (Subscript[G, R12])(u) = (u^2) and (Subscript[F, 2])(x)  =  (((Subscript[F, 1])(x))^2) and we recognize that the ditribution of (Subscript[F, 2]) corresponds to that of the maximum of two independent copies of the (Subscript[F, 1])distribution. Correspondingly λ = +1 generates the distribution of the minimum. We could use higher powers and distributions of max/min of k. But the quadratic is tractable and natural, and is easily seen to give us a natural skewing mechanism.

A Skew-Uniform example

H[z_] := z; h[z_]  :=  1;
H2[
          z_,  λ_]  :=  (1 + λ) 
            H[z] - λ H[z]^2;
h2a[x_,  λ_] =  
        D[H2[x,  λ],  x];
h2[x_,  λ_]  :=  h2a[x,  λ]

h2[x,  p]

Do[Plot[h2[x, p],  {x,  0,  1},  
        PlotRange  ->  {0, 2}],  {p,  -1,  1,  
        0.5}]

This is a very natural interpretation of a "skew-uniform".

A Skew-Normal example

Φ[z_] := 1/2 (1 + Erf[z/Sqrt[2]]);

f2[x_,  λ_]  :=  
      1/Sqrt[2 Pi] Exp[-x^2/2] (1 + λ - 
            2  λ Φ[x])

The first 4 moments are then

fourmoments = Integrate[{x, x^2,  x^3,  
            x^4} f2[x,  λ],  {x,  -Infinity,  
          Infinity}]

so that the variance is

fourmoments[[2]] - fourmoments[[1]]^2

and the skewness is

Integrate[(x - fourmoments[[1]])^3 ftwo[
          x,  λ],  {x,  -Infinity,  
        Infinity}]

Plot[
        Evaluate[
          Table[ftwo[x,  λ],  {λ,  0,  1,  
              0.5}]],  {x,  -4,  4}];

If we want to standardize this distribution to have zero mean and unit variance we form the following from an obvious linear transformation:

f3[x_,  λ_]  :=  
      Sqrt[1 - λ^2/Pi]*
        f2[x Sqrt[1 - λ^2/Pi] - λ/
              Sqrt[Pi], λ]

We could, if we wanted to, do an series expansion of this around the standard Normal distribution, in powers of both x and λ, in order to see the relation to an Gram-Charlier-type expansion, but, as with the Cornish-Fisher expansion on the Quantile, this is now unnecessary!

We could also just as easily do skew-Student or skew anything else, and the Monte Carlo sampling is precisely as tractable as it is for the base distribution. These ideas should of course be compared with the work of many others, in particular A. Azzalini and co-workers, who use a different mechanism. See his web site at http://tango.stat.unipd.it/SN/ for details of their approach, which is based on multiplying a density by a CDF. Our approach coincides with his in certain special cases (e.g. distributions of max/min are also contained) but appears to be different in general.

Summary

Computer algebra is a powerful technique for:

(a) generating quantile functions in the first place.
(b) turning samples from one distribution into another;
(c) turning the ranks of one distribution into the ranks of another, e,g, to introduce skewness in a universal way;

These techniques are well adapted for Quasi-Monte-Carlo and Copula simulation methods.