  The purpose of this paper is to present a recursive algorithm and its
implementation in Maple and Mathematica for the computation of joint moments
and cumulants of Hawkes processes with exponential kernels. Numerical results
and computation times are also discussed. Obtaining closed form expressions can
be computationally intensive, as joint fifth cumulant and moment formulas can
be respectively expanded into up to 3,288 and 27,116 summands.
