qfrm {qfratio}

Moment of ratio of quadratic forms in normal variables

Description

qfrm() is a front-end function to obtain the (compound) moment of a ratio of quadratic forms in normal variables, i.e., \mathrm{E} \left( \frac{(\mathbf{x^\mathit{T} A x})^p }{(\mathbf{x^\mathit{T} B x})^q} \right) , where \mathbf{x} \sim N_n(\bm{\mu}, \mathbf{\Sigma}). Internally, qfrm() calls one of the following functions which does the actual calculation, depending on \mathbf{A}, \mathbf{B}, and p. Usually the best one is automatically selected.

qfrm_ApIq_int(): For \mathbf{B} = \mathbf{I}_n and positive-integral p.

qfrm_ApIq_npi(): For \mathbf{B} = \mathbf{I}_n and non-positive-integral p (fraction or negative).

qfrm_ApBq_int(): For general \mathbf{B} and positive-integral p.

qfrm_ApBq_npi(): For general \mathbf{B} and non-integral p.

Usage

qfrm(
  A,
  B,
  p = 1,
  q = p,
  m = 100L,
  mu = rep.int(0, n),
  Sigma = diag(n),
  tol_zero = .Machine$double.eps * 100,
  tol_sing = tol_zero,
  ...
)

qfrm_ApIq_int(
  A,
  p = 1,
  q = p,
  m = 100L,
  mu = rep.int(0, n),
  use_cpp = TRUE,
  exact_method = TRUE,
  cpp_method = "double",
  nthreads = 1,
  tol_zero = .Machine$double.eps * 100,
  thr_margin = 100
)

qfrm_ApIq_npi(
  A,
  p = 1,
  q = p,
  m = 100L,
  mu = rep.int(0, n),
  error_bound = TRUE,
  check_convergence = c("relative", "strict_relative", "absolute", "none"),
  use_cpp = TRUE,
  cpp_method = c("double", "long_double", "coef_wise"),
  nthreads = 1,
  alphaA = 1,
  tol_conv = .Machine$double.eps^(1/4),
  tol_zero = .Machine$double.eps * 100,
  tol_sing = tol_zero,
  thr_margin = 100
)

qfrm_ApBq_int(
  A,
  B,
  p = 1,
  q = p,
  m = 100L,
  mu = rep.int(0, n),
  error_bound = TRUE,
  check_convergence = c("relative", "strict_relative", "absolute", "none"),
  use_cpp = TRUE,
  cpp_method = "double",
  nthreads = 1,
  alphaB = 1,
  tol_conv = .Machine$double.eps^(1/4),
  tol_zero = .Machine$double.eps * 100,
  tol_sing = tol_zero,
  thr_margin = 100
)

qfrm_ApBq_npi(
  A,
  B,
  p = 1,
  q = p,
  m = 100L,
  mu = rep.int(0, n),
  check_convergence = c("relative", "strict_relative", "absolute", "none"),
  use_cpp = TRUE,
  cpp_method = c("double", "long_double", "coef_wise"),
  nthreads = 0,
  alphaA = 1,
  alphaB = 1,
  tol_conv = .Machine$double.eps^(1/4),
  tol_zero = .Machine$double.eps * 100,
  tol_sing = tol_zero,
  thr_margin = 100
)

Arguments

A, B

Argument matrices. Should be square. Will be automatically symmetrized.

p, q

Exponents corresponding to \mathbf{A} and \mathbf{B}, respectively. When only one is provided, the other is set to the same value. Should be length-one numeric (see “Details” for further conditions).

m

Order of polynomials at which the series expression is truncated. M in Hillier et al. (2009, 2014).

mu

Mean vector \bm{\mu} for \mathbf{x}

Sigma

Covariance matrix \mathbf{\Sigma} for \mathbf{x}. Accommodated only by the front-end qfrm(). See “Details”.

tol_zero

Tolerance against which numerical zero is determined. Used to determine, e.g., whether mu is a zero vector, A or B equals the identity matrix, etc.

tol_sing

Tolerance against which matrix singularity and rank are determined. The eigenvalues smaller than this are considered zero.

...

Additional arguments in the front-end qfrm() will be passed to the appropriate “internal” function.

use_cpp

Logical to specify whether the calculation is done with C++ functions via Rcpp. TRUE by default.

exact_method

Logical to specify whether the exact method is used in qfrm_ApIq_int() (see “Details”).

cpp_method

Method used in C++ calculations to avoid numerical overflow/underflow (see “Details”). Options:

"double"

default; fastest but prone to underflow in some conditions

"long_double"

same algorithm but using the long double variable type; robust but slow and memory-inefficient

"coef_wise"

coefficient-wise scaling algorithm; most robust but variably slow

nthreads

Number of threads used in OpenMP-enabled C++ functions. 0 or any negative value is special and means one-half of the number of processors detected. See “Multithreading” in “Details”.

thr_margin

Optional argument to adjust the threshold for scaling (see “Scaling” in d1_i). Passed to internal functions (d1_i, d2_ij, d3_ijk) or their C++ equivalents.

error_bound

Logical to specify whether an error bound is returned (if available).

check_convergence

Specifies how numerical convergence is checked (see “Details”). Options:

"relative"

default; magnitude of the last term of the series relative to the sum is compared with tol_conv

"strict_relative" or TRUE

same, but stricter than default by setting tol_conv = .Machine$double.eps (unless a smaller value is specified by the user)

"absolute"

absolute magnitude of the last term is compared with tol_conv

"none" or FALSE

skips convergence check

alphaA, alphaB

Factors for the scaling constants for \mathbf{A} and \mathbf{B}, respectively. See “Details”.

tol_conv

Tolerance against which numerical convergence of series is checked. Used with check_convergence.

Details

These functions use infinite series expressions based on the joint moment-generating function (with the top-order zonal/invariant polynomials) (see Smith 1989, Hillier et al. 2009, 2014; Bao and Kan 2013), and the results are typically partial (truncated) sums from these infinite series, which necessarily involve truncation errors. (An exception is when \mathbf{B} = \mathbf{I}_n and p is a positive integer, the case handled by qfrm_ApIq_int().)

The returned value is a list consisting of the truncated sequence up to the order specified by m, its sum, and error bounds corresponding to these (see “Values”). The print method only displays the terminal partial sum and its error bound (when available). Use plot() for visual inspection, or the ordinary list element access as required.

In most cases, p and q must be nonnegative (in addition, p must be an integer in qfrm_ApIq_int() and qfrm_ApBq_int() when used directly), and an error is thrown otherwise. The only exception is qfrm_ApIq_npi() which accepts negative exponents to accommodate \frac{(\mathbf{x^\mathit{T} x})^q }{(\mathbf{x^\mathit{T} A x})^p} . Even in the latter case, the exponents must have the same sign. (Technically, not all of these conditions are necessary for the mathematical results to hold, but they are enforced for simplicity).

When error_bound = TRUE (default), qfrm_ApBq_int() evaluates a truncation error bound following Hillier et al. (2009: theorem 6) or Hillier et al. (2014: theorem 7) (for zero and nonzero means, respectively). qfrm_ApIq_npi() implements similar error bounds. No error bound is known for qfrm_ApBq_npi() to the author's knowledge.

For situations when the error bound is unavailable, a very rough check of numerical convergence is also conducted; a warning is thrown if the magnitude of the last term does not look small enough. By default, its relative magnitude to the sum is compared with the tolerance controlled by tol_conv, whose default is .Machine$double.eps^(1/4) (= ~1.2e-04) (see check_convergence).

When Sigma is provided, the quadratic forms are transformed into a canonical form; that is, using the decomposition \mathbf{\Sigma} = \mathbf{K} \mathbf{K}^T, where the number of columns m of \mathbf{K} equals the rank of \mathbf{\Sigma}, \mathbf{A}_\mathrm{new} = \mathbf{K^\mathit{T} A K}, \mathbf{B}_\mathrm{new} = \mathbf{K^\mathit{T} B K}, and \mathbf{x}_\mathrm{new} = \mathbf{K}^{-} \mathbf{x} \sim N_m(\mathbf{K}^{-} \bm{\mu}, \mathbf{I}_m) . qfrm() handles this by transforming A, B, and mu and calling itself recursively with these new arguments. Note that the “internal” functions do not accommodate Sigma (the error for unused arguments will happen). For singular \mathbf{\Sigma}, one of the following conditions must be met for the above transformation to be valid: 1) \bm{\mu} is in the range of \mathbf{\Sigma}; 2) \mathbf{A} and \mathbf{B} are in the range of \mathbf{\Sigma}; or 3) \mathbf{A} \bm{\mu} = \mathbf{B} \bm{\mu} = \mathbf{0}_n . An error is thrown if none is met with a singular Sigma.

The existence of the moment is assessed by the eigenstructures of \mathbf{A} and \mathbf{B}, p, and q, according to Bao and Kan (2013: proposition 1). An error will result if the conditions are not met.

Straightforward implementation of the original recursive algorithms can suffer from numerical overflow when the problem is large. Internal functions (d1_i, d2_ij, d3_ijk) are designed to avoid overflow by order-wise scaling. However, when evaluation of multiple series is required (qfrm_ApIq_npi() with nonzero mu and qfrm_ApBq_npi()), the scaling occasionally yields underflow/diminishing of some terms to numerical 0, causing inaccuracy. A warning is thrown in this case. (See also “Scaling” in d1_i.) To avoid this problem, the C++ versions of these functions have two workarounds, as controlled by cpp_method. 1) The "long_double" option uses the long double variable type instead of the regular double. This is generally slow and most memory-inefficient. 2) The "coef_wise" option uses a coefficient-wise scaling algorithm with the double variable type. This is generally robust against underflow issues. Computational time varies a lot with conditions; generally only modestly slower than the "double" option, but can be the slowest in some extreme conditions.

For the sake of completeness (only), the scaling parameters \beta (see the package vignette) can be modified via the arguments alphaA and alphaB. These are the factors for the inverses of the largest eigenvalues of \mathbf{A} and \mathbf{B}, respectively, and must be between 0 and 2. The default is 1, which should suffice for most purposes. Values larger than 1 often yield faster convergence, but are not recommended as the error bound will not strictly hold (see Hillier et al. 2009, 2014).

Multithreading

All these functions use C++ versions to speed up computation by default. Furthermore, some of the C++ functions, in particular those using more than one matrix arguments, are parallelized with OpenMP (when available). Use the argument nthreads to control the number of OpenMP threads. By default (nthreads = 0), one-half of the processors detected with omp_get_num_procs() are used. This is except when all the argument matrices share the same eigenvectors and hence the calculation only involves element-wise operations of eigenvalues. In that case, the calculation is typically fast without parallelization, so nthreads is automatically set to 1 unless explicitly specified otherwise; the user can still specify a larger value or 0 for (typically marginal) speed gains in large problems.

Exact method for qfrm_ApIq_int()

An exact expression of the moment is available when p is integer and \mathbf{B} = \mathbf{I}_n (handled by qfrm_ApIq_int()), whose expression involves a confluent hypergeometric function when \bm{\mu} is nonzero (Hillier et al. 2014: theorem 4). There is an option (exact_method = FALSE) to use the ordinary infinite series expression (Hillier et al. 2009), which is less accurate and slow.

Value

A qfrm object consisting of the following:

$statistic

evaluation result (sum(terms))

$terms

vector of 0th to mth order terms

$error_bound

error bound of statistic

$seq_error

vector of error bounds corresponding to partial sums (cumsum(terms))

References

Bao, Y. and Kan, R. (2013) On the moments of ratios of quadratic forms in normal random variables. Journal of Multivariate Analysis, 117, 229–245. doi:10.1016/j.jmva.2013.03.002.

Hillier, G., Kan, R. and Wang, X. (2009) Computationally efficient recursions for top-order invariant polynomials with applications. Econometric Theory, 25, 211–242. doi:10.1017/S0266466608090075.

Hillier, G., Kan, R. and Wang, X. (2014) Generating functions and short recursions, with applications to the moments of quadratic forms in noncentral normal vectors. Econometric Theory, 30, 436–473. doi:10.1017/S0266466613000364.

Smith, M. D. (1989) On the expectation of a ratio of quadratic forms in normal variables. Journal of Multivariate Analysis, 31, 244–257. doi:10.1016/0047-259X(89)90065-1.

Smith, M. D. (1993) Expectations of ratios of quadratic forms in normal variables: evaluating some top-order invariant polynomials. Australian Journal of Statistics, 35, 271–282. doi:10.1111/j.1467-842X.1993.tb01335.x.

See Also

qfmrm for multiple ratio

Examples

## Some symmetric matrices and parameters
nv <- 4
A <- diag(nv:1)
B <- diag(sqrt(1:nv))
mu <- nv:1 / nv
Sigma <- matrix(0.5, nv, nv)
diag(Sigma) <- 1

## Expectation of (x^T A x)^2 / (x^T x)^2 where x ~ N(0, I)
## An exact expression is available
(res1 <- qfrm(A, p = 2))

# The above internally calls the following:
qfrm_ApIq_int(A, p = 2) ## The same

# Similar result with different expression
# This is a suboptimal option and throws a warning
qfrm_ApIq_npi(A, p = 2)

## Expectation of (x^T A x)^1/2 / (x^T x)^1/2 where x ~ N(0, I)
## Note how quickly the series converges in this case
(res2 <- qfrm(A, p = 1/2))
plot(res2)

# The above calls:
qfrm_ApIq_npi(A, p = 0.5)

# This is not allowed (throws an error):
try(qfrm_ApIq_int(A, p = 0.5))

## (x^T A x)^2 / (x^T B x)^3 where x ~ N(0, I)
(res3 <- qfrm(A, B, 2, 3))
plot(res3)

## (x^T A x)^2 / (x^T B x)^2 where x ~ N(mu, I)
## Note the two-sided error bound
(res4 <- qfrm(A, B, 2, 2, mu = mu))
plot(res4)

## (x^T A x)^2 / (x^T B x)^2 where x ~ N(mu, Sigma)
(res5 <- qfrm(A, B, p = 2, q = 2, mu = mu, Sigma = Sigma))
plot(res5)

# Sigma is not allowed in the "internal" functions:
try(qfrm_ApBq_int(A, B, p = 2, q = 2, Sigma = Sigma))

# In res5 above, the error bound didn't converge
# Use larger m to evaluate higher-order terms
plot(print(qfrm(A, B, p = 2, q = 2, mu = mu, Sigma = Sigma, m = 300)))