These are programs for my Econometric Theory paper "Computationally Efficient 
Recursions for Top-Order Invariant Polynomials with Applications" (with 
Grant Hillier and Xiaolu Wang)
If you have questions, comments, or bug reports, please send them to 
kan@chass.utoronto.ca

Raymond Kan
Rotman School of Management
University of Toronto

Note 1: The Fortran subdirectory contains the source codes for all the Fortran 
programs.  For cross-platform flexibility (i.e., for both 32 bit and 64 bit Windows,
I use Matlab's preprocessor macro by calling a Fortran header file using the line 

#include "fintrf.h"

at the beginning.  If your Fortran compiler does not support preprocessor macro,
then remove this line and replace mwpointer and mwsize in the declaration by
INTEGER*4 (for 32 bit Windows) or INTEGER*8 (for 64 bit Windows).  The 32 bit
mex files are included in the subdirectory 32 bit and the 64 bit ones are
included in the subdirectory 64 bit.

Note 2: If you receive an error message when running some of the programs,  
you may need to first install Microsoft Visual C++ 2008 SP1 Redistributable
Package (x86) (or x64 if you use 64 bit Windows).  For the 32 bit one, it is 
available at

http://www.microsoft.com/downloads/details.aspx?familyid=A5C84275-3B97-4AB7-A40D-3802B2AF5FC2&displaylang=en

For the 64 bit one, it is available at

http://www.microsoft.com/downloads/details.aspx?familyid=BA9257CA-337F-4B40-8C14-157CFDFFEE4E&displaylang=en

Note 3: The Notes subdirectory contains two notes that discuss the 
        implementation issues for mdks.m and mdkl.m.  Due to the constraints
        on k, we need to use a modification of the algorithm in the paper to 
        compute the coefficients of p and e. 


Version 1.0: 2/15/2008, initial release
Version 1.1: 3/25/2008, bug fix for qratiomoml.m and qratmoml.m, add
             Fortran version of qratmoml.m and qratmoms.m
Version 1.2: 6/17/2008, bug fix for mdkl2.m, mdkl2q.m, mdks2.m, mdks2q.m
Version 1.3: 12/31/2008, bug fix for dkrecurs.m and mhyperpfqevens.m
Version 1.4: 5/15/2009, qmomlf.f is now quad precision, used to be just
             double precision.  hyperg.m is updated so that Matlab's
             own hypergeom.m is used for version R2008b or after. 
             tau2quad.f is fixed to properly deal with matrix A that
             is nonsymmetric.
Version 1.5: 5/22/2009, speed up qratmoml.m, qratiomoml.m, qmoml.m by
             (1) using a matrix to compute p_{i,j}, and (2) only 
             computing the necessary p_{i,j}.  Fix a minor bug in
             qmoml.f and qmoms.f.
Version 1.6: 11/14/2009, change the Fortran programs to work with
             both 32 bit and 64 bit Windows.  In addition, the mex1.m
             now links the BLAS and LAPACK from Matalb, reducing the
             need for external libraries.
Version 1.7: 3/20/2010, change of taum.m and taumc.m so that they work
             also for asymmetric matrices.  Thanks to Yong Bao (of Purdue
             University) for pointing out to me that the previous version 
             does not work properly for asymmetric matrices.  ecoeff2s.m,
             ecoeff.m and ecoeffc.m are also updated to work with asymmetric 
             matrices.
Version 1.8: 1/10/2012, minor change in qratiomoms.m and qratiomoml.m to
             improve efficiency of error checking.
Version 1.9: 3/5/2012, minor update of mdks.m, mdkl.m, mdksquad.f and
             and mdklquad.f to speed up the programs.

             
Single matrix case:
dkrecur.m: Generates d_k(A) for k=0,...,M, where a is the vector
           of eigenvalues of A, and n is the vector of 
           multiplicities
dkr.m: Quadruple precision version of dkrecur.m
dkrecurs.m: Symbolic version of dkrecur.m
dkeven.m: Generates d_k(A) for k=0,...,M when all the eigenvalues
          of A have even multiplicities.
dkevens.m: Symbolic version of dkeven.m
hyperg.m: Matlab hypergeometric function replacement.
          Note that Matlab's hypergeom.m has a bug (up to R2008a) 
          that makes it slow and sometimes return incorrect
          answer.  Example: hypergeom([10 30],20,0.5) = 3.5253e15
          whereas the correct answer is 7.4080e04.  
          If you use Matlab R2008b or later, then you can simply 
          replace hyperg with hypergeom.
mhyperpfqeven.m: {}_1F_1(1/2,alpha;beta;X) when eigenvalues of X 
                 have even multiplicities
mhyperpfqevens.m: Symbolic version of mhyperpfqeven.m


Two matrix case: (A and B are nxn symmetric matrices)
tau2.m/tau2q.m: Computes p_{i,j} by using an efficient algorithm.
                tau2.m also works for A and B that are not
                symmetric.
ecoeff2.m/ecoeff2q.m: Computes the coefficeints of |I_n-t_1*A-t_2*B|
                      using Newton-Girard formula.  These two programs
                      also work for A and B that are not symmetric.
ecoeff2cq.m : Computes the coefficients of |I_n-t_1*A-t_2*diag(B)|
              but with a constraint on the power of t_1 to be 
              less than or equal to p.
ecoeff2s.m: Symbolic version of ecoeff2.m, also works for A and B
            that are not symmetric.
mdks2.m/mdks2q.m: Computes d_{i,j} by using the short recursive
                  algorithm.
mdkl2.m/mdkl2q.m: Computes d_{i,j} by using the long recursive 
                  algorithm.

Multiple matrix case: (A_1 to A_r are symmetric nxn matrices)
ecoeff.m: Generates the coefficients in |I_n-t_1A_1-...-t_rA_r|.
          This program works for asymmetric A_i's.
ecoeffc.m: Generates the coefficients in |I_n-t_1A_1-...-t_rA_r|
           but with the constraint of nu<=kappa.  This program
           works for asymmetric A_i's.
taum.m: Computes p_nu(A_1,...,A_r) by using an efficient algorithm,
        also works for asymmetric A_i's.
taumc.m: Computes p_nu(A_1,...,A_r) by using an efficient algorithm 
         but with the constraint of nu<=kappa, also works for
         asymmetric A_i's.
mdk_kan.m: Computes d_{k_1,...,k_r}(A_1,...,A_r) based on 
           Kan's (2008) identity. 
invpq.m: Computes d_{k_1,...,k_r}(A_1,...,A_r) based on
         Kan's (2008) identity, quadruple precision version
invpoly_new.m: Computes d_{k_1,...,k_r}(A_1,...,A_r) using the
               identity in Kan (2008).  This program generates the
               entire table of d_{k_1,...,k_r}.
mdkl.m/mdklq.m: Computes d_{k_1,...,k_r}(A_1,...,A_r) using the long
                recursive algorithm
mdks.m/mdksq.m: Computes d_{k_1,...,k_r}(A_1,...,A_r) using the short
                recursive algorithm

Moment of ratio of two quadratic forms: 
(A and B are nxn symmetric matrices)
qratmoml.m: Computes E[(z'Az)^p/(z'Bz)^q] using the long recursive 
            algorithm to generate d_{p,j}, where z ~ N(0_n,I_n).
qratmoms.m: Computes E[(z'Az)^p/(z'Bz)^q] using the short recursive 
            algorithm to generate d_{p,j}, where z ~ N(0_n,I_n).
qmoml.m: Fortran version of qratmoml.m, uses quadruple precision.
qmoms.m: Fortran version of qratmoms.m, uses quadruple precision.
qratiomoml.m: Same as qratmoml.m but generates the entire table.
qratiomoms.m: Same as qratmoms.m but generates the entire table.


Testing programs:
test1.m: Tests mhyperpfqeven.m and mhyperpfqevens.m
test2.m: Tests the timing and accuracy of various versions of
         computing d_{i,j}(A,B).  Short recursive relation can
         be numerically unstable when we use double precision.
test3.m: Tests the timing and accuracy of various versions of
         computing d_nu(A_1,...,A_r).
example1.m: Generates Table 1 in the paper.  Note that the result
            based on qratiomoms.m is off a little bit due to
            error propagation in the short recursive algorithm.
example1a.m: Generates Table 1 in the paper based on Fortran
             programs that use quadruple precision.
