These are programs for my Econometric Theory paper "Generating Functions 
and Short Recursions, with Applications to the Moments of Quadratic Forms 
in Noncentral Normal Vectors" (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).  

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: You should use the super-short recursion programs for computation.
        The short recursion can be numerically unstable, and the long recursion
        can be slow.

Version 1.0: 7/1/2008, initial release
Version 1.1: 5/24/2009, speed up the long recursive programs by only
             computing the necessary p_{i,j}.  
Version 1.2: 11/15/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.3: 1/10/2012, fix a bug in mdklnc, mdklnc_fast, mdklnc2, mdklnc2_fast,
             mdksnc, mdksnc2 when only one element of kappa is nonzero.
             Fix a bug in qratmomcl, qratmomcl1, qratmomcs, qratmomcs1 when q 
             is not an integer.
Version 2.0: 2/26/2012.  This is a major update.  We introduce a new super-short
             recursion in this version.  
Version 2.1: 1/15/2017, fix a minor bug in mdk.m and mdk2.m for the central case. 


Single matrix case:
dknc.m: Generates \tilde{d}_k(a,delta,n) for k=0,...,M, where a is the vector
        of eigenvalues of A, delta is the vector of noncentrality, and n is the 
        vector of multiplicities.  It uses a super-short recursive algorithm.
dkncs.m: Generates \tilde{d}_k(a,delta,n) for k=0,...,M, where a is the vector
         of eigenvalues of A, delta is the vector of noncentrality, and n is the 
         vector of multiplicities.  It uses a short recursive algorithm.
dkncl.m: Generates \tilde{d}_k(a,delta,n) for k=0,...,M, where a is the vector
         of eigenvalues of A, delta is the vector of noncentrality, and n is the 
         vector of multiplicities.  It uses a long recursive algorithm.

Two matrix case:
mdknc2.m: Computes \tilde{d}_{i,j}(A,B,mu) using the super-short recursive algorithm.
mdkncs2.m: Computes \tilde{d}_{i,j}(A,B,mu) using the short recursive algorithm.
mdkncl2.m: Computes \tilde{d}_{i,j}(A,B,mu) using the long recursive algorithm.
mdkncl2a.m: Computes \tilde{d}_{i,j}(A,B,mu) using a different version of long 
            recursive algorithm.

Multiple matrix case:
mdknc.m: Computes \tilde{d}_{k_1,...,k_r}(A_1,...,A_r,mu) using the super-short 
         recursive algorithm.
mdkncs.m: Computes \tilde{d}_{k_1,...,k_r}(A_1,...,A_r,mu) using the short 
          recursive algorithm.
mdkncl.m: Computes \tilde{d}_{k_1,...,k_r}(A_1,...,A_r,mu) using the long
          recursive algorithm.
mdkncla.m: Computes \tilde{d}_{k_1,...,k_r}(A_1,...,A_r,mu) using a different
           version of the long recursive algorithm.
mdknc_kan.m: Computes \tilde{d}_{k_1,...,k_r}(A_1,...,A_r,mu) based on
             Kan's (2008) identity. 

Moment of ratio of two quadratic forms:
B = I_n:
ratiomom.m: Computes E[(z'Az)^p/(z'z)^q] using the super-short recursive
            algorithm, where z ~ N(mu,I_n).
ratiomoms.m: Computes E[(z'Az)^p/(z'z)^q] using the short recursive
             algorithm, where z ~ N(mu,I_n).
ratiomoml.m: Computes E[(z'Az)^p/(z'z)^q] using the long recursive
             algorithm, where z ~ N(mu,I_n).

General B:
qratmomnc.m: Computes E[(z'Az)^p/(z'Bz)^q] using the super-short recursive 
             algorithm to generate h_{p,j}, where z ~ N(mu,I_n).
qratmomncs.m: Computes E[(z'Az)^p/(z'Bz)^q] using the short recursive 
              algorithm to generate h_{p,j}, where z ~ N(mu,I_n).
qratmomncl.m: Computes E[(z'Az)^p/(z'Bz)^q] using the long recursive 
              algorithm to generate h_{p,j}, where z ~ N(mu,I_n).
qmomnc.m: Fortran version of qratmomnc.m, uses quadruple precision.
qmomncs.m: Fortran version of qratmomncs.m, uses quadruple precision.
qmomncl.m: Fortran version of qratmomncl.m, uses quadruple precision.
qratiomomnc.m: Same as qratmomnc.m but generates the entire table.
qratiomomncl.m: Same as qratmomncl.m but generates the entire table.

Testing programs:
test1.m: Test the speed and accuracy of mdknc2, mdkncs2, mdkncl2, mdkncl2a.
test2.m: Test the speed and accuracy of mdknc, mdkncs, mdkncl, mdkncla,
         and mdknc_kan.
example1.m: Generates Table 2 in the paper using qratiomomnc.m and
            qratiomomncl.m.  Note that the number of steps required 
            based on qratiomomncl.m is off a little bit for p=1, q=10.
eg1.m: Generates Table 2 in the paper using qratmomnc.m, qratmomncs.m 
       and qratmomncl.m. 
       Note that the results based on qratmomncs.m are all wrong due
       to error propagation of the short recursive algorithm.
eg2.m: Generates Table 2 in the paper using qmomnc.m, qmomncs.m and 
       qmomncl.m, which calls external Fortran programs to do the
       calculations in quadrupe precision.
       Note that the number of steps required based on qmomncs.m is
       off a little bit for p=0, q=9 and p=1, q=10.

Subdirectory "central" contain the programs for the special case
of mu=0_n.  All the programs in this subdirectory are based on
super-short recursions.
dk.m: a program to compute normalized top-order zonal polynomials.
mdk2.m: a program to compute normalized top-order invariant 
        polynomials of two matrix.
mdk.m: a program to compute normalized top-order invariant polynomials
       of multiple matrix.
qratmom.m: Computes E[(z'Az)^p/(z'Bz)^q] where z ~ N(0_n,I_n)
qratiomom.m: same as qratmom.m but generates the entire table
qmom.m: Fortran verion of qratmom.m, uses quadruple precision.
