A Class of Markovian Self-Reinforcing Processes with Power-Law Distributions

TL;DR

This work introduces a local Markovian self-reinforcing point process in which the instantaneous intensity decays quadratically between events via $\dot{\lambda} = -a\lambda^2$ and is reset after each event by a function $\lambda^+ = f(\lambda^-)$. It shows that constant-reset and linear-reset variants generate power-law tails in inter-event intervals, with $\alpha_{IEI} = 1 + \frac{1}{a}$, and that, under linear resets, the intensity tail obeys a transcendental constraint $a\gamma = 1 - k^{-{\gamma}}$, with $\alpha_{\lambda} = 2 - \gamma$. The framework avoid non-local memory kernels like Hawkes processes and provides analytically tractable expressions for stationary distributions, mean density, and correlations, including fractal-like burst structures and log-correlations between consecutive IEIs. It also offers practical methods for simulation, parameter estimation via maximum likelihood, and model validation, enabling application to bursty phenomena across natural and social systems.

Abstract

Solar flares, email exchanges, and many natural or social systems exhibit bursty dynamics, with periods of intense activity separated by long inactivity. These patterns often follow power- law distributions in inter-event intervals or event rates. Existing models typically capture only one of these features and rely on non-local memory, which complicates analysis and mechanistic interpretation. We introduce a novel self-reinforcing point process whose event rates are governed by local, Markovian nonlinear dynamics and post-event resets. The model generates power-law tails for both inter-event intervals and event rates over a broad range of exponents observed empirically across natural and human phenomena. Compared to non-local models such as Hawkes processes, our approach is mechanistically simpler, highly analytically tractable, and also easier to simulate. We provide methods for model fitting and validation, establishing this framework as a versatile foundation for the study of bursty phenomena.

Figures (10)

• Figure 1: Characteristics of the stochastic process with constant reset, defined by Eq. \ref{['eq:baseline']}. a) Illustration of a simulation run with parameters $a=1.5$ and $\lambda_0 = 1$. b) Event plots for different values of the decay parameter $a$ and total number of simulated events. c) Simulated survival function for the inter-event interval (IEI) distribution. Dashed lines correspond to the theoretical stationary distributions. d) Empirical probability density distribution of the intensity $\lambda$ for decay factor $a = 1.5$ and different simulation durations. The probability density follows the theoretical power-law of $\lambda^{-4/3}$, but since the distribution is non-normalizable, the curve shifts downwards as the simulation length increases. e) Empirical probability density distribution of the pre-event intensity $\lambda^-$ for decay factor $a = 1.5$ and different numbers of simulated events. f) Mean event density ($\rho$) as a function of the decay parameter $a$. The dashed line shows the exact analytical solution; other lines indicate the simulated median $\rho$ for different total numbers of events. Shaded areas represent the 10th–90th percentile range. g) The conditional intensity of observing an event at time $t+\Delta t$ given an event at $t$, $p(\text{event at } t+\Delta t| \text{event at }t)$, as a function of the time lag $\Delta t$. h) The mean duration of a continuous burst of events (as seen in panel Fig 1b) vs. the minimal gap resolution $\tau_{\text{res}}$ used to define it.
• Figure 2: Characteristics of the process with linear reset function. a) Example simulations of the intensity evolution for $a=0.5, c=1$, and three different values of the gain parameter $k$, corresponding to three dynamical regimes. Top: Bounded, self-excited ($k=2/3$). Middle: Unbounded ($k=1.5$). Bottom: Bounded, mixed-feedback ($k=-0.8$). b) An illustration of the event-to-event intensity map, $\lambda^+ = k\lambda^- + c$. The cobweb plot shows how the pre-event ($\lambda^-$) and post-event ($\lambda^+$) intensities evolve over successive events. c) The parameter space of the model, showing the unbounded, bounded self-excited, and bounded mixed-feedback regimes as a function of the gain parameter $k$ d) The survival function of inter-event intervals for $a=1$ and different values of $k$. The identical power-law tail demonstrates that the IEI exponent is independent of $k$ e) The power-law exponent of the stationary intensity distribution in the limit $\lambda \rightarrow \infty$. The red dashed line delineates the stability threshold beyond which the model becomes unstable. f) Empirical survival function of the pre-event intensity $\lambda^-$ for different values of the gain parameter $k$ and $a = 1$. The dashed lines correspond to the theoretical power-law trend for the stationary distribution. g) Empirical survival function of the pre-event intensity for different simulation lengths, $a = 1$ and $k = 2$. h) The mean event density, $\rho$, as a function of $a$ and $k$, calculated from simulations. i) The conditional intensity $h(\Delta t)\equiv p(\textrm{event at } t+\Delta t| \textrm{event at } t)$, for $a=1$ and different values of $k$. j) Covariance of the logarithm for neighboring inter-event intervals as a function of the model parameters $a$ and $k$. k) Covariance of the logarithm for inter-event intervals with $j$ events between them for $a = 1$.
• Figure 3: Behavior of the scale free-process near the critical regime $k=e^a$. a) The evolution of the pre-event intensity, $\lambda^-$, as a function of the event number for three values of $k$: subcritical ($k=2.4$), near-critical ($k=2.7$), and supercritical ($k=3.0$). For $k < e^a$, the process tends to decay, while for $k > e^a$, it tends to grow. b) A schematic illustrating the biased random walk performed by the logarithmic intensity, $\ln(\lambda^-)$. The top panel shows the distribution of the step size, $\Delta(\ln\lambda^-)$, for each of the three regimes. The mean of this distribution (the drift) is negative for the subcritical case (red), positive for the supercritical case (blue), and near zero at the critical threshold (black), explaining the behaviors shown in (a).
• Figure 4: Dynamics and forms of non-linear reset functions. a) Dynamics of a process governed by a general non-linear reset function, $f(\lambda^-)$ (solid black curve), shown in the log-log plane. The intersection of the reset function with the critical line $\lambda^+ = e^a \lambda^-$ (dashed line) defines a reversal point. The process experiences a net drift towards this point, causing the stationary distribution of the pre-event intensity, $p(\lambda^-)$, to concentrate around it (gray histogram). This distribution is well-approximated by the corresponding canonical process (red line, red histogram). b) A family of analytically tractable reset functions, $f(\lambda)=(k\lambda^q+c)^{1/q}$, for several values of the non-linearity parameter $q$. The parameter $k$ is set to $2^q$ for illustration.
• Figure 5: An example of the fitting of a model with $a = 1$ and $k = 1.5$ (indicated by the red dot) for $N_{\text{events}} = 10^3$. The plots show the negative log-likelihood of the model, as a function of the corresponding pair of parameters. The white line indicates $95\%$ MLE confidence estimate.