  Aftershocks of aftershocks - and their aftershock cascades - substantially
contribute to the increased seismicity rate and the associated elevated seismic
hazard after the occurrence of a large earthquake. Current state-of-the-art
earthquake forecasting models therefore describe earthquake occurrence using
self-exciting point processes, where events can recursively trigger more events
according to empirical laws. To estimate earthquake probabilities within future
time horizons of interest, a large number of possible realizations of a process
are simulated, which is typically associated with long computation times that
increase with the desired resolution of the forecast in space, time, or
magnitude range. We here propose a machine learning approach to estimate the
temporal evolution of the rate of higher-order aftershocks. For this, we train
a deep neural network to predict the output of the simulation-based approach,
given a parametric description of the rate of direct aftershocks. A comparison
of the two approaches reveals that they perform very similarly in describing
synthetic datasets generated with the simulation-based approach. Our method has
two major benefits over the traditional approach. It is faster by several
orders of magnitude, and it is not susceptible to being influenced by the
presence or absence of individual `extreme' realizations of the process, and
thus enables accurate earthquake forecasting in near-real-time.
