Simulation of the fractional noncentral Wishart distribution

Posted on December 9, 2017 by Stéphane Laurent

Tags: R, maths, statistics

It is well known how to simulate the noncentral Wishart distribution when the number of degrees of freedom \(\nu\) and the dimension \(d\) satisfy \(\nu > 2d-1\), or when \(\nu \geq d\) is an integer. In their paper Exact and high-order discretization schemes for Wishart processes and their affine extensions, Ahdida & Alfonsi provide a method that allows to simulate the Wishart process of dimension \(d\) for any number of degrees of freedom \(\nu \geq d-1\) and without restrictions on the other parameters. This method allows to simulate the noncentral Wishart distribution, in the way we will expose now.

Two properties of the noncentral Wishart distribution

We will need the two following properties.

Recall that the characteristic function of the noncentral Wishart distribution \(\mathcal{W}(\nu, \Sigma, \Theta)\) at \(Z\) is \[ \phi_{\nu,\Sigma,\Theta}(Z) = \frac{\exp\Bigl(i\textrm{tr}\bigl({(I_d - 2i Z\Sigma)}^{-1}Z\Theta\bigr)\Bigr)}{{\det(I_d - 2i Z\Sigma)}^{\frac{\nu}{2}}}. \]

\(\bullet\) First property of the Wishart distribution. Using the characteristic function, it is easy to check that \(A W A' \sim \mathcal{W}(\nu, A\Sigma A', A\Theta A')\) when \(W \sim \mathcal{W}(\nu, \Sigma, \Theta)\).

\(\bullet\) Second property of the Wishart distribution. Using the characteristic function, it is not hard to check that if \(W_1\) and \(W_2\) are two random matrices such that \(W_1 \sim \mathcal{W}(\nu, \Sigma_1, \Theta)\) and \((W_2 \mid W_1) \sim \mathcal{W}(\nu, \Sigma_2, W_1)\), then \(W_2 \sim \mathcal{W}(\nu, \Sigma_1+\Sigma_2, \Theta)\). This result is proved in A&A’s paper only for the covariance matrices \(J_d^i\) we will see later, by means of another method.

Let’s prove this result with the characteristic function. The conditional characteristic function of \(W_2\) given \(W_1\) at \(Z\) is \[ \frac{\exp\Bigl(i\textrm{tr}\bigl({(I_d - 2i Z\Sigma_2)}^{-1}ZW_1\bigr)\Bigr)}{{\det(I_d - 2i Z\Sigma_2)}^{\frac{\nu}{2}}}. \] The characteristic function of \(W_2\) is obtained by taking the expectation of this expression, and doing so we get \[ \frac{\phi_{\nu,\Sigma_1,\Theta}\bigl({(I_d - 2i Z\Sigma_2)}^{-1}Z\bigr)}{{\det(I_d - 2i Z\Sigma_2)}^{\frac{\nu}{2}}} = \frac{\exp\left(i\textrm{tr}\Bigl({\bigr(I_d - 2i {(I_d - 2i Z\Sigma_2)}^{-1}Z\Sigma_1\bigl)}^{-1}{(I_d - 2i Z\Sigma_2)}^{-1}Z\Theta\Bigr)\right)}{{\det(I_d - 2i Z\Sigma_2)}^{\frac{\nu}{2}} {\det(I_d - 2i {\bigl(I_d - 2i Z\Sigma_2)}^{-1}Z\Sigma_1\bigr)}^{\frac{\nu}{2}}}. \]

It is easy to check that the denominator is \({\det\bigl(I_d - 2iZ(\Sigma_1+\Sigma_2)\bigr)}^{\frac{\nu}{2}}\). The expression inside \(\textrm{tr}(\ldots)\) at the numerator is \[ {\Bigl((I_d - 2i Z\Sigma_2)\bigr(I_d - 2i {(I_d - 2i Z\Sigma_2)}^{-1}Z\Sigma_1\bigl)\Bigr)}^{-1}Z\Theta = {\bigr(I_d - 2iZ(\Sigma_1+\Sigma_2)\bigl)}^{-1}Z\Theta. \]

Thus we find that the the characteristic function of \(W_2\) is \(\phi_{\nu, \Sigma_1+\Sigma_2,\Theta}\), that is to say \(W_2 \sim \mathcal{W}(\nu, \Sigma_1+\Sigma_2, \Theta)\).

A&A’s simulation method

A&A’s simulation method has three steps:

it firstly gives an algorithm to simulate \(\mathcal{W}(\nu, J_d^1, \Theta)\), denoting by \(J_d^i\) the \(d \times d\) covariance matrix whose all entries are equal to zero except the \((i,i)\)-entry which is equal to one;
using the first step and the two properties of the Wishart distribution that we have seen, it provides a way to simulate \(\mathcal{W}(\nu, I_d^n, \Theta)\) where \(I_d^n = J_d^1 + \ldots + J_d^n\);
using the second step and the first property of the Wishart distribution that we have seen, it provides a way to simulate \(\mathcal{W}(\nu, \Sigma, \Theta)\) for any covariance matrix \(\Sigma\).

\(\bullet\) Simulation of \(\mathcal{W}(\nu, J_d^1, \Theta)\). This algorithm runs as follows. Let \((L,M,P)\) be an extended Cholesky decomposition of \(\Theta_{2:d,2:d}\). Set \(Q = \begin{pmatrix} 1 & 0 \\ 0 & P \end{pmatrix}\) and \(\widetilde{\Theta} = Q\Theta Q'\), then set \(u = L^{-1}\widetilde{\Theta}_{1, 2:(r+1)}'\) and \(v = \widetilde{\Theta}_{1,1} - \sum_{i=1}^r u_i^2\). Take \(Z_1, \ldots, Z_r \sim_{\text{iid}} \mathcal{N}(0,1)\) and set \(G_i = u_i + Z_i\). Finally, take \(X \sim \chi^2_{\nu-r, v}\) (noncentral chi-squared distribution) independent of the \(Z_i\), and set \[ W = Q' \begin{pmatrix} 1 & 0 & 0 \\ 0 & L & 0 \\ 0 & M & I_{d-r-1} \end{pmatrix} \begin{pmatrix} X + \sum_{i=1}^r G_i^2 & G' & 0 \\ G & I_r & 0 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} 1 & 0 & 0 \\ 0 & L' & M' \\ 0 & 0 & I_{d-r-1} \end{pmatrix} Q. \] Then A&A have shown that \(W \sim \mathcal{W}(\nu, J_d^1, \Theta)\).

\(\bullet\) Simulation of \(\mathcal{W}(\nu, I_d^n, \Theta)\). Let \(P\) be the permutation matrix exchanging rows \(1\) and \(2\). Using the previous algorithm, simulate \(W_1 \sim \mathcal{W}(\nu, J_d^1, \Theta)\). By the first property of \(\mathcal{W}\) we have seen, \(P W_1 P \sim \mathcal{W}(\nu, J_d^2, P\Theta P)\). Then, still using the previous algorithm, simulate \((W_2 \mid W_1) \sim \mathcal{W}(\nu, J_d^1, P W_1 P)\). By the second property of \(\mathcal{W}\) we have seen, \(W_2 \sim \mathcal{W}(\nu, I_d^2, P \Theta P)\). And by the first property, \(P W_2 P \sim \mathcal{W}(\nu, I_d^2, \Theta)\). Continuing so on, we can simulate \(\mathcal{W}(\nu, I_d^n, \Theta)\) for any \(n \leq d\).

\(\bullet\) Simulation of \(\mathcal{W}(\nu, \Sigma, \Theta)\). Finally, given any covariance matrix \(\Sigma\) of rank \(n\), take the \(\widetilde{C}\) matrix of an extended Cholesky decomposition of \(\Sigma\) with permutation matrix \(P\), and set \(A = P'\widetilde{C}\). Simulate \(Y \sim \mathcal{W}\bigl(\nu, I_d^n, A^{-1}\Theta{(A^{-1})}'\bigr)\) with the previous algorithm and finally set \(W = AYA'\), so that \(W \sim \mathcal{W}(\nu, \Sigma, \Theta)\) by the first property and by the property of the Cholesky decomposition.

Checking

The algorithm is implemented in my package matrixsampling. Let’s try it.

library(matrixsampling)
p <- 6
nu <- 6.3
Sigma <- toeplitz(p:1)
Theta <- matrix(1, p, p)
nsims <- 100000
W <- rwishart(nsims, nu, Sigma, Theta)

As expected, the average simulated matrix is close to the theoretical mean \(\nu \Sigma + \Theta\):

round((nu*Sigma + Theta) - apply(W, 1:2, mean), 2)
##       [,1]  [,2]  [,3]  [,4]  [,5]  [,6]
## [1,] -0.05 -0.04  0.00  0.00  0.02  0.04
## [2,] -0.04 -0.06 -0.03 -0.03  0.00  0.01
## [3,]  0.00 -0.03 -0.04 -0.06 -0.05 -0.05
## [4,]  0.00 -0.03 -0.06 -0.07 -0.04 -0.07
## [5,]  0.02  0.00 -0.05 -0.04 -0.03 -0.05
## [6,]  0.04  0.01 -0.05 -0.07 -0.05 -0.09

Let’s compare the theoretical characteristic function to its approximation obtained from the simulations:

z <- seq(0.001, 0.004, length.out = 20)
Z <- sapply(z, function(z){
  z*diag(p) + matrix(z, p, p)
}, simplify=FALSE)
tr <- function(A) sum(diag(A))
Phi_theoretical <- sapply(Z, function(Z){
             complexplus::Det(diag(p) - 2*1i*Z%*%Sigma)^(-nu/2) * 
             exp(1i*tr(solve(diag(p) - 2*1i*Z%*%Sigma) %*% Z %*% Theta))
})
Phi_sims <- sapply(Z, function(Z){
  mean(apply(W, 3, function(W){
    exp(1i*tr(Z%*%W))
  }))
})
layout(t(1:2))
plot(z, Re(Phi_theoretical), type="o", pch=19)
lines(z, Re(Phi_sims), type="o", pch=19, col="red")
plot(z, Im(Phi_theoretical), type="o", pch=19)
lines(z, Im(Phi_sims), type="o", pch=19, col="red")