Efficient Generation of Random Variates from a Copula Distribution

Question

I have a 7-asset portfolio for which I want to simulate daily log-deltas using a Student T copula. The marginal distributions are all either Stable or TsallisQGaussian. Using NMaximize, I have already determined that the copula dof = 9.

The asset correlation matrix, corrmat, is:

corrmat = {{1., 0.491789, 0.382652, -0.580449, 0.34607, -0.0887926, 
  0.292343}, {0.491789, 1., 0.475475, -0.713583, 0.364887, -0.112926, 
  0.372419}, {0.382652, 0.475475, 1., -0.695327, 0.370684, -0.0894264,
   0.359768}, {-0.580449, -0.713583, -0.695327, 1., -0.498612, 
  0.1529, -0.482753}, {0.34607, 0.364887, 0.370684, -0.498612, 
  1., -0.0433323, 0.291971}, {-0.0887926, -0.112926, -0.0894264, 
  0.1529, -0.0433323, 1., -0.0954981}, {0.292343, 0.372419, 
  0.359768, -0.482753, 0.291971, -0.0954981, 1.}}

from which the conditioned VCV matrix, vcvmat, is:

vcvmat = {{0.9586172917, -0.03410852909, -0.2933963066, 0.1208800096, 
  0.1706880212, 0.1239885803, -0.06784705084}, {-0.03410852909, 
  0.6979797397, 
  0.02009209742, -0.001634295474, -0.02837172149, -0.02564790389, 
  0.01303223508}, {-0.2933963066, 0.02009209742, 
  0.7314875880, -0.1494059087, -0.3716208429, -0.2204546813, 
  0.1000107891}, {0.1208800096, -0.001634295474, -0.1494059087, 
  0.8597301696, 0.1366021203, 
  0.07930474269, -0.05334637064}, {0.1706880212, -0.02837172149, \
-0.3716208429, 0.1366021203, 1.136447366, 
  0.2592871862, -0.06289877472}, {0.1239885803, -0.02564790389, \
-0.2204546813, 0.07930474269, 0.2592871862, 
  1.154068011, -0.02358813328}, {-0.06784705084, 0.01303223508, 
  0.1000107891, -0.05334637064, -0.06289877472, -0.02358813328, 
  0.8678362987}}

The respective marginal distributions are:

mda = StableDistribution[1, 1.66576, -0.205451, -0.00149914, 0.00932388]
mdb = TsallisQGaussianDistribution[0.000198498, 0.0124003, 1.43148]
mdc = StableDistribution[1, 1.84146, -0.453362, 0.000669985, 0.00802551]
mdd = TsallisQGaussianDistribution[-0.00129644,0.0112156, 1.52237]
mde = TsallisQGaussianDistribution[0.000625683, 0.0147539, 1.46668]
mdf = StableDistribution[1, 1.45877, 0.253482, 0.000708238, 0.00111962]
mdg = TsallisQGaussianDistribution[0.00064191, 0.0121773, 1.52216]

Given these values, the entire copula expression is:

tCopula = CopulaDistribution[{"MultivariateT", vcvmat, 9}, 
              {mda, mdb, mdc, mdd, mde, mdf, mdg}];

So far, so good.

The problem is the requisite calculation time. Here are the outputs:

In: Timing[x = RandomVariate[tCopula, 100];]
Out: {8.917131, Null}

In: Timing[y = RandomVariate[tCopula, 1000];]
Out: {92.197472, Null}

In: Timing[z = RandomVariate[tCopula, 10000];]
Out: {921.512002, Null}

The computation times are roughly linear, and since 10,000 samples is a bare minimum for meaningful analysis, simulating a mere 7-asset portfolio is at least a fifteen minute proposition.

Is there a way to call this simulation more efficiently within Mathematica 9? FWIW, I'm using Mathematica on a Mac, so Finance Platform isn't an option.

---

Sasha,

Thanks for taking the time to unbundle all that functionality.

For the benefit of anyone else who wants to run the simulation, be aware that Mma threw a non-Hermitian error when I evaluated the data argument the way Sasha initially set it up. I went back and adjusted the machine precision of my VCV matrix, and I created a separate variable, v2mat, for Sasha's expression

vcvmat/KronekerProduct[Sqrt[Diagonal[vcvmat]],Sqrt[Diagonal[vcvmat]]]

Then I evaluated

RandomVariate[MultivariateTDistribution[v2mat, dof], 10^4]

and it worked fine. For purposes of comparison, I reran the original tCopula argument for a 10^4 x 7 set of observations:

In: Timing[data1 = RandomVariate[tCopula, 10000];]
Out: {888.577837, Null}

versus a 10^4 x 7 set of values from Sasha's method:

In: Timing[data2 = Transpose[MapThread[#1[#2]&, {qfs, Transpose[data]}]];]
Out: {4.322776, Null}

Quite a difference.

Here's a graphical juxtaposition after running the same Mahalanobis distance multivariate outlier routine on both data sets.

Not only is Sasha's method much faster, if anything it appears to yield more representative random variates.

Answer 1

A nice question.

Sampling from tCopula is done in stages. First a sample is generated from the copula with uniform marginal distributions, and then quantiles of appropriate marginal distributions are applied to the respective slots.

Most of the time goes into evaluation of these quantiles, and they are expensive to compute. Being interested in $\geqslant 10^4$ sample points, the performance can be enhanced by interpolating quantile functions of marginals using monotonic interpolation:

FritschCarlsonInterpolation[data_] := 
 Module[{del, slopes, tau}, del = #2/#1 & @@@ Differences[data];
  slopes = Flatten[{del[[1]], MovingAverage[del, 2], del[[-1]]}];
  tau = MapThread[Min, Through[{Append, Prepend}[
      Min[#, 1] & /@ (3 del/(Norm /@ Partition[slopes, 2, 1])), 1]]];
  Interpolation[ Transpose[{List /@ data[[All, 1]], data[[All, -1]], slopes tau}], 
   InterpolationOrder -> 3, Method -> "Hermite"]]

Let's build interpolants in the bulk of the quantile space.

qgrid = With[{n = 100}, 
   Join[{0.001, 0.002, 0.005, 0.007}, 
    Range[1./n, (n - 1.)/n, 1./n], {0.993, 0.995, 0.998, 0.999}]];

AbsoluteTiming[
 interpolants = 
   Table[FritschCarlsonInterpolation[
     Table[{q, Quantile[di, q]}, {q, qgrid}]], {di, {mda, mdb, mdc, 
      mdd, mde, mdf, mdg}}];]

(* {20.770029, Null} *)

Now build quantile functions by using interpolants in the bulk, and the original Quantile[dist,q] for tails:

AbsoluteTiming[
 qfs = MapThread[
    Function[{di, qf}, 
     Function @@ {q, 
       Piecewise[{{qf[q], Min[qgrid] <= q <= Max[qgrid]}}, Quantile[di, q]], 
       Listable}], {{mda, mdb, mdc, mdd, mde, mdf, mdg}, 
     interpolants}];]

(* {0.010000, Null} *)

Unfortunately it seems one can beat multivariate T CopulaDistribution at random number generation:

AbsoluteTiming[
 data = RandomVariate[
   MultivariateTDistribution[
    vcvmat/KroneckerProduct[Sqrt[Diagonal[vcvmat]], 
      Sqrt[Diagonal[vcvmat]]], dof], 10^5];
 data = ArrayReshape[CDF[StudentTDistribution[dof], Flatten[data]], 
   Dimensions[data]];]

(* {0.250000, Null} *)

The variable data now contains $10^5$ sample data points from the multivariate-T copula with 9 degrees of freedom and correlation structure determined by vcvmat and uniform marginals.

It now remains to apply the constructed quantile functions to these sample vectors:

In[49]:= AbsoluteTiming[
 data = Transpose[MapThread[#1[#2] &, {qfs, Transpose[data]}]];]

Out[49]= {62.410087, Null}

Thus we got around to generate $10^5$ 7-tuples in under 2 minutes of the total time. You can adjust the interpolation grid to improve the accuracy as needed.

Answer 2

With the help of comments on this and other stackexchange pages I managed to solve the problem of how to use custom distributions in things like CopulaDistribution (and other functions like RandomVariate, Expectation, etc), and given that it took me a couple of days’ hard slog I thought I’d share my discoveries with this community. Please excuse the flippant nature of the following – it started life as an email to a friend and colleague of mine… :)

(Also, please note that Random`DistributionVector, mentioned by Sasha in one of the linked posts above, appears to have been changed to Random`Private`DistributionVector in the production release of v9.0.1)

I needed a Mathematica version of Tukey’s “g-and-h” distribution, which is a transform of Z ~ N[0,1] ...

Xgh = A + B/g(Exp[g*Z] - 1)Exp[(h*Z^2)/2]

OK, you’d think just use

ghd = TransformedDistribution[Xgh, Distributed[Z, NormalDistribution[0,1]]

would you not?

Turns out that, even though the above is valid, it takes a very (VERY) long time to compute CDF[ghd, X], or even Quantile[ghd, p]. I tried a number of tricks including Compile, taking a Taylor series, etc, but I’ve concluded (reluctantly) it’s not viable to use TransformedDistribution for g-and-h.

OK then, you’d think just use the g-and-h cdf prob in

ghd = ProbabilityDistribution[{“CDF”, prob}, {x, -Infinity, Infinity}] 

would you not?

Turns out that a symbolic form of the pdf and/or cdf is problematic, and using

prob = CDF[NormalDistribution[0, 1], (Z /. FindRoot[Xgh == x, {Z,0}])] 

doesn’t cut it with Mathematica. It wants something symbolic inside ProbabilityDistribution.

So, I reverted back to my “hand-rolled” (custom) distribution that I’d written before the new fancy enhanced Distributions came out in v8 & v9. It already had all the std Mathematica things like CDF (implemented as the FindRoot of Xgh), PDF, InverseCDF, Random, Quantile, Mean, Variance, etc, via the TagSet feature (ie, defining “up values” of the distribution name). Eg...

GandHDistribution/: InverseCDF[GandHDistribution[A_, B_, g_, h_], fraction_] :=
  Xgh[A, B, g, h, Quantile[NormalDistribution[0, 1], fraction]] /; 
    DistributionParameterQ[GandHDistribution[A, B, g, h]]

But it turns out that there are a lot of undocumented test functions used by CopulaDistribution that you also need to define. Also, some of the pre-v8 internals had changed, e.g., ParameterQ is now DistributionParameterQ. I found out what was needed by Unprotect[CopulaDistribution] and then removing the ReadProtected attribute so that I can see all the required up-value definitions. They appear to be (aside from the “well-known” things like InverseCDF and DistributionParameterQ mentioned above)...

GandHDistribution/: DistributionParameterQ[GandHDistribution[A_, B_, g_, h_]] := And[
    If[FreeQ[N[A], Complex], True,
       Message[GandHDistribution::realparm, A]; False],
    If[FreeQ[N[B], Complex], True,
       Message[GandHDistribution::realparm, B]; False],
    If[FreeQ[N[g], Complex], True,
       Message[GandHDistribution::realparm, g]; False],
    If[FreeQ[N[h], Complex], True,
       Message[GandHDistribution::realparm, h]; False]
];

GandHDistribution::posparm =
  "Parameter `1` is expected to be positive."

GandHDistribution::realparm = "Parameter `1` is expected to be real."

GandHDistribution/: 
  DistributionParameterAssumptions[GandHDistribution[A_, B_, g_, h_]]:=
Element[{A,B,g,h},Reals]   

GandHDistribution/: 
  RandomVariate[GandHDistribution[A_, B_, g_, h_], dim_] :=
    Module[
      {dimv=Flatten[{dim}] (*if dim is single int, convert to single element list*)},
      Map[Xgh[A, B, g, h, #]&, RandomVariate[NormalDistribution[0,1],dimv], {Length@dimv}]
    ] /; (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&];

GandHDistribution/: 
  Random`Private`
    DistributionVector[GandHDistribution[A_, B_, g_, h_],
                       n_Integer, prec_?Positive] :=
      Xgh[A, B, g, h, RandomVariate[NormalDistribution[0, 1], n, 
        WorkingPrecision -> prec]];

GandHDistribution/: 
 Statistics`CopulaDistributionDump` 
   UnivariateDistributionListQ[GandHDistribution[A_, B_, g_, h_]] := True;

GandHDistribution/: 
  Statistics`Library`
    ContinuousUnivariateDistributionQ[GandHDistribution[A_, B_, g_, h_]] := True;

GandHDistribution/: 
  Statistics`Library`
    DiscreteUnivariateDistributionQ[GandHDistribution[A_, B_, g_, h_]] := False;

GandHDistribution/: 
  Statistics`Library`
    ContinuousMultivariateDistributionQ[GandHDistribution[A_, B_, g_, h_]] := False;

GandHDistribution/: 
  Statistics`Library`
    DiscreteMultivariateDistributionQ[GandHDistribution[A_, B_, g_, h_]] := False;

GandHDistribution/: 
  Statistics`Library`
    DistributionNParameterQ[GandHDistribution[A_, B_, g_, h_]]:=
     DistributionParameterQ[GandHDistribution[A, B, g, h]];

Now it works like a dream – the CopulaDistribution function is fantastic.