These are programs for my paper "Explicit Expression of Expectation 
of Product of Quadratic Forms in Multivariate Normal Random Variables"
(with Jiening Pan).  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: 8/15/2025, initial release
Version 1.1: 9/2/2025, qprods.m, qprodsc.m, qprods_comp.m, and qprodsc_comp.m
             were updated to reduce the number of computations for the 
             coefficients
Version 1.2: 9/.8/2025, added example5.m

Main programs:
qcount.m: a function that computes the number of unique terms in the 
          expression for E[\prod_{i=1}^k z'*A_i*z] when z~N(\mu,\Sigma)
qcountc.m: a function that computes the number of unique terms in the 
           expression for E[\prod_{i=1}^k z'*A_i*z] when z~N(0,\Sigma)
qcountc_rec.m: a function that computes the number of unique terms in the 
               experssion for E[\prod_{i=1}^k z'*A_i*z] when z~N(0,\Sigma).
               This function uses a recursive relation.
qcounts.m: a function that computes the number of unique terms in the 
           expression for E[\prod_{i=1}^k (z'*A_i*z)^{s_i}] when z~N(\mu,\Sigma)
qcountsc.m: a function that computes the number of unique terms in the 
            expression for E[\prod_{i=1}^k (z'*A_i*z)^{s_i}] when z~N(0,\Sigma)

qmom.m: provide the coefficients of the expression for E[(z'*A*z)^k] when z~N(\mu,\Sigma)
qmomc.m: provide the coefficients of the expression for E[(z'*A*z)^k] when z~N(0,\Sigma)
qmom_rec.m: provide the coefficients of the expression for E[(z'*A*z)^k] when z~N(\mu,\Sigma),
            it uses a recursive algorithm
qmomc_rec.m: provide the coefficients of the expression for E[(z'*A*z)^k] when z~N(0,\Sigma),
             it uses a recursive algorithm

qprodastr.m: a function that returns the string representation of the terms
            in the expression for E[\prod_{i=1}^k z'*A_i*z] when z~N(\mu,\Sigma)
qprodastr_comp.m: a function that returns the string representation of the terms
                  in the expression for E[\prod_{i=1}^k z'*A_i*z] when z~CN(\mu,\Sigma)
qprodastrc.m: a function that returns the string representation of the terms
              in the expression for E[\prod_{i=1}^k z'*A_i*z] when z~N(0,\Sigma)
qprodastrc_comp.m: a function that returns the string representation of the terms
                   in the expression for E[\prod_{i=1}^k z'*A_i*z] when z~CN(0,\Sigma)
qprodsstr.m: a function that returns the string representation of the terms
             in the expression for E[\prod_{i=1}^k (z'*A_i*z)^{s_i}] when z~N(\mu,\Sigma)
qprodsstr_comp.m: a function that returns the string representation of the terms
                 in the expression for E[\prod_{i=1}^k (z'*A_i*z)^{s_i}] when z~CN(\mu,\Sigma)
qprodsstrc.m: a function that returns the string representation of the terms
              in the expression for E[\prod_{i=1}^k (z'*A_i*z)^{s_i}] when z~N(0,\Sigma)
qprodsstrc_comp.m: a function that returns the string representation of the terms
                   in the expression for E[\prod_{i=1}^k (z'*A_i*z)^{s_i}] when z~CN(0,\Sigma)

qproda.m: a function that computes the value of E[\prod_{i=1}^k z'*A_i*z] when 
          z~N(\mu,\Sigma) using the efficient explicit expression          
qproda_slow.m: a function that computes the value of E[\prod_{i=1}^k z'*A_i*z] when 
               z~N(\mu,\Sigma) using the less efficient explicit expression          
qproda_lm.m: a function that computes the value of E[\prod_{i=1}^k z'*A_i*z] when 
             z~N(\mu,\Sigma) using the expression from Letac and Massam 
             (2008, JMVA)
qprods.m: a function that computes the value of E[\prod_{i=1}^k (z'*A_i*z)^{s_i}] when 
          z~N(\mu,\Sigma) using the efficient explicit expression

qprodac.m: a function that computes the value of E[\prod_{i=1}^k z'*A_i*z] when 
           z~N(0,\Sigma) using the efficient explicit expression          
qprodac_slow.m: a function that computes the value of E[\prod_{i=1}^k z'*A_i*z] when 
                z~N(0,\Sigma) using the less efficient explicit expression          
qprodac_glm.m: a function that computes the value of E[\prod_{i=1}^k z'*A_i*z] when 
               z~N(0,\Sigma) using the expression from Graczyk, Letac and Massam 
               (2005, JTP)
qprodsc.m: a function that computes the value of E[\prod_{i=1}^k (z'*A_i*z)^{s_i}] when 
           z~N(0,\Sigma) using the efficient explicit expression

qproda_comp.m: a function that computes the value of E[\prod_{i=1}^k z'*A_i*z] when 
               z~CN(\mu,\Sigma) using the efficient explicit expression
qproda_lm_comp.m: a function that computes the value of E[\prod_{i=1}^k z'*A_i*z] when 
                  z~CN(\mu,\Sigma) using the expression from Letac and Massam 
                  (2008, JMVA)
qprods_comp.m: a function that computes the value of E[\prod_{i=1}^k (z'*A_i*z)^{s_i}] when 
               z~CN(\mu,\Sigma) using the efficient explicit expression

qprodac_comp.m: a function that computes the value of E[\prod_{i=1}^k z'*A_i*z] when 
                z~CN(0,\Sigma) using the efficient explicit expression
qprodac_glm_comp.m: a function that computes the value of E[\prod_{i=1}^k z'*A_i*z] when 
                    z~CN(0,\Sigma) using the expression from Graczyk, Letac and Massam 
                    (2003, AS)
qprodsc_comp.m: a function that computes the value of E[\prod_{i=1}^k (z'*A_i*z)^{s_i}] when 
                z~CN(0,\Sigma) using the efficient explicit expression

example1.m: a program that demonstrates the use of qproda.m, qprodac.m, qproda_slow.m,
            qprodac_slow.m, qproda_lm.m, and qprodac_glm.m
example2.m: a program that demonstrates the use of qproda_comp.m, qproda_lm_comp.m,
            qprodac_comp.m, and qprodac_glm_comp.m
example3.m: a program that demonstrates the use of qprods.m and qprodsc.m
example4.m: a program that demonstrates the use of qprods_comp.m and qprodsc_comp.m
example5.m: a program that shows how to use qprodastr.m and qprodsstr_comp.m to
            get the correct expressions for E[tr(W*A1)*tr(W*A2)*tr(W*A3)] for
            real Wishart and E[tr(W*A1)*tr(W*A2)^2] for complex Wishat.
            
Auxiliary programs:
bell.m: compute the Bell number B_k
fixbrace.m: A function that enumerates the bracelets with
            a fixed content
fixbracef.m: A function to enumerate bracelets with fixed
             content of (1,2,...,k).  The output is the
             factoradic of the bracelets
ip_desc.m: A function to generate integer partition in descending 
           lexicographical order
nthprem.m: A function that returns the n-th lexicogrphical permutation of 
           a vector, it was downloaded from Matlab Central FileExchange
bracelet.m: A function that computes the number of bracelets with fixed
            content
necklace.m: A function that computes the number of necklaces with fixed
            content
istheta.m: A function that tests if an index is valid for unique \theta term
totient.m: A function that computes the Euler totient function
divisor.m: A function that returns the list of integer divisors of an integer
gcdv.m: A function that computes the greatest common divisor of a vector of
        integers.
nmp.m: A function that generates the next multiset permutation
setpart.m: A function that generates all set partitions of {1,...,k}
