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\).
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
)
A qfrm object consisting of the following:
$statisticevaluation result (sum(terms))
$termsvector of \(0\)th to \(m\)th order terms
$error_bounderror bound of statistic
$seq_errorvector of error bounds corresponding to
partial sums (cumsum(terms))
Argument matrices. Should be square. Will be automatically symmetrized.
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).
Order of polynomials at which the series expression is truncated. \(M\) in Hillier et al. (2009, 2014).
Mean vector \(\bm{\mu}\) for \(\mathbf{x}\)
Covariance matrix \(\mathbf{\Sigma}\) for
\(\mathbf{x}\). Accommodated only by the front-end
qfrm(). See “Details”.
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.
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.
Logical to specify whether the calculation is done with C++
functions via Rcpp. TRUE by default.
Logical to specify whether the exact method is used in
qfrm_ApIq_int() (see “Details”).
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
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”.
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.
Logical to specify whether an error bound is returned (if available).
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 TRUEsame, 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 FALSEskips convergence check
Factors for the scaling constants for \(\mathbf{A}\) and \(\mathbf{B}\), respectively. See “Details”.
Tolerance against which numerical convergence of series is checked. Used
with check_convergence.
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).
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.
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.
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. tools:::Rd_expr_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. tools:::Rd_expr_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. tools:::Rd_expr_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. tools:::Rd_expr_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. tools:::Rd_expr_doi("10.1111/j.1467-842X.1993.tb01335.x").
qfmrm for multiple ratio
## 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)))
Run the code above in your browser using DataLab