These are programs for my JMVA paper "From Moments of Sum to Moments of Product."
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 zip file 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: In the Fortran programs, I use allocatable arrays.  If you use automatic
        array (like B(N,N)), then the stack will overflow when the size of 
        the automatic array is large.  This will lead to Matlab crashing.

Version 1.0: 2/15/2008, initial release
Version 1.1: 4/30/2009, get rid of the function in qmom.m, qmom_new.m,
             qmomnc.m, and qmomnc_new.m
Version 1.2: 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.3: 6/20/2010, revise qmom_new.m and qmomnc_new.m based on a
             recursive algorithm of Brown (1986).
Version 1.4: 4/19/2015, fix a bug in prodmomd.m and prodmomq.m
             revise prodmom.m and significantly improve prodmom_rec.m
             add prodmomnc.m and prodmomnc_rec.m

prodmom.m: Moment of products of normal random variables z, with z distributed
           as N(0,V).  It makes use of the new identity.
prodmomd.m/prodmomq.m: Same as prodmom.m but calls an external Fortran program,
                       prodmomd.m is double precision, prodmomq.m is quadruple
                       precision.
prodmom_haf.m: Moment of products of normal random variables z, with z 
               distributed as N(0,V).  It uses Isserlis formula.
prodmom_rec.m: Moment of products of normal random variables z, with z 
               distributed as N(0,V).  It uses a recursive formula, and
               it can be faster than prodmom_haf.m when some of the variables
               in the product are repeated.
bivmom.m: Computes the product moment of x_1^{p_1}x_2^{p_2} when x_1 and x_2
          have standard bivariate normal distribution with correlation rho.

prodmomnc.m: Moment of products of normal random variables z, with z distributed
             as N(mu,V).  It makes use of the new identity.
prodmomnc_rec.m: Moment of products of normal random variables z, with z 
                 distributed as N(mu,V).  It uses a recursive formula.

qmom.m: Moment of products of quadratic forms of N(0,V), using the new
        identity
qmom_new.m: Moment of products of quadratic forms of N(0,V), using the new
            identity and also a new recursive algorithm to compute E[(z'Bz)^s].
qmomq.m: Same as qmom.m but calls an external Fortran program.  qmomq.m is
         quadruple precision.

qmomnc.m: Moment of products of quadratic forms of N(mu,V), using the new
          identity
qmomnc_new.m: Moment of products of quadratic forms of N(mu,V), using the new
              identity and also a new recursive algorithm to compute E[(z'Bz)^s].
              This one is often slower than qmomnc.m because it needs to compute
              the eigenvectors in each update whereas qmomnc.m only needs to
              compute the eigenvalues.
qmomncq.m: Same as qmomnc.m but calls an external Fortran program.  qmomncq.m is
           quadruple precision.

test1.m: Compares the speed of various versions of prodmom.m 
test2.m: Compares the speed of prodmomd.m and prodmomq.m
test3.m: Compares the speed of various versions of qmom.m
test4.m: Compares the speed of various versions of qmomnc.m
