These are programs for my Journal of Multivariate Analysis paper 
"On the Moments of Ratios of Quadratic Forms in Normal Random 
Variables" (with Yong Bao).  If you have questions, comments, or bug 
reports, please send them to Raymond.kan@rotman.utoronto.ca

Raymond Kan
Rotman School of Management
University of Toronto

Version 1.0: 2/5/2012, initial release
Version 1.1: 8/15/2016, vectorize qratint.m and qratintnc.m to improve
             speed.
Version 1.2: 1/1/2017, improve qratint.m and qratintnc.m when p is an
             integer.  No longer needs to do eigenvalue decomposition.
Version 1.3: 7/22/2026, bug fixes and robustness improvements:
             - qratser.m/qratsernc.m: fix a floating-point bug in a colon
               expression that crashed the non-integer p branch for some
               values of q (e.g., q=1.8).
             - qratser.m/qratsernc.m: when B is rank deficient, the series
               terms decay only polynomially and the maximum number of
               terms can be reached; an estimated tail is now added to
               the partial sum and reported.
             - qratsernc.m: fix a problem reported by Alexey Kaplan: the
               series coefficients alternate in sign and can be as large
               as exp(mu'mu/2), so for large mu'mu the series suffers
               catastrophic cancellation and can return badly wrong
               values without warning (e.g., n=2, A=B=I_2, p=0.5, q=0,
               mu=7*ones(2,1)).  The program now estimates the
               cancellation error from the sum of absolute terms and
               automatically switches to the integration method
               (qratintnc.m, which must be on the path) when the series
               cannot deliver the requested accuracy.
             - qratser.m: with mu = 0_n the above failure mode does not
               arise (the relevant coefficients are nonnegative); the
               only delicate case is an indefinite A with integer p and
               a moment close to zero, where the relative accuracy is
               intrinsically limited and the integration method would do
               no better.  The program therefore only issues a warning
               when cancellation limits the achievable accuracy.
             - qratintnc.m/qratsernc.m: the n=1 case with non-integer
               p-q<0 no longer fails; the confluent hypergeometric
               function is computed by a direct series, which also
               removes the Statistics Toolbox (ncfpdf) dependency.
             - qratint.m/qratintnc.m: the relative tolerance of the
               numerical integration is now set explicitly (1e-9).
             - Zero A or B matrices are now handled gracefully.
             - test3.m added: regression tests for the above.

Integration approach
qratint.m: Compute E[(z'Az)^p/(z'Bz)^q] for z ~ N(0_n,I_n) using an
           integration approach
qratintnc.m: Compute E[(z'Az)^p/(z'Bz)^q] for z ~ N(mu,I_n) using an
             integration approach

Infinite series approach
qratser.m: Compute E[(z'Az)^p/(z'Bz)^q] for z ~ N(0_n,I_n) using an
           infinite series approach
qratsernc.m: Compute E[(z'Az)^p/(z'Bz)^q] for z ~ N(mu,I_n) using an
             infinite series approach
             
Testing programs:
test1.m: test the relative speed of the two different approaches 
         for computing E[(z'Az)^p/(z'Bz)^q] for a set of p and q's.
test2.m: test the relative speed of the two different approaches 
         for computing E[(z'Az)^p/(z'Bz)^q] for different values of n.
test3.m: regression tests for the fixes introduced in Version 1.3.
           