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.