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.
