**Appendix A**

The estimation of the MAP parameter set *θ* = {*<sup>λ</sup>i*, *qij*} is based on the maximization of its likelihood function, *L*(*θ*|*T*) = *<sup>π</sup>Tarre***<sup>D</sup>**0*τ*1**D**1 ...*e***D**0*τN***D**1**1***K*, with interevent times, *τi* = *ti*+<sup>1</sup> − *ti*, *i* = 1, ... , *N*, comprising the trace *T* = {*<sup>τ</sup>*1, ... , *<sup>τ</sup>N*} of the data, *K* hidden states and *N* +1 number of events. The EM algorithm [73] is used for the optimization of the likelihood function, which is a common procedure for applications with hidden data, i.e., the change points of seismicity rate. At each iteration of the algorithm, the log-likelihood (LL) function is computed through the forward and backward vectors, which describe the evolution of the process recursively. The *i*-th element of the forward vector **f**[*k*] = { *fi*(*k*), *i* = 1, ... , *K*} gives the probability to be in state *i* by taking into account the history of occurrences up to time *tk*+1. Their values are obtained recursively through **<sup>f</sup>**[*k*]*j* = ∑*<sup>K</sup> i*=1 **f**[*k* − <sup>1</sup>]*ie*<sup>−</sup>*λiτ<sup>k</sup> qij*(1), with **f**[0] = *πarr*. Similarly, the backward vectors **b**[*k*] = {**b***i*(*k*), *i* = 1, . . . , *K*} are defined giving the likelihood function *L*(*θ*|*T*) = **f**[*k*]**b**[*k* + 1]. Additionally, the forward and backward equations are used for the evaluation of the transitions among the states of the Markov process *Jt*. This is crucial for the implementation of our method, since it allows the detection of changes in the seismicity rate. The state probabilities of the hidden process *Jt* at a given time *tk* are obtained by

$$p\_{l\_k}(i) = P(I\_{l\_k} = i) = \frac{p(\tau\_{l'}, \dots, \tau\_{N'} | I\_{l\_k} = i)}{p(\tau\_{l'}, \dots, \tau\_N)} = \frac{\mathbf{f}[k-1]\_i \mathbf{b}(k | l\_i)}{L(\boldsymbol{\theta} | T)},\tag{A1}$$

For the determination of the hidden states number that is appropriate for capturing the seismicity rate changes, the Bayesian Information Criterion (BIC) [74] is used, which is a metric based on the maximum log-likelihood of each model. It is expressed through *BIC* = −2 × *LL* + *log*(*N*∗) × *k*, where *k* is the number of estimated parameters and *N*∗ corresponds to the number of observations.
