<figure><img src="image_1.png" /></figure>

HAL
open science

# Bayesian inference of a spatially dependent semi-Markovian model with application to Madagascar Covid'19 data

Angelo Fulgence Raherinirina, Stefana Tabera Tsilefa, Tsidikaina Nirilanto, Solym Manou-Abi

► To cite this version:

Angelo Fulgence Raherinirina, Stefana Tabera Tsilefa, Tsidikaina Nirilanto, Solym Manou-Abi. Bayesian inference of a spatially dependent semi-Markovian model with application to Madagascar Covid'19 data. PLoS ONE, 2025, 20 (7), pp.e0326264. 10.1371/journal.pone.0326264. hal-05166340

HAL Id: hal-05166340

https://hal.science/hal-05166340v1

Submitted on 16 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Bayesian inference of a spatially dependent semi-Markovian model with application to Madagascar Covid'19 data

Angelo Raherinirina¹,²✉, Stefana Tabera Tsilefa¹✉, Tsidikaina Nirilanto¹, Solym M. Manou-Abi³,⁴

1 Centre de Recherche sur l'Enseignement des Mathématiques, Ecole Normale Supérieure Fianarantsoa Madagascar

2 Laboratoire Informatique et Mathématiques Appliquées pour le Développement, University of Fianarantsoa Madagascar

3 Institut Montpelliérain Alexander Grothendieck, Montpellier, UMR CNRS 5149, France

4 Laboratoire de Mathématiques et Applications, Université de Poitiers, UMR CNRS 7348, France

✉These authors contributed equally to this work.

* angelo.raherinirina@univ-fianarantsoa.mg

## Abstract

This paper presents a Markovian analysis of the spatio-temporal dynamics of an epidemic. We develop an explicit semi-Markovian model that accounts for spatial dependence, operating in discrete time over a finite state space. The model allowed us to have a propagation model conditioned by neighboring states and quantifies two key characteristics : spatial propagation timescales and propagation law in an region dependent on neighbor states. The model is inferred from data collected on the spread of Covid'19 in Madagascar's 22 regions, using a bayesian approach to take into account the knowledge of disease specialists. The result has demonstrated the effect of neighborhoods on the propagation dynamics of the deseases. We conclude with a discussion of potential future theoretical developments.

## Introduction

Following the emergence of Covid'19 and its spread around the world, many countries have implemented multiple interventions, including non-pharmaceutical and pharmaceutical interventions. The Malagasy government has launched awareness, protection and intervention programmes and campaigns. The observed epidemic data studied in this paper relate to the analysis of the spread of the Covid'19 epidemic during the first wave of the epidemic, which began in March 2020 and ended in October 2020, and was reported in all the *R* = 22 regions of the island of Madagascar. Many mathematical models of the Covid'19 coronavirus epidemic have been developed over many years [1,2]. To model phenomena that evolve over time, such as the evolution of a small epidemic, it is often appropriate to represent it by a stochastic process. The Markov model is widely used in epidemiology when working in a discrete state space [3-5] and also in environmental and control analysis for land use dynamics [6-11]. The power and efficiency of the Markovian approach in analysing and modelling complex systems has been known for decades [12, 13]. In its various forms and extensions, the use of Markov models has become intuitive and natural for dealing with
sequence data, including in epidemiology [1, 14], the work of Linda J.S. Allen being one of the best references on this point [15, 16]. The most natural model is that of a classical Markov chain, which imposes that the sojourn times in each state are geometrically or exponential distributed [10, 11, 17]. The adequacy tests performed in this study indicate that this assumption does not always hold for the observed data. A semi-Markov model may be appropriate when the sojourn times laws in different states are not geometric or exponential [7, 18, 19]. The results of preliminary tests on sojourn times allow us to consider other alternatives such as discrete Weibull and Poisson distributions. In addition to suitability tests on the sojourn times, we also perform spatial dependence tests to check the spatial proximity of the observed state transitions [20]. In this study, we will operate on this aspect by integrating a spatial parameter for the state transition according to the neighbouring observations. These models were introduced by [21] in order to propose new strategies for analysing the evolution of the spatial (joint) distributions of revenues in time and space, using this spatialised Markovian extension to measure the spatial dependence of a given dynamic. The semi-Markov hypothesis makes the models more flexible and tractable [22, 23]. In this paper, we try to model the dynamics of Covid'19 from its arrival in Antananarivo (the capital) in March 2020 to its appearance in the other regions of Madagascar. By plotting the intensity rates of the appearance of Covid'19 in each region, we can study the dynamics of its spread in a given region according to the neighbouring regions. In this approach, we propose a classification of the intensity rates into three levels: low, medium and high. Based on the adequacy tests conducted for these three states, the classical Markov assumption appears not to always hold. This motivates the use of semi-Markov models, which provide a natural extension of classical Markov models by allowing greater flexibility in the modeling of sojourn times. We will also address a spatial extension of the semi-Markov model. The remainder of the paper is organized as follows. First, the data-driven model selection process by performing statistical tests on the fit of the data to the sojourn time laws and the spatial dependence of the observations. The results of the statistical tests then lead us to the spatialisation of the semi-Markovian model. Secondly, we perform statistical inference on the model parameters. To do this, we use a Bayesian approach to overcome the data size constraint using expert knowledge. Finally, we discuss the estimation and simulation results of the model and highlight future theoretical and future work.

# 1 Material and Methods

The methods outlined below are based on the daily number of Covid'19 cases during the first wave of the pandemic, spanning March to October 2020, across the 22 regions of Madagascar presented in Figure 1.

## 1.1 Data

The daily number of new cases reported during the above study period is shown in Figure 2.
<figure><img src="image_2.png" /><figcaption>**Fig 1.** Map of the 22 regions of Madagascar</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig 2.** Daily number of cases during the first wave in Madagascar ("https://ourworldindata.org/coronavirus-data").</figcaption></figure>

The spatio-temporal evolution of confirmed cases across these regions is presented in Figure 3. To compute the frequency distribution of the cases for each region, we will also consider the population size as well as the total number of cases cumulated.

57

58

59
<figure><img src="image_4.png" /><figcaption>Fig 3. Daily confirmed cases in Madagascar’s 22 regions</figcaption></figure>

## 1.2 Markov models and spatial dependency

Markov models are the most popular form of analytical framework for statistical models used to describe systems transition from one state to another in a state space(finite). These models are widely used in various fields, such as epidemiology, engineering, decision-making processes and so on. They are simple, well-suited to systems with clear state transitions and computationally efficient, with sound theoretical foundations and well-established analysis methodology. Preliminary results presented in [20] used transition matrices with spatial dependency and neighborhood-conditioned sojourn time distribution to test whether there is a significant difference in spatial lag, determined by the states of neighboring regions. This earlier study highlighted the importance of incorporating spatial dependence to accurately model propagation dynamics, with the difference of sojourn time laws and independence of neighborhood in sojourn times concluding the relevance of a non-Markovian model. On this basis, the present paper adopts semi-Markovian models with spatial dependency in state transition. Maximum likelihood estimation (MLE) is a point estimator used to estimate the appropriate parameters given the observed data. In contrast, the Bayesian method [24] incorporates prior information about the parameters by estimating their distribution. It is widely recognized that the Bayesian approach often provides more accurate results than MLE. The following outlines some key methodological considerations.

For a given region $r \in \{1, \dots, 22\}$, let $n_{rj}$ denote the number of Covid’19 cases recorded during the $j$-th observation week ($j \in \{1, \dots, 21\}$, with 21 weeks of observation in our dataset). The intensity of Covid’19, denoted by $Ic_{rj}$, for region $r$ during week $j$ is defined as (see [20, 25]):

$$Ic_{rj} = \frac{n_{rj}}{\sum_{j=1}^{21} n_{rj}/21}.$$

We then highlight three distinct state profile observations as follows. The choice of the threshold, as described in Table 1, corresponds to the statistical position parameter,
ensuring an approximately homogeneous distribution of states. This approach allows for the classification of regions into comparable categories based on their epidemic intensity.

<table><caption>**Table 1.** The proposed intensity classification of Covid-19</caption><thead><tr><th>Covid-intensity (lc) range</th><th>Covid-intensity level</th><th>Covid-intensity code or state</th></tr></thead><tbody><tr><td>lc ≤ 0.5</td><td>Low</td><td>L (Class 1)</td></tr><tr><td>0.5 &lt; lc ≤ 1.5</td><td>Medium</td><td>M (Class 2)</td></tr><tr><td>lc &gt; 1.5</td><td>High</td><td>H (Class 3)</td></tr></tbody></table>

Consider a stochastic process $(Z_n)_{n \in \mathbb{N}}$ with value in this state space, denoted by E. For a state $i \in E$, $T_i$ defines the sojourn time in state $i$ and is given by $T_i = \inf\{n \in \mathbb{N}; Z_n \neq i | Z_0 = i\}$. After identifying the different types of sojourn time distributions, as detailed in the next section, we assume that the distribution of sojourn times in the LowLow state depends on the nature of the transition originating from that state. More precisely, we assume that for transitions from the LowLow state to the HighHigh state, the sojourn times follow a discrete Weibull distribution, whereas for transitions from the LowLow state to the MediumMedium state, the sojourn times follow a Poisson distribution.

### 1.2.1 Spatial dependency analysis

We are concerned here with spatial homogeneity tests, which consist in identifying some geographical areas of Madagascar according to the intensities of their spatial neighbors. The given $k$ nearest neighbors of the structure in terms of intensity [20, 26] is given in Figure 4 for $k = 3$. The Markov transition matrix can be modified so that the transition probabilities of a given region depend on the intensity class of its spatial neighborhood. Rey [21] developed a theoretical framework that expanded the conventional Markov chain to include the influence of neighbors (spatial dependence) on the progression of the distribution and the probability of transition, incorporating a spatial delay. The first column of Table 2 displays the state of the spatial lag (space neighbor), i.e. the state in which the average Covid'19 intensity per capita of the neighbors of a region belongs to it. The transition probability from the neighboring class $i$ to the neighbouring class $j$, given that the neighbors are in class $v$, is denoted by $p_{ij}/(v)$. This results in a three-dimensional probability hypermatrix $\hat{P}^{(v)}$ corresponding to the transition matrices of the process conditioned by the neighborhood states [21].

<table><caption>**Table 2.** The spatial Markov matrix of Rey [21].</caption><thead><tr><th rowspan="2">Space neighbour</th><th rowspan="2">$t_i$</th><th colspan="3">$t_{i+1}$</th></tr><tr><th>class 1</th><th>class 2</th><th>class 3</th></tr></thead><tbody><tr><td rowspan="3">Class 1</td><td>Class 1</td><td>$p_{11}/1$</td><td>$p_{12}/1$</td><td>$p_{13}/1$</td></tr><tr><td>Class 2</td><td>$p_{21}/1$</td><td>$p_{22}/1$</td><td>$p_{23}/1$</td></tr><tr><td>Class 3</td><td>$p_{31}/1$</td><td>$p_{32}/1$</td><td>$p_{33}/1$</td></tr><tr><td rowspan="3">Class 2</td><td>Class 1</td><td>$p_{11}/2$</td><td>$p_{12}/2$</td><td>$p_{13}/2$</td></tr><tr><td>Class 2</td><td>$p_{21}/2$</td><td>$p_{22}/2$</td><td>$p_{23}/2$</td></tr><tr><td>Class 3</td><td>$p_{31}/2$</td><td>$p_{32}/2$</td><td>$p_{33}/2$</td></tr><tr><td rowspan="3">Class 3</td><td>Class 1</td><td>$p_{11}/3$</td><td>$p_{12}/3$</td><td>$p_{13}/3$</td></tr><tr><td>Class 2</td><td>$p_{21}/3$</td><td>$p_{22}/3$</td><td>$p_{23}/3$</td></tr><tr><td>Class 3</td><td>$p_{31}/3$</td><td>$p_{32}/3$</td><td>$p_{33}/3$</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>Fig 4. Neighborhood identification by nearest neighbor method (Closest neighbor of order k=3)</figcaption></figure>

The spatial homogeneity test is conducted using the previous $k = 3$ nearest neighbours, with spatial autocorrelation indices measuring the spatial dependence between values at different locations in space. The Moran diagram facilitates the interpretation of the spatial structure see [27]. The Local Moran's I statistic quantifies the degree of spatial autocorrelation around each spatial unit within the study area. It serves as a Local Indicator of Spatial Association (LISA), capturing localized clustering patterns of similar or dissimilar values. Formally, it is defined as:

$$I_i = \frac{R}{\sum_{j=1}^{n} (A_j - \bar{A})^2} (A_i - \bar{A}) \sum_{j=1}^{n} w_{ij} (A_j - \bar{A})$$

where:

* R is the total number of spatial units (regions),

* i and j are indices corresponding to two different spatial units ($i, j \in \{1, ..., R\}$),

* $A_i$ and $A_j$ denote the observed values at locations i and j, respectively,

* $\bar{A}$ is the mean of the observed values across all locations,

* We transform the link between points i and j into the element $w_{ij}$ of the weight matrix W.

$$w_{ij} = \begin{cases} 1 & \text{if i and j are considered neighbours} \\ 0 & \text{otherwise} \end{cases}$$

* The weight matrix quantifies the spatial correlation between geographical areas through the neighbourhood matrix denoted W.

The numerator measures the spatial cross-product of deviations from the mean, weighted by spatial proximity, while the denominator normalizes this by the overall variance of the observations. The value of the global Moran's I statistic lies between -1 and 1. A value close to 1 indicates strong positive spatial autocorrelation (similar values cluster together). A value close to -1 indicates strong negative spatial autocorrelation (dissimilar values are neighbours) and a value near 0 suggests spatial randomness.
### 1.2.2 A semi-Markovian model with spatially dependence

We denote by $\mathcal{M}_E$ the set of real matrices on $E \times E$ and by $\mathcal{M}_E(\mathbb{N})$ the set of matrix-valued functions defined on $\mathbb{N}$, with values in $\mathcal{M}_E$. Consider a discrete-time stochastic process $(Z_n)_{n \in \mathbb{N}}$ on the probability space $(\Omega, \mathcal{F}, \mathbb{P})$ with value in $E$ that describes the state of the system through the following tools:

1. The discrete process $J = (J_n)_{n \in \mathbb{N}}$ with state space $E$ measures the n-th jump time.

2. The discrete process $S = (S_n)_{n \in \mathbb{N}}$ represents the successive time points at which state transitions occur in $(Z_n)_{n \in \mathbb{N}}$.

3. The chain $X = (X_n)_{n \in \mathbb{N}}$ with state space $\mathbb{N}^*$, $X_n := S_n - S_{n-1}$ describes the sojourn time in state $J_{n-1}$.

The matrix-valued function $q=(q_{ij}(k)) \in \mathcal{M}_E(\mathbb{N})$ is a discrete-time semi-Markov kernel define by:

$$q_{ij}(k) := \mathbb{P}(J_n = j, S_n - S_{n-1} = k | J_{n-1} = i).$$

It denotes the probability of transitioning from state $i$ to another state $j$ after a stay of $k$ time units in $i$. We refer to [28, 29] for more details on the construction of Semi-Markov chain. The discrete time process $(Z_n)_{n \in \mathbb{N}}$ is the semi-Markovian Chain associated with the Markovian renewal chain $(J_n, S_n)_{n \in \mathbb{N}}$ if:

$$\mathbb{P}(J_n = i_n, S_n - S_{n-1} = X_n | J_{0:n-1} = i_{0:n-1}, X_{0:n-1} = X_{0:n-1}) = \mathbb{P}(J_n = i_n, S_n - S_{n-1} = X_n | J_{n-1} = i_{n-1}),$$

where $X_{0:n-1}$ denotes $(X_0, X_1, \dots, X_{n-1})$. When Equation is independent of $n$, the Markov renewal chain $(J, S)$ is said to be homogeneous. Generally, semi-Markov chains are not Markov processes, as they do not satisfy the Markov property when the sojourn time distribution deviates from a geometric distribution. Furthermore, a semi-Markov process allows arbitrarily distributed sojourn times within any state while preserving the Markov property for the embedded (discrete-time) Markov chain. For more details on the construction of semi-Markov chains, refer to [28, 29]. We denote by $P = (p_{ij/(v)})_{i,j,v \in E}$ the transition matrix of $(J_n)$, incorporating the spatial conditionality as defined in [20, 30]:

$$P_{ij/(v)} := \mathbb{P}(J_{n+1} = j | J_n = i, V_n = v), i, j, v \in E, n \in \mathbb{N}$$

We will consider the following assumption for all $(i, j) \in E \times E$:

1. The conditional distribution $f_{ij}(\cdot)$ of $X_{n+1} = S_{n+1} - S_n, n \in \mathbb{N}$ is independent of the neighborhood state:

$$f_{ij}(k) := \mathbb{P}(X_{n+1} = k | J_n = i, J_{n+1} = j), k \in \mathbb{N}$$

2. The conditional cumulative distribution $F_{ij}(\cdot)$ of $X_{n+1} = S_{n+1} - S_n, n \in \mathbb{N}$ is given by

$$F_{ij}(k) := \mathbb{P}(X_{n+1} \le k | J_n = i, J_{n+1} = j) = \sum_{l=0}^{k} f_{ij}(l), k \in \mathbb{N}.$$

Note that, for all $i, j, v \in E$ and for all $k \in \mathbb{N}$, the following relation holds:

$$f_{ij}(k) = \begin{cases} \frac{q_{ijv}(k)}{P_{ij(v)}}, & n, k \in \mathbb{N} \quad \text{si } P_{ijv} \neq 0 \\ 1_{\{k=\infty\}} & \text{si } P_{ij(v)} = 0 \end{cases} \quad (1)$$
3. The sojourn distribution time, $h_i(\cdot)$, in state $i$ is given by:

$$h_i(k) := \mathbb{P}(T_{n+1} - T_n = k | Y_n = i) = \sum_{j \in E} q_{ij}(k), k \in \mathbb{N}$$

4. The sojourn cumulative distribution time, $H_i(\cdot)$, in state $i$ is given by the following formula:

$$H_i(k) := \mathbb{P}(T_{n+1} - T_n \le k | Y_n = i) = \sum_{l=1}^{k} h_i(l), k \in \mathbb{N}.$$

Note that, from Equation (1), the semi-Markov kernel satisfies the following relation:

$$q_{ij}(k) = P_{ij/(v)} f_{ij}(k)$$

for all $i, j \in E$ and $k \in \mathbb{N}$ such that $p_{ij} \neq 0$. It is clear that the changing state is controlled by the following hypermatrix, according to the dependence of the neighborhood region.

$$P_{ij/(v)} = \begin{pmatrix} 0 & \theta_1^{(v)} & 1 - \theta_1^{(v)} \\ \theta_2^{(v)} & 0 & 1 - \theta_2^{(v)} \\ \theta_3^{(v)} & 1 - \theta_3^{(v)} & 0 \end{pmatrix}, \quad (2)$$

for $i, j \in \{\mathbb{L}, \mathbb{M}, \mathbb{H}\}$ and $v \in \{1, 2, 3\}$. and then, the sojourn time distribution is as follows:

$$f(k) = \begin{pmatrix} 0 & \mathcal{P}(\lambda, k) & \mathcal{W}(\eta, \beta, k) \\ \mathcal{G}(\gamma_1, k) & 0 & \mathcal{G}(\gamma_2, k) \\ \mathcal{G}(\gamma_3, k) & \mathcal{G}(\gamma_4, k) & 0 \end{pmatrix}, k \in \mathbb{N} \quad (3)$$

where $\lambda$ is the parameter of the Poisson distribution, $(\eta, \beta)$ represents the parameters of the discrete Weibull distribution, and $\gamma = (\gamma_1, \gamma_2, \gamma_3, \gamma_4)$ denotes the parameters of the geometric distributions considered.

## 1.3 Bayesian statistical inference for parameter estimation

We consider a Bayesian approach to estimate the parameters described above. The posterior distributions will be computed using the Monte Carlo Markov Chain (MCMC) approximation method, employing algorithms that have proven effective in similar contexts [18]. Specifically, we aim to estimate the elements characterizing the semi-Markov process with spatial effects, including the transition hypermatrix governing state transitions and the conditional sojourn time distributions. To this end, we employ statistical estimation within a Bayesian framework [31, 32]. As previously mentioned, the transition dynamics of a semi-Markov process are fully characterized by the transition matrix of the embedded chain $P$ and the conditional sojourn time distribution $f$ [11]. Let $n_r$ denote the length of the observations associated with region $r$ and $e_k^{(r)}$ the observations of the Covid-19 intensity in region $r$. The neighborhood effect can be represented using a hypermatrix as described above. After recording the sojourn times $(X_k^{(r)})$ for a given region $r$, we can compute the times of state changes and the censored sojourn time for the last state visited, denoted by $u_{n_r} = n_r - T_{n_r}$. Now, denote by $k^{(r)}$ the number of state changes observed for region $r$. We introduce the following notations:

* The value $n_{ij}(k)$ or

$$n_{ij} \stackrel{\text{def}}{=} \sum_{l=0}^{k} 1_{\{J_l=i, J_{l+1}=j\}}$$
is defined as the transition number from i to j (assuming that there is a total of k transitions), where $J_l$ is the included chain that govern the succession states of X.

• The value

$$n_{ij/(v)}^{(r)} \stackrel{\text{def}}{=} \sum_{l=0}^{k^{(r)}} 1_{\{J_l^{(r)}=i, J_{l+1}^{(r)}=j, V_l^{(r)}=v\}}$$

is defined as the transition number from i to j, knowing that the region in the neighbourhood is v, where $J_l^{(r)}$ is the included chain that govern the succession states of X and $V_l^{(r)}$ the states of the process in the neighboring regions of r.

• The quantity

$$n_{i/v}^{(r)} \stackrel{\text{def}}{=} \sum_{j} n_{ij/v}^{(r)}$$

is the number of transitions from i to any other state, knowing that the neighbouring region is v and set $n_{i/v} \stackrel{\text{def}}{=} \sum_{r=1}^{R} n_{i/v}^{(r)}$.

The empirical estimator of the transition matrix $P_{ij}/(v)$ and the conditional sojourn time distribution $f_{ij}$ are defined as follows:

$$\hat{P}_{ij}/(v) = \hat{\theta}^{(v)} = \frac{n_{ij}/(v)}{n_i/(v)}; \quad \hat{f}_{ij}(k) = \frac{n_{ij}(k)}{n_i} \quad n_i = \sum_j n_{ij}.$$

We have

$$\begin{aligned}
\mathbb{P}(J_{0:k^{(r)}}^{(r)} = e_{0:k^{(r)}}^{(r)}, X_{1:k^{(r)}}^{(r)} = x_{1:k^{(r)}}^{(r)}, V_{0:k^{(r)}}^{(r)} = v_{0:k^{(r)}}^{(r)}) &= \mathbb{P}(X_{k^{(r)}+1}^{(r)} \ge u_{n_r} | J_{k^{(r)}}^{(r)} = e_{k^{(r)}}^{(r)}) \\
&= \prod_{r=1}^{R} \mathbb{P}(J_{1:k^{(r)}}^{(r)} = e_{1:k^{(r)}}^{(r)}, X_{1:k^{(r)}}^{(r)} = x_{1:k^{(r)}}^{(r)}, V_{0:k^{(r)}}^{(r)} = v_{0:k^{(r)}}^{(r)}) \\
&= \prod_{r=1}^{R} \prod_{l=1}^{k^{(r)}} \delta_L(e_0^{(r)}) P_{e_{l-1}^{(r)}, e_l^{(r)}, v_l^{(r)}/(r)} \bar{H}_{e_{k^{(r)}}^{(r)}}(u_{n_r})
\end{aligned}$$

where $\delta_L$ is the initial law of $(J_l)_{l \in \mathbb{N}}$ and

$$\bar{H}_i^{(r)}(k) = \mathbb{P}(X_{k^{(r)}+1}^{(r)} > k | J_{k^{(r)}} = i), i \in E, k \in \mathbb{N},$$

corresponds to the survival function associated with the time spent in the last visited state [11]. We assume in the following lines that such survival quantities will be neglected (censored data) [33]. The likelihood function associated with the model is given by

$$L(P, f, r) = \prod_{k=1}^{k^{(r)}} \prod_{v \in E} \prod_{i,j \in E} [P_{ij/v}]^{n_{ij/v}^{(r)}} [f_{ij}(k)]^{n_{ij/v}^{(r)}(k)}. \quad (4)$$

Using the expressions for $P$ and $f$ in the equations (2) and (3), the likelihood can be written as :

$$L(\theta^{(v)}, \lambda, \eta, \beta, \gamma) = \prod_{k=1}^{k^{(r)}} \prod_{v \in E} (\theta_1^{(v)})^{n_{LM/v}} (1 - \theta_1^{(v)})^{n_{LH/v}} (\theta_2^{(v)})^{n_{ML/v}} (1 - \theta_2^{(v)})^{n_{MH/v}} (\theta_3^{(v)})^{n_{HL/v}} (1 - \theta_3^{(v)})^{n_{HM/v}} \\ \left[ (\mathcal{P}_{\lambda}(k))^{n_{LM}(k)} (\mathcal{W}_{(\eta, \beta)}(k))^{n_{LH}(k)} (\mathcal{G}_{\gamma_1}(k))^{n_{ML}(k)} (\mathcal{G}_{\gamma_2}(k))^{n_{MH}(k)} (\mathcal{G}_{\gamma_3}(k))^{n_{HL}(k)} (\mathcal{G}_{\gamma_4}(k))^{n_{HM}(k)} \right] \quad (5)$$

where $\theta^{(v)} = (\theta_1^{(v)}, \theta_2^{(v)}, \theta_3^{(v)})$ is the parameter value of the transition hypermatrix.
### 1.3.1 Bayesian Approach with Informative Prior Distributions for Estimating Sojourn Time Laws

213

214

Now, let us proceed to estimate the parameters described above using a Bayesian approach. We assume that the conditional sojourn time distributions are spatially homogeneous but characterized by different parameters, as follows:

215

216

217

$$f_{ij}(k) = \begin{pmatrix} 0 & \mathcal{P}(\lambda, k) & \mathcal{W}(\eta, \beta, k) \\ \mathcal{G}(\gamma_1, k) & 0 & \mathcal{G}(\gamma_2, k) \\ \mathcal{G}(\gamma_3, k) & \mathcal{G}(\gamma_4, k) & 0 \end{pmatrix}, k \in \mathbb{N}$$

where $\lambda$ represents the parameter of the Poisson distribution, $(\eta, \beta)$ denotes the parameters of the discrete Weibull distribution, and $\gamma = (\gamma_1, \gamma_2, \gamma_3, \gamma_4)$ corresponds to the parameters of the geometric distributions under consideration. In the following, we choose to estimate only the parameters $\lambda$ and $\gamma$ using a Bayesian approach. This decision stems from the fact that estimating the parameters of the discrete Weibull distribution within this framework would incur significant computational costs. For the geometric distribution parameter, we will use a Beta distribution as an informative prior, since it is defined over the interval $[0, 1]$. For the Poisson distribution parameter, we will adopt a Gaussian distribution as the prior. We can then decompose (4) as follows:

$$L(\theta^{(v)}, \lambda, \eta, \beta, \gamma) = L_1(\theta^{(v)}) L_2(\lambda) L_3(\eta, \beta) L_4(\gamma),$$

where

218

$$L_1(\theta^{(v)}) = \prod_{v \in E} (\theta_1^{(v)})^{n_{LM/v}} (1 - \theta_1^{(v)})^{n_{LR/v}} (\theta_2^{(v)})^{n_{ML/v}} (1 - \theta_2^{(v)})^{n_{MH/v}} (\theta_3^{(v)})^{n_{HL/v}} (1 - \theta_3^{(v)})^{n_{HM/v}},$$

$$L_2(\lambda) = \prod_{k=1}^{N} (\mathcal{P}_{\lambda}(k))^{n_{LM}(k)}, \quad L_3(\eta, \beta) = \prod_{k=1}^{N} (\mathcal{W}_{(\eta, \beta)}(k))^{n_{LR}(k)},$$

$$L_4(\gamma) = \prod_{k=1}^{N} (\mathcal{G}_{\gamma_1}(k))^{n_{ML}(k)} (\mathcal{G}_{\gamma_2}(k))^{n_{MH}(k)} (\mathcal{G}_{\gamma_3}(k))^{n_{HL}(k)} (\mathcal{G}_{\gamma_4}(k))^{n_{HM}(k)}.$$

Under the classical framework, the Maximum Likelihood estimators of the parameters are as follows [6]:

219

220

$$\hat{\eta} = 1 - \frac{\sum_{k=1}^{M} k n_{ij}(k)}{\sum_{k=1}^{M} n_{ij}(k)}, \quad \hat{\beta} = \frac{1}{\log(2)} \log \left( \frac{\log(\hat{\eta}) - \frac{\sum_{k=1}^{M} n_{ij}(k+1)}{n_{ij}}}{\log(\hat{\eta})} \right).$$

$$\hat{\gamma} = \frac{\sum_{k=1}^{M} n_{ij}(k)}{\sum_{k=1}^{M} k n_{ij}(k)}, \quad \hat{\lambda} = \frac{\sum_{k=1}^{M} k n_{ij}(k)}{\sum_{k=1}^{M} n_{ij}(k)}.$$

221

Now, based on the Bayesian estimation framework, the distribution of an entire parameter vector $\Sigma$, conditional on the observations, referred to as the posterior distribution $\pi_{posterior}$ is derived using Bayes' theorem:

222

223

224

$$\pi_{posterior} \propto L(\Sigma) \pi_{prior}(\Sigma)$$

where $L(\Sigma)$ is the likelihood function, which is simply the distribution of the observations conditional on $\Sigma$.

225

226
### 1.3.2 Bayesian Approach with Non-Informative Prior Distributions for Estimating the Transition Matrix

This approach leverages Bayesian inference to incorporate uncertainty in disease transmission dynamics, particularly useful when prior knowledge is limited or data is sparse. For the parameters of the transition hypermatrix of the embedded chain, denoted as $\theta^{(v)}$, $v \in \{L, M, H\}$, we adopt the non-informative Jeffreys prior [31]. This choice ensures invariance under reparameterization and avoids introducing subjective bias, which is critical in epidemiological applications where unbiased estimation is paramount [34]. The Jeffreys prior is derived from the Fisher information matrix $\mathbb{I}(\theta^{(v)})$, which quantifies the expected curvature of the log-likelihood function:

$$ \pi_{\text{prior}}^{1}(\theta^{(v)}) \propto \sqrt{\det[\mathbb{I}(\theta^{(v)})]}, \quad \text{where} \quad \mathbb{I}(\theta^{(v)}) = \left[ \mathbb{E}_{\theta^{(v)}} \left( -\frac{\partial^2 l_1(\theta^{(v)})}{\partial \theta_k^{(v)} \partial \theta_l^{(v)}} \right) \right]_{1 \le k, l \le 3}, $$

and $l_1(\theta^{(v)}) = \log L_1(\theta^{(v)})$ is the log-likelihood of the observed transitions. This prior is particularly advantageous for multinomial transition models, as it automatically adapts to the scale of the parameters.

Applying Bayes' theorem, the posterior distribution combines the likelihood $L_1(\theta^{(v)})$ with the Jeffreys prior:

$$ \pi_{\text{posterior}}^{1}(\theta^{(v)}) \propto L_1(\theta^{(v)}) \pi_{\text{prior}}^{1}(\theta^{(v)}). $$

The Bayesian estimator $\hat{\theta}^{(v)}$ is the posterior mean, which minimizes the expected quadratic loss:

$$ \begin{aligned} \hat{\theta}^{(v)} &= \int_{\Theta^{(v)}} \theta^{(v)} \pi_{\text{posterior}}^{1}(\theta^{(v)}) d\theta^{(v)} \\ &= \frac{\int_{\Theta^{(v)}} \theta^{(v)} L_1(\theta^{(v)}) \pi_{\text{prior}}^{1}(\theta^{(v)}) d\theta^{(v)}}{\int_{\Theta^{(v)}} L_1(\theta^{(v)}) \pi_{\text{prior}}^{1}(\theta^{(v)}) d\theta^{(v)}}. \end{aligned} $$

## 2 Results and Discussion

This section begins with a statistical analysis of the Covid'19 epidemic curve on the island of Madagascar using a simple Markov model. The above data seems to reveal that the epidemic exhibited distinct dynamic behaviors across Madagascar's 22 regions. In our study of COVID-19 intensity, each region is classified into one of three distinct states reflecting increasing levels of epidemic severity: Low ($L = 1$), Medium ($M = 2$), and High ($H = 3$). At any given time $t$, the Covid'19 intensity level variable takes one of these three values, corresponding to a particular state, see Figure 5 for the observed weekly intensity sequence.
<figure><img src="image_6.png" /><figcaption>**Fig 5. Weekly intensity sequence tables**</figcaption></figure>

Since the transition between two successive states depends only on the Covid'19 intensity level at the start time, a first-order time-dependent 3-state Markov chain can be constructed. We now turn to justify the choice of our modelling approach based on Markovian testing hypotheses [35].

253

254

255

256
<figure><img src="image_7.png" /><figcaption>Fig 6. Regional trends recorded for the profiles of the above mentioned three states</figcaption></figure>

Note that in the classical Markov assumption it remains to assume that the sojourn-time laws are geometrically distributed [18]. To evaluate such hypotheses, we propose some Statistical test on the sojourn times laws. For instance, one can use the classical Kolmogorov-Smirnov test [36] to check the adequacy of the geometric law for the sojourn times. The geometric distribution on $N^*$ of parameter $\gamma$, $0 < \gamma < 1$ is defined by:

$$G_{\gamma}(k) = \gamma(1 - \gamma)^{k-1}, k \in N^*$$

Figure 7 shows the superposition with a geometric distribution and the results of the goodness-of-fit test in Table 3.
<figure><img src="image_8.png" /><figcaption>**Fig 7.** Superposition of the geometric probability distribution (red colour) on the observed sojourn times (black colour)</figcaption></figure>

<table><caption>**Table 3.** Goodness-of-fit test</caption><thead><tr><td>State</td><td>p-value</td><td>Goodness-of-fit test</td></tr></thead><tbody><tr><td>Low</td><td>0.1233</td><td>hypothesis rejected</td></tr><tr><td>Medium</td><td>0.001</td><td>hypothesis accepted</td></tr><tr><td>High</td><td>0.007</td><td>hypothesis accepted</td></tr></tbody></table>

The above result leads us to examine the following laws for the sejour times laws of the intensity rate. In the case of the state Low, we studied the sejourn time on the states separately according to the type of transition. For the transitions Low → Medium and Low → High, we suggest performing goodness-of-fit tests with other distributions, such as the Poisson distribution and the discrete Weibull distribution as shown in Figure 8. We also present the associated goodness-of-fit tests in Table 4. Recall that the discrete-time Weibull distribution with parameter vector $(\eta, \beta)$ is defined as follows:

$$ \begin{cases} W_{\eta,\beta}(0) = 0 \\ W_{\eta,\beta}(k) = \eta^{(k-1)\beta} - \eta^{k\beta} \quad \forall k \ge 1 \end{cases} $$

where $\beta > 0$ and $\eta \in [0, 1]$.

<figure><img src="image_9.png" /><figcaption>**Fig 8.** Test on other sejour times laws</figcaption></figure>

In the following lines, we'll analyse the spatial effect in the transition probability. Building on the previous statement, it is also insightful to consider the use of non-geometric laws for the sojourn times within the states. Therefore, for the rest we will consider a semi-Markovian model in which the transition probabilities of a region
<table><caption>**Table 4.** Adequacy test result with p-value 0.05</caption><thead><tr><th>State</th><th>Geometric law</th><th>Poisson law</th><th>Discrete Weibull law</th></tr></thead><tbody><tr><td>Low to Medium</td><td>hypothesis rejected</td><td>hypothesis accepted</td><td>hypothesis rejected</td></tr><tr><td>Low to High</td><td>hypothesis rejected</td><td>hypothesis rejected</td><td>hypothesis accepted</td></tr></tbody></table>

depend on its spatial lag. Such spatial dependence will also be analysed. The Moran statistical test, as outlined in Figure 9 for the global spatial autocorrelation and 10 for the Local indicators of spatial association (LISA) has been employed to validate the neighbourhood zones. It demonstrates a p-value of 0.1 when evaluated against the criterion of local spatial autocorrelation.

277

278

279

280

281

<figure><img src="image_10.png" /><figcaption>**Fig 9.** The global indicators of spatial autocorrelation</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig 10.** The Local indicators of spatial association map</figcaption></figure>

We can also assess the spatial effect by comparing the elements in Table 2. The influence of neighboring regions becomes evident, for instance, when $p_{23} > p_{23}/(1)$. This indicates that the probability of regions in class 2 (medium intensity) transitioning

282

283

284
to class 3 (high intensity) is higher without considering their neighbors is higher than when these regions have low-intensity neighbors. Similarly, if $p_{13} > p_{13}/(1)$, the probability of regions transitioning from class 1 (low intensity) to class 2 (medium intensity) is higher when they have neighbors in class 3 (high intensity) compared to the average transition probability. Conversely, if neighbors exert no influence, the transition probabilities would satisfy the condition $p_{ij}/(1) = p_{ij}/(2) = p_{ij}/(3)$, indicating uniformity of the neighboring regions' classes [20,37]. The above data analysis show the following. Firstly, the law of the sojourn times on a given state does not allow a geometric distribution, depending on the nature of the transitions observed. Secondly, the evolution of the intensity of the Covid'19 in a given region depends on the values of the Covid'19 intensity in neighbouring regions.

This means that a semi-Markovian model that includes a spatial effect on the transition matrix could be insightful. Applying log-likelihood maximization of (5) like that [12, 18, 38], we obtain:

$$ \hat{f}_{MLE}(k) = \begin{pmatrix} 0 & \mathcal{P}(\lambda = 5.37, k) & \mathcal{W}(\eta_1 = 0.95, \beta_1 = 5.09, k) \\ \mathcal{G}(\gamma_1 = 0.75, k) & 0 & \mathcal{G}(\gamma_2 = 0.41, k) \\ \mathcal{G}(\gamma_3 = 1, k) & \mathcal{G}(\gamma_4 = 0.31, k) & 0 \end{pmatrix}, k \in \mathbb{N}. $$

The Bayesian estimator $\hat{\Sigma}$ of $\Sigma$ is the mean of the posterior distribution, and it is given by:

$$ \hat{\Sigma} = \int \Sigma \pi_{\text{posterior}}(\Sigma) d\Sigma. $$

We will use the Markov Chain Monte Carlo (MCMC) method [31, 32]. This methods consist in simulating a target probability distribution from an ergodic Markov chain. These techniques were introduced by Metropolis et al [39], they enable numerical sampling according to the a posteriori law. This approach has been previously applied to agrarian dynamics problems [18].

**Algorithm 1:** Metropolis-Hastings algorithm with Gaussian kernel: the target distribution is $\pi_{post}(\theta)$ and the proposal distribution is $g(\cdot, -\theta)$, where g is the density associated with the centered normal distribution $N(0, \sigma^2)$.

Choose Iteration-number;
Sample $\theta \sim U[0, 1]$;
$\tilde{\theta} \leftarrow \theta$;
for $k = 1, 2, \dots, \text{Iteration - number do}$
  $\epsilon \sim N(0, \sigma^2)$;
  $\theta^{prop} \leftarrow \theta + \epsilon$;
  sample $u \sim U[0, 1]$;
  $\alpha \leftarrow \min \left\{ 1, \frac{\pi_{post}(\theta^{prop})g(\theta - \theta^{prop})}{\pi_{post}(\theta)g(\theta^{prop} - \theta)} \right\}$
  if $u \le \alpha$ then
    $\lfloor \theta \leftarrow \theta^{prop} \rfloor$
  $\tilde{\theta} \leftarrow \frac{k-1}{k}\tilde{\theta} + \frac{1}{k}\theta$
<figure><img src="image_12.png" /><figcaption>**Fig 11.** Bayesian estimator of sojourn time parameters: sampling from the posterior distribution using the MCMC method "in black", the associated mean in "red", the maximum likelihood estimator in "blue" and the upper graph is the sampling histogram.</figcaption></figure>

Applying the Bayesian estimation method, the following estimations are obtained: 307

$$ \hat{f}_{\text{Bayesian}}(k) = \begin{pmatrix} 0 & \mathcal{P}(\lambda = 5.02, k) & \mathcal{W}(\eta_1 = 0.95, \beta_1 = 5.09, k) \\ \mathcal{G}(\gamma_1 = 0.654, k) & 0 & \mathcal{G}(\gamma_2 = 0.444, k) \\ \mathcal{G}(\gamma_3 = 0.62, k) & \mathcal{G}(\gamma_4 = 0.33, k) & 0 \end{pmatrix}, k \in \mathbb{N} $$

The count of transitions in the dataset is presented in Table 5.

<table><caption>**Table 5.** Number of transition in the data</caption><thead><tr><th rowspan="2">Space neighbour</th><th rowspan="2">t<sub>i</sub></th><th colspan="3">t<sub>i+1</sub></th></tr><tr><th>class 1</th><th>class 2</th><th>class 3</th></tr></thead><tbody><tr><td rowspan="3">Class 1</td><td>Class 1</td><td>0</td><td>8</td><td>2</td></tr><tr><td>Class 2</td><td>6</td><td>0</td><td>1</td></tr><tr><td>Class 3</td><td>1</td><td>1</td><td>0</td></tr><tr><td rowspan="3">Class 2</td><td>Class 1</td><td>0</td><td>15</td><td>3</td></tr><tr><td>Class 2</td><td>8</td><td>0</td><td>7</td></tr><tr><td>Class 3</td><td>1</td><td>2</td><td>0</td></tr><tr><td rowspan="3">Class 3</td><td>Class 1</td><td>0</td><td>6</td><td>10</td></tr><tr><td>Class 2</td><td>7</td><td>0</td><td>13</td></tr><tr><td>Class 3</td><td>4</td><td>15</td><td>0</td></tr></tbody></table>

We displayed in Figure 12, the Bayesian estimation $\hat{\theta}^{(v)}$ using Markov Chain Monte Carlo (MCMC) method as well as the Maximum Likelihood estimation.
<figure><img src="image_13.png" /><figcaption>Fig 12. Bayesian estimator and maximum likelihood estimator of the parameter $\theta$ of the transition matrix of the spatialized embedded chain: sample of the posterior distribution obtained by MCMC in black, the associated mean in "red", the maximum likelihood estimator in "blue" and the upper graph is the sampling histogram.</figcaption></figure>

Applying the Bayesian approach with Non-Informative Prior Distributions with neighbor dependency, we obtain the following estimation.

<figure><img src="image_14.png" /><figcaption>Fig 13. Transition matrix with neighbor dependency</figcaption></figure>

## 2.1 Limit laws, simulation and time scale of dynamics

We will first examine the law of limits as conditioned by neighbourhood. With the above estimators, the included Markov chains are all irreducible, and the state sojourn times have finite means. Therefore, the spatial semi-Markov chain admits limit laws conditionally on the neighbors.

$$ \Pi_j = \frac{\pi_j \epsilon_j}{\sum_{i \in E} \pi_j \epsilon_i}, j \in E $$
where $\pi$ is the invariant law of the included Markov chain $(Y_k)_{k \in \mathbb{N}}$ and $\epsilon$ the vector containing the mean of the sojourn times on the states [33, 38, 40]. After calculation, we have determined that:

$$ \Pi_{low} = (0.094 \ 0.716 \ 0.19) \quad \Pi_{medium} = (0.057 \ 0.649 \ 0.293) $$

and

$$ \Pi_{high} = (0.126 \ 0.498 \ 0.376) $$

The limit laws represent the equilibrium state of the pandemic, conditional on the neighboring regions.

<figure><img src="image_15.png" /><figcaption>**Fig 14.** Equilibrium state of the pandemic, conditional on the neighborhood states (Low, Medium, High)</figcaption></figure>

From such results, we observe that once the transition laws stabilize, the Covid'19 cases will predominantly remain in the Medium state, regardless of the state of neighboring regions

### 2.1.1 Numerical simulation

Using the numerical simulation method, we can simulate the model over a long trajectory. As the simulated trajectory progresses, the time-scale laws will converge towards a distribution. The simulation method is based on the following algorithm.

**Algorithm 2:** Simulation algorithm for the semi-Markov model with spatial dependence

**Input** : Horizon, initial state $\mu_{01:R}$;

**Output** : $X_{0:Horizon}$;

Initialisation $X_{01:R} \sim \mu_{01:R}$;

for n in Horizon do

$$ \left\{ \begin{array}{l} v \sim f(X_{n1:R}); \\ X_{n+1} \sim q(X_n|v); \\ S \sim f_{ij}(k); \end{array} \right. $$

return $X_{0:Horizon}$

This algorithm simulates semi-Markov dynamics with spatial dependency by generating a state trajectory $X_{0:Horizon}$ over a given horizon for the different regions. Starting from a state i, it identifies the state of influence neighbors v depending on the neighboring states $X_{n1:R}$, then updates the next state j according to a transition law q conditioned by v; and finally it simulates the sojourn time S on i towards j according to a law $f_{ij}(k)$. This approach makes it possible to simulate several model trajectories for numerical study after a few evolution time distributions.
<figure><img src="image_16.png" /><figcaption>**Fig 15.** Example for simulation of the proposed model</figcaption></figure>

Now, we can compute the law of the expected return time for instance to the low state, i.e. the recovery time, and the time to return to the high state, knowing the state of the neighbors based on several simulation samples.

<figure><img src="image_17.png" /><figcaption>**Fig 16.** Law on the expected return time to a low state after contamination</figcaption></figure>

For the recovery time following contamination in a region with "low" intensity neighbors, the expected return time is 5.45 weeks. In contrast, for a region with "Medium" intensity neighbors, the expected recovery time is 10.95 weeks, and for a region with "High" intensity neighbors, the expected recovery time is 19.66 weeks.

<figure><img src="image_18.png" /><figcaption>**Fig 17.** Expected return time to High state after Covid'19 recovery</figcaption></figure>
For the return to epidemic state after recovery in a region with "low" intensity neighbors, the expected time is 62.95 weeks. For a region with "Medium" intensity neighbors, the expected return time is 19.47 weeks, while for a region neighboring "High" intensity areas, the return time is 9.03 weeks. To conclude this section we conducted a simulation of the spread of the epidemic in space.

<figure><img src="image_19.png" /><figcaption>**Fig 18.** The simulation both in time and space. The line represents the kernel estimate of the empirical distribution</figcaption></figure>

The Covid'19 began in Antananarivo, the capital city and from the international airport and then spread across Madagascar. Using simulation methods, we estimated the expected duration for each phase of propagation. It was observed that, after the epidemic reached the capital, the province of Atsinanana would transition to the "High" state most rapidly, with an expected duration of 6.7 weeks. In contrast, the duration in the other provinces showed little variation: Antsiranana reached this state in 7.87 weeks, Mahajanga in 7.53 weeks, Matsiatra Ambony in 7.96 weeks, and Toliary in 7.81 weeks.

# 3 Conclusion and perspectives

This study focuses on spatializing a semi-Markovian dynamic, with particular attention to the intensity conditions of neighboring regions. In our model, the spatial component is integrated by considering the states of neighboring areas, with a specific application to the dynamics of Covid'19 in Madagascar. We began by conducting statistical tests to assess our hypotheses. The semi-Markovian nature of the dynamics was confirmed through a goodness-of-fit test for the geometric distribution, which was rejected for certain state residence times. Regarding the influence of neighboring regions, statistical tests based on the global and local Moran's index confirmed the dependence of several regions on their three nearest neighbors. The parameters of the spatially dependent semi-Markovian model were estimated using both maximum likelihood maximization and a Bayesian approach. These methodologies provided precise estimates, reinforcing the reliability of our model in reflecting the epidemiological reality in Madagascar. By estimating the spatial transition matrix conditioned on neighbors—where the state of a region is determined by the maximum intensity of its three nearest neighbors—and
confirming this spatial dependence through a homogeneity test, we constructed a model characterized by the hypermatrix of the embedded chain and the laws of residence times. The model's analysis led to key insights regarding the spread of Covid'19 in Madagascar, particularly the temporal scales of propagation and spatial diffusion across the country. In conclusion, our study introduces an innovative approach to modeling dynamic systems using a spatially dependent semi-Markovian model. The statistical tests confirmed the validity of our model from both the semi-Markovian and spatial perspectives. The results offer valuable insights into the spread of Covid'19, with both temporal and spatial applications in Madagascar. These findings can guide public health policies and help implement targeted strategies in regions that are highly influenced by their neighbors. Looking forward, we plan to investigate a deep understanding of a kind of spatial markov models including a mathematical analysis of their asymptotic behavior. Also, more realistic distance measures to improve the accuracy of our modeling approach. These future developments could contribute to a deeper understanding of spatial Markov models in decision-making for managing transitional dynamics in various contexts.

## Data Availability Statement

**Availability of data and materials** All code and data-sets accompanying this manuscript are available on GitHub at https://github.com/Jlingz14/Markov_chain_with_spatial_dependency.

**Funding** This research has not received any type of funding

**Conflict of Interests** The authors declare that there is no conflict of interest regarding the publication of this paper.

## Acknowledgments

We certify the submission is an original work and is not under review at any other publication. The corresponding authors thank the University of Mayotte for providing best conditions during his visit that allows the finalisation of this paper.

## References

1. Dietz K, Heesterbeek JAP. Daniel Bernoulli's epidemiological model revisited. 2012;180(1):1-21. doi:10.1016/S0025-5564(02)00122-0.

2. Kermack WO, McKendrick AG, Walker GT. A contribution to the mathematical theory of epidemics. 1927;115(772):700-721. doi:10.1098/rspa.1927.0118.

3. Bouzalmat I, de Saporta B, Manou-Abi SM. Parameter estimation for a hidden linear birth and death process with immigration. arXiv preprint arXiv:230300531v2. 2023;

4. Tsilefa SF, Manou-Abi S, Raherinirina A. Numerical convergence of some stochastic compartmental models in epidemiology; 2022. Available from: https://inria.hal.science/hal-03715840.
5. Raherinirina A, Fandresena TS, Hajalalaina AR, Rabetafika H, Rakotoarivelo RA, Rafamatanantsoa F. Probabilistic Modelling of COVID-19 Dynamic in the Context of Madagascar. 2021;9(3):211-230. doi:10.4236/ojmsi.2021.93014.

6. Campillo F, Hervé D, Raherinirina A, Rakotozafy R. Markov analysis of land use dynamics - A Case Study in Madagascar. Revue Africaine de Recherche en Informatique et Mathématiques Appliquées. 2014;Volume 17 - 2014 - Special issue CARI'12:1-22. doi:10.46298/arima.1964.

7. Kharchenko V, Ponochovnyi Y, Ivanchenko O, Fesenko H, Illiashenko O. Combining Markov and Semi-Markov Modelling for Assessing Availability and Cybersecurity of Cloud and IoT Systems. 2022;6(3):44. doi:10.3390/cryptography6030044.

8. Bickenbach F, Bode E. Evaluating the Markov Property in Studies of Economic Convergence. International Regional Science Review. 2003;26(3):363-392. doi:10.1177/0160017603253789.

9. Claudio D, Moyce S, Albano T, Ibe E, Miller N, O'Leary M. A Markov Chain Model for Mental Health Interventions. 2023;20(4):3525. doi:10.3390/ijerph20043525.

10. Privault N. Discrete-Time Markov Chains. In: Privault N, editor. Understanding Markov Chains: Examples and Applications. Springer Undergraduate Mathematics Series. Springer; 2018. p. 89-113. Available from: https://doi.org/10.1007/978-981-13-0659-4_4.

11. Raherinirina A. Modélisation markovienne des dynamiques d'usage des sols. Cas des parcelles situées sur le bord du corridor forestier Ranomafana-Andringitra [phdthesis]; 2018. Available from: https://tel.archives-ouvertes.fr/tel-00870679.

12. Anderson TW, Goodman LA. Statistical Inference about Markov Chains. The Annals of Mathematical Statistics. 1957;28(1):89-110. doi:10.1214/aoms/1177707039.

13. Yao Q, Wang Y, Xiong X, Wang P, Li Y. Adversarial Decision-Making for Moving Target Defense: A Multi-Agent Markov Game and Reinforcement Learning Approach. Entropy. 2023;25(4). doi:10.3390/e25040605.

14. Gani J, Jerwood D. Markov Chain Methods in Chain Binomial Epidemic Models. Biometrics. 1971;27(3):591-603.

15. Allen LJS. An Introduction to Stochastic Epidemic Models. In: Brauer F, van den Driessche P, Wu J, editors. Mathematical Epidemiology. Lecture Notes in Mathematics. Berlin, Heidelberg: Springer; 2008. p. 81-130. Available from: https://doi.org/10.1007/978-3-540-78911-6_3.

16. Allen EJ, Allen LJS, Arciniega A, Greenwood PE. Construction of Equivalent Stochastic Differential Equation Models. Stochastic Analysis and Applications. 2008;26(2):274-297. doi:10.1080/07362990701857129.

17. Smith WL, Stoneley R. Regenerative stochastic processes. 1997;232(1188):6-31. doi:10.1098/rspa.1955.0198.

18. Raherinirina A, Hervé D, Campillo F. Inferring land use dynamics by semi-Markov model; 2014. p. 1. Available from: https://hal.inria.fr/hal-01100087.
19. Gibbon J. Scalar timing and semi-markov chains in free-operant avoidance. 1971;8(1):109-138. doi:10.1016/0022-2496(71)90025-3.

20. Tabera Tsilefa S, Raherinirina A. Spatial Markov matrices for measuring the spatial dependencies of an epidemiological spread : case Covid'19 Madagascar. BMC Public Health. 2024;24(1):2243. doi:10.1186/s12889-024-19654-9.

21. Rey SJ. Spatial Empirics for Economic Growth and Convergence. Geographical Analysis. 2001;33(3):195-214. doi:10.1111/j.1538-4632.2001.tb00444.x.

22. Janssen J, Limnios N. Semi-Markov Models and Applications; 2011. Available from: https://api.semanticscholar.org/CorpusID:117813955.

23. Vassiliou PCG, Georgiou AC. Markov and Semi-Markov Chains, Processes, Systems, and Emerging Related Fields. Mathematics. 2021;9(19). doi:10.3390/math9192490.

24. Christian R. The Bayesian Choice. Springer Texts in Statistics. New York, NY: Springer; 2007. Available from: http://link.springer.com/10.1007/0-387-71599-1.

25. Wang Y, Zhao C, Liu Z, Gao D. Spatiotemporal Analysis of AIDS Incidence and Its Influencing Factors on the Chinese Mainland, 2005-2017. International Journal of Environmental Research and Public Health. 2021;18(3):1043. doi:10.3390/ijerph18031043.

26. Can A. The Measurement of Neighborhood Dynamics in Urban House Prices. Economic Geography. 1990;66(3):254-272. doi:10.2307/143400.

27. Getis A. Spatial autocorrelation. In: Handbook of applied spatial analysis: Software tools, methods and applications. Springer; 2009. p. 255-278.

28. Limnios N, Oprisan G. Ch. 14. An introduction to semi-markov processes with application to reliability. In: Shanbhag DN, Rao CR, editors. Stochastic Processes: Modelling and Simulation. vol. 21. Elsevier; 2003. p. 515-556.

29. Barbu V, Limnios N. Semi-Markov Chains. In: Semi-Markov Chains and Hidden Semi-Markov Models toward Applications: Their use in Reliability and DNA Analysis. New York, NY: Springer; 2008. p. 1-32. Available from: https://doi.org/10.1007/978-0-387-73173-5_3.

30. Silver D, Silva TH. A Markov model of urban evolution: Neighbourhood change as a complex process. PLOS ONE. 2021;16(1):e0245357. doi:10.1371/journal.pone.0245357.

31. Marin JM, Robert C. Bayesian Core: A Practical Approach to Computational Bayesian Statistics. Springer Science & Business Media; 2007.

32. Alyousifi Y, Ibrahim K, Kang W, Zin WZW. Robust empirical Bayes approach for Markov chain modeling of air pollution index. 2021;19(1):343-356. doi:10.1007/s40201-020-00607-4.

33. Barbu VS, Limnios N. Reliability of Semi-Markov Systems in Discrete Time: Modeling and Estimation. In: Misra KB, editor. Handbook of Performability Engineering. Springer; 2008. p. 369-380. Available from: https://doi.org/10.1007/978-1-84800-131-2_24.
34. Kass RE, Wasserman L. The Selection of Prior Distributions by Formal Rules. Journal of the American Statistical Association. 1996;91(435):1343-1370. doi:10.1080/01621459.1996.10477003.

35. SOLOW AR, SMITH WK. Using Markov chain successional models backwards. Journal of Applied Ecology. 2006;43(1):185-188. doi:https://doi.org/10.1111/j.1365-2664.2005.01127.x.

36. Massey Jr FJ. The Kolmogorov-Smirnov test for goodness of fit. Journal of the American statistical Association. 1951;46(253):68-78.

37. Shabani ZD, Shahnazi R. Spatial distribution dynamics and prediction of COVID-19 in Asian countries : spatial Markov chain approach. Regional Science Policy & Practice. 2020;doi:10.1111/rsp3.12372.

38. Barbu VS, Bulla J, Maruotti A. ESTIMATION OF THE STATIONARY DISTRIBUTION OF A SEMI-MARKOV CHAIN. 2012; p. 15-26.

39. Metropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH, Teller E. Equation of State Calculations by Fast Computing Machines. Journal of Chemical Physics. 1953;21(6):1087-1092. doi:10.1063/1.1699114.

40. Moore AD. The semi-Markov process: a useful tool in the analysis of vegetation dynamics for management. 1990;30(2):111-130. doi:10.1016/0301-4797(90)90009-L.