# Mean field analysis of stochastic networks with reservation

Christine Fricker, Hanene Mohamed

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

To cite this version:

Christine Fricker, Hanene Mohamed. Mean field analysis of stochastic networks with reservation.
2023. hal-03539104v2

HAL Id: hal-03539104

https://hal.science/hal-03539104v2

Preprint submitted on 4 Feb 2024

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.

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

Distributed under a Creative Commons Attribution 4.0 International License
MEAN-FIELD ANALYSIS OF STOCHASTIC NETWORKS WITH
RESERVATION

CHRISTINE FRICKER AND HANENE MOHAMED

ABSTRACT. The problem of reservation in a large distributed system is analyzed via a new mathematical model. The target application is car-sharing systems. This model is precisely motivated by the large station-based car-sharing system in France, called Autolib'. This system can be described as a closed stochastic network where the nodes are the stations and the customers are the cars. The user can reserve the car and the parking space. In the paper, we study the evolution of the system when the reservation of parking spaces and cars is effective for all users. The asymptotic behavior of the underlying stochastic network is given when the number N of stations and the fleet size M increase at the same rate. The analysis involves a Markov process on a state space with dimension of order N². It is quite remarkable that the state process describing the evolution of the stations, whose dimension is of order N, converges in distribution, although not Markov, to an non-homogeneous Markov process. We prove this mean-field convergence. We also prove, using combinatorial arguments, that the mean-field limit has a unique equilibrium measure when the time between reserving and picking up the car is sufficiently small. This result extends the case where only the parking space can be reserved.

CONTENTS

<table><tr><td>1. Introduction</td><td>1</td></tr><tr><td>1.1. Simple reservation model</td><td>2</td></tr><tr><td>1.2. Double reservation model</td><td>2</td></tr><tr><td>1.3. Discussion of the model</td><td>2</td></tr><tr><td>1.4. Main results and contribution</td><td>3</td></tr><tr><td>1.5. Related works</td><td>4</td></tr><tr><td>Outline of the paper</td><td>5</td></tr><tr><td>2. The model</td><td>5</td></tr><tr><td>2.1. Transitions of the Markov process</td><td>5</td></tr><tr><td>2.2. Stochastic evolution equations</td><td>6</td></tr><tr><td>3. The asymptotic process</td><td>7</td></tr><tr><td>3.1. Introduction to the asymptotic process</td><td>7</td></tr><tr><td>3.2. Heuristic computation of the asymptotic process</td><td>8</td></tr><tr><td>3.3. A first result</td><td>10</td></tr><tr><td>4. Mean-field limit</td><td>14</td></tr><tr><td>4.1. Evolution equations of the empirical measure</td><td>14</td></tr><tr><td>4.2. The martingale term</td><td>16</td></tr><tr><td>4.3. Convergence of the empirical measure process</td><td>17</td></tr><tr><td>4.4. Probabilistic interpretation of the asymptotic process</td><td>19</td></tr><tr><td>5. Equilibrium of the asymptotic process</td><td>19</td></tr><tr><td>Acknowledgement</td><td>25</td></tr><tr><td>References</td><td>25</td></tr></table>

<footnote n="2020">2020 Mathematics Subject Classification. 60J27, 60K25, 60K30.</footnote>

<footnote>Key words and phrases. Stochastic Networks, Reservation, Mean-Field, Non-Linear Markov Process, Finite Capacity Queues, Product-Form Distribution.</footnote>
# MEAN-FIELD ANALYSIS OF STOCHASTIC NETWORKS WITH RESERVATION

CHRISTINE FRICKER AND HANENE MOHAMED

**ABSTRACT.** The problem of reservation in a large distributed system is analyzed via a new mathematical model. The target application is car-sharing systems. This model is precisely motivated by the large station-based car-sharing system in France, called Autolib'. This system can be described as a closed stochastic network where the nodes are the stations and the customers are the cars. The user can reserve the car and the parking space. In the paper, we study the evolution of the system when the reservation of parking spaces and cars is effective for all users. The asymptotic behavior of the underlying stochastic network is given when the number *N* of stations and the fleet size *M* increase at the same rate. The analysis involves a Markov process on a state space with dimension of order *N*². It is quite remarkable that the state process describing the evolution of the stations, whose dimension is of order *N*, converges in distribution, although not Markov, to an non-homogeneous Markov process. We prove this mean-field convergence. We also prove, using combinatorial arguments, that the mean-field limit has a unique equilibrium measure when the time between reserving and picking up the car is sufficiently small. This result extends the case where only the parking space can be reserved.

## CONTENTS

<table><tr><td>1. Introduction</td><td>1</td></tr><tr><td>1.1. Simple reservation model</td><td>2</td></tr><tr><td>1.2. Double reservation model</td><td>2</td></tr><tr><td>1.3. Discussion of the model</td><td>2</td></tr><tr><td>1.4. Main results and contribution</td><td>3</td></tr><tr><td>1.5. Related works</td><td>4</td></tr><tr><td>Outline of the paper</td><td>5</td></tr><tr><td>2. The model</td><td>5</td></tr><tr><td>2.1. Transitions of the Markov process</td><td>5</td></tr><tr><td>2.2. Stochastic evolution equations</td><td>6</td></tr><tr><td>3. The asymptotic process</td><td>7</td></tr><tr><td>3.1. Introduction to the asymptotic process</td><td>7</td></tr><tr><td>3.2. Heuristic computation of the asymptotic process</td><td>8</td></tr><tr><td>3.3. A first result</td><td>10</td></tr><tr><td>4. Mean-field limit</td><td>14</td></tr><tr><td>4.1. Evolution equations of the empirical measure</td><td>14</td></tr><tr><td>4.2. The martingale term</td><td>16</td></tr><tr><td>4.3. Convergence of the empirical measure process</td><td>17</td></tr><tr><td>4.4. Probabilistic interpretation of the asymptotic process</td><td>19</td></tr><tr><td>5. Equilibrium of the asymptotic process</td><td>19</td></tr><tr><td>Acknowledgement</td><td>25</td></tr><tr><td>References</td><td>25</td></tr></table>

<footnote n="2020">2020 Mathematics Subject Classification. 60J27, 60K25, 60K30.</footnote>

<footnote>Key words and phrases. Stochastic Networks, Reservation, Mean-Field, Non-Linear Markov Process, Finite Capacity Queues, Product-Form Distribution.</footnote>
# 1. INTRODUCTION

The paper deals with a problem of reservation in a large distributed system. Our motivation is a car-sharing system in which a fleet of cars move around and is parked in a set of stations, mainly for electric issues. A crucial problem is the presence of empty and full stations. In an empty station, users cannot pick up a car, while in a full station, also called saturated, users cannot park the car. Reservation could help the user to find both a car at the departure and a parking space at the destination. In this paper, we focus on a reservation policy called *double reservation* which is to reserve both the car and the parking space at the same time, a moment before picking up the car. Cars and parking spaces reservation are proposed by Autolib', the car-sharing system that existed in Paris from 2011 to 2018. Its fleet was composed, in July 2016, of 3980 electric vehicles, called *Bluecars*, distributed in 1084 stations in Paris metropolitan area with 5935 charging points. More than 126900 subscribers had been registered for the service. See [22] for more details about Autolib'.

Note that if the time between the reservation and the pick-up is zero, the policy is to reserve only the parking space when the car is picked up. This policy is called *simple reservation*. In free-floating systems, users cannot reserve the parking spaces as the cars are parked in the public space. These systems are outside the scope of the paper.

## 1.1. Simple reservation model.
The simple reservation model can be described as follows. It consists of $N$ stations of finite capacity $K$ and $M_N$ cars. Users arrive at rate $\lambda$ in each station. A user reserves a parking space at a randomly chosen destination station when he or she picks up a car. If it is not possible, the unhappy user leaves the system. Otherwise, after a trip with exponential distribution of parameter $\mu$, the car is parked at its reserved parking space in the destination station.

The system is said to be large when $N$ and $M_N$ tend to infinity at the same order. Only the fleet size $M_N$ is of the order of $N$. The other variables such as the arrival rate per station $\lambda$ or the mean trip duration $1/\mu$ are bounded, so of the order of 1. Indeed, the increase in the number of $N$ stations can be considered as a densification of the service area, and not as an extension. There is no reason for the mean trip time to increase. The aim of studying the model when $N$ tends to infinity is to obtain an approximation of a car-sharing system with a large fleet and a large number of stations. The aim is not to study the physical extension of the service area or the number of stations.

Let us denote by $s_N$ the ratio $M_N/N$, tending to a constant $s$ which is the *average number of cars per station*. This sizing parameter $s$ is a key parameter of the system. In [6], this policy is studied in a homogeneous framework by a mean-field approach, using a large-scale analysis similar to that of bike-sharing systems in [10].

Due to the the parking space reservation, the state of each station is described in [6] by a vector with two components: the number of cars and the number of reserved parking spaces, which is the main difference with [10]. The state process is $2N$-dimensional and Markovian. With standard arguments, the mean-field convergence is established. Beyond this, the analysis of the equilibrium point is much more delicate with parking space reservation.

The aim of the paper is to prove these results extended to the double reservation model introduced as follows. In particular, our papers gives the proof of [6, Theorem 2] omitted in [6] on the simple reservation model.

## 1.2. Double reservation model.
Station-based car-sharing systems such Autolib' offer the possibility of reserving a car and a parking space in the desired station online, a moment before actually picking up the car. In this paper, we consider a model that meets this user demand. In our model, the user reserves a car at a given departure station and, at the same time, a parking space at a destination station. If there is no car or no parking space available, the user leaves the system. Otherwise, it takes a time called *reservation time* between the car reservation and the car
# 1. INTRODUCTION

The paper deals with a problem of reservation in a large distributed system. Our motivation is a car-sharing system in which a fleet of cars move around and is parked in a set of stations, mainly for electric issues. A crucial problem is the presence of empty and full stations. In an empty station, users cannot pick up a car, while in a full station, also called saturated, users cannot park the car. Reservation could help the user to find both a car at the departure and a parking space at the destination. In this paper, we focus on a reservation policy called *double reservation* which is to reserve both the car and the parking space at the same time, a moment before picking up the car. Cars and parking spaces reservation are proposed by Autolib', the car-sharing system that existed in Paris from 2011 to 2018. Its fleet was composed, in July 2016, of 3980 electric vehicles, called *Bluecars*, distributed in 1084 stations in Paris metropolitan area with 5935 charging points. More than 126900 subscribers had been registered for the service. See [22] for more details about Autolib'.

Note that if the time between the reservation and the pick-up is zero, the policy is to reserve only the parking space when the car is picked up. This policy is called *simple reservation*. In free-floating systems, users cannot reserve the parking spaces as the cars are parked in the public space. These systems are outside the scope of the paper.

## 1.1. Simple reservation model.
The simple reservation model can be described as follows. It consists of $N$ stations of finite capacity $K$ and $M_N$ cars. Users arrive at rate $\lambda$ in each station. A user reserves a parking space at a randomly chosen destination station when he or she picks up a car. If it is not possible, the unhappy user leaves the system. Otherwise, after a trip with exponential distribution of parameter $\mu$, the car is parked at its reserved parking space in the destination station.

The system is said to be large when $N$ and $M_N$ tend to infinity at the same order. Only the fleet size $M_N$ is of the order of $N$. The other variables such as the arrival rate per station $\lambda$ or the mean trip duration $1/\mu$ are bounded, so of the order of 1. Indeed, the increase in the number of $N$ stations can be considered as a densification of the service area, and not as an extension. There is no reason for the mean trip time to increase. The aim of studying the model when $N$ tends to infinity is to obtain an approximation of a car-sharing system with a large fleet and a large number of stations. The aim is not to study the physical extension of the service area or the number of stations.

Let us denote by $s_N$ the ratio $M_N/N$, tending to a constant $s$ which is the *average number of cars per station*. This sizing parameter $s$ is a key parameter of the system. In [6], this policy is studied in a homogeneous framework by a mean-field approach, using a large-scale analysis similar to that of bike-sharing systems in [10].

Due to the the parking space reservation, the state of each station is described in [6] by a vector with two components: the number of cars and the number of reserved parking spaces, which is the main difference with [10]. The state process is $2N$-dimensional and Markovian. With standard arguments, the mean-field convergence is established. Beyond this, the analysis of the equilibrium point is much more delicate with parking space reservation.

The aim of the paper is to prove these results extended to the double reservation model introduced as follows. In particular, our papers gives the proof of [6, Theorem 2] omitted in [6] on the simple reservation model.

## 1.2. Double reservation model.
Station-based car-sharing systems such Autolib' offer the possibility of reserving a car and a parking space in the desired station online, a moment before actually picking up the car. In this paper, we consider a model that meets this user demand. In our model, the user reserves a car at a given departure station and, at the same time, a parking space at a destination station. If there is no car or no parking space available, the user leaves the system. Otherwise, it takes a time called *reservation time* between the car reservation and the car
pick-up. The reservation time has an exponential distribution with mean $1/\nu$. Then follows the trip whose duration is assumed to have an exponential distribution with mean $1/\mu$. Eventually the user parks the car at the reserved space in the destination station and leaves the system.

## 1.3. Discussion of the model.

**Double Vs simple reservation:** For an order of magnitude of the variables of our model, note that Autolib' offered 30 minutes as the maximum reservation duration for the car and 1 hour 30 minutes for the maximum reservation duration of the parking space at the destination. Moreover, the mean trip duration was around 38 minutes (see [19] for 2013). Thus, the mean trip time $1/\mu$ and the mean reservation time $1/\nu$ are comparable. This fully justifies the motivation of studying the double reservation policy.

**Extension to an heterogeneous model:** For sake of simplicity, our choice of a homogeneous network is motivated by the mean-field approach used in our study. A homogeneous framework is generally presented as the simplest, allowing the difficulties of the model to be highlighted and simple explicit results to be obtained. It is still possible to extend this study to a heterogeneous model using clusters by grouping stations with similar parameters in the real system. Since the trip times are random with an exponential distribution of parameter $\mu$, the heterogeneity of the trips is carried by the randomness. Finally, network heterogeneity is carried by the arrival rate depending on the cluster and the probability to reserve a parking space in a station of a given cluster (for the homogeneous model, this probability is $1/N$). See [9] for details. For models with state process having a product-form invariant measure, the heterogeneous framework is quite natural. In the context of bike and car-sharing systems, see [11, 12]. However, our model does not fit into this framework. For such models, the mean-field approach remains effective, hence our choice of a homogeneous model.

**Extension to heterogeneous users:** To take account of different reservation behaviors of users, two classes of users can be considered: users who make reserve early and users who reserve at the last minute. This yields to introduce two different parameters $\nu_1$ and $\nu_2$ for the exponential distribution of the reservation time. This extension remains within a Markovian framework.

Furthermore, the real-world trip time distribution can be fitted by an Erlang distribution. Replacing the exponential distribution of the trip time by an Erlang distribution, i.e. the distribution of the sum of i.i.d. random variables with exponential distribution, a state process capturing the phases of the Erlang distribution is still Markovian.

## 1.4. Main results and contribution.

### 1.4.1. Mean-field approach.
A three-dimensional state space (reserved cars, reserved spaces and available cars) for each station does not allow this process to be fully described as a Markov process, because the parking space reservation time is the sum of two exponentially distributed variables (car pick-up time and trip time). This leads to the introduction of a fourth variable to distinguish between parking spaces reserved by *users not yet travelling* and *users travelling*. But even the associated state process is not Markov. As we will see on the transitions described in Section 2, the Markov process associated to this model, denoted by $(X^N(t))$, is more complicated and especially of dimension $(N+1)^2$, where $N$ is the total number of stations. Therefore, for a given station, the state space is of dimension $N$, which is not suitable for an asymptotic analysis when $N$ tends to infinity. To solve this problem, we introduce the process, denoted by $(Z^N(t))$, which describes the state of each station $i$, $1 \le i \le N$, as a function of the Markov process $(X^N(t))$. Indeed $Z_i^N(t)$ has four components,

- $R_i^{r,N}(t)$ the number of parking places reserved by non-driving users at station $i$ at time $t$,

- $R_i^N(t)$ the number of parking places reserved by users driving at station $i$ at time $t$,

- $V_i^N(t)$ the number of cars available at station $i$ at time $t$ and

- $V_i^{r,N}(t)$ the number of reserved cars at station $i$ at time $t$.
pick-up. The reservation time has an exponential distribution with mean $1/\nu$. Then follows the trip whose duration is assumed to have an exponential distribution with mean $1/\mu$. Eventually the user parks the car at the reserved space in the destination station and leaves the system.

## 1.3. Discussion of the model.

**Double Vs simple reservation:** For an order of magnitude of the variables of our model, note that Autolib' offered 30 minutes as the maximum reservation duration for the car and 1 hour 30 minutes for the maximum reservation duration of the parking space at the destination. Moreover, the mean trip duration was around 38 minutes (see [19] for 2013). Thus, the mean trip time $1/\mu$ and the mean reservation time $1/\nu$ are comparable. This fully justifies the motivation of studying the double reservation policy.

**Extension to an heterogeneous model:** For sake of simplicity, our choice of a homogeneous network is motivated by the mean-field approach used in our study. A homogeneous framework is generally presented as the simplest, allowing the difficulties of the model to be highlighted and simple explicit results to be obtained. It is still possible to extend this study to a heterogeneous model using clusters by grouping stations with similar parameters in the real system. Since the trip times are random with an exponential distribution of parameter $\mu$, the heterogeneity of the trips is carried by the randomness. Finally, network heterogeneity is carried by the arrival rate depending on the cluster and the probability to reserve a parking space in a station of a given cluster (for the homogeneous model, this probability is $1/N$). See [9] for details. For models with state process having a product-form invariant measure, the heterogeneous framework is quite natural. In the context of bike and car-sharing systems, see [11, 12]. However, our model does not fit into this framework. For such models, the mean-field approach remains effective, hence our choice of a homogeneous model.

**Extension to heterogeneous users:** To take account of different reservation behaviors of users, two classes of users can be considered: users who make reserve early and users who reserve at the last minute. This yields to introduce two different parameters $\nu_1$ and $\nu_2$ for the exponential distribution of the reservation time. This extension remains within a Markovian framework.

Furthermore, the real-world trip time distribution can be fitted by an Erlang distribution. Replacing the exponential distribution of the trip time by an Erlang distribution, i.e. the distribution of the sum of i.i.d. random variables with exponential distribution, a state process capturing the phases of the Erlang distribution is still Markovian.

## 1.4. Main results and contribution.

### 1.4.1. Mean-field approach.
A three-dimensional state space (reserved cars, reserved spaces and available cars) for each station does not allow this process to be fully described as a Markov process, because the parking space reservation time is the sum of two exponentially distributed variables (car pick-up time and trip time). This leads to the introduction of a fourth variable to distinguish between parking spaces reserved by *users not yet travelling* and *users travelling*. But even the associated state process is not Markov. As we will see on the transitions described in Section 2, the Markov process associated to this model, denoted by $(X^N(t))$, is more complicated and especially of dimension $(N+1)^2$, where $N$ is the total number of stations. Therefore, for a given station, the state space is of dimension $N$, which is not suitable for an asymptotic analysis when $N$ tends to infinity. To solve this problem, we introduce the process, denoted by $(Z^N(t))$, which describes the state of each station $i$, $1 \le i \le N$, as a function of the Markov process $(X^N(t))$. Indeed $Z_i^N(t)$ has four components,

- $R_i^{r,N}(t)$ the number of parking places reserved by non-driving users at station $i$ at time $t$,

- $R_i^N(t)$ the number of parking places reserved by users driving at station $i$ at time $t$,

- $V_i^N(t)$ the number of cars available at station $i$ at time $t$ and

- $V_i^{r,N}(t)$ the number of reserved cars at station $i$ at time $t$.
As previously mentioned, process $(Z^N(t)) = (R^r, R^N(t), V^N(t), V^r, V^N(t))$ is of dimension $4N$ but is not a Markov process. The loss of the Markov property is the price to pay to this dimension reduction. Nevertheless, it is remarkable and we prove that, despite this non-Markovian description, the state $(Z_i^N(t))$ of a given station $i$ ($1 \le i \le N$) converges in distribution, as $N$ goes to infinity, to a non-homogeneous Markov process $(\bar{Z}(t)) = (\bar{R}^r(t), \bar{R}(t), \bar{V}(t), \bar{V}^r(t))$ satisfying the following Fokker-Planck equation

$$ (1) \qquad \begin{aligned} \frac{d}{dt} \mathbb{E} (f(\bar{Z}(t))) = & \lambda \mathbb{P}(\bar{V}(t) > 0) \mathbb{E} \left( (f(\bar{Z}(t) + e_1) - f(\bar{Z}(t))) 1_{\{\bar{S}(t) < K\}} \right) \\ & + \nu \mathbb{E} \left( (f(\bar{Z}(t) + e_2 - e_1) - f(\bar{Z}(t))) 1_{\{\bar{R}^r(t) > 0\}} \right) \\ & + \mu \mathbb{E} \left( (f(\bar{Z}(t) + e_3 - e_2) - f(\bar{Z}(t))) 1_{\{\bar{R}(t) > 0\}} \right) \\ & + \lambda \mathbb{P}(\bar{S}(t) < K) \mathbb{E} \left( (f(\bar{Z}(t) + e_4 - e_3) - f(\bar{Z}(t))) 1_{\{\bar{V}(t) > 0\}} \right) \\ & + \nu \mathbb{E} \left( (f(\bar{Z}(t) - e_4) - f(\bar{Z}(t))) 1_{\{\bar{V}^r(t) > 0\}} \right) \end{aligned} $$

with $e_1 = (1, 0, 0, 0)$, $e_2 = (0, 1, 0, 0)$, $e_3 = (0, 0, 1, 0)$, $e_4 = (0, 0, 0, 1)$, $f$ a function with finite support on $\mathbb{N}^4$ and $\bar{S}(t) = \bar{R}^r(t) + \bar{R}(t) + \bar{V}(t) + \bar{V}^r(t)$ the limiting number of unavailable parking spaces at station $i$ at time $t$. The asymptotic process $(\bar{Z}(t)) = ((\bar{R}^r(t), \bar{R}(t), \bar{V}(t), \bar{V}^r(t))$ is a jump process with time-dependent rates, the so-called McKean-Vlasov process. This mean-field convergence theorem is not standard at all. Indeed, the mean-field limit is usually obtained for a Markov state process.

### 1.4.2. The equilibrium.
For non-homogeneous Markov processes, recall that there can be several invariant measures. We prove that there is a unique invariant measure in a restricted framework. The proof is based on three main arguments. First, by applying queueing theory, the McKean-Vlasov process on the basic state space is identified with a tandem of four queues with an invariant measure of explicit product-form. Then, by simplifications, the problem of existence and uniqueness of the invariant measure of the non-homogeneous Markov process amounts to the same problem for a fixed point equation in dimension 2. Finally, the global inversion theorem and a monotonicity property allow us to conclude. The last two arguments are based on combinatorial calculations. Monotonicity is just proved when the mean reservation time is sufficiently small. We are convinced that this assumption is technical but the proof needs other tools.

### 1.5. Related works.
For sharing systems, despite the need of analysis of these stochastic systems, most of the literature concerns operation research and data analysis. Probabilistic models have been proposed and analyzed for bike-sharing systems [10, 11], where usage does not allow reservation. As far as we know, there is no stochastic analysis for car-sharing systems with reservation before [6].

For car-sharing systems, a first part of the literature investigates the location problem. In Boyaci et al. [7], OR optimization is used to plan an efficient car-sharing system in terms of the number, location and capacity of stations and fleet size, applied to the case of Nice, France.

Vehicle redistribution and staff rebalancing is a big issue in car-sharing systems (see for example [18]). For simple reservation, called complete parking reservation (CPR), Kaspi et al. show by simulation in [15] and then with OR techniques in [16] that CPR outperforms NR (no reservation) for a specific user-oriented metric. This metric is the excess time users spend in the system due to lack of cars or parking spaces, i.e. the difference between actual time and ideal trip time, including walking times to stations. Paper [8] does not investigate the reservation impact, but the charging problem. A queueing analysis for dimensioning the fleet size is proposed for a closed network taking account car charging.
As previously mentioned, process $(Z^N(t)) = (R^r, R^N(t), V^N(t), V^r, V^N(t))$ is of dimension $4N$ but is not a Markov process. The loss of the Markov property is the price to pay to this dimension reduction. Nevertheless, it is remarkable and we prove that, despite this non-Markovian description, the state $(Z_i^N(t))$ of a given station $i$ ($1 \le i \le N$) converges in distribution, as $N$ goes to infinity, to a non-homogeneous Markov process $(\bar{Z}(t)) = (\bar{R}^r(t), \bar{R}(t), \bar{V}(t), \bar{V}^r(t))$ satisfying the following Fokker-Planck equation

$$ (1) \qquad \begin{aligned} \frac{d}{dt} \mathbb{E} (f(\bar{Z}(t))) = & \lambda \mathbb{P}(\bar{V}(t) > 0) \mathbb{E} \left( (f(\bar{Z}(t) + e_1) - f(\bar{Z}(t))) 1_{\{\bar{S}(t) < K\}} \right) \\ & + \nu \mathbb{E} \left( (f(\bar{Z}(t) + e_2 - e_1) - f(\bar{Z}(t))) 1_{\{\bar{R}^r(t) > 0\}} \right) \\ & + \mu \mathbb{E} \left( (f(\bar{Z}(t) + e_3 - e_2) - f(\bar{Z}(t))) 1_{\{\bar{R}(t) > 0\}} \right) \\ & + \lambda \mathbb{P}(\bar{S}(t) < K) \mathbb{E} \left( (f(\bar{Z}(t) + e_4 - e_3) - f(\bar{Z}(t))) 1_{\{\bar{V}(t) > 0\}} \right) \\ & + \nu \mathbb{E} \left( (f(\bar{Z}(t) - e_4) - f(\bar{Z}(t))) 1_{\{\bar{V}^r(t) > 0\}} \right) \end{aligned} $$

with $e_1 = (1, 0, 0, 0)$, $e_2 = (0, 1, 0, 0)$, $e_3 = (0, 0, 1, 0)$, $e_4 = (0, 0, 0, 1)$, $f$ a function with finite support on $\mathbb{N}^4$ and $\bar{S}(t) = \bar{R}^r(t) + \bar{R}(t) + \bar{V}(t) + \bar{V}^r(t)$ the limiting number of unavailable parking spaces at station $i$ at time $t$. The asymptotic process $(\bar{Z}(t)) = ((\bar{R}^r(t), \bar{R}(t), \bar{V}(t), \bar{V}^r(t))$ is a jump process with time-dependent rates, the so-called McKean-Vlasov process. This mean-field convergence theorem is not standard at all. Indeed, the mean-field limit is usually obtained for a Markov state process.

### 1.4.2. The equilibrium.
For non-homogeneous Markov processes, recall that there can be several invariant measures. We prove that there is a unique invariant measure in a restricted framework. The proof is based on three main arguments. First, by applying queueing theory, the McKean-Vlasov process on the basic state space is identified with a tandem of four queues with an invariant measure of explicit product-form. Then, by simplifications, the problem of existence and uniqueness of the invariant measure of the non-homogeneous Markov process amounts to the same problem for a fixed point equation in dimension 2. Finally, the global inversion theorem and a monotonicity property allow us to conclude. The last two arguments are based on combinatorial calculations. Monotonicity is just proved when the mean reservation time is sufficiently small. We are convinced that this assumption is technical but the proof needs other tools.

### 1.5. Related works.
For sharing systems, despite the need of analysis of these stochastic systems, most of the literature concerns operation research and data analysis. Probabilistic models have been proposed and analyzed for bike-sharing systems [10, 11], where usage does not allow reservation. As far as we know, there is no stochastic analysis for car-sharing systems with reservation before [6].

For car-sharing systems, a first part of the literature investigates the location problem. In Boyaci et al. [7], OR optimization is used to plan an efficient car-sharing system in terms of the number, location and capacity of stations and fleet size, applied to the case of Nice, France.

Vehicle redistribution and staff rebalancing is a big issue in car-sharing systems (see for example [18]). For simple reservation, called complete parking reservation (CPR), Kaspi et al. show by simulation in [15] and then with OR techniques in [16] that CPR outperforms NR (no reservation) for a specific user-oriented metric. This metric is the excess time users spend in the system due to lack of cars or parking spaces, i.e. the difference between actual time and ideal trip time, including walking times to stations. Paper [8] does not investigate the reservation impact, but the charging problem. A queueing analysis for dimensioning the fleet size is proposed for a closed network taking account car charging.
Despite the potential of car sharing, even data analysis remains largely unexplored. This is also due to the lack of data provided by the operators. In [4], Boldrini et al. exploit one month (April 2015) of publicly available data from Autolib', in Paris, France, with 960 stations and 2700 electric cars at this time. We get an idea of the average car pickup rate and the availability of cars at a station. Furthermore, a dichotomy between Paris and the suburbs is highlighted.

Mean-field is an efficient tool for studying the behavior of large distributed systems or interacting systems in many different domains for applications. These systems have both a large number of nodes and customers (particles, cars, etc). To the best of our knowledge, our large-scale stochastic analysis is the first for a stochastic model of car sharing systems with double reservation. Our first main result is that the non-Markovian state process for a given station converges to a non-homogeneous Markov process which is not standard. Furthermore, we note that this limiting Markov process can be described using a simple queueing system. This result was first obtained by simulation in [5, Section 5]. Indeed, in [5], an artificial Markov model on the state of the stations, called the approximate model process, is introduced for the double reservation. Its mean-field limit at equilibrium fits with the real dynamics at equilibrium, which allows the authors to guess such an outcome. The same phenomenon is proved in an entirely different framework, for a model of a network with failures in [1]. Note that the framework in [1] induces simultaneous jumps, which make the proofs more technical. Our paper discusses a simpler framework which focuses on the main arguments and it provides the proof of the result expected in [5]. Like our model, the celebrated Gibbens-Hunt-Kelly model exhibits strong interactions ([13]). In [14], it is proved that these interactions disappear for the mean-field limit. Because of these three models, we believe that this methodology can be useful in many other contexts. A main contribution of the paper relies on the result of uniqueness of the equilibrium point in high dimension (dimension 4), thanks to a nice interpretation in terms of queues. As far as we know, there is no such result in the literature.

**Outline of the paper.** In Section 2, the Markov process ($X^N(t)$) associated to our model is defined and the stochastic evolution equations are given. In Section 3, the second process ($Z^N(t)$) describing the state of the stations is introduced. A heuristic computation of its McKean-Vlasov asymptotic process is derived. Theorem 1 which establishes the existence and uniqueness of this stochastic process is proved. Section 4 is devoted to Theorem 2 giving the mean-field convergence for ($Z^N(t)$) and highlights the probabilistic interpretation of the McKean-Vlasov process. Section 5 analyses of its invariant distribution.

## 2. THE MODEL

In this section, we describe the dynamics of our stochastic model. We recall that the system has $N$ stations of capacity $K$. The Markov process that gives the dynamics of the system is the following,

$$X(t) = (X_{i,j}^N(t), 0 \le i, j \le N)$$

where, for $1 \le i, j \le N$, at time $t$,

- $X_{i,j}^N(t)$ is the number of cars reserved at $i$ with parking space reserved at $j$,

- $X_{0,j}^N(t)$ is the number of parking spaces reserved at $j$ by users driving

- and $X_{i,0}^N(t)$ is the number of cars available at station $i$.

The total number of cars is $M_N$ and the fleet size parameter, defined as the mean number of cars per station if they are all parked, is

$$s \stackrel{\text{def.}}{=} \lim_{N \to +\infty} \frac{M_N}{N}.$$
Despite the potential of car sharing, even data analysis remains largely unexplored. This is also due to the lack of data provided by the operators. In [4], Boldrini et al. exploit one month (April 2015) of publicly available data from Autolib', in Paris, France, with 960 stations and 2 700 electric cars at this time. We get an idea of the average car pickup rate and the availability of cars at a station. Furthermore, a dichotomy between Paris and the suburbs is highlighted.

Mean-field is an efficient tool for studying the behavior of large distributed systems or interacting systems in many different domains for applications. These systems have both a large number of nodes and customers (particles, cars, etc). To the best of our knowledge, our large-scale stochastic analysis is the first for a stochastic model of car sharing systems with double reservation. Our first main result is that the non-Markovian state process for a given station converges to a non-homogeneous Markov process which is not standard. Furthermore, we note that this limiting Markov process can be described using a simple queueing system. This result was first obtained by simulation in [5, Section 5]. Indeed, in [5], an artificial Markov model on the state of the stations, called the approximate model process, is introduced for the double reservation. Its mean-field limit at equilibrium fits with the real dynamics at equilibrium, which allows the authors to guess such an outcome. The same phenomenon is proved in an entirely different framework, for a model of a network with failures in [1]. Note that the framework in [1] induces simultaneous jumps, which make the proofs more technical. Our paper discusses a simpler framework which focuses on the main arguments and it provides the proof of the result expected in [5]. Like our model, the celebrated Gibbens-Hunt-Kelly model exhibits strong interactions ([13]). In [14], it is proved that these interactions disappear for the mean-field limit. Because of these three models, we believe that this methodology can be useful in many other contexts. A main contribution of the paper relies on the result of uniqueness of the equilibrium point in high dimension (dimension 4), thanks to a nice interpretation in terms of queues. As far as we know, there is no such result in the literature.

**Outline of the paper.** In Section 2, the Markov process ($X^N(t)$) associated to our model is defined and the stochastic evolution equations are given. In Section 3, the second process ($Z^N(t)$) describing the state of the stations is introduced. A heuristic computation of its McKean-Vlasov asymptotic process is derived. Theorem 1 which establishes the existence and uniqueness of this stochastic process is proved. Section 4 is devoted to Theorem 2 giving the mean-field convergence for ($Z^N(t)$) and highlights the probabilistic interpretation of the McKean-Vlasov process. Section 5 analyses of its invariant distribution.

## 2. THE MODEL

In this section, we describe the dynamics of our stochastic model. We recall that the system has $N$ stations of capacity $K$. The Markov process that gives the dynamics of the system is the following,

$$X(t) = (X_{i,j}^N(t), 0 \le i, j \le N)$$

where, for $1 \le i, j \le N$, at time $t$,

- $X_{i,j}^N(t)$ is the number of cars reserved at $i$ with parking space reserved at $j$,

- $X_{0,j}^N(t)$ is the number of parking spaces reserved at $j$ by users driving

- and $X_{i,0}^N(t)$ is the number of cars available at station $i$.

The total number of cars is $M_N$ and the fleet size parameter, defined as the mean number of cars per station if they are all parked, is

$$s \stackrel{\text{def.}}{=} \lim_{N \to +\infty} \frac{M_N}{N}.$$
2.1. **Transitions of the Markov process.** The transitions of the Markov process ($X^N(t)$) are reservations, car picks up and car returns.

- Reservations. At station $i$, at rate $\lambda$, an available car is replaced by a reserved car and an available parking space is reserved at the same time at a random station, say $j$. If there is either no available car at $i$ or no available parking space at $j$, the reservation fails. If a reservation is made at time $t$,

$$ \begin{cases} X_{i,0}^N(t) = X_{i,0}^N(t^-) - 1 & \text{for } 1 \le i \le N, \\ X_{i,j}^N(t) = X_{i,j}^N(t^-) + 1 & \text{for } 1 \le i, j \le N. \end{cases} $$

where the limit from the left of function $f$ at $t$ is denoted by $f(t^-)$.

- Car picks up. After a reservation, the user takes a time with an exponential distribution with parameter $\nu$ to come and pick the car. So, at rate $\nu$, each car reserved at station $i$ disappears and the associated parking space reserved at station $j$ moves to a parking space reserved by a driving user. When taking such a car at time $t$,

$$ \begin{cases} X_{i,j}^N(t) = X_{i,j}^N(t^-) - 1, & \text{for } 1 \le i, j \le N, \\ X_{0,j}^N(t) = X_{0,j}^N(t^-) + 1 & \text{for } 1 \le j \le N, \end{cases} $$

- Car returns. After a trip with exponential distribution with parameter $\mu$, a user driving returns his car. Thus, at rate $\mu$, each parking space at station $j$ reserved by a user driving his car becomes an available car. When a car is returned at station $j$ at time $t$,

$$ \begin{cases} X_{0,j}^N(t) = X_{0,j}^N(t^-) - 1, \\ X_{j,0}^N(t) = X_{j,0}^N(t^-) + 1. \end{cases} $$

Note that $(X^N(t))$ is an irreducible Markov process on the finite state space

$$ (2) \qquad \{x = (x_{i,j}) \in \mathbb{N}^{(N+1)^2}, \sum_{j=0}^{N} x_{i,j} + \sum_{j=0}^{N} x_{j,i} \le K, \sum_{0 \le i,j \le N}^{N} x_{i,j} = M_N\}, $$

where $K$ is the finite capacity of each station, thus $(X^N(t))$ is ergodic.

2.2. **Stochastic evolution equations.** The dynamics of $(X^N(t))$ can be given in terms of stochastic integrals with respect to Poisson processes. Let us introduce the following notations.

- A Poisson process on $\mathbb{R}_+$ with parameter $\xi$ is denoted by $N_\xi$. A sequence of such i.i.d. processes is denoted by $(N_{\xi,i}, i \in \mathbb{N})$.

- A Poisson process on $\mathbb{R}_+^2$ with intensity $\xi dt dh$ is denoted by $\bar{N}_\xi$. A sequence of such i.i.d. processes is denoted by $(\bar{N}_{\xi,i}, i \in \mathbb{N})$.

- A marked Poisson process $(t_n, U_n)$, where $(t_n, n \in \mathbb{N})$ is a Poisson process on $\mathbb{R}_+$ with parameter $\xi$ and $(U_n, n \in \mathbb{N})$ is a sequence of i.i.d random variables with uniform distribution on $\{1, \dots, N\}$, is denoted by $N_\xi^{U,N}$. Note that, for $1 \le i \le N$, $N_\xi^{U,N}(\cdot, \{i\})$ is a Poisson process on $\mathbb{R}_+$ with parameter $\xi/N$ and $N_\xi^{U,N}(\cdot, \mathbb{N})$ is a Poisson process on $\mathbb{R}_+$ with parameter $\xi$.

Let us introduce the following independent point processes. A reservation of a car in station $i$ is a point of a Poisson process $N_{\lambda,i}^{U,N}$ on $\mathbb{R}_+^2$. For reservations of both a car at station $i$ and a parking space at station $j$, the times from the moment the user makes a reservation to the moment the car is picked up form a Poisson process $N_{\nu,i,j}$ on $\mathbb{R}_+$. We need a sequence of such i.i.d processes $(N_{\nu,i,j,l}, l \in \mathbb{N})$ as cars are picked up independently. And the same for the trip times of cars returned at station $j$, associated to a sequence $(N_{\mu,j,l}, l \in \mathbb{N})$ of independent Poisson processes.
2.1. **Transitions of the Markov process.** The transitions of the Markov process ($X^N(t)$) are reservations, car picks up and car returns.

- Reservations. At station $i$, at rate $\lambda$, an available car is replaced by a reserved car and an available parking space is reserved at the same time at a random station, say $j$. If there is either no available car at $i$ or no available parking space at $j$, the reservation fails. If a reservation is made at time $t$,

$$ \begin{cases} X_{i,0}^N(t) = X_{i,0}^N(t^-) - 1 & \text{for } 1 \le i \le N, \\ X_{i,j}^N(t) = X_{i,j}^N(t^-) + 1 & \text{for } 1 \le i, j \le N. \end{cases} $$

where the limit from the left of function $f$ at $t$ is denoted by $f(t^-)$.

- Car picks up. After a reservation, the user takes a time with an exponential distribution with parameter $\nu$ to come and pick the car. So, at rate $\nu$, each car reserved at station $i$ disappears and the associated parking space reserved at station $j$ moves to a parking space reserved by a driving user. When taking such a car at time $t$,

$$ \begin{cases} X_{i,j}^N(t) = X_{i,j}^N(t^-) - 1, & \text{for } 1 \le i, j \le N, \\ X_{0,j}^N(t) = X_{0,j}^N(t^-) + 1 & \text{for } 1 \le j \le N, \end{cases} $$

- Car returns. After a trip with exponential distribution with parameter $\mu$, a user driving returns his car. Thus, at rate $\mu$, each parking space at station $j$ reserved by a user driving his car becomes an available car. When a car is returned at station $j$ at time $t$,

$$ \begin{cases} X_{0,j}^N(t) = X_{0,j}^N(t^-) - 1, \\ X_{j,0}^N(t) = X_{j,0}^N(t^-) + 1. \end{cases} $$

Note that $(X^N(t))$ is an irreducible Markov process on the finite state space

$$ (2) \qquad \{x = (x_{i,j}) \in \mathbb{N}^{(N+1)^2}, \sum_{j=0}^{N} x_{i,j} + \sum_{j=0}^{N} x_{j,i} \le K, \sum_{0 \le i,j \le N} x_{i,j} = M_N\}, $$

where $K$ is the finite capacity of each station, thus $(X^N(t))$ is ergodic.

2.2. **Stochastic evolution equations.** The dynamics of $(X^N(t))$ can be given in terms of stochastic integrals with respect to Poisson processes. Let us introduce the following notations.

- A Poisson process on $\mathbb{R}_+$ with parameter $\xi$ is denoted by $N_\xi$. A sequence of such i.i.d. processes is denoted by $(N_{\xi,i}, i \in \mathbb{N})$.

- A Poisson process on $\mathbb{R}_+^2$ with intensity $\xi dt dh$ is denoted by $\bar{N}_\xi$. A sequence of such i.i.d. processes is denoted by $(\bar{N}_{\xi,i}, i \in \mathbb{N})$.

- A marked Poisson process $(t_n, U_n)$, where $(t_n, n \in \mathbb{N})$ is a Poisson process on $\mathbb{R}_+$ with parameter $\xi$ and $(U_n, n \in \mathbb{N})$ is a sequence of i.i.d random variables with uniform distribution on $\{1, \dots, N\}$, is denoted by $N_\xi^{U,N}$. Note that, for $1 \le i \le N$, $N_\xi^{U,N}(\cdot, \{i\})$ is a Poisson process on $\mathbb{R}_+$ with parameter $\xi/N$ and $N_\xi^{U,N}(\cdot, \mathbb{N})$ is a Poisson process on $\mathbb{R}_+$ with parameter $\xi$.

Let us introduce the following independent point processes. A reservation of a car in station $i$ is a point of a Poisson process $N_{\lambda,i}^{U,N}$ on $\mathbb{R}_+^2$. For reservations of both a car at station $i$ and a parking space at station $j$, the times from the moment the user makes a reservation to the moment the car is picked up form a Poisson process $N_{\nu,i,j}$ on $\mathbb{R}_+$. We need a sequence of such i.i.d processes $(N_{\nu,i,j,l}, l \in \mathbb{N})$ as cars are picked up independently. And the same for the trip times of cars returned at station $j$, associated to a sequence $(N_{\mu,j,l}, l \in \mathbb{N})$ of independent Poisson processes.
Using the previous notations, process $(X^N(t))$ is given by the following stochastic differential equations. For $1 \le i, j \le N$ and $t \ge 0$,

$$ (3) \qquad dX_{i,j}^N(t) = - \sum_{l=1}^{\infty} \mathbf{1}_{\{l \le X_{i,j}^N(t^{-})\}} \mathcal{N}_{\nu,i,j,l}(dt) + \mathbf{1}_{\{X_{i,0}^N(t^{-}) > 0, \sum_{k=0}^N (X_{k,j}^N(t^{-}) + X_{j,k}^N(t^{-})) < K\}} \mathcal{N}_{\lambda,i}^{U,N}(dt, \{j\}), $$

$$ (4) \qquad dX_{0,j}^N(t) = \sum_{l=1}^{\infty} \sum_{i=1}^{N} \mathbf{1}_{\{l \le X_{i,j}^N(t^{-})\}} \mathcal{N}_{\nu,i,j,l}(dt) - \sum_{l=1}^{\infty} \mathbf{1}_{\{l \le X_{0,j}^N(t^{-})\}} \mathcal{N}_{\mu,j,l}(dt) $$

and

$$ (5) \qquad dX_{i,0}^N(t) = - \sum_{j=1}^{N} \mathbf{1}_{\{X_{i,0}^N(t^{-}) > 0, \sum_{k=0}^N (X_{k,j}^N(t^{-}) + X_{j,k}^N(t^{-})) < K\}} \mathcal{N}_{\lambda,i}^{U,N}(dt, \{j\}) \\ + \sum_{l=1}^{\infty} \mathbf{1}_{\{l \le X_{0,i}^N(t^{-})\}} \mathcal{N}_{\mu,i,l}(dt). $$

### 3. THE ASYMPTOTIC PROCESS

3.1. **Introduction to the asymptotic process.** Some other notations are needed for that. The state of each station $i$ is given by the quadruplet

$$ Z_i^N(t) = (R_i^{r,N}(t), R_i^N(t), V_i^N(t), V_i^{r,N}(t)) $$

where, at node $i$ at time $t$,

- $R_i^{r,N}(t)$ is the number of parking spaces reserved by non-driving users,

- $R_i^N(t)$ is the number of reserved parking spaces by users driving,

- $V_i^N(t)$ is the number of available cars,

- $V_i^{r,N}(t)$ is the number of reserved cars.

Process $(Z^N(t))$ takes values on

$$ S^N = \{(w_i, x_i, y_i, z_i)_{1 \le i \le N}, (w_i, x_i, y_i, z_i) \in \Sigma_K, \sum_{i=1}^N w_i + x_i + y_i + z_i = M_N\} $$

where

$$ \Sigma_K = \{(w, x, y, z) \in \mathbb{N}^4, w + x + y + z \le K\}. $$

Moreover, let us denote by

$$ S_i^N(t) = R_i^{r,N}(t) + R_i^N(t) + V_i^N(t) + V_i^{r,N}(t) $$

the number of unavailable parking spaces at station $i$ at time $t$. Note that $S_i^N(t)$ is a function of $Z_i^N(t)$. This process $(Z^N(t))$ gives some refined state of the stations and can be obtained as a function of the Markov process $(X^N(t))$ by

$$ R_i^{r,N}(t) = \sum_{j=1}^{N} X_{j,i}^N(t), \qquad R_i^N(t) = X_{0,i}^N(t), $$

$$ V_i^N(t) = X_{i,0}^N(t), \qquad V_i^{r,N}(t) = \sum_{j=1}^{N} X_{i,j}^N(t). $$
Using the previous notations, process $(X^N(t))$ is given by the following stochastic differential equations. For $1 \le i, j \le N$ and $t \ge 0$,

$$ (3) \qquad dX_{i,j}^N(t) = -\sum_{l=1}^{\infty} \mathbf{1}_{\{l \le X_{i,j}^N(t^{-})\}} \mathcal{N}_{\nu,i,j,l}(dt) + \mathbf{1}_{\{X_{i,0}^N(t^{-}) > 0, \sum_{k=0}^N (X_{k,j}^N(t^{-}) + X_{j,k}^N(t^{-})) < K\}} \mathcal{N}_{\lambda,i}^{U,N}(dt, \{j\}), $$

$$ (4) \qquad dX_{0,j}^N(t) = \sum_{l=1}^{\infty} \sum_{i=1}^{N} \mathbf{1}_{\{l \le X_{i,j}^N(t^{-})\}} \mathcal{N}_{\nu,i,j,l}(dt) - \sum_{l=1}^{\infty} \mathbf{1}_{\{l \le X_{0,j}^N(t^{-})\}} \mathcal{N}_{\mu,j,l}(dt) $$

and

$$ (5) \qquad dX_{i,0}^N(t) = -\sum_{j=1}^{N} \mathbf{1}_{\{X_{i,0}^N(t^{-}) > 0, \sum_{k=0}^N (X_{k,j}^N(t^{-}) + X_{j,k}^N(t^{-})) < K\}} \mathcal{N}_{\lambda,i}^{U,N}(dt, \{j\}) \\ \qquad \qquad + \sum_{l=1}^{\infty} \mathbf{1}_{\{l \le X_{0,i}^N(t^{-})\}} \mathcal{N}_{\mu,i,l}(dt). $$

### 3. THE ASYMPTOTIC PROCESS

3.1. **Introduction to the asymptotic process.** Some other notations are needed for that. The state of each station $i$ is given by the quadruplet

$$ Z_i^N(t) = (R_i^{r,N}(t), R_i^N(t), V_i^N(t), V_i^{r,N}(t)) $$

where, at node $i$ at time $t$,

- $R_i^{r,N}(t)$ is the number of parking spaces reserved by non-driving users,

- $R_i^N(t)$ is the number of reserved parking spaces by users driving,

- $V_i^N(t)$ is the number of available cars,

- $V_i^{r,N}(t)$ is the number of reserved cars.

Process $(Z^N(t))$ takes values on

$$ S^N = \{(w_i, x_i, y_i, z_i)_{1 \le i \le N}, (w_i, x_i, y_i, z_i) \in \Sigma_K, \sum_{i=1}^N w_i + x_i + y_i + z_i = M_N\} $$

where

$$ \Sigma_K = \{(w, x, y, z) \in \mathbb{N}^4, w + x + y + z \le K\}. $$

Moreover, let us denote by

$$ S_i^N(t) = R_i^{r,N}(t) + R_i^N(t) + V_i^N(t) + V_i^{r,N}(t) $$

the number of unavailable parking spaces at station $i$ at time $t$. Note that $S_i^N(t)$ is a function of $Z_i^N(t)$. This process $(Z^N(t))$ gives some refined state of the stations and can be obtained as a function of the Markov process $(X^N(t))$ by

$$ R_i^{r,N}(t) = \sum_{j=1}^{N} X_{j,i}^N(t), \qquad R_i^N(t) = X_{0,i}^N(t), $$

$$ V_i^N(t) = X_{i,0}^N(t), \qquad V_i^{r,N}(t) = \sum_{j=1}^{N} X_{i,j}^N(t). $$
So the evolution equations of $(Z^N(t))$ can be obtained from equations (3)–(5) as follows

$$ (6) \qquad dR_i^{r,N}(t) = \sum_{j=1}^{N} \mathbf{1}_{\{V_j^N(t^-)>0, S_i^N(t^-)<K\}} \mathcal{N}_{\lambda,j}^{U,N}(dt, \{i\}) \\ \qquad \qquad - \sum_{j=1}^{N} \left( \sum_{l=1}^{\infty} \mathbf{1}_{\{l \le X_{j,i}^N(t^-)\}} \mathcal{N}_{\nu,j,i,l}(dt) \right), $$

$$ (7) \qquad dR_i^N(t) = \sum_{l=1}^{\infty} \sum_{j=1}^{N} \mathbf{1}_{\{l \le X_{j,i}^N(t^-)\}} \mathcal{N}_{\nu,j,i,l}(dt) - \sum_{l=1}^{\infty} \mathbf{1}_{\{l \le R_i^N(t^-)\}} \mathcal{N}_{\mu,i,l}(dt), $$

$$ (8) \qquad dV_i^N(t) = \sum_{l=1}^{\infty} \mathbf{1}_{\{l \le R_i^N(t^-)\}} \mathcal{N}_{\mu,i,l}(dt) \\ \qquad \qquad - \sum_{j=1}^{N} \mathbf{1}_{\{V_i^N(t^-)>0, S_j^N(t^-)<K\}} \mathcal{N}_{\lambda,i}^{U,N}(dt, \{j\}), $$

$$ (9) \qquad dV_i^{r,N}(t) = \sum_{j=1}^{N} \mathbf{1}_{\{V_j^N(t^-)>0, S_j^N(t^-)<K\}} \mathcal{N}_{\lambda,i}^{U,N}(dt, \{j\}) \\ \qquad \qquad - \sum_{j=1}^{N} \left( \sum_{l=1}^{\infty} \mathbf{1}_{\{l \le X_{i,j}^N(t^-)\}} \mathcal{N}_{\nu,i,j,l}(dt) \right). $$

This process $(Z^N(t))$, in dimension $4N$, is not a Markov process because the evolution equations for $(Z^N(t))$ are not autonomous. They depend on the Markov process $(X^N(t))$ which lives in dimension $(N+1)^2$. See equation (2). We introduce process $(Z^N(t))$ because it is sufficient to capture the performance of the system, such as the fact that a station is empty or full. The quite remarkable property is that the limit as $N$ gets large of $(Z^N(t))$ is a non-linear Markov process. We will present this process in this section and prove the convergence result in the next section.

### 3.2. Heuristic computation of the asymptotic process.
The proof of the convergence will be given in the next section. Here we show how we can guess the asymptotic process. Suppose that, for $1 \le i \le N$, $(Z_i^N(t))$ converges in distribution to some process

$$ (\bar{Z}(t)) = (\bar{R}^r(t), \bar{R}(t), \bar{V}(t), \bar{V}^r(t)). $$

Let $\mathcal{P}_i^N$ be the random measure on $\mathbb{R}_+$ defined by

$$ (10) \qquad \mathcal{P}_i^N([0, t]) = \int_0^t \sum_{j=1}^N \left( \sum_{l=1}^\infty \mathbf{1}_{\{l \le X_{j,i}^N(s^-)\}} \mathcal{N}_{\nu,j,i,l}(ds) \right). $$

It shows that $\mathcal{P}_i^N$ is a counting process (with jump size 1) on $\mathbb{R}_+$ and its compensator is

$$ \nu \int_0^t \sum_{j=1}^N X_{j,i}^N(s) \, ds = \nu \int_0^t R_i^{r,N}(s) \, ds. $$

See Robert [20, Proposition A.9] for example. Thus, due to the convergence in distribution of $(Z_i^N(t))$ and the standard results on convergence of point processes, $\mathcal{P}_i^N$ converges to an non-homogeneous Poisson process $\mathcal{P}^\infty$ with intensity $(\nu \bar{R}^r(t))$. It can be written as follows

$$ \mathcal{P}^\infty(dt) = \int_{\mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \bar{R}^r(t^-)\}} \bar{\mathcal{N}}_{\nu,1}(dt, dh) $$

where, with our notations, $\bar{\mathcal{N}}_{\nu,1}$ is a Poisson process with intensity $\nu dhdt$, using the characterization of a Poisson process by the martingale of its stochastic integral.
So the evolution equations of $(Z^N(t))$ can be obtained from equations (3)–(5) as follows

$$ (6) \qquad dR_i^{r,N}(t) = \sum_{j=1}^{N} \mathbf{1}_{\{V_j^N(t^-)>0, S_i^N(t^-)<K\}} \mathcal{N}_{\lambda,j}^{U,N}(dt, \{i\}) \\ \qquad \qquad - \sum_{j=1}^{N} \left( \sum_{l=1}^{\infty} \mathbf{1}_{\{l \le X_{j,i}^N(t^-)\}} \mathcal{N}_{\nu,j,i,l}(dt) \right), $$

$$ (7) \qquad dR_i^N(t) = \sum_{l=1}^{\infty} \sum_{j=1}^{N} \mathbf{1}_{\{l \le X_{j,i}^N(t^-)\}} \mathcal{N}_{\nu,j,i,l}(dt) - \sum_{l=1}^{\infty} \mathbf{1}_{\{l \le R_i^N(t^-)\}} \mathcal{N}_{\mu,i,l}(dt), $$

$$ (8) \qquad dV_i^N(t) = \sum_{l=1}^{\infty} \mathbf{1}_{\{l \le R_i^N(t^-)\}} \mathcal{N}_{\mu,i,l}(dt) \\ \qquad \qquad - \sum_{j=1}^{N} \mathbf{1}_{\{V_i^N(t^-)>0, S_j^N(t^-)<K\}} \mathcal{N}_{\lambda,i}^{U,N}(dt, \{j\}), $$

$$ (9) \qquad dV_i^{r,N}(t) = \sum_{j=1}^{N} \mathbf{1}_{\{V_j^N(t^-)>0, S_j^N(t^-)<K\}} \mathcal{N}_{\lambda,i}^{U,N}(dt, \{j\}) \\ \qquad \qquad - \sum_{j=1}^{N} \left( \sum_{l=1}^{\infty} \mathbf{1}_{\{l \le X_{i,j}^N(t^-)\}} \mathcal{N}_{\nu,i,j,l}(dt) \right). $$

This process $(Z^N(t))$, in dimension $4N$, is not a Markov process because the evolution equations for $(Z^N(t))$ are not autonomous. They depend on the Markov process $(X^N(t))$ which lives in dimension $(N+1)^2$. See equation (2). We introduce process $(Z^N(t))$ because it is sufficient to capture the performance of the system, such as the fact that a station is empty or full. The quite remarkable property is that the limit as $N$ gets large of $(Z^N(t))$ is a non-linear Markov process. We will present this process in this section and prove the convergence result in the next section.

### 3.2. Heuristic computation of the asymptotic process.
The proof of the convergence will be given in the next section. Here we show how we can guess the asymptotic process. Suppose that, for $1 \le i \le N$, $(Z_i^N(t))$ converges in distribution to some process

$$ (\bar{Z}(t)) = (\bar{R}^r(t), \bar{R}(t), \bar{V}(t), \bar{V}^r(t)). $$

Let $\mathcal{P}_i^N$ be the random measure on $\mathbb{R}_+$ defined by

$$ (10) \qquad \mathcal{P}_i^N([0, t]) = \int_0^t \sum_{j=1}^N \left( \sum_{l=1}^\infty \mathbf{1}_{\{l \le X_{j,i}^N(s^-)\}} \mathcal{N}_{\nu,j,i,l}(ds) \right). $$

It shows that $\mathcal{P}_i^N$ is a counting process (with jump size 1) on $\mathbb{R}_+$ and its compensator is

$$ \nu \int_0^t \sum_{j=1}^N X_{j,i}^N(s) \, ds = \nu \int_0^t R_i^{r,N}(s) \, ds. $$

See Robert [20, Proposition A.9] for example. Thus, due to the convergence in distribution of $(Z_i^N(t))$ and the standard results on convergence of point processes, $\mathcal{P}_i^N$ converges to an non-homogeneous Poisson process $\mathcal{P}^\infty$ with intensity $(\nu \bar{R}^r(t))$. It can be written as follows

$$ \mathcal{P}^\infty(dt) = \int_{\mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \bar{R}^r(t^-)\}} \bar{\mathcal{N}}_{\nu,1}(dt, dh) $$

where, with our notations, $\bar{\mathcal{N}}_{\nu,1}$ is a Poisson process with intensity $\nu dhdt$, using the characterization of a Poisson process by the martingale of its stochastic integral.
Following the same lines, it gives also that $\tilde{P}_i^N$, the random measure on $\mathbb{R}_+$ defined by

$$ (11) \qquad \tilde{P}_i^N([0, t]) = \int_0^t \sum_{j=1}^N \left( \sum_{l=1}^\infty \mathbf{1}_{\{l \le X_{i,j}^N(s^-)\}} \mathcal{N}_{\nu,i,j,l}(ds) \right) $$

is a counting process (with jump size 1) on $\mathbb{R}_+$ with compensator

$$ \nu \int_0^t \sum_{j=1}^N X_{i,j}^N(s) \, ds = \nu \int_0^t V_i^{r,N}(s) \, ds. $$

It converges to an non-homogeneous Poisson process $\tilde{P}^\infty$ with intensity $(\nu \overline{V^r}(t))$, i.e.

$$ \tilde{P}^\infty(dt) = \int_{\mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \overline{V^r}(t^-)\}} \overline{N}_{\nu,2}(dt, dh). $$

**Remark 1.** Point processes $\overline{N}_{\nu,1}$ and $\overline{N}_{\nu,2}$ are independent because $P_i^N$ (respectively $\tilde{P}_i^N$) is a function of $(N_{\nu,i,j,\cdot}, j)$ (respectively $(N_{\nu,j,i,\cdot}, j)$), where $(N_{\nu,i,j,\cdot}, j \neq i)$ and $(N_{\nu,j,i,\cdot}, j \neq i)$ are independent.

Then let us consider the random measure $Q_i^N$ on $\mathbb{R}_+$ defined by

$$ (12) \qquad Q_i^N([0, t]) = \int_0^t \sum_{j=1}^N \mathbf{1}_{\{V_j^N(s^-) > 0, S_i^N(s^-) < K\}} \mathcal{N}_{\lambda,j}^{U,N}(ds, \{i\}) $$

with compensator

$$ \lambda \int_0^t \mathbf{1}_{\{S_i^N(s^-) < K\}} \frac{1}{N} \sum_{j=1}^N \mathbf{1}_{\{V_j^N(s^-) > 0\}} ds. $$

Heuristically, by asymptotic independence of stations and the law of large numbers, as $N$ tends to infinity,

$$ \frac{1}{N} \sum_{j=1}^{N} \mathbf{1}_{\{V_j^N(t)>0\}} \to \mathbb{P}(\overline{V}(t)>0). $$

Thus formally, as $N$ tends to $+\infty$, $Q_i^N$ converges to an non-homogeneous Poisson process $Q^\infty$ with intensity $(\lambda \mathbf{1}_{\{\overline{S}(t)<K\}} \mathbb{P}(\overline{V}(t)>0))$. In other words,

$$ Q^\infty(dt) = \int_{\mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \mathbf{1}_{\{\overline{S}(t^-)<K\}} \mathbb{P}(\overline{V}(t^-)>0)\}} \overline{N}_{\lambda,1}(dt, dh). $$

With exactly the same work, $\tilde{Q}_i^N$ on $\mathbb{R}_+$ defined by

$$ (13) \qquad \tilde{Q}_i^N([0, t]) = \int_0^t \sum_{j=1}^N \mathbf{1}_{\{V_i^N(s^-) > 0, S_j^N(s^-) < K\}} \mathcal{N}_{\lambda,i}^{U,N}(ds, \{j\}) $$

formally converges in distribution when $N$ tends to $+\infty$ to the non-homogeneous Poisson process $\tilde{Q}^\infty$ with intensity $(\lambda \mathbf{1}_{\{\overline{V}(t)>0\}} \mathbb{P}(\overline{S}(t)<K))$ also defined by

$$ \tilde{Q}^\infty(dt) = \int_{\mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \mathbf{1}_{\{\overline{V}(t^-)>0\}} \mathbb{P}(\overline{S}(t^-)<K)\}} \overline{N}_{\lambda,2}(dt, dh). $$
Following the same lines, it gives also that $\tilde{P}_i^N$, the random measure on $\mathbb{R}_+$ defined by

$$ (11) \qquad \tilde{P}_i^N([0, t]) = \int_0^t \sum_{j=1}^N \left( \sum_{l=1}^\infty \mathbf{1}_{\{l \le X_{i,j}^N(s^{-})\}} \mathcal{N}_{\nu,i,j,l}(ds) \right) $$

is a counting process (with jump size 1) on $\mathbb{R}_+$ with compensator

$$ \nu \int_0^t \sum_{j=1}^N X_{i,j}^N(s) ds = \nu \int_0^t V_i^{r,N}(s) ds. $$

It converges to an non-homogeneous Poisson process $\tilde{P}^\infty$ with intensity $(\nu \overline{V^r}(t))$, i.e.

$$ \tilde{P}^\infty(dt) = \int_{\mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \overline{V^r}(t^{-})\}} \overline{N}_{\nu,2}(dt, dh). $$

**Remark 1.** Point processes $\overline{N}_{\nu,1}$ and $\overline{N}_{\nu,2}$ are independent because $P_i^N$ (respectively $\tilde{P}_i^N$) is a function of $(N_{\nu,i,j,\cdot}, j)$ (respectively $(N_{\nu,j,i,\cdot}, j)$), where $(N_{\nu,i,j,\cdot}, j \neq i)$ and $(N_{\nu,j,i,\cdot}, j \neq i)$ are independent.

Then let us consider the random measure $Q_i^N$ on $\mathbb{R}_+$ defined by

$$ (12) \qquad Q_i^N([0, t]) = \int_0^t \sum_{j=1}^N \mathbf{1}_{\{V_j^N(s^{-}) > 0, S_i^N(s^{-}) < K\}} \mathcal{N}_{\lambda,j}^{U,N}(ds, \{i\}) $$

with compensator

$$ \lambda \int_0^t \mathbf{1}_{\{S_i^N(s) < K\}} \frac{1}{N} \sum_{j=1}^N \mathbf{1}_{\{V_j^N(s) > 0\}} ds. $$

Heuristically, by asymptotic independence of stations and the law of large numbers, as $N$ tends to infinity,

$$ \frac{1}{N} \sum_{j=1}^{N} \mathbf{1}_{\{V_j^N(t)>0\}} \to \mathbb{P}(\overline{V}(t)>0). $$

Thus formally, as $N$ tends to $+\infty$, $Q_i^N$ converges to an non-homogeneous Poisson process $Q^\infty$ with intensity $(\lambda \mathbf{1}_{\{\overline{S}(t)<K\}} \mathbb{P}(\overline{V}(t)>0))$. In other words,

$$ Q^\infty(dt) = \int_{\mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \mathbf{1}_{\{\overline{S}(t^{-})<K\}} \mathbb{P}(\overline{V}(t^{-})>0)\}} \overline{N}_{\lambda,1}(dt, dh). $$

With exactly the same work, $\tilde{Q}_i^N$ on $\mathbb{R}_+$ defined by

$$ (13) \qquad \tilde{Q}_i^N([0, t]) = \int_0^t \sum_{j=1}^N \mathbf{1}_{\{V_i^N(s^{-}) > 0, S_j^N(s^{-}) < K\}} \mathcal{N}_{\lambda,i}^{U,N}(ds, \{j\}) $$

formally converges in distribution when $N$ tends to $+\infty$ to the non-homogeneous Poisson process $\tilde{Q}^\infty$ with intensity $(\lambda \mathbf{1}_{\{\overline{V}(t)>0\}} \mathbb{P}(\overline{S}(t)<K))$ also defined by

$$ \tilde{Q}^\infty(dt) = \int_{\mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \mathbf{1}_{\{\overline{V}(t^{-})>0\}} \mathbb{P}(\overline{S}(t^{-})<K)\}} \overline{N}_{\lambda,2}(dt, dh). $$
Note that for the same reason as previously, $\overline{\mathcal{N}}_{\lambda,2}$ is independent of $\overline{\mathcal{N}}_{\lambda,1}$. By formally taking the limit in equations (6)–(9), it leads to

$$ (14) \quad \left\{ \begin{aligned} d\overline{R}^r(t) &= \int_{\mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \mathbf{1}_{\{\overline{S}(t^-) < K\}} \mathbb{P}(\overline{V}(t^-) > 0)\}} \overline{\mathcal{N}}_{\lambda,1}(dt, dh) - \int_{\mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \overline{R}^r(t^-)\}} \overline{\mathcal{N}}_{\nu,1}(dt, dh), \\ d\overline{R}(t) &= \int_{\mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \overline{R}^r(t^-)\}} \overline{\mathcal{N}}_{\nu,1}(dt, dh) - \sum_{l=1}^{+\infty} \mathbf{1}_{\{l \le \overline{R}(t^-)\}} \mathcal{N}_{\mu,l}(dt), \\ d\overline{V}(t) &= \sum_{l=1}^{+\infty} \mathbf{1}_{\{l \le \overline{R}(t^-)\}} \mathcal{N}_{\mu,l}(dt) - \int_{\mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \mathbf{1}_{\{\overline{V}(t^-) > 0\}} \mathbb{P}(\overline{S}(t^-) < K)\}} \overline{\mathcal{N}}_{\lambda,2}(dt, dh), \\ d\overline{V}^r(t) &= \int_{\mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \mathbf{1}_{\{\overline{V}(t^-) > 0\}} \mathbb{P}(\overline{S}(t^-) < K)\}} \overline{\mathcal{N}}_{\lambda,2}(dt, dh) - \int_{\mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \overline{V}^r(t^-)\}} \overline{\mathcal{N}}_{\nu,2}(dt, dh). \end{aligned} \right. $$

**3.3. A first result.** Then the first result gives the existence and uniqueness of a stochastic process solution of the system of SDEs (14). Let $T > 0$ be fixed. Let $\mathcal{D}_T = \mathcal{D}([0, T], \mathcal{P}(\Sigma_K))$ be the set of cad-lag functions from $[0, T]$ to $\mathcal{P}(\Sigma_K)$.

**Theorem 1 (McKean–Vlasov process).** For every $(w, x, y, z) \in \Sigma_K$, the system of equations

$$ (15) \quad \left\{ \begin{aligned} \overline{R}^r(t) &= w + \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \mathbf{1}_{\{\overline{S}(s^-) < K\}} \mathbb{P}(\overline{V}(s^-) > 0)\}} \overline{\mathcal{N}}_{\lambda,1}(ds, dh) \\ &\quad - \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \overline{R}^r(s^-)\}} \overline{\mathcal{N}}_{\nu,1}(ds, dh), \\ \overline{R}(t) &= x + \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \overline{R}^r(s^-)\}} \overline{\mathcal{N}}_{\nu,1}(ds, dh) \\ &\quad - \int_0^t \sum_{l=1}^{+\infty} \mathbf{1}_{\{l \le \overline{R}(s^-)\}} \mathcal{N}_{\mu,l}(ds), \\ \overline{V}(t) &= y + \int_0^t \sum_{l=1}^{+\infty} \mathbf{1}_{\{l \le \overline{R}(s^-)\}} \mathcal{N}_{\mu,l}(ds) \\ &\quad - \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \mathbf{1}_{\{\overline{V}(s^-) > 0\}} \mathbb{P}(\overline{S}(s^-) < K)\}} \overline{\mathcal{N}}_{\lambda,2}(ds, dh), \\ \overline{V}^r(t) &= z + \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \mathbf{1}_{\{\overline{V}(s^-) > 0\}} \mathbb{P}(\overline{S}(s^-) < K)\}} \overline{\mathcal{N}}_{\lambda,2}(ds, dh) \\ &\quad - \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \overline{V}^r(s^-)\}} \overline{\mathcal{N}}_{\nu,2}(ds, dh), \end{aligned} \right. $$

has a unique solution $(\overline{R}^r(t), \overline{R}(t), \overline{V}(t), \overline{V}^r(t))$ in $\mathcal{D}_T$.

Note that the solution of equation (15) satisfies the Fokker–Planck equation (1) in the introduction.

**Proof.** The proof is standard and quite technical. Let us introduce the Wasserstein distances on $\mathcal{P}(\mathcal{D}_T)$. For $\pi_1, \pi_2 \in \mathcal{P}(\mathcal{D}_T)$,

$$ W_T(\pi_1, \pi_2) = \inf_{\pi \in \mathcal{C}_T(\pi_1, \pi_2)} \int_{\omega=(\omega_1, \omega_2) \in \mathcal{D}_T^2} (d_T(\omega_1, \omega_2) \wedge 1) d\pi(\omega), \\ \rho_T(\pi_1, \pi_2) = \inf_{\pi \in \mathcal{C}_T(\pi_1, \pi_2)} \int_{\omega=(\omega_1, \omega_2) \in \mathcal{D}_T^2} (\|\omega_1 - \omega_2\|_{\infty, T} \wedge 1) d\pi(\omega) $$

where, for $f \in \mathcal{D}_T$, $\|f\|_{\infty,T} = \sup\{\|f\|, 0 \le t \le T\} = \sup\{\sum_{i=1}^4 |f_i(t)|, 0 \le t \le T\}$, $C_T(\pi_1, \pi_2)$ is the set of couplings of $\pi_1$ and $\pi_2$, i.e. the subset of $\mathcal{P}(\mathcal{D}_T^2)$ with first marginal is $\pi_1$ and the
Note that for the same reason as previously, $\overline{\mathcal{N}}_{\lambda,2}$ is independent of $\overline{\mathcal{N}}_{\lambda,1}$. By formally taking the limit in equations (6)–(9), it leads to

$$ (14) \quad \left\{ \begin{aligned} d\overline{R}^r(t) &= \int_{\mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \mathbf{1}_{\{\overline{S}(t^-) < K\}} \mathbb{P}(\overline{V}(t^-) > 0)\}} \overline{\mathcal{N}}_{\lambda,1}(dt, dh) - \int_{\mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \overline{R}^r(t^-)\}} \overline{\mathcal{N}}_{\nu,1}(dt, dh), \\ d\overline{R}(t) &= \int_{\mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \overline{R}^r(t^-)\}} \overline{\mathcal{N}}_{\nu,1}(dt, dh) - \sum_{l=1}^{+\infty} \mathbf{1}_{\{l \le \overline{R}(t^-)\}} \mathcal{N}_{\mu,l}(dt), \\ d\overline{V}(t) &= \sum_{l=1}^{+\infty} \mathbf{1}_{\{l \le \overline{R}(t^-)\}} \mathcal{N}_{\mu,l}(dt) - \int_{\mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \mathbf{1}_{\{\overline{V}(t^-) > 0\}} \mathbb{P}(\overline{S}(t^-) < K)\}} \overline{\mathcal{N}}_{\lambda,2}(dt, dh), \\ d\overline{V}^r(t) &= \int_{\mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \mathbf{1}_{\{\overline{V}(t^-) > 0\}} \mathbb{P}(\overline{S}(t^-) < K)\}} \overline{\mathcal{N}}_{\lambda,2}(dt, dh) - \int_{\mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \overline{V}^r(t^-)\}} \overline{\mathcal{N}}_{\nu,2}(dt, dh). \end{aligned} \right. $$

**3.3. A first result.** Then the first result gives the existence and uniqueness of a stochastic process solution of the system of SDEs (14). Let $T > 0$ be fixed. Let $\mathcal{D}_T = \mathcal{D}([0, T], \mathcal{P}(\Sigma_K))$ be the set of cad-lag functions from $[0, T]$ to $\mathcal{P}(\Sigma_K)$.

**Theorem 1 (McKean–Vlasov process).** For every $(w, x, y, z) \in \Sigma_K$, the system of equations

$$ (15) \quad \left\{ \begin{aligned} \overline{R}^r(t) &= w + \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \mathbf{1}_{\{\overline{S}(s^-) < K\}} \mathbb{P}(\overline{V}(s^-) > 0)\}} \overline{\mathcal{N}}_{\lambda,1}(ds, dh) \\ &\quad - \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \overline{R}^r(s^-)\}} \overline{\mathcal{N}}_{\nu,1}(ds, dh), \\ \overline{R}(t) &= x + \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \overline{R}^r(s^-)\}} \overline{\mathcal{N}}_{\nu,1}(ds, dh) \\ &\quad - \int_0^t \sum_{l=1}^{+\infty} \mathbf{1}_{\{l \le \overline{R}(s^-)\}} \mathcal{N}_{\mu,l}(ds), \\ \overline{V}(t) &= y + \int_0^t \sum_{l=1}^{+\infty} \mathbf{1}_{\{l \le \overline{R}(s^-)\}} \mathcal{N}_{\mu,l}(ds) \\ &\quad - \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \mathbf{1}_{\{\overline{V}(s^-) > 0\}} \mathbb{P}(\overline{S}(s^-) < K)\}} \overline{\mathcal{N}}_{\lambda,2}(ds, dh), \\ \overline{V}^r(t) &= z + \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \mathbf{1}_{\{\overline{V}(s^-) > 0\}} \mathbb{P}(\overline{S}(s^-) < K)\}} \overline{\mathcal{N}}_{\lambda,2}(ds, dh) \\ &\quad - \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \overline{V}^r(s^-)\}} \overline{\mathcal{N}}_{\nu,2}(ds, dh), \end{aligned} \right. $$

has a unique solution $(\overline{R}^r(t), \overline{R}(t), \overline{V}(t), \overline{V}^r(t))$ in $\mathcal{D}_T$.

Note that the solution of equation (15) satisfies the Fokker–Planck equation (1) in the introduction.

**Proof.** The proof is standard and quite technical. Let us introduce the Wasserstein distances on $\mathcal{P}(\mathcal{D}_T)$. For $\pi_1, \pi_2 \in \mathcal{P}(\mathcal{D}_T)$,

$$ \begin{aligned} W_T(\pi_1, \pi_2) &= \inf_{\pi \in \mathcal{C}_T(\pi_1, \pi_2)} \int_{\omega=(\omega_1, \omega_2) \in \mathcal{D}_T^2} (d_T(\omega_1, \omega_2) \wedge 1) d\pi(\omega), \\ \rho_T(\pi_1, \pi_2) &= \inf_{\pi \in \mathcal{C}_T(\pi_1, \pi_2)} \int_{\omega=(\omega_1, \omega_2) \in \mathcal{D}_T^2} (\|\omega_1 - \omega_2\|_{\infty, T} \wedge 1) d\pi(\omega) \end{aligned} $$

where, for $f \in \mathcal{D}_T$, $\|f\|_{\infty,T} = \sup\{\|f\|, 0 \le t \le T\} = \sup\{\sum_{i=1}^4 |f_i(t)|, 0 \le t \le T\}$, $C_T(\pi_1, \pi_2)$ is the set of couplings of $\pi_1$ and $\pi_2$, i.e. the subset of $\mathcal{P}(\mathcal{D}_T^2)$ with first marginal is $\pi_1$ and the
second is $\pi_2$, ($\mathcal{D}_T$, $d_T$) with $d_T$ the distance associated to the Skorohod topology is complete and separable thus ($\mathcal{P}(\mathcal{D}_T)$, $W_T$) is complete and separable and $W_T(\pi_1, \pi_2) \le \rho_T(\pi_1, \pi_2)$.

Let us define

$$ \Phi : (\mathcal{P}(\mathcal{D}_T), W_T) \to (\mathcal{P}(\mathcal{D}_T), W_T) \\ \pi \mapsto \Phi(\pi) $$

where $\Phi(\pi)$ is the distribution of $(Z_\pi(t)) = (R_\pi^r(t), R_\pi(t), V_\pi(t), V_\pi^r(t))$, unique solution of the

$$ \text{SDEs} \quad \left\{ \begin{aligned} R_{\pi}^{r}(t) &= w + \iint_{[0,t] \times \mathbb{R}_{+}} \mathbf{1}_{\{0 \le h \le \mathbf{1}_{\{\|Z_{\pi}(s^{-})\| < K\}} \pi(v(s^{-}) > 0)\}} \overline{\mathcal{N}}_{\lambda,1}(ds, dh) \\ & \quad - \iint_{[0,t] \times \mathbb{R}_{+}} \mathbf{1}_{\{0 \le h \le R_{\pi}^{r}(s^{-})\}} \overline{\mathcal{N}}_{\nu,1}(ds, dh), \\ R_{\pi}(t) &= x + \iint_{[0,t] \times \mathbb{R}_{+}} \mathbf{1}_{\{0 \le h \le R_{\pi}^{r}(s^{-})\}} \overline{\mathcal{N}}_{\nu,1}(ds, dh) \\ & \quad - \int_{0}^{t} \sum_{l=1}^{+\infty} \mathbf{1}_{\{l \le R_{\pi}(s^{-})\}} \mathcal{N}_{\mu,l}(ds), \\ V_{\pi}(t) &= y + \int_{0}^{t} \sum_{l=1}^{+\infty} \mathbf{1}_{\{l \le R_{\pi}(s^{-})\}} \mathcal{N}_{\mu,l}(ds) \\ & \quad - \iint_{[0,t] \times \mathbb{R}_{+}} \mathbf{1}_{\{0 \le h \le \mathbf{1}_{\{V_{\pi}(s^{-}) > 0\}} \pi(\|z(s^{-})\| < K)\}} \overline{\mathcal{N}}_{\lambda,2}(ds, dh) \\ V_{\pi}^{r}(t) &= z + \iint_{[0,t] \times \mathbb{R}_{+}} \mathbf{1}_{\{0 \le h \le \mathbf{1}_{\{V_{\pi}(s^{-}) > 0\}} \pi(\|z(s^{-})\| < K)\}} \overline{\mathcal{N}}_{\lambda,2}(ds, dh) \\ & \quad - \iint_{[0,t] \times \mathbb{R}_{+}} \mathbf{1}_{\{0 \le h \le V_{\pi}^{r}(s^{-})\}} \overline{\mathcal{N}}_{\nu,2}(ds, dh). \end{aligned} \right. \quad (16) $$

Note that $\pi(v(t) > 0) = \int_{z=(r^r, r, v, v^r) \in \mathcal{D}_T} \mathbf{1}_{\{v(t) > 0\}} d\pi(z)$. The existence and uniqueness of a solution to (15) is equivalent to the existence and uniqueness of a fixed point $\pi = \Phi(\pi)$.

Let us prove that $\pi = \Phi(\pi)$ has a unique fixed point. For $\pi_1, \pi_2 \in \mathcal{P}(\mathcal{D}_T)$, let $Z_{\pi_1}$ and $Z_{\pi_2}$ be solutions of (16). Thus $(Z_{\pi_1}, Z_{\pi_2})$ is a coupling of $\Phi(\pi_1)$ and $\Phi(\pi_2)$ and, for $t \le T$, $\rho_t(\Phi(\pi_1), \Phi(\pi_2)) \le \mathbb{E}(\|Z_{\pi_1} - Z_{\pi_2}\|_{\infty,t})$.

For $t \le T$, using the definition (16) of $Z_{\pi_1}$ and $Z_{\pi_2}$,

$$ \begin{aligned} \|Z_{\pi_1} - Z_{\pi_2}\|_{\infty,t} &= \sup_{0 \le s \le t} (|R_{\pi_1}^r(s) - R_{\pi_2}^r(s)| + |R_{\pi_1}(s) - R_{\pi_2}(s)| + |V_{\pi_1}(s) - V_{\pi_2}(s)| + |V_{\pi_1}^r(s) - V_{\pi_2}^r(s)|) \\ &\le \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{A_{\pi_1}(s^{-}) \wedge A_{\pi_2}(s^{-}) \le h \le A_{\pi_1}(s^{-}) \vee A_{\pi_2}(s^{-})\}} \overline{\mathcal{N}}_{\lambda,1}(ds, dh) \\ & \quad + 2 \int_0^t \sum_{l=1}^{+\infty} |\mathbf{1}_{\{l \le R_{\pi_1}(s^{-})\}} - \mathbf{1}_{\{l \le R_{\pi_2}(s^{-})\}}| \mathcal{N}_{\mu,l}(ds) \\ & \quad + 2 \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{R_{\pi_1}^r(s^{-}) \wedge R_{\pi_2}^r(s^{-}) \le h \le R_{\pi_1}^r(s^{-}) \vee R_{\pi_2}^r(s^{-})\}} \overline{\mathcal{N}}_{\nu,1}(ds, dh) \\ & \quad + 2 \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{B_{\pi_1}(s^{-}) \wedge B_{\pi_2}(s^{-}) \le h \le B_{\pi_1}(s^{-}) \vee B_{\pi_2}(s^{-})\}} \overline{\mathcal{N}}_{\lambda,2}(ds, dh) \\ & \quad + \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{V_{\pi_1}^r(s^{-}) \wedge V_{\pi_2}^r(s^{-}) \le h \le V_{\pi_1}^r(s^{-}) \vee V_{\pi_2}^r(s^{-})\}} \overline{\mathcal{N}}_{\nu,2}(ds, dh). \end{aligned} \quad (17) $$

where

$$ A_{\pi}(t) = \mathbf{1}_{\{\|Z_{\pi}(t)\| < K\}} \pi(v(t) > 0) \quad \text{and} \quad B_{\pi}(t) = \mathbf{1}_{\{V_{\pi}(t) > 0\}} \pi(\|z(t)\| < K). $$
second is $\pi_2$, ($\mathcal{D}_T$, $d_T$) with $d_T$ the distance associated to the Skorohod topology is complete and separable thus ($\mathcal{P}(\mathcal{D}_T)$, $W_T$) is complete and separable and $W_T(\pi_1, \pi_2) \le \rho_T(\pi_1, \pi_2)$.

Let us define

$$ \Phi : (\mathcal{P}(\mathcal{D}_T), W_T) \to (\mathcal{P}(\mathcal{D}_T), W_T) \\ \pi \mapsto \Phi(\pi) $$

where $\Phi(\pi)$ is the distribution of $(Z_\pi(t)) = (R_\pi^r(t), R_\pi(t), V_\pi(t), V_\pi^r(t))$, unique solution of the

$$ \text{SDEs} \quad \left\{ \begin{aligned} R_{\pi}^{r}(t) &= w + \iint_{[0,t] \times \mathbb{R}_{+}} \mathbf{1}_{\{0 \le h \le \mathbf{1}_{\{\|Z_{\pi}(s^{-})\| < K\}} \pi(v(s^{-}) > 0)\}} \overline{\mathcal{N}}_{\lambda,1}(ds, dh) \\ & \quad - \iint_{[0,t] \times \mathbb{R}_{+}} \mathbf{1}_{\{0 \le h \le R_{\pi}^{r}(s^{-})\}} \overline{\mathcal{N}}_{\nu,1}(ds, dh), \\ R_{\pi}(t) &= x + \iint_{[0,t] \times \mathbb{R}_{+}} \mathbf{1}_{\{0 \le h \le R_{\pi}^{r}(s^{-})\}} \overline{\mathcal{N}}_{\nu,1}(ds, dh) \\ & \quad - \int_{0}^{t} \sum_{l=1}^{+\infty} \mathbf{1}_{\{l \le R_{\pi}(s^{-})\}} \mathcal{N}_{\mu,l}(ds), \\ V_{\pi}(t) &= y + \int_{0}^{t} \sum_{l=1}^{+\infty} \mathbf{1}_{\{l \le R_{\pi}(s^{-})\}} \mathcal{N}_{\mu,l}(ds) \\ & \quad - \iint_{[0,t] \times \mathbb{R}_{+}} \mathbf{1}_{\{0 \le h \le \mathbf{1}_{\{V_{\pi}(s^{-}) > 0\}} \pi(\|z(s^{-})\| < K)\}} \overline{\mathcal{N}}_{\lambda,2}(ds, dh) \\ V_{\pi}^{r}(t) &= z + \iint_{[0,t] \times \mathbb{R}_{+}} \mathbf{1}_{\{0 \le h \le \mathbf{1}_{\{V_{\pi}(s^{-}) > 0\}} \pi(\|z(s^{-})\| < K)\}} \overline{\mathcal{N}}_{\lambda,2}(ds, dh) \\ & \quad - \iint_{[0,t] \times \mathbb{R}_{+}} \mathbf{1}_{\{0 \le h \le V_{\pi}^{r}(s^{-})\}} \overline{\mathcal{N}}_{\nu,2}(ds, dh). \end{aligned} \right. \quad (16) $$

Note that $\pi(v(t) > 0) = \int_{z=(r^r, r, v, v^r) \in \mathcal{D}_T} \mathbf{1}_{\{v(t) > 0\}} d\pi(z)$. The existence and uniqueness of a solution to (15) is equivalent to the existence and uniqueness of a fixed point $\pi = \Phi(\pi)$.

Let us prove that $\pi = \Phi(\pi)$ has a unique fixed point. For $\pi_1, \pi_2 \in \mathcal{P}(\mathcal{D}_T)$, let $Z_{\pi_1}$ and $Z_{\pi_2}$ be solutions of (16). Thus $(Z_{\pi_1}, Z_{\pi_2})$ is a coupling of $\Phi(\pi_1)$ and $\Phi(\pi_2)$ and, for $t \le T$, $\rho_t(\Phi(\pi_1), \Phi(\pi_2)) \le \mathbb{E}(\|Z_{\pi_1} - Z_{\pi_2}\|_{\infty,t})$.

For $t \le T$, using the definition (16) of $Z_{\pi_1}$ and $Z_{\pi_2}$,

$$ \begin{aligned} \|Z_{\pi_1} - Z_{\pi_2}\|_{\infty,t} &= \sup_{0 \le s \le t} (|R_{\pi_1}^r(s) - R_{\pi_2}^r(s)| + |R_{\pi_1}(s) - R_{\pi_2}(s)| + |V_{\pi_1}(s) - V_{\pi_2}(s)| + |V_{\pi_1}^r(s) - V_{\pi_2}^r(s)|) \\ &\le \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{A_{\pi_1}(s^{-}) \wedge A_{\pi_2}(s^{-}) \le h \le A_{\pi_1}(s^{-}) \vee A_{\pi_2}(s^{-})\}} \overline{\mathcal{N}}_{\lambda,1}(ds, dh) \\ & \quad + 2 \int_0^t \sum_{l=1}^{+\infty} |\mathbf{1}_{\{l \le R_{\pi_1}(s^{-})\}} - \mathbf{1}_{\{l \le R_{\pi_2}(s^{-})\}}| \mathcal{N}_{\mu,l}(ds) \\ & \quad + 2 \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{R_{\pi_1}^r(s^{-}) \wedge R_{\pi_2}^r(s^{-}) \le h \le R_{\pi_1}^r(s^{-}) \vee R_{\pi_2}^r(s^{-})\}} \overline{\mathcal{N}}_{\nu,1}(ds, dh) \\ & \quad + 2 \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{B_{\pi_1}(s^{-}) \wedge B_{\pi_2}(s^{-}) \le h \le B_{\pi_1}(s^{-}) \vee B_{\pi_2}(s^{-})\}} \overline{\mathcal{N}}_{\lambda,2}(ds, dh) \\ & \quad + \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{V_{\pi_1}^r(s^{-}) \wedge V_{\pi_2}^r(s^{-}) \le h \le V_{\pi_1}^r(s^{-}) \vee V_{\pi_2}^r(s^{-})\}} \overline{\mathcal{N}}_{\nu,2}(ds, dh). \end{aligned} \quad (17) $$

where

$$ A_{\pi}(t) = \mathbf{1}_{\{\|Z_{\pi}(t)\| < K\}} \pi(v(t) > 0) \quad \text{and} \quad B_{\pi}(t) = \mathbf{1}_{\{V_{\pi}(t) > 0\}} \pi(\|z(t)\| < K). $$
The mean of each term of the right-hand side of the previous equation is bounded as follows. For the first term,

$$
\begin{aligned}
& \mathbb{E} \left( \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{A_{\pi_1}(s^-) \wedge A_{\pi_2}(s^-) \le h \le A_{\pi_1}(s^-) \vee A_{\pi_2}(s^-)\}} \overline{\mathcal{N}}_{\lambda,1}(ds, dh) \right) \\
& \le \lambda \int_0^t |\pi_1(v(s) > 0) - \pi_2(v(s) > 0)| ds \\
& = \lambda \int_0^t |\pi(z, v_1(s) > 0) - \pi(z, v_2(s) > 0)| ds \\
& = \lambda \int_0^t \int_{\omega=(z_1,z_2) \in \mathcal{D}_T^2} |\mathbf{1}_{\{v_1(s)>0\}} - \mathbf{1}_{\{v_2(s)>0\}}| \pi(d\omega) ds \\
& \le \lambda \int_0^t \int_{\omega=(z_1,z_2) \in \mathcal{D}_T^2} |v_1(s) - v_2(s)| \wedge 1 \ \pi(d\omega) ds \\
& \le \lambda \int_0^t \rho_s(\pi_1, \pi_2) ds.
\end{aligned}
$$

For the second term of the right-hand side of equation (17),

$$
\begin{aligned}
& \mathbb{E} \left( \int_0^t \sum_{l=1}^{+\infty} \left| \mathbf{1}_{\{l \le R_{\pi_1}(s^-)\}} - \mathbf{1}_{\{l \le R_{\pi_2}(s^-)\}} \right| \mathcal{N}_{\mu,l}(ds) \right) \\
& \le \mu \int_0^t \mathbb{E} \left( \sum_{l=1}^{+\infty} \left| \mathbf{1}_{\{l \le R_{\pi_1}(s)\}} - \mathbf{1}_{\{l \le R_{\pi_2}(s)\}} \right| \right) ds \\
& \le \mu \int_0^t \mathbb{E} \left( |R_{\pi_1}(s) - R_{\pi_2}(s)| \right) ds \\
& \le \mu \int_0^t \mathbb{E} \left( \| Z_{\pi_1} - Z_{\pi_2} \|_{\infty,s} \right) ds.
\end{aligned}
$$

For the third term of the right-hand side of (17),

$$
\begin{aligned}
& \mathbb{E} \left( \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{R_{\pi_1}^r(s^-) \wedge R_{\pi_2}^r(s^-) \le h \le R_{\pi_1}^r(s^-) \vee R_{\pi_2}^r(s^-)\}} \overline{\mathcal{N}}_{\nu,1}(ds, dh) \right) \\
& \le \nu \int_0^t \mathbb{E}(|R_{\pi_1}^r(s) - R_{\pi_2}^r(s)|) ds \\
& \le \nu \int_0^t \mathbb{E}(\|Z_{\pi_1} - Z_{\pi_2}\|_{\infty,s}) ds.
\end{aligned}
$$
The mean of each term of the right-hand side of the previous equation is bounded as follows. For the first term,

$$
\begin{aligned}
& \mathbb{E} \left( \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{A_{\pi_1}(s^-) \wedge A_{\pi_2}(s^-) \le h \le A_{\pi_1}(s^-) \vee A_{\pi_2}(s^-)\}} \overline{\mathcal{N}}_{\lambda,1}(ds, dh) \right) \\
& \le \lambda \int_0^t |\pi_1(v(s) > 0) - \pi_2(v(s) > 0)| ds \\
& = \lambda \int_0^t |\pi(z, v_1(s) > 0) - \pi(z, v_2(s) > 0)| ds \\
& = \lambda \int_0^t \int_{\omega=(z_1,z_2) \in \mathcal{D}_T^2} |\mathbf{1}_{\{v_1(s)>0\}} - \mathbf{1}_{\{v_2(s)>0\}}| \pi(d\omega) ds \\
& \le \lambda \int_0^t \int_{\omega=(z_1,z_2) \in \mathcal{D}_T^2} |v_1(s) - v_2(s)| \wedge 1 \ \pi(d\omega) ds \\
& \le \lambda \int_0^t \rho_s(\pi_1, \pi_2) ds.
\end{aligned}
$$

For the second term of the right-hand side of equation (17),

$$
\begin{aligned}
& \mathbb{E} \left( \int_0^t \sum_{l=1}^{+\infty} \left| \mathbf{1}_{\{l \le R_{\pi_1}(s^-)\}} - \mathbf{1}_{\{l \le R_{\pi_2}(s^-)\}} \right| \mathcal{N}_{\mu,l}(ds) \right) \\
& \le \mu \int_0^t \mathbb{E} \left( \sum_{l=1}^{+\infty} \left| \mathbf{1}_{\{l \le R_{\pi_1}(s)\}} - \mathbf{1}_{\{l \le R_{\pi_2}(s)\}} \right| \right) ds \\
& \le \mu \int_0^t \mathbb{E} \left( |R_{\pi_1}(s) - R_{\pi_2}(s)| \right) ds \\
& \le \mu \int_0^t \mathbb{E} \left( \|Z_{\pi_1} - Z_{\pi_2}\|_{\infty,s} \right) ds.
\end{aligned}
$$

For the third term of the right-hand side of (17),

$$
\begin{aligned}
& \mathbb{E} \left( \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{R_{\pi_1}^r(s^-) \wedge R_{\pi_2}^r(s^-) \le h \le R_{\pi_1}^r(s^-) \vee R_{\pi_2}^r(s^-)\}} \overline{\mathcal{N}}_{\nu,1}(ds, dh) \right) \\
& \le \nu \int_0^t \mathbb{E}(|R_{\pi_1}^r(s) - R_{\pi_2}^r(s)|) ds \\
& \le \nu \int_0^t \mathbb{E}(\|Z_{\pi_1} - Z_{\pi_2}\|_{\infty,s}) ds.
\end{aligned}
$$
For the fourth term,

$$
\begin{align*}
& \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{B_{\pi_1}(s^-) \wedge B_{\pi_2}(s^-) \le h \le B_{\pi_1}(s^-) \vee B_{\pi_2}(s^-)\}} \overline{\mathcal{N}}_{\lambda,2}(ds, dh) \\
& \le \lambda \int_0^t \mathbb{E} \left( \left| \mathbf{1}_{\{V_{\pi_1}(s)>0\}} \pi_1(\|z(s)\| < K) - \mathbf{1}_{\{V_{\pi_2}(s)>0\}} \pi_2(\|z(s)\| < K) \right| \right) ds \\
& \le \lambda \int_0^t \mathbb{E} \left( \left| \mathbf{1}_{\{V_{\pi_1}(s)>0\}} - \mathbf{1}_{\{V_{\pi_2}(s)>0\}} \right| \right) ds \\
& \le \lambda \int_0^t \mathbb{E} \left( |V_{\pi_1}(s) - V_{\pi_2}(s)| \right) ds \\
& \le \lambda \int_0^t \mathbb{E} \left( \|Z_{\pi_1} - V_{\pi_2}\|_{\infty,s} \right) ds.
\end{align*}
$$

As for the third term, for the fifth term,

$$
\begin{align*}
& \mathbb{E} \left( \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{V_{\pi_1}^r(s^-) \wedge V_{\pi_2}^r(s^-) \le h \le V_{\pi_1}^r(s^-) \vee V_{\pi_2}^r(s^-)\}} \overline{\mathcal{N}}_{\nu,2}(ds, dh) \right) \\
& \le \nu \int_0^t \mathbb{E}(\|Z_{\pi_1} - Z_{\pi_2}\|_{\infty,s}) ds.
\end{align*}
$$

Thus

$$
\mathbb{E}(\|Z_{\pi_1} - Z_{\pi_2}\|_{\infty,t}) \leq (2\mu + 3\nu + 2\lambda) \int_0^t \mathbb{E}(\|Z_{\pi_1} - Z_{\pi_2}\|_{\infty,s}) ds + \lambda \int_0^t \rho_s(\pi_1, \pi_2) ds.
$$

By Grönwall's inequality,

$$
\mathbb{E}(\|Z_{\pi_1} - Z_{\pi_2}\|_{\infty,t}) \le C_t \int_0^t \rho_s(\pi_1, \pi_2) ds
$$

with $C_t = \lambda \exp((2\mu + 3\nu + 2\lambda)t)$. For $t \le T$,

$$
(18) \qquad \rho_t(\Phi(\pi_1), \Phi(\pi_2)) \le C_T \int_0^t \rho_s(\pi_1, \pi_2) ds.
$$

It leads to the uniqueness of the solution of $\Phi(\pi) = \pi$. Indeed, if $\pi_1$ and $\pi_2$ are fixed points of $\Phi$, then equation (18) is rewritten

$$
\rho_t(\pi_1, \pi_2) \le C_T \int_0^t \rho_s(\pi_1, \pi_2) ds.
$$

By Grönwall's inequality again, for each $t \le T$, $\rho_t(\pi_1, \pi_2) = 0$ thus $\pi_1 = \pi_2$. The existence is proved by an iteration argument. Let $\pi_0 \in \mathcal{P}(\mathcal{D}_T)$ and $\pi_{n+1} = \Phi(\pi_n)$. By equation (18),

$$
\begin{align*}
W_T(\pi_{n+1}, \pi_n) &\le \rho_T(\pi_{n+1}, \pi_n) \\
&\le C_T^n \rho_T(\pi_1, \pi_0) \int_{0\le s_1 \le s_2 \dots s_n \le T} ds_1 \dots ds_n \\
&\le \frac{(C_T T)^n}{n!} \rho_T(\pi_1, \pi_0).
\end{align*}
$$

Thus $(\pi_n)$ converges since $(\mathcal{P}(\mathcal{D}_T), W_T)$ is complete. Because of equation (18), $\Phi$ is continuous for the Skorohod topology and its limit is a fixed point of $\Phi$. $\square$
For the fourth term,

$$
\begin{align*}
& \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{B_{\pi_1}(s^-) \wedge B_{\pi_2}(s^-) \le h \le B_{\pi_1}(s^-) \vee B_{\pi_2}(s^-)\}} \overline{\mathcal{N}}_{\lambda,2}(ds, dh) \\
& \le \lambda \int_0^t \mathbb{E} \left( \left| \mathbf{1}_{\{V_{\pi_1}(s)>0\}} \pi_1(\|z(s)\| < K) - \mathbf{1}_{\{V_{\pi_2}(s)>0\}} \pi_2(\|z(s)\| < K) \right| \right) ds \\
& \le \lambda \int_0^t \mathbb{E} \left( \left| \mathbf{1}_{\{V_{\pi_1}(s)>0\}} - \mathbf{1}_{\{V_{\pi_2}(s)>0\}} \right| \right) ds \\
& \le \lambda \int_0^t \mathbb{E} \left( |V_{\pi_1}(s) - V_{\pi_2}(s)| \right) ds \\
& \le \lambda \int_0^t \mathbb{E} \left( \|Z_{\pi_1} - V_{\pi_2}\|_{\infty,s} \right) ds.
\end{align*}
$$

As for the third term, for the fifth term,

$$
\begin{align*}
& \mathbb{E} \left( \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{V_{\pi_1}^r(s^-) \wedge V_{\pi_2}^r(s^-) \le h \le V_{\pi_1}^r(s^-) \vee V_{\pi_2}^r(s^-)\}} \overline{\mathcal{N}}_{\nu,2}(ds, dh) \right) \\
& \le \nu \int_0^t \mathbb{E}(\|Z_{\pi_1} - Z_{\pi_2}\|_{\infty,s}) ds.
\end{align*}
$$

Thus

$$
\mathbb{E}(\|Z_{\pi_1} - Z_{\pi_2}\|_{\infty,t}) \leq (2\mu + 3\nu + 2\lambda) \int_0^t \mathbb{E}(\|Z_{\pi_1} - Z_{\pi_2}\|_{\infty,s}) ds + \lambda \int_0^t \rho_s(\pi_1, \pi_2) ds.
$$

By Grönwall's inequality,

$$
\mathbb{E}(\|Z_{\pi_1} - Z_{\pi_2}\|_{\infty,t}) \le C_t \int_0^t \rho_s(\pi_1, \pi_2) ds
$$

with $C_t = \lambda \exp((2\mu + 3\nu + 2\lambda)t)$. For $t \le T$,

$$
(18) \qquad \rho_t(\Phi(\pi_1), \Phi(\pi_2)) \le C_T \int_0^t \rho_s(\pi_1, \pi_2) ds.
$$

It leads to the uniqueness of the solution of $\Phi(\pi) = \pi$. Indeed, if $\pi_1$ and $\pi_2$ are fixed points of $\Phi$,
then equation (18) is rewritten

$$
\rho_t(\pi_1, \pi_2) \le C_T \int_0^t \rho_s(\pi_1, \pi_2) ds.
$$

By Grönwall's inequality again, for each $t \le T$, $\rho_t(\pi_1, \pi_2) = 0$ thus $\pi_1 = \pi_2$. The existence is proved by an iteration argument. Let $\pi_0 \in \mathcal{P}(\mathcal{D}_T)$ and $\pi_{n+1} = \Phi(\pi_n)$. By equation (18),

$$
\begin{align*}
W_T(\pi_{n+1}, \pi_n) &\le \rho_T(\pi_{n+1}, \pi_n) \\
&\le C_T^n \rho_T(\pi_1, \pi_0) \int_{0\le s_1 \le s_2 \dots s_n \le T} ds_1 \dots ds_n \\
&\le \frac{(C_T T)^n}{n!} \rho_T(\pi_1, \pi_0).
\end{align*}
$$

Thus $(\pi_n)$ converges since $(\mathcal{P}(\mathcal{D}_T), W_T)$ is complete. Because of equation (18), $\Phi$ is continuous for the Skorohod topology and its limit is a fixed point of $\Phi$. $\square$
## 4. MEAN-FIELD LIMIT

Recall that, by Theorem 1, the so-called McKean-Vlasov process ($\bar{Z}(t)$) is the unique solution of the system of equations (15). The empirical distribution $\Lambda^N(t)$ of ($Z_i^N(t)$, $1 \le i \le N$) is defined, for $f$ on $\Sigma_K$, by

$$ \Lambda^N(t)(f) = \frac{1}{N} \sum_{i=1}^{N} f(Z_i^N(t)) = \frac{1}{N} \sum_{i=1}^{N} f(R_i^{r,N}(t), R_i^N(t), V_i^N(t), V_i^{r,N}(t)) $$

Process ($\Lambda^N(t)$) takes values in

$$ (19) \qquad \mathcal{Y}^N = \{ \Lambda \in \mathcal{P}(\Sigma_K), \Lambda_{(w,x,y,z)} \in \frac{\mathbb{N}}{N}, (w,x,y,z) \in \Sigma_K, \sum_{(w,x,y,z) \in \Sigma_K} (x+y+z)\Lambda_{(w,x,y,z)}N = M_N \}. $$

As process $Z^N$ is not Markov, process ($\Lambda^N(t)$) is not Markov. The aim of this section is to prove the mean-field convergence for ($Z^N(t)$), i.e. that the sequence of processes ($\Lambda^N(t)$) converges in distribution to ($\bar{Z}(t)$). It means that, for any function $f$ with finite support, the sequence of processes ($\Lambda^N(t)(f)$) converges in distribution to ($\mathbb{E}(f(\bar{Z}(t)))$).

**Theorem 2 (Mean-field convergence).** The sequence of empirical distribution process ($\Lambda^N(t)$) converges in distribution to a process ($\Lambda(t)$) $\in \mathcal{D}([0,T], \mathcal{P}(\Sigma_K))$ defined by, for $f$ with finite support on $\Sigma_K$,

$$ \Lambda(t)(f) = \mathbb{E}(f(\bar{Z}(t))) $$

with ($\bar{Z}(t)$) the unique solution of equation (15). Moreover, for any $k \ge 1$ and for $1 \le i_1 < \dots < i_k \le N$, the sequence of finite marginals ($Z_{i_1}^N(t), \dots, Z_{i_k}^N(t)$) converges in distribution to ($\bar{Z}_{i_1}(t), \dots, \bar{Z}_{i_k}(t)$), where ($\bar{Z}_{i_1}(t), \dots, \bar{Z}_{i_k}(t)$) are independent random variables with the same distribution as ($\bar{Z}(t)$).

The last property is the propagation of chaos. The proof is presented in Section 4.3.

### 4.1. Evolution equations of the empirical measure.
Let us introduce the following notations. For $z \in \Sigma_K$, $f: \Sigma_K \to \mathbb{R}_+$ and ($e_i$, $1 \le i \le 4$) the canonical basis of $\mathbb{R}^4$ i.e. $e_1 = (1, 0, 0, 0), \dots, e_4 = (0, 0, 0, 1)$,

$$ \begin{aligned} \Delta_{i,i+1}(f)(z) &= f(z - e_i + e_{i+1}) - f(z), & 1 \le i \le 3, \\ \Delta_1^+(f)(z) &= f(z + e_1) - f(z), \\ \Delta_4^-(f)(z) &= f(z - e_4) - f(z). \end{aligned} $$
## 4. MEAN-FIELD LIMIT

Recall that, by Theorem 1, the so-called McKean-Vlasov process ($\bar{Z}(t)$) is the unique solution of the system of equations (15). The empirical distribution $\Lambda^N(t)$ of ($Z_i^N(t)$, $1 \le i \le N$) is defined, for $f$ on $\Sigma_K$, by

$$ \Lambda^N(t)(f) = \frac{1}{N} \sum_{i=1}^{N} f(Z_i^N(t)) = \frac{1}{N} \sum_{i=1}^{N} f(R_i^{r,N}(t), R_i^N(t), V_i^N(t), V_i^{r,N}(t)) $$

Process ($\Lambda^N(t)$) takes values in

$$ (19) \qquad \mathcal{Y}^N = \{ \Lambda \in \mathcal{P}(\Sigma_K), \Lambda_{(w,x,y,z)} \in \frac{\mathbb{N}}{N}, (w,x,y,z) \in \Sigma_K, \sum_{(w,x,y,z) \in \Sigma_K} (x+y+z) \Lambda_{(w,x,y,z)} N = M_N \}. $$

As process $Z^N$ is not Markov, process ($\Lambda^N(t)$) is not Markov. The aim of this section is to prove the mean-field convergence for ($Z^N(t)$), i.e. that the sequence of processes ($\Lambda^N(t)$) converges in distribution to ($\bar{Z}(t)$). It means that, for any function $f$ with finite support, the sequence of processes ($\Lambda^N(t)(f)$) converges in distribution to ($\mathbb{E}(f(\bar{Z}(t)))$).

**Theorem 2 (Mean-field convergence).** The sequence of empirical distribution process ($\Lambda^N(t)$) converges in distribution to a process ($\Lambda(t)$) $\in \mathcal{D}([0, T], \mathcal{P}(\Sigma_K))$ defined by, for $f$ with finite support on $\Sigma_K$,

$$ \Lambda(t)(f) = \mathbb{E}(f(\bar{Z}(t))) $$

with ($\bar{Z}(t)$) the unique solution of equation (15). Moreover, for any $k \ge 1$ and for $1 \le i_1 < \dots < i_k \le N$, the sequence of finite marginals ($Z_{i_1}^N(t), \dots, Z_{i_k}^N(t)$) converges in distribution to ($\bar{Z}_{i_1}(t), \dots, \bar{Z}_{i_k}(t)$), where ($\bar{Z}_{i_1}(t), \dots, \bar{Z}_{i_k}(t)$) are independent random variables with the same distribution as ($\bar{Z}(t)$).

The last property is the propagation of chaos. The proof is presented in Section 4.3.

### 4.1. Evolution equations of the empirical measure.
Let us introduce the following notations. For $z \in \Sigma_K$, $f: \Sigma_K \to \mathbb{R}_+$ and $(e_i, 1 \le i \le 4)$ the canonical basis of $\mathbb{R}^4$ i.e. $e_1 = (1, 0, 0, 0), \dots, e_4 = (0, 0, 0, 1)$,

$$ \begin{aligned} \Delta_{i,i+1}(f)(z) &= f(z - e_i + e_{i+1}) - f(z), & 1 \le i \le 3, \\ \Delta_1^+(f)(z) &= f(z + e_1) - f(z), \\ \Delta_4^-(f)(z) &= f(z - e_4) - f(z). \end{aligned} $$
Let $\Sigma_K = \{(w, x, y, z) \in \mathbb{N}^4, w + x + y + z \le K\}$. For $f: \Sigma_K \to \mathbb{R}_+$, straightforwardly by equations (6)-(9),

$$
\begin{align*}
df(Z_i^N(t)) ={}& \Delta_1^+(f)(Z_i^N(t)) \sum_{j=1}^N \mathbf{1}_{\{V_j^N(t^-)>0, S_i^N(t^-)<K\}} \mathcal{N}_{\lambda,j}^{U,N}(dt, \{i\}) \\
& + \Delta_{1,2}(f)(Z_i^N(t)) \sum_{j=1}^N \left( \sum_{l=1}^\infty \mathbf{1}_{\{l \le X_{j,i}^N(t^-)\}} \mathcal{N}_{\nu,j,i,l}(dt) \right) \\
& + \Delta_{2,3}(f)(Z_i^N(t)) \sum_{l=1}^\infty \mathbf{1}_{\{l \le R_i^N(t^-)\}} \mathcal{N}_{\mu,i,l}(dt) \\
& + \Delta_{3,4}(f)(Z_i^N(t)) \sum_{j=1}^N \mathbf{1}_{\{V_i^N(t^-)>0, S_j^N(t^-)<K\}} \mathcal{N}_{\lambda,i}^{U,N}(dt, \{j\}) \\
& + \Delta_4^-(f)(Z_i^N(t)) \sum_{j=1}^N \left( \sum_{l=1}^\infty \mathbf{1}_{\{l \le X_{i,j}^N(t^-)\}} \mathcal{N}_{\nu,i,j,l}(dt) \right)
\end{align*}
$$

Thus, using the martingale decomposition for Poisson processes,

$$
\begin{align}
(20) \quad df(Z_i^N(t)) ={}& \Delta_1^+(f)(Z_i^N(t)) \mathbf{1}_{\{S_i^N(t)<K\}} \sum_{j=1}^N \mathbf{1}_{\{V_j^N(t)>0\}} \frac{\lambda}{N} dt \nonumber \\
& + \Delta_{1,2}(f)(Z_i^N(t)) R_i^{r,N}(t) \nu dt + \Delta_{2,3}(f)(Z_i^N(t)) R_i^N(t) \mu dt \nonumber \\
& + \Delta_{3,4}(f)(Z_i^N(t)) \mathbf{1}_{\{V_i^N(t)>0\}} \sum_{j=1}^N \mathbf{1}_{\{S_j^N(t)<K\}} \frac{\lambda}{N} dt \nonumber \\
& + \Delta_4^-(f)(Z_i^N(t)) V_i^{r,N}(t) \nu dt + d\mathcal{M}_{f,i}^N(t) \tag{20}
\end{align}
$$

where $(\mathcal{M}_{f,i}^N(t))$ is a martingale which will be detailed in Section 4.2. By summing equation (20) for $i$ from 1 to $N$ and dividing by $N$,

$$
\begin{align}
\Lambda^N(t)(f) ={}& \Lambda^N(0)(f) + \mathcal{M}_f^N(t) \nonumber \\
& + \lambda \int_0^t \Lambda^N(s) (\Sigma_K \cap \{y > 0\}) \Lambda^N(s) (\Delta_1^+(f) \mathbf{1}_{\Sigma_K}) ds \tag{21} \\
& + \nu \int_0^t \Lambda^N(s) (\Delta_{1,2}(f) p_1) ds + \mu \int_0^t \Lambda^N(s) (\Delta_{2,3}(f) p_2) ds \nonumber \\
& + \lambda \int_0^t \Lambda^N(s) (\Sigma_{<K}) \Lambda^N(s) (\Delta_{3,4}(f) \mathbf{1}_{\{y>0\}}) ds + \nu \int_0^t \Lambda^N(s) (\Delta_4^-(f) p_4) ds \nonumber
\end{align}
$$

where

$$
(22) \qquad \mathcal{M}_f^N(t) = \frac{1}{N} (\mathcal{M}_{f,1}^N(t) + \mathcal{M}_{f,2}^N(t) \dots + \mathcal{M}_{f,N}^N(t)),
$$

$$
(23) \qquad \Sigma_{<K} = \{(w, x, y, z) \in \mathbb{N}^4, w + x + y + z < K\}
$$

and

$$
(24) \quad p_i : \mathbb{N}^4 \to \mathbb{N} \text{ is the } i \text{-th projection (for example } p_1(w, x, y, z) = w).
$$
Let $\Sigma_K = \{(w, x, y, z) \in \mathbb{N}^4, w + x + y + z \le K\}$. For $f: \Sigma_K \to \mathbb{R}_+$, straightforwardly by equations (6)-(9),

$$ 
\begin{aligned}
df(Z_i^N(t)) = & \Delta_1^+(f)(Z_i^N(t)) \sum_{j=1}^N \mathbf{1}_{\{V_j^N(t^-)>0, S_i^N(t^-)<K\}} \mathcal{N}_{\lambda,j}^{U,N}(dt, \{i\}) \\
& + \Delta_{1,2}(f)(Z_i^N(t)) \sum_{j=1}^N \left( \sum_{l=1}^\infty \mathbf{1}_{\{l \le X_{j,i}^N(t^-)\}} \mathcal{N}_{\nu,j,i,l}(dt) \right) \\
& + \Delta_{2,3}(f)(Z_i^N(t)) \sum_{l=1}^\infty \mathbf{1}_{\{l \le R_i^N(t^-)\}} \mathcal{N}_{\mu,i,l}(dt) \\
& + \Delta_{3,4}(f)(Z_i^N(t)) \sum_{j=1}^N \mathbf{1}_{\{V_i^N(t^-)>0, S_j^N(t^-)<K\}} \mathcal{N}_{\lambda,i}^{U,N}(dt, \{j\}) \\
& + \Delta_4^-(f)(Z_i^N(t)) \sum_{j=1}^N \left( \sum_{l=1}^\infty \mathbf{1}_{\{l \le X_{i,j}^N(t^-)\}} \mathcal{N}_{\nu,i,j,l}(dt) \right)
\end{aligned}
$$

Thus, using the martingale decomposition for Poisson processes,

$$ 
\begin{aligned}
(20) \quad df(Z_i^N(t)) = & \Delta_1^+(f)(Z_i^N(t)) \mathbf{1}_{\{S_i^N(t)<K\}} \sum_{j=1}^N \mathbf{1}_{\{V_j^N(t)>0\}} \frac{\lambda}{N} dt \\
& + \Delta_{1,2}(f)(Z_i^N(t)) R_i^{r,N}(t) \nu dt + \Delta_{2,3}(f)(Z_i^N(t)) R_i^N(t) \mu dt \\
& + \Delta_{3,4}(f)(Z_i^N(t)) \mathbf{1}_{\{V_i^N(t)>0\}} \sum_{j=1}^N \mathbf{1}_{\{S_j^N(t)<K\}} \frac{\lambda}{N} dt \\
& + \Delta_4^-(f)(Z_i^N(t)) V_i^{r,N}(t) \nu dt + d\mathcal{M}_{f,i}^N(t)
\end{aligned}
$$

where $(\mathcal{M}_{f,i}^N(t))$ is a martingale which will be detailed in Section 4.2. By summing equation (20) for $i$ from 1 to $N$ and dividing by $N$,

$$ \Lambda^N(t)(f) = \Lambda^N(0)(f) + \mathcal{M}_f^N(t) $$

$$ + \lambda \int_0^t \Lambda^N(s) (\Sigma_K \cap \{y > 0\}) \Lambda^N(s) (\Delta_1^+(f) \mathbf{1}_{\Sigma_K}) ds $$

$$ (21) \quad + \nu \int_0^t \Lambda^N(s) (\Delta_{1,2}(f) p_1) ds + \mu \int_0^t \Lambda^N(s) (\Delta_{2,3}(f) p_2) ds $$

$$ + \lambda \int_0^t \Lambda^N(s) (\Sigma_{<K}) \Lambda^N(s) (\Delta_{3,4}(f) \mathbf{1}_{\{y>0\}}) ds + \nu \int_0^t \Lambda^N(s) (\Delta_4^-(f) p_4) ds $$

where

$$ (22) \qquad \mathcal{M}_f^N(t) = \frac{1}{N} (\mathcal{M}_{f,1}^N(t) + \mathcal{M}_{f,2}^N(t) \dots + \mathcal{M}_{f,N}^N(t)), $$

$$ (23) \qquad \Sigma_{<K} = \{(w, x, y, z) \in \mathbb{N}^4, w + x + y + z < K\} $$

and

$$ (24) \qquad p_i: \mathbb{N}^4 \to \mathbb{N} \text{ is the } i\text{-th projection (for example } p_1(w, x, y, z) = w). $$
4.2. **The martingale term.** The martingale term is $M_f^N(t)$ given by equation (22) where, for $1 \le i \le N$,

$$dM_{f,i}^N(t) = \Delta_1^+(f)(Z_i^N(t)) \sum_{j=1}^N \mathbf{1}_{\{V_j^N(t^-)>0, S_i^N(t^-)<K\}} (\mathcal{N}_{\lambda,j}^{U,N}(dt, \{i\}) - \frac{\lambda}{N}dt) \\ + \Delta_{1,2}(f)(Z_i^N(t)) \sum_{j=1}^N \sum_{l=1}^\infty \mathbf{1}_{\{l \le X_{j,i}^N(t^-)\}} (\mathcal{N}_{\nu,j,i,l}(dt) - \nu dt) \\ + \Delta_{2,3}(f)(Z_i^N(t)) \sum_{l=1}^\infty \mathbf{1}_{\{l \le R_i^N(t^-)\}} (\mathcal{N}_{\mu,i,l}(dt) - \mu dt) \\ + \Delta_{3,4}(f)(Z_i^N(t)) \sum_{j=1}^N \mathbf{1}_{\{V_j^N(t^-)>0, S_j^N(t^-)<K\}} (\mathcal{N}_{\lambda,i}^{U,N}(dt, \{j\}) - \frac{\lambda}{N}dt) \\ + \Delta_4^-(f)(Z_i^N(t)) \sum_{j=1}^N \sum_{l=1}^\infty \mathbf{1}_{\{l \le X_{i,j}^N(t^-)\}} (\mathcal{N}_{\nu,i,j,l}(dt) - \nu dt)$$

whose increasing process is expressed as follows

$$\langle M_f^N \rangle (t) = \frac{1}{N^2} (\lambda I_1^N(t) + \mu I_2^N(t) + \nu I_3^N(t))$$

where, from careful calculations,

$$I_1^N(t) = \sum_{i=1}^{N} \int_{0}^{t} \left( \begin{aligned} & (\Delta_1^+(f)(Z_i^N(s)))^2 \mathbf{1}_{\{S_i^N(s)<K\}} \Lambda^N(s) (\Sigma_K \cap \{y>0\}) \\ & + (\Delta_{3,4}(f)(Z_i^N(s)))^2 \mathbf{1}_{\{V_i^N(s)>0\}} \Lambda^N(s) (\Sigma_{<K}) \\ & + \frac{2}{N} \Delta_1^+(f)(Z_i^N(s)) \Delta_{3,4}(f)(Z_i^N(s)) \mathbf{1}_{\{V_i^N(s)>0, S_i^N(s)<K\}} \end{aligned} \right) ds$$

$$I_2^N(t) = \sum_{i=1}^{N} \int_{0}^{t} (\Delta_{2,3}(f)(Z_i^N(s)))^2 R_i^N(s) ds$$

$$I_3^N(t) = \sum_{i=1}^{N} \int_{0}^{t} \left( \begin{aligned} & (\Delta_{1,2}(f)(Z_i^N(s)))^2 R_i^{r,N}(s) + (\Delta_4^-(f)(Z_i^N(s)))^2 V_i^{r,N}(s) \\ & + \frac{2}{N} \Delta_{1,2}(f)(Z_i^N(s)) \Delta_4^-(f)(Z_i^N(s)) X_{i,i}(s) \end{aligned} \right) ds.$$

The term with $X_{i,i}^N(s)$ comes from the fact that only sequences $(\mathcal{N}_{\nu,i,j,.,j \neq i})$ and $(\mathcal{N}_{\nu,j,i,.,j \neq i})$ are independent. See Remark 1. It yields then straightforwardly that there exist $C_0, C_1 > 0$ such that, for $1 \le i \le 3$,

$$\|I_i^N\|_{\infty,T} \le (C_0 N + C_1) T \|f\|_{\infty}^2.$$

Thus, for

$$(25) \qquad \|\langle M_f^N \rangle\|_{\infty,T} \le (\lambda + \mu + \nu) \left(\frac{C_0}{N} + \frac{C_1}{N^2}\right) T \|f\|_{\infty}^2.$$
4.2. The martingale term. The martingale term is $M_f^N(t)$ given by equation (22) where, for $1 \le i \le N$,

$$
\begin{aligned}
d\mathcal{M}_{f,i}^N(t) = & \Delta_1^+(f)(Z_i^N(t)) \sum_{j=1}^N \mathbf{1}_{\{V_j^N(t^-)>0, S_i^N(t^-)<K\}} (\mathcal{N}_{\lambda,j}^{U,N}(dt, \{i\}) - \frac{\lambda}{N}dt) \\
& + \Delta_{1,2}(f)(Z_i^N(t)) \sum_{j=1}^N \sum_{l=1}^\infty \mathbf{1}_{\{l \le X_{j,i}^N(t^-)\}} (\mathcal{N}_{\nu,j,i,l}(dt) - \nu dt) \\
& + \Delta_{2,3}(f)(Z_i^N(t)) \sum_{l=1}^\infty \mathbf{1}_{\{l \le R_i^N(t^-)\}} (\mathcal{N}_{\mu,i,l}(dt) - \mu dt) \\
& + \Delta_{3,4}(f)(Z_i^N(t)) \sum_{j=1}^N \mathbf{1}_{\{V_j^N(t^-)>0, S_j^N(t^-)<K\}} (\mathcal{N}_{\lambda,i}^{U,N}(dt, \{j\}) - \frac{\lambda}{N}dt) \\
& + \Delta_4^-(f)(Z_i^N(t)) \sum_{j=1}^N \sum_{l=1}^\infty \mathbf{1}_{\{l \le X_{i,j}^N(t^-)\}} (\mathcal{N}_{\nu,i,j,l}(dt) - \nu dt)
\end{aligned}
$$

whose increasing process is expressed as follows

$$
\langle \mathcal{M}_f^N \rangle (t) = \frac{1}{N^2} (\lambda I_1^N(t) + \mu I_2^N(t) + \nu I_3^N(t))
$$

where, from careful calculations,

$$
I_1^N(t) = \sum_{i=1}^{N} \int_{0}^{t} \left( \begin{aligned}[t]
& (\Delta_1^+(f)(Z_i^N(s)))^2 \mathbf{1}_{\{S_i^N(s)<K\}} \Lambda^N(s) (\Sigma_K \cap \{y>0\}) \\
& + (\Delta_{3,4}(f)(Z_i^N(s)))^2 \mathbf{1}_{\{V_i^N(s)>0\}} \Lambda^N(s) (\Sigma_{<K}) \\
& + \frac{2}{N} \Delta_1^+(f)(Z_i^N(s)) \Delta_{3,4}(f)(Z_i^N(s)) \mathbf{1}_{\{V_i^N(s)>0, S_i^N(s)<K\}}
\end{aligned} \right) ds
$$

$$
I_2^N(t) = \sum_{i=1}^{N} \int_{0}^{t} \left( \Delta_{2,3}(f)(Z_i^N(s)) \right)^2 R_i^N(s) ds
$$

$$
I_3^N(t) = \sum_{i=1}^{N} \int_{0}^{t} \left( \begin{aligned}
& (\Delta_{1,2}(f)(Z_i^N(s)))^2 R_i^{r,N}(s) + (\Delta_4^-(f)(Z_i^N(s)))^2 V_i^{r,N}(s) \\
& + \frac{2}{N} \Delta_{1,2}(f)(Z_i^N(s)) \Delta_4^-(f)(Z_i^N(s)) X_{i,i}(s)
\end{aligned} \right) ds.
$$

The term with $X_{i,i}^N(s)$ comes from the fact that only sequences $(\mathcal{N}_{\nu,i,j,.,j \neq i})$ and $(\mathcal{N}_{\nu,j,i,.,j \neq i})$ are independent. See Remark 1. It yields then straightforwardly that there exist $C_0$, $C_1 > 0$ such that, for $1 \le i \le 3$,

$$
\| I_i^N \|_{\infty, T} \le (C_0 N + C_1) T \| f \|_{\infty}^2.
$$

Thus, for

$$
(25) \qquad \| \langle \mathcal{M}_f^N \rangle \|_{\infty, T} \le (\lambda + \mu + \nu) \left( \frac{C_0}{N} + \frac{C_1}{N^2} \right) T \|f\|_{\infty}^2.
$$
Thus, applying Cauchy-Schwarz then Doob's inequalities,

$$ 
\begin{aligned}
\left(\mathbb{E}\left(\sup_{0 \le s \le T} |\mathcal{M}_f^N(s)|\right)\right)^2 &\le \mathbb{E}\left(\sup_{0 \le s \le T} |\mathcal{M}_f^N(s)|^2\right) \\
&\le 4\mathbb{E}((\mathcal{M}_f^N)^2(T)) \\
&= 4\mathbb{E}(\langle \mathcal{M}_f^N \rangle(T)),
\end{aligned}
$$

and the martingale $(\mathcal{M}_f^N(t))$ converges in distribution to 0 when $N$ tends to $\infty$.

**4.3. Convergence of the empirical measure process.** This section is devoted to the proof of the mean-field convergence theorem (Theorem 2). The proof is based on tightness and uniqueness standard arguments using stochastic calculus and martingale theory. To be self-contained, the paper presents the detailed proof, via Propositions 1 and 2.

**Proposition 1 (Tightness of the empirical measure process).** The sequence $(\Lambda^N(t))$ is tight with respect to the convergence in distribution in $\mathcal{D}(\mathbb{R}_+, \mathcal{Y}^N)$, $\mathcal{Y}^N$ defined by equation (19). Any limiting point $(\Lambda(t))$ is a continuous process with values in

$$ (26) \qquad \mathcal{Y} = \{\Lambda \in \mathcal{P}(\Sigma_K), \sum_{(w,x,y,z) \in \Sigma_K} (x+y+z)\Lambda_{(w,x,y,z)} = s\}, $$

solution of

$$ (27) \qquad \Lambda(t)(f) = \Lambda(0)(f) + \lambda \int_0^t \Lambda(s)(\Sigma_K \cap \{y > 0\})\Lambda(s)(\Delta_1^+(f)\mathbf{1}_{\Sigma_{<K}})ds \\ 
+ \nu \int_0^t \Lambda(s)(\Delta_{1,2}(f)p_1)ds + \mu \int_0^t \Lambda(s)(\Delta_{2,3}(f)p_2)ds \\ 
+ \lambda \int_0^t \Lambda(s)(\Sigma_{<K})\Lambda(s)(\Delta_{3,4}(f)\mathbf{1}_{\{y>0\}})ds + \nu \int_0^t \Lambda(s)(\Delta_4^-(f)p_4)ds $$

for any function $f$ on $\Sigma_K = \{(w, x, y, z) \in \mathbb{N}^4, w+x+y+z \le K\}$ and $\Sigma_{<K}$ and the $p_i$'s defined by equations (23) and (24).

**Proof.** This amounts to proving that, for any function $f$ on $\Sigma_K$, the sequence of processes $(\Lambda^N(t)(f))$ is tight with respect to the topology of the uniform norm on compact sets. For this, using the modulus of continuity criterion (see Billinsley [3]), it suffices to prove that, for $T$, $\varepsilon$, $\eta > 0$, there exist $\delta_0 > 0$ and $N_0 \in \mathbb{N}$ such that, for all $\delta < \delta_0$ and all $N \ge N_0$,

$$ (28) \qquad \mathbb{P} \left( \sup_{\substack{0 \le s \le t \le T \\ |s-t| < \delta}} |\Lambda^N(t)(f) - \Lambda^N(s)(f)| > \eta \right) < \varepsilon. $$

Let $T > 0$, $\varepsilon > 0$, $\eta > 0$, $\delta > 0$ and $s, t \in [0, T]$ such that $|s-t| < \delta$ be fixed. For the third term of the right-hand side of equation (21), there exists $C_2 > 0$ such that

$$ \left| \int_s^t \Lambda^N(u)(\Sigma_K \cap \{y > 0\}) \Lambda^N(u)(\Delta_1^+(f) \mathbf{1}_{\{\Sigma_{<K}\}}) du \right| \le \delta C_2 \|f\|_{\infty} $$

and the same holds for the other terms. Thus there exists $C_3 > 0$ such that

$$ |\Lambda^N(t)(f) - \Lambda^N(s)(f)| \le \delta C_3 \|f\|_{\infty} + |\mathcal{M}_f^N(t) - \mathcal{M}_f^N(s)|. $$

Using equation (25) for the martingale term, it holds that there exists $C_4 > 0$ such that

$$ \mathbb{E} \left( \sup_{\substack{0 \le s \le t \le T \\ |s-t| < \delta}} |\Lambda^N(t)(f) - \Lambda^N(s)(f)| \right) \le \delta C_4 \|f\|_{\infty} + 2 \mathbb{E} \left( \sup_{0 \le t \le T} |\mathcal{M}_f^N(t)| \right). $$
Thus, applying Cauchy-Schwarz then Doob's inequalities,

$$ \left( \mathbb{E} \left( \sup_{0 \le s \le T} |\mathcal{M}_f^N(s)| \right) \right)^2 \le \mathbb{E} \left( \sup_{0 \le s \le T} |\mathcal{M}_f^N(s)|^2 \right) \\ \le 4 \mathbb{E} \left( (\mathcal{M}_f^N)^2(T) \right) \\ = 4 \mathbb{E} \left( \langle \mathcal{M}_f^N \rangle(T) \right), $$

and the martingale $(\mathcal{M}_f^N(t))$ converges in distribution to 0 when $N$ tends to $\infty$.

**4.3. Convergence of the empirical measure process.** This section is devoted to the proof of the mean-field convergence theorem (Theorem 2). The proof is based on tightness and uniqueness standard arguments using stochastic calculus and martingale theory. To be self-contained, the paper presents the detailed proof, via Propositions 1 and 2.

**Proposition 1 (Tightness of the empirical measure process).** The sequence $(\Lambda^N(t))$ is tight with respect to the convergence in distribution in $\mathcal{D}(\mathbb{R}_+, \mathcal{Y}^N)$, $\mathcal{Y}^N$ defined by equation (19). Any limiting point $(\Lambda(t))$ is a continuous process with values in

$$ (26) \qquad \mathcal{Y} = \{ \Lambda \in \mathcal{P}(\Sigma_K), \sum_{(w,x,y,z) \in \Sigma_K} (x+y+z) \Lambda_{(w,x,y,z)} = s \}, $$

solution of

$$ \Lambda(t)(f) = \Lambda(0)(f) + \lambda \int_0^t \Lambda(s)(\Sigma_K \cap \{y > 0\})\Lambda(s)(\Delta_1^+(f)\mathbf{1}_{\Sigma_{<K}})ds \\ + \nu \int_0^t \Lambda(s)(\Delta_{1,2}(f)p_1)ds + \mu \int_0^t \Lambda(s)(\Delta_{2,3}(f)p_2)ds \\ (27) \qquad + \lambda \int_0^t \Lambda(s)(\Sigma_{<K})\Lambda(s)(\Delta_{3,4}(f)\mathbf{1}_{\{y>0\}})ds + \nu \int_0^t \Lambda(s)(\Delta_4^-(f)p_4)ds $$

for any function $f$ on $\Sigma_K = \{(w, x, y, z) \in \mathbb{N}^4, w + x + y + z \le K\}$ and $\Sigma_{<K}$ and the $p_i$'s defined by equations (23) and (24).

**Proof.** This amounts to proving that, for any function $f$ on $\Sigma_K$, the sequence of processes $(\Lambda^N(t)(f))$ is tight with respect to the topology of the uniform norm on compact sets. For this, using the modulus of continuity criterion (see Billinsley [3]), it suffices to prove that, for $T$, $\varepsilon$, $\eta > 0$, there exist $\delta_0 > 0$ and $N_0 \in \mathbb{N}$ such that, for all $\delta < \delta_0$ and all $N \ge N_0$,

$$ (28) \qquad \mathbb{P} \left( \sup_{\substack{0 \le s \le t \le T \\ |s-t| < \delta}} |\Lambda^N(t)(f) - \Lambda^N(s)(f)| > \eta \right) < \varepsilon. $$

Let $T > 0$, $\varepsilon > 0$, $\eta > 0$, $\delta > 0$ and $s, t \in [0, T]$ such that $|s - t| < \delta$ be fixed. For the third term of the right-hand side of equation (21), there exists $C_2 > 0$ such that

$$ \left| \int_{s}^{t} \Lambda^{N}(u)(\Sigma_{K} \cap \{y > 0\}) \Lambda^{N}(u)(\Delta_{1}^{+}(f) \mathbf{1}_{\{\Sigma_{<K}\}}) du \right| \le \delta C_{2} \|f\|_{\infty} $$

and the same holds for the other terms. Thus there exists $C_3 > 0$ such that

$$ |\Lambda^N(t)(f) - \Lambda^N(s)(f)| \le \delta C_3 \|f\|_{\infty} + |\mathcal{M}_f^N(t) - \mathcal{M}_f^N(s)|. $$

Using equation (25) for the martingale term, it holds that there exists $C_4 > 0$ such that

$$ \mathbb{E} \left( \sup_{\substack{0 \le s \le t \le T \\ |s-t| < \delta}} |\Lambda^N(t)(f) - \Lambda^N(s)(f)| \right) \le \delta C_4 \|f\|_{\infty} + 2 \mathbb{E} \left( \sup_{0 \le t \le T} |\mathcal{M}_f^N(t)| \right). $$
Thus, as the martingale term converges in distribution to 0, there exist $\delta_0 > 0$ and $N_0 \in \mathbb{N}$ such that, for all $\delta < \delta_0$ and all $N \ge N_0$, the left-hand side of the previous equation is less than $\varepsilon$. Then, using Markov's inequality, the sequence of processes $(\Lambda^N(t)(f))$ satisfies equation (28) thus is tight. Therefore, if $\Lambda$ is a limiting point of $\Lambda^N$, using again equation (21), as $(M_f^N(t))$ converges in distribution to 0, equation (27) holds. As function $f$ has finite support, all the terms of the right-hand side are straightforwardly continuous on $t$. $\square$

The following proposition gives the uniqueness of a limiting point of $(\Lambda^N(t))$.

**Proposition 2 (Uniqueness).** For every probability $\Lambda_0$ on $\Sigma_K$, equation (27) has at most one solution $(\Lambda(t))$ in $\mathcal{D}(\mathbb{R}_+, \mathcal{P}(\Sigma_K))$ with initial condition $\Lambda_0$.

*Proof.* Let $\Lambda^1(t))$ and $\Lambda^2(t))$ in $\mathcal{D}(\mathbb{R}_+, \mathcal{P}(\Sigma_K))$ be two solutions of equation (27) with initial condition $\Lambda_0$. For $f$ function on $\Sigma_K$ and $t \ge 0$,

$$ 
\begin{aligned}
& \Lambda^1(t)(f) - \Lambda^2(t)(f) \\
& = \lambda \int_0^t (\Lambda^1(s) - \Lambda^2(s))(\Sigma_K \cap \{y > 0\})\Lambda^1(s)(\Delta_1^+(f)\mathbf{1}_{\Sigma_{<K}})ds \\
& + \lambda \int_0^t \Lambda^2(s)(\Sigma_K \cap \{y > 0\})(\Lambda^1(s) - \Lambda^2(s))(\Delta_1^+(f)\mathbf{1}_{\Sigma_{<K}})ds \\
& + \nu \int_0^t (\Lambda^1(s) - \Lambda^2(s))(\Delta_{1,2}(f)p_1 + \Delta_4^-(f)p_4)ds \\
& + \mu \int_0^t (\Lambda^1(s) - \Lambda^2(s))(\Delta_{2,3}(f)p_2)ds \\
& + \lambda \int_0^t (\Lambda^1(s) - \Lambda^2(s))(\Sigma_{<K})\Lambda^1(s)(\Delta_{3,4}(f)\mathbf{1}_{\{y>0\}})ds \\
& + \lambda \int_0^t \Lambda^2(s)(\Sigma_{<K})(\Lambda^1(s) - \Lambda^2(s))(\Delta_{3,4}(f)\mathbf{1}_{\{y>0\}})ds.
\end{aligned}
$$

Recall that, for a signed measure $\pi$ on $\Sigma_K$,

$$ \|\pi\|_{TV} = \sup\{|\pi(f)|, f : \Sigma_K \to \mathbb{R}, \|f\|_{\infty} \le 1\}. $$

From the previous equation, it holds that

$$ \|\Lambda^1(t) - \Lambda^2(t)\|_{TV} \le (8\lambda + 4\nu + 2\mu) \int_0^t \|\Lambda^1(s) - \Lambda^2(s)\|_{TV} ds. $$

Applying Grönwall's lemma,

$$ \|\Lambda^1(t) - \Lambda^2(t)\|_{TV} = 0. $$

It ends the proof.

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

**Proof of Theorem 2.** Let $(w, x, y, z) \in \Sigma_K$ and $\Lambda_0 = \delta_{(w,x,y,z)}$. If $(\bar{Z}(t))$ is the unique solution of equation (15) and the measure valued process $(\Lambda(t))$ is defined, for $f$ function on $\Sigma_K$, by

$$ \Lambda(t)(f) = \mathbb{E}(f(\bar{Z}(t))) $$

then it is easy to check that $(\Lambda(t))$ is a solution of equation (27). The convergence of $(\Lambda^N(t))$ follows from Propositions 1 and 2. See Sznitman [21, Proposition 2.2] for the propagation of chaos property.

**4.4. Probabilistic interpretation of the asymptotic process.** Note that the Fokker-Planck equation (1) is the functional form of the stochastic equation (27). Recall that, in equation (1), equalities in distribution hold, as

$$ \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \overline{R^r}(s^-)\}} \overline{N}_\nu(ds, dh) = \int_0^t \sum_{l=0}^\infty \mathbf{1}_{\{l \le \overline{R^r}(s^-)\}} \mathcal{N}_{\nu,l}(ds). $$
Thus, as the martingale term converges in distribution to 0, there exist $\delta_0 > 0$ and $N_0 \in \mathbb{N}$ such that, for all $\delta < \delta_0$ and all $N \ge N_0$, the left-hand side of the previous equation is less than $\varepsilon$. Then, using Markov's inequality, the sequence of processes $(\Lambda^N(t)(f))$ satisfies equation (28) thus is tight. Therefore, if $\Lambda$ is a limiting point of $\Lambda^N$, using again equation (21), as $(M_f^N(t))$ converges in distribution to 0, equation (27) holds. As function $f$ has finite support, all the terms of the right-hand side are straightforwardly continuous on $t$. $\square$

The following proposition gives the uniqueness of a limiting point of $(\Lambda^N(t))$.

**Proposition 2 (Uniqueness).** For every probability $\Lambda_0$ on $\Sigma_K$, equation (27) has at most one solution $(\Lambda(t))$ in $D(\mathbb{R}_+, P(\Sigma_K))$ with initial condition $\Lambda_0$.

*Proof.* Let $\Lambda^1(t))$ and $\Lambda^2(t))$ in $D(\mathbb{R}_+, P(\Sigma_K))$ be two solutions of equation (27) with initial condition $\Lambda_0$. For $f$ function on $\Sigma_K$ and $t \ge 0$,

$$ \begin{aligned} & \Lambda^1(t)(f) - \Lambda^2(t)(f) \\ & = \lambda \int_0^t (\Lambda^1(s) - \Lambda^2(s))(\Sigma_K \cap \{y > 0\})\Lambda^1(s)(\Delta_1^+(f)\mathbf{1}_{\Sigma_{<K}})ds \\ & + \lambda \int_0^t \Lambda^2(s)(\Sigma_K \cap \{y > 0\})(\Lambda^1(s) - \Lambda^2(s))(\Delta_1^+(f)\mathbf{1}_{\Sigma_{<K}})ds \\ & + \nu \int_0^t (\Lambda^1(s) - \Lambda^2(s))(\Delta_{1,2}(f)p_1 + \Delta_4^-(f)p_4)ds + \mu \int_0^t (\Lambda^1(s) - \Lambda^2(s))(\Delta_{2,3}(f)p_2)ds \\ & + \lambda \int_0^t (\Lambda^1(s) - \Lambda^2(s))(\Sigma_{<K})\Lambda^1(s)(\Delta_{3,4}(f)\mathbf{1}_{\{y>0\}})ds \\ & + \lambda \int_0^t \Lambda^2(s)(\Sigma_{<K})(\Lambda^1(s) - \Lambda^2(s))(\Delta_{3,4}(f)\mathbf{1}_{\{y>0\}})ds. \end{aligned} $$

Recall that, for a signed measure $\pi$ on $\Sigma_K$,

$$ \|\pi\|_{TV} = \sup\{|\pi(f)|, f : \Sigma_K \to \mathbb{R}, \|f\|_{\infty} \le 1\}. $$

From the previous equation, it holds that

$$ \|\Lambda^1(t) - \Lambda^2(t)\|_{TV} \le (8\lambda + 4\nu + 2\mu) \int_0^t \|\Lambda^1(s) - \Lambda^2(s)\|_{TV} ds. $$

Applying Grönwall's lemma,

$$ \|\Lambda^1(t) - \Lambda^2(t)\|_{TV} = 0. $$

It ends the proof. $\square$

**Proof of Theorem 2.** Let $(w, x, y, z) \in \Sigma_K$ and $\Lambda_0 = \delta_{(w,x,y,z)}$. If $(\bar{Z}(t))$ is the unique solution of equation (15) and the measure valued process $(\Lambda(t))$ is defined, for $f$ function on $\Sigma_K$, by

$$ \Lambda(t)(f) = \mathbb{E}(f(\bar{Z}(t))) $$

then it is easy to check that $(\Lambda(t))$ is a solution of equation (27). The convergence of $(\Lambda^N(t))$ follows from Propositions 1 and 2. See Sznitman [21, Proposition 2.2] for the propagation of chaos property.

**4.4. Probabilistic interpretation of the asymptotic process.** Note that the Fokker-Planck equation (1) is the functional form of the stochastic equation (27). Recall that, in equation (1), equalities in distribution hold, as

$$ \iint_{[0,t] \times \mathbb{R}_+} \mathbf{1}_{\{0 \le h \le \overline{R^r}(s^-)\}} \overline{N}_\nu(ds, dh) = \int_0^t \sum_{l=0}^\infty \mathbf{1}_{\{l \le \overline{R^r}(s^-)\}} \mathcal{N}_{\nu,l}(ds). $$
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Dynamics of ($\bar{Z}(t)$) as a tandem of four queues. The *vertical* queues are $M/M/\infty$ queues, the *horizontal* is a $M/M/1$ queue. The overall capacity is $K$.</figcaption></figure>

Thus the non-homogeneous Markov process ($\bar{Z}(t)$) can be seen as the state process of four queues in tandem, with overall capacity $K$. This means that $\bar{R}^r(t)$, $\bar{R}(t)$, $\bar{V}(t)$ and $\bar{V}^r(t)$ are the numbers of customers in respectively

- the first queue, an infinite-server queue with service rate $\nu$,

- the second one, an infinite-server queue with service rate $\mu$,

- the third one, a one-server queue with variable service rate $\lambda \mathbb{P}(\bar{S}(t) < K)$ at time $t$

- and the last one, an infinite-server queue with service rate $\nu$,

while the arrival process is a non-homogeneous Poisson process with intensity $\lambda \mathbb{P}(\bar{V}(t) > 0)dt$.

Let $\eta_1, \rho_1, \rho_2$ and $\eta_2$ be the arrival-to-service rate ratios for the four queues from the left to the right (as shown on Figure 1). By definition, one gets

$$ (29) \qquad \begin{cases} \eta_1(t) = \frac{\lambda}{\nu} \mathbb{P}(\bar{V}(t) > 0), & \rho_1(t) = \frac{\lambda}{\mu} \mathbb{P}(\bar{V}(t) > 0), \\ \rho_2(t) = \frac{\mathbb{P}(\bar{V}(t) > 0)}{\mathbb{P}(\bar{S}(t) < K)}, & \eta_2(t) = \frac{\lambda}{\nu} \mathbb{P}(\bar{V}(t) > 0). \end{cases} $$

## 5. EQUILIBRIUM OF THE ASYMPTOTIC PROCESS

The quadruplet of the numbers of customers in such four queues in tandem with fixed arrival-to-service rate ratios $\eta_1, \rho_1, \rho_2$ and $\eta_2$ and finite overall capacity $K$ is an ergodic Markov process as an irreducible Markov process on a finite state space. Moreover, the unique invariant probability measure has a well-known product form given by

$$ (30) \qquad \pi_{j,k,l,m}(\rho) = \frac{1}{Z(\rho)} \frac{\eta_1^j \rho_1^k}{j!k!} \rho_2^l \frac{\eta_2^m}{m!} $$

where, to shorten the notations, $(\eta_1, \rho_1, \rho_2, \eta_2)$ is denoted by $\rho$ and the normalizing constant is

$$ (31) \qquad Z(\rho) = \sum_{j+k+l+m \le K} \frac{\eta_1^j \rho_1^k}{j!k!} \rho_2^l \frac{\eta_2^m}{m!}. $$

If the process ($\bar{Z}(t)$) of the number of customers in the four queues in tandem is at equilibrium then, denoting by $L_{\rho(t)}$ its generator, it holds that

$$ 0 = \pi(\rho(t))L_{\rho(t)} $$
<figure><img src="image_3.png" /><figcaption>FIGURE 1. Dynamics of $(\bar{Z}(t))$ as a tandem of four queues. The *vertical* queues are $M/M/\infty$ queues, the *horizontal* is a $M/M/1$ queue. The overall capacity is $K$.</figcaption></figure>

Thus the non-homogeneous Markov process $(\bar{Z}(t))$ can be seen as the state process of four queues in tandem, with overall capacity $K$. This means that $\bar{R}^r(t)$, $\bar{R}(t)$, $\bar{V}(t)$ and $\bar{V}^r(t)$ are the numbers of customers in respectively

- the first queue, an infinite-server queue with service rate $\nu$,

- the second one, an infinite-server queue with service rate $\mu$,

- the third one, a one-server queue with variable service rate $\lambda \mathbb{P}(\bar{S}(t) < K)$ at time $t$

- and the last one, an infinite-server queue with service rate $\nu$,

while the arrival process is a non-homogeneous Poisson process with intensity $\lambda \mathbb{P}(\bar{V}(t) > 0)dt$.

Let $\eta_1, \rho_1, \rho_2$ and $\eta_2$ be the arrival-to-service rate ratios for the four queues from the left to the right (as shown on Figure 1). By definition, one gets

$$ (29) \qquad \begin{cases} \eta_1(t) = \frac{\lambda}{\nu} \mathbb{P}(\bar{V}(t) > 0), & \rho_1(t) = \frac{\lambda}{\mu} \mathbb{P}(\bar{V}(t) > 0), \\ \rho_2(t) = \frac{\mathbb{P}(\bar{V}(t) > 0)}{\mathbb{P}(\bar{S}(t) < K)}, & \eta_2(t) = \frac{\lambda}{\nu} \mathbb{P}(\bar{V}(t) > 0). \end{cases} $$

## 5. EQUILIBRIUM OF THE ASYMPTOTIC PROCESS

The quadruplet of the numbers of customers in such four queues in tandem with fixed arrival-to-service rate ratios $\eta_1, \rho_1, \rho_2$ and $\eta_2$ and finite overall capacity $K$ is an ergodic Markov process as an irreducible Markov process on a finite state space. Moreover, the unique invariant probability measure has a well-known product form given by

$$ (30) \qquad \pi_{j,k,l,m}(\rho) = \frac{1}{Z(\rho)} \frac{\eta_1^j \rho_1^k}{j!k!} \rho_2^l \frac{\eta_2^m}{m!} $$

where, to shorten the notations, $(\eta_1, \rho_1, \rho_2, \eta_2)$ is denoted by $\rho$ and the normalizing constant is

$$ (31) \qquad Z(\rho) = \sum_{j+k+l+m \le K} \frac{\eta_1^j \rho_1^k}{j!k!} \rho_2^l \frac{\eta_2^m}{m!}. $$

If the process $(\bar{Z}(t))$ of the number of customers in the four queues in tandem is at equilibrium then, denoting by $L_{\rho(t)}$ its generator, it holds that

$$ 0 = \pi(\rho(t)) L_{\rho(t)} $$
where $\rho(t) = (\eta_1, \rho_1, \rho_2, \eta_2)(t)$ is given by equation (29). It means that the equilibrium point is the probability measure $\pi(\rho)$ on $\mathcal{Y}$ defined by (26) where $\rho = (\eta_1, \rho_1, \rho_2, \eta_2)$ satisfies

$$ (32) \qquad \eta_1 = \frac{\lambda}{\nu}(1 - \pi_{0V}(\rho)), $$

$$ (33) \qquad \rho_1 = \frac{\lambda}{\mu}(1 - \pi_{0V}(\rho)), $$

$$ (34) \qquad \rho_2 = \frac{1 - \pi_{0V}(\rho)}{1 - \pi_S(\rho)}, $$

$$ (35) \qquad \eta_2 = \frac{\lambda}{\nu}(1 - \pi_{0V}(\rho)) $$

with $\pi_{0V}(\rho) = \sum_{j+l+m \le K} \pi_{j,0,l,m}(\rho)$ and $\pi_S(\rho) = \sum_{j+k+l+m=K} \pi_{j,k,l,m}(\rho)$. Viewing the tandem of four queues as a station, $\pi_{0V}(\rho)$ is the probability that there is no car available and $\pi_S(\rho)$ the probability that the station is saturated. Because $\pi(\rho)$ has support on $\mathcal{Y}$, it holds that

$$ (36) \qquad s = \sum_{j+k+l+m \le K} (k+l+m)\pi_{j,k,l,m}(\rho). $$

**Theorem 3** (Uniqueness and characterization of the equilibrium point). If $\nu$ is enough large then there exists a unique equilibrium point for the Fokker-Planck equation (27) which is $\pi(\rho)$ defined by equation (30) where $\rho = (\frac{\mu}{\nu}\rho_1, \rho_1, \rho_2, \frac{\mu}{\nu}\rho_1)$ and $(\rho_1, \rho_2)$ is the unique solution of

$$ (37) \qquad \rho_1 = \frac{\lambda}{\mu}(1 - \pi_{0V}(\rho)), $$

$$ (38) \qquad s = \sum_{j+k+l+m \le K} (k+l+m)\pi_{j,k,l,m}(\rho). $$

**Remark 2.** The uniqueness of the equilibrium point is the main issue. Moreover, its characterization given by Theorem 3 is a major contribution which allows to derive quantitative results. Proving the uniqueness of the equilibrium point in dimension 2 is quite rare in the literature. We can cite two papers. For the celebrated model by Gibbens, Hunt and Kelly [13], simulations highlight in this paper a range of parameters with two stable equilibrium points, called metastability phenomena. This fact was proved some 30 years later in [17]. In [2], the same question is solved for a migration-contagion model using a convexity argument. The arguments used in our paper are totally different.

**Remark 3.** The assumption that $\nu$ is enough large is a technical assumption for the proof. It seems that the result is true for all $\nu > 0$ but the proof is for the moment out of reach.

Let us prove Theorem 3 in five steps. Steps 2, 3 and 4 are devoted to the special case where $\nu$ tends to infinity. This case is called the simple reservation case. Indeed, when $\nu$ gets large, it turns out that the model corresponds to the case where the car is not reserved in advance, and the parking space is just reserved when the user takes the car. This model is studied in [6]. Step 1 is here to present the framework for the simple reservation case. Steps 2, 3 and 4 exhibit the proof of [6, Theorem 2] omitted in [6] for the simple reservation model.
where $\rho(t) = (\eta_1, \rho_1, \rho_2, \eta_2)(t)$ is given by equation (29). It means that the equilibrium point is the probability measure $\pi(\rho)$ on $\mathcal{Y}$ defined by (26) where $\rho = (\eta_1, \rho_1, \rho_2, \eta_2)$ satisfies

$$ (32) \qquad \eta_1 = \frac{\lambda}{\nu}(1 - \pi_{0V}(\rho)), $$

$$ (33) \qquad \rho_1 = \frac{\lambda}{\mu}(1 - \pi_{0V}(\rho)), $$

$$ (34) \qquad \rho_2 = \frac{1 - \pi_{0V}(\rho)}{1 - \pi_S(\rho)}, $$

$$ (35) \qquad \eta_2 = \frac{\lambda}{\nu}(1 - \pi_{0V}(\rho)) $$

with $\pi_{0V}(\rho) = \sum_{j+l+m \le K} \pi_{j,0,l,m}(\rho)$ and $\pi_S(\rho) = \sum_{j+k+l+m=K} \pi_{j,k,l,m}(\rho)$. Viewing the tandem of four queues as a station, $\pi_{0V}(\rho)$ is the probability that there is no car available and $\pi_S(\rho)$ the probability that the station is saturated. Because $\pi(\rho)$ has support on $\mathcal{Y}$, it holds that

$$ (36) \qquad s = \sum_{j+k+l+m \le K} (k+l+m)\pi_{j,k,l,m}(\rho). $$

**Theorem 3** (Uniqueness and characterization of the equilibrium point). If $\nu$ is enough large then there exists a unique equilibrium point for the Fokker-Planck equation (27) which is $\pi(\rho)$ defined by equation (30) where $\rho = (\frac{\mu}{\nu}\rho_1, \rho_1, \rho_2, \frac{\mu}{\nu}\rho_1)$ and $(\rho_1, \rho_2)$ is the unique solution of

$$ (37) \qquad \rho_1 = \frac{\lambda}{\mu}(1 - \pi_{0V}(\rho)), $$

$$ (38) \qquad s = \sum_{j+k+l+m \le K} (k+l+m)\pi_{j,k,l,m}(\rho). $$

**Remark 2.** The uniqueness of the equilibrium point is the main issue. Moreover, its characterization given by Theorem 3 is a major contribution which allows to derive quantitative results. Proving the uniqueness of the equilibrium point in dimension 2 is quite rare in the literature. We can cite two papers. For the celebrated model by Gibbens, Hunt and Kelly [13], simulations highlight in this paper a range of parameters with two stable equilibrium points, called metastability phenomena. This fact was proved some 30 years later in [17]. In [2], the same question is solved for a migration-contagion model using a convexity argument. The arguments used in our paper are totally different.

**Remark 3.** The assumption that $\nu$ is enough large is a technical assumption for the proof. It seems that the result is true for all $\nu > 0$ but the proof is for the moment out of reach.

Let us prove Theorem 3 in five steps. Steps 2, 3 and 4 are devoted to the special case where $\nu$ tends to infinity. This case is called the simple reservation case. Indeed, when $\nu$ gets large, it turns out that the model corresponds to the case where the car is not reserved in advance, and the parking space is just reserved when the user takes the car. This model is studied in [6]. Step 1 is here to present the framework for the simple reservation case. Steps 2, 3 and 4 exhibit the proof of [6, Theorem 2] omitted in [6] for the simple reservation model.
**STEP 1: THE SIMPLE RESERVATION CASE.** For the model with simple reservation, the problem of existence and uniqueness of an equilibrium point amounts to finding $(\rho_1, \rho_2)$ such that

$$ (39) \qquad \rho_1 = \frac{\lambda}{\mu}(1 - \pi_{.,0}(\rho_1, \rho_2)), $$

$$ (40) \qquad \rho_2 = \frac{1 - \pi_{.,0}(\rho_1, \rho_2)}{1 - \pi_S(\rho_1, \rho_2)}, $$

$$ (41) \qquad s = \sum_{i+j \le K} (i+j)\pi_{i,j}(\rho_1, \rho_2). $$

where the invariant probability measure is now

$$ \pi_{i,j}(\rho_1, \rho_2) = \frac{1}{Z(\rho_1, \rho_2)} \frac{\rho_1^i}{i!} \rho_2^j $$

with $Z(\rho_1, \rho_2) = \sum_{i+j \le K} \pi_{i,j}(\rho_1, \rho_2)$, $\pi_{.,0} = \sum_{i=0}^K \pi_{i,0}$ and $\pi_S = \sum_{i+j=K} \pi_{i,j}$.

**STEP 2: EQUATIONS (39)-(41) AMOUNT TO EQUATIONS (39) and (41).** Note that any $\rho$ solution of (39) and (40) lies in $\Gamma = [0, \lambda/\mu] \times \mathbb{R}^+$. For any $\rho \in \Gamma$, $\rho$ is solution of equation (40) because

$$ \rho_2(1 - \pi_S(\rho_1, \rho_2)) = \rho_2 \frac{1}{Z(\rho_1, \rho_2)} \sum_{i+j<K} \frac{\rho_1^i}{i!} \rho_2^j = \frac{1}{Z(\rho_1, \rho_2)} \sum_{j>0, i+j \le K} \frac{\rho_1^i}{i!} \rho_2^j = 1 - \pi_{.,0}(\rho_1, \rho_2). $$

**STEP 3: DIFFEOMORPHISM FROM EQUATION (39).**

**Lemma 1 (Diffeomorphism).** There exists a strictly increasing diffeomorphism $\phi : ]0, \lambda/\mu[ \to ]0, \infty[, \text{ and } \psi = \phi^{-1}, \text{ such that } (\rho_1, \rho_2) \text{ is a solution of equation (39) if and only if } \rho_2 = \phi(\rho_1).$

*Proof.* First, $(\rho_1, \rho_2) \in \Gamma$ is solution of equation (39) if and only if $f(\rho_1, \rho_2) = 0$ where $f$ is the $C^\infty$ function defined by

$$ (42) \qquad f(\rho_1, \rho_2) = \left(\frac{\lambda}{\mu} - \rho_1\right) Z(\rho_1, \rho_2) - \frac{\lambda}{\mu} \sum_{i=0}^{K} \frac{\rho_1^i}{i!}. $$

To prove the existence of a strictly increasing diffeomorphism $\phi$ which maps $\rho_1$ to $\rho_2$, we apply the global inverse function theorem to auxiliary function $h$ defined on $\Gamma$ by

$$ h(\rho_1, \rho_2) = (\rho_1, f(\rho_1, \rho_2)). $$

Indeed, $h$ is injective if and only if, for any $\rho_1 \in ]0, \lambda/\mu[ \text{ and } \rho_2, \rho'_2 \in ]0, +\infty[, f(\rho_1, \rho_2) = f(\rho_1, \rho'_2)$ implies that $\rho_2 = \rho'_2$. As $f$ is $C^1$ on $]0, \lambda/\mu[ \times ]0, +\infty[$, it is sufficient to prove that $\partial f / \partial \rho_2 \neq 0$ on $]0, \lambda/\mu[ \times ]0, +\infty[$ to obtain that $\rho_2 \mapsto f(\rho_1, \rho_2)$ is strictly monotone thus injective.

Note that it is easy to see that $\partial f / \partial \rho_2$ is positive on $]0, \lambda/\mu[ \times ]0, +\infty[$. Indeed, since $\rho_1 < \lambda/\mu$ and $\rho_2 \mapsto Z(\rho_1, \rho_2)$ is non decreasing,

$$ \frac{\partial f}{\partial \rho_2}(\rho_1, \rho_2) = \left(\frac{\lambda}{\mu} - \rho_1\right) \frac{\partial Z}{\partial \rho_2}(\rho_1, \rho_2) > 0. $$

As a consequence, $h$ is injective and $C^1$ on $]0, \lambda/\mu[ \times ]0, +\infty[$. By global inversion function theorem, $h^{-1}$ is $C^1$ on $h([0, \lambda/\mu[ \times ]0, +\infty[$) and it exists $\phi: C^1$ defined on $]0, \lambda/\mu[$ by

$$ h^{-1}(\rho_1, 0) = (\rho_1, \phi(\rho_1)). $$

Thus, to prove that diffeomorphism $\phi$ is strictly increasing on $]0, \lambda/\mu[$, it amounts to showing that, for all $\rho_1 \in ]0, \lambda/\mu[$, one gets

$$ \frac{\partial f}{\partial \rho_1}(\rho_1, \phi(\rho_1)) < 0, $$
**STEP 1: THE SIMPLE RESERVATION CASE.** For the model with simple reservation, the problem of existence and uniqueness of an equilibrium point amounts to finding $(\rho_1, \rho_2)$ such that

$$ (39) \qquad \rho_1 = \frac{\lambda}{\mu}(1 - \pi_{.,0}(\rho_1, \rho_2)), $$

$$ (40) \qquad \rho_2 = \frac{1 - \pi_{.,0}(\rho_1, \rho_2)}{1 - \pi_S(\rho_1, \rho_2)}, $$

$$ (41) \qquad s = \sum_{i+j \le K} (i+j)\pi_{i,j}(\rho_1, \rho_2). $$

where the invariant probability measure is now

$$ \pi_{i,j}(\rho_1, \rho_2) = \frac{1}{Z(\rho_1, \rho_2)} \frac{\rho_1^i}{i!} \rho_2^j $$

with $Z(\rho_1, \rho_2) = \sum_{i+j \le K} \pi_{i,j}(\rho_1, \rho_2)$, $\pi_{.,0} = \sum_{i=0}^K \pi_{i,0}$ and $\pi_S = \sum_{i+j=K} \pi_{i,j}$.

**STEP 2: EQUATIONS (39)-(41) AMOUNT TO EQUATIONS (39) and (41).** Note that any $\rho$ solution of (39) and (40) lies in $\Gamma = [0, \lambda/\mu] \times \mathbb{R}^+$. For any $\rho \in \Gamma$, $\rho$ is solution of equation (40) because

$$ \rho_2(1 - \pi_S(\rho_1, \rho_2)) = \rho_2 \frac{1}{Z(\rho_1, \rho_2)} \sum_{i+j<K} \frac{\rho_1^i}{i!} \rho_2^j = \frac{1}{Z(\rho_1, \rho_2)} \sum_{j>0, i+j \le K} \frac{\rho_1^i}{i!} \rho_2^j = 1 - \pi_{.,0}(\rho_1, \rho_2). $$

**STEP 3: DIFFEOMORPHISM FROM EQUATION (39).**

**Lemma 1 (Diffeomorphism).** There exists a strictly increasing diffeomorphism $\phi : ]0, \lambda/\mu[ \to ]0, \infty[, \text{ and } \psi = \phi^{-1}, \text{ such that } (\rho_1, \rho_2) \text{ is a solution of equation (39) if and only if } \rho_2 = \phi(\rho_1).$

*Proof.* First, $(\rho_1, \rho_2) \in \Gamma$ is solution of equation (39) if and only if $f(\rho_1, \rho_2) = 0$ where $f$ is the $C^\infty$ function defined by

$$ (42) \qquad f(\rho_1, \rho_2) = \left(\frac{\lambda}{\mu} - \rho_1\right) Z(\rho_1, \rho_2) - \frac{\lambda}{\mu} \sum_{i=0}^{K} \frac{\rho_1^i}{i!}. $$

To prove the existence of a strictly increasing diffeomorphism $\phi$ which maps $\rho_1$ to $\rho_2$, we apply the global inverse function theorem to auxiliary function $h$ defined on $\Gamma$ by

$$ h(\rho_1, \rho_2) = (\rho_1, f(\rho_1, \rho_2)). $$

Indeed, $h$ is injective if and only if, for any $\rho_1 \in ]0, \lambda/\mu[ \text{ and } \rho_2, \rho'_2 \in ]0, +\infty[, f(\rho_1, \rho_2) = f(\rho_1, \rho'_2)$ implies that $\rho_2 = \rho'_2$. As $f$ is $C^1$ on $]0, \lambda/\mu[ \times ]0, +\infty[$, it is sufficient to prove that $\partial f / \partial \rho_2 \neq 0$ on $]0, \lambda/\mu[ \times ]0, +\infty[$ to obtain that $\rho_2 \mapsto f(\rho_1, \rho_2)$ is strictly monotone thus injective.

Note that it is easy to see that $\partial f / \partial \rho_2$ is positive on $]0, \lambda/\mu[ \times ]0, +\infty[$. Indeed, since $\rho_1 < \lambda/\mu$ and $\rho_2 \mapsto Z(\rho_1, \rho_2)$ is non decreasing,

$$ \frac{\partial f}{\partial \rho_2}(\rho_1, \rho_2) = \left(\frac{\lambda}{\mu} - \rho_1\right) \frac{\partial Z}{\partial \rho_2}(\rho_1, \rho_2) > 0. $$

As a consequence, $h$ is injective and $C^1$ on $]0, \lambda/\mu[ \times ]0, +\infty[$. By global inversion function theorem, $h^{-1}$ is $C^1$ on $h(]0, \lambda/\mu[ \times ]0, +\infty[$) and it exists $\phi: C^1$ defined on $]0, \lambda/\mu[$ by

$$ h^{-1}(\rho_1, 0) = (\rho_1, \phi(\rho_1)). $$

Thus, to prove that diffeomorphism $\phi$ is strictly increasing on $]0, \lambda/\mu[$, it amounts to showing that, for all $\rho_1 \in ]0, \lambda/\mu[$, one gets

$$ \frac{\partial f}{\partial \rho_1}(\rho_1, \phi(\rho_1)) < 0, $$
or equivalently that, for all $(\rho_1, \rho_2) \in ]0, \lambda/\mu[\times]0, +\infty[$ such that $f(\rho_1, \rho_2) = 0$,

$$ \frac{\partial f}{\partial \rho_1}(\rho_1, \rho_2) < 0. $$

First, from equation (42), one obtains the first partial derivative of function $f$,

$$ (43) \qquad \frac{\partial f}{\partial \rho_1}(\rho_1, \rho_2) = -Z(\rho_1, \rho_2) + \left(\frac{\lambda}{\mu} - \rho_1\right) \sum_{i+j \le K-1} \frac{\rho_1^i}{i!} \rho_2^j - \frac{\lambda}{\mu} \sum_{i=0}^{K-1} \frac{\rho_1^i}{i!} $$

thus, subtracting equation (43) to equation (42) yields

$$ (44) \qquad f(\rho_1, \rho_2) - \frac{\partial f}{\partial \rho_1}(\rho_1, \rho_2) = Z(\rho_1, \rho_2) + \left(\frac{\lambda}{\mu} - \rho_1\right) \sum_{j=0}^{K} \frac{\rho_2^j \rho_1^{K-j}}{(K-j)!} - \frac{\lambda}{\mu} \frac{\rho_1^K}{K!}. $$

By equation (42), $f(\rho_1, \rho_2) = 0$ can be rewritten

$$ \frac{\lambda}{\mu} - \rho_1 = \frac{\lambda}{\mu Z(\rho_1, \rho_2)} \sum_{i=0}^{K} \frac{\rho_1^i}{i!}. $$

Thus, the second term of the right-hand side of equation (44) is

$$ (45) \qquad \left(\frac{\lambda}{\mu} - \rho_1\right) \sum_{j=0}^{K} \frac{\rho_2^j \rho_1^{K-j}}{(K-j)!} = \frac{\lambda}{\mu} \frac{1}{Z(\rho_1, \rho_2)} \sum_{i=0}^{K} \frac{\rho_1^i}{i!} \sum_{j=0}^{K} \frac{\rho_2^j \rho_1^{K-j}}{(K-j)!} \\ = \frac{\lambda}{\mu} \rho_1^K \frac{1}{Z(\rho_1, \rho_2)} \sum_{i,j=0}^{K} \frac{\rho_1^{i-j} \rho_2^j}{i!(K-j)!} \ge \frac{\lambda}{\mu} \rho_1^K \frac{1}{Z(\rho_1, \rho_2)} \sum_{0 \le j \le i \le K} \frac{\rho_1^{i-j} \rho_2^j}{i!(K-j)!} $$

using that all the terms are positive. For the sum of the right-hand side of (45), it holds that

$$ \sum_{0 \le j \le i \le K} \frac{\rho_1^{i-j} \rho_2^j}{i!(K-j)!} = \frac{1}{K!} \sum_{j+k \le K} \frac{K! \rho_1^k \rho_2^j}{(j+k)!(K-j)!} \ge \frac{Z(\rho_1, \rho_2)}{K!} $$

using that, for any $j, k \in \mathbb{N}, j + k \le K$,

$$ \frac{K!}{(j+k)!(K-j)!} = \frac{1}{k!} \frac{K!}{(K-j)!} \frac{k!}{(j+k)!} = \frac{1}{k!} \prod_{i=0}^{j-1} \frac{K-i}{j+k-i} \ge \frac{1}{k!}. $$

In conclusion, equation (45) gives that

$$ (46) \qquad \left(\frac{\lambda}{\mu} - \rho_1\right) \sum_{k=0}^{K} \frac{\rho_2^K \rho_1^{K-k}}{(K-k)!} \ge \frac{\lambda}{\mu} \frac{\rho_1^K}{K!}. $$

Plugging equation (46) in equation (44) and using that $Z(\rho_1, \rho_2) > 0$, it turns out that

$$ f(\rho_1, \rho_2) - \frac{\partial f}{\partial \rho_1}(\rho_1, \rho_2) > 0. $$

Therefore, as $f(\rho_1, \rho_2) = 0$,

$$ \frac{\partial f}{\partial \rho_1}(\rho_1, \rho_2) < 0 $$

for all $(\rho_1, \rho_2) \in ]0, \lambda/\mu[\times]0, +\infty[$ such that $f(\rho_1, \rho_2) = 0$.

□
or equivalently that, for all $(\rho_1, \rho_2) \in ]0, \lambda/\mu[\times]0, +\infty[$ such that $f(\rho_1, \rho_2) = 0$,

$$ \frac{\partial f}{\partial \rho_1}(\rho_1, \rho_2) < 0. $$

First, from equation (42), one obtains the first partial derivative of function $f$,

$$ (43) \qquad \frac{\partial f}{\partial \rho_1}(\rho_1, \rho_2) = -Z(\rho_1, \rho_2) + \left(\frac{\lambda}{\mu} - \rho_1\right) \sum_{i+j \le K-1} \frac{\rho_1^i}{i!} \rho_2^j - \frac{\lambda}{\mu} \sum_{i=0}^{K-1} \frac{\rho_1^i}{i!} $$

thus, subtracting equation (43) to equation (42) yields

$$ (44) \qquad f(\rho_1, \rho_2) - \frac{\partial f}{\partial \rho_1}(\rho_1, \rho_2) = Z(\rho_1, \rho_2) + \left(\frac{\lambda}{\mu} - \rho_1\right) \sum_{j=0}^{K} \frac{\rho_2^j \rho_1^{K-j}}{(K-j)!} - \frac{\lambda}{\mu} \frac{\rho_1^K}{K!}. $$

By equation (42), $f(\rho_1, \rho_2) = 0$ can be rewritten

$$ \frac{\lambda}{\mu} - \rho_1 = \frac{\lambda}{\mu Z(\rho_1, \rho_2)} \sum_{i=0}^{K} \frac{\rho_1^i}{i!}. $$

Thus, the second term of the right-hand side of equation (44) is

$$ (45) \qquad \left(\frac{\lambda}{\mu} - \rho_1\right) \sum_{j=0}^{K} \frac{\rho_2^j \rho_1^{K-j}}{(K-j)!} = \frac{\lambda}{\mu} \frac{1}{Z(\rho_1, \rho_2)} \sum_{i=0}^{K} \frac{\rho_1^i}{i!} \sum_{j=0}^{K} \frac{\rho_2^j \rho_1^{K-j}}{(K-j)!} \\ = \frac{\lambda}{\mu} \rho_1^K \frac{1}{Z(\rho_1, \rho_2)} \sum_{i,j=0}^{K} \frac{\rho_1^{i-j} \rho_2^j}{i!(K-j)!} \ge \frac{\lambda}{\mu} \rho_1^K \frac{1}{Z(\rho_1, \rho_2)} \sum_{0 \le j \le i \le K} \frac{\rho_1^{i-j} \rho_2^j}{i!(K-j)!} $$

using that all the terms are positive. For the sum of the right-hand side of (45), it holds that

$$ \sum_{0 \le j \le i \le K} \frac{\rho_1^{i-j} \rho_2^j}{i!(K-j)!} = \frac{1}{K!} \sum_{j+k \le K} \frac{K! \rho_1^k \rho_2^j}{(j+k)!(K-j)!} \ge \frac{Z(\rho_1, \rho_2)}{K!} $$

using that, for any $j, k \in \mathbb{N}, j + k \le K$,

$$ \frac{K!}{(j+k)!(K-j)!} = \frac{1}{k!} \frac{K!}{(K-j)!} \frac{k!}{(j+k)!} = \frac{1}{k!} \prod_{i=0}^{j-1} \frac{K-i}{j+k-i} \ge \frac{1}{k!}. $$

In conclusion, equation (45) gives that

$$ (46) \qquad \left(\frac{\lambda}{\mu} - \rho_1\right) \sum_{k=0}^{K} \frac{\rho_2^K \rho_1^{K-k}}{(K-k)!} \ge \frac{\lambda}{\mu} \frac{\rho_1^K}{K!}. $$

Plugging equation (46) in equation (44) and using that $Z(\rho_1, \rho_2) > 0$, it turns out that

$$ f(\rho_1, \rho_2) - \frac{\partial f}{\partial \rho_1}(\rho_1, \rho_2) > 0. $$

Therefore, as $f(\rho_1, \rho_2) = 0$,

$$ \frac{\partial f}{\partial \rho_1}(\rho_1, \rho_2) < 0 $$

for all $(\rho_1, \rho_2) \in ]0, \lambda/\mu[\times]0, +\infty[$ such that $f(\rho_1, \rho_2) = 0$.

□
Lemma 1 means that the solutions of equation (39) can be expressed with only one parameter. Note that it will be useful for the study of the equilibrium, especially in calculations to obtain asymptotics. This concludes the third step of the proof.

STEP 4: A MONOTONICITY ARGUMENT TO CONCLUDE SIMPLE RESERVATION CASE.
After equations (39) and (40), let us focus now on equation (41). The idea is to prove that the mean number in a tandem of two queues with total capacity K is a strictly increasing function of both arrival-to-service rates. It generalises the monotonicity argument in the similar proof in [10, Section 3.1]. Then using Lemma 1, it concludes to the existence and uniqueness of the equilibrium point.

**Lemma 2 (Monotonicity).** The average number $\mathbb{E}(R+V)$ of vehicles and reserved spaces per station, where $(R, V)$ is a random variable with distribution $\pi(\rho_2, \rho_1)$, is a strictly increasing function of both $\rho_2$ and $\rho_1$.

*Proof.* Let $(\rho_1, \rho_2) \mapsto \mathbb{E}(R+V)$ be denoted by $g_K$. It is sufficient to prove that, for all $(\rho_1, \rho_2) \in \Gamma$,

$$ \frac{\partial g_K}{\partial \rho_1}(\rho_1, \rho_2) > 0 \text{ and } \frac{\partial g_K}{\partial \rho_2}(\rho_1, \rho_2) > 0 $$

by induction on K.

By a change of indexes, $g_K$ can be rewritten

$$ g_K = \frac{\sum_{i+j \le K} (i+j) \frac{\rho_1^i}{i!} \rho_2^j}{\sum_{i+j \le K} \frac{\rho_1^i}{i!} \rho_2^j} = \frac{\sum_{k=0}^{K} k p_k}{\sum_{k=0}^{K} p_k} $$

where, by definition,

$$ p_k = \sum_{j=0}^{k} \rho_2^j \frac{\rho_1^{k-j}}{(k-j)!} $$

Define also, for $(k, l) \in \mathbb{N}^2$,

$$ r_{l,k} = \frac{p_l}{p_k}. $$

Let $k > 0$ be fixed. We first show that $r_{k,k-1}$ is an increasing function of both $\rho_1$ and $\rho_2$. Indeed

$$ r_{k,k-1} = \frac{p_k}{p_{k-1}} = \frac{\rho_2 p_{k-1} + \rho_1^k / k!}{p_{k-1}} = \rho_2 + \frac{\rho_1^k}{k!} \frac{1}{p_{k-1}}, $$

thus

$$ (47) \qquad \frac{\partial r_{k,k-1}}{\partial \rho_2} = 1 - \frac{\rho_1^k}{k!} \frac{\partial p_{k-1}}{\partial \rho_2} \frac{1}{p_{k-1}^2}. $$

But, by definition of $p_k$,

$$ (48) \qquad \rho_1^k \frac{\partial p_{k-1}}{\partial \rho_2} = \sum_{j=1}^{k-1} j \rho_2^{j-1} \frac{\rho_1^{2k-j-1}}{(k-j-1)!} = \sum_{j=0}^{k-2} (j+1) \rho_2^j \frac{\rho_1^{2k-j-2}}{(k-j-2)!}. $$

And, using that all terms of the sum in the following equation are positive,

$$ (49) \qquad k! p_{k-1}^2 = k! \sum_{1 \le u, v \le k} \rho_2^{u+v-2} \frac{\rho_1^{2k-u-v}}{(k-u)!(k-v)!} \\ > \sum_{i=0}^{k-2} \rho_2^i \rho_1^{2k-i-2} \sum_{j=1}^{i+1} \frac{k!}{(k-2-i+j)!(k-j)!}. $$

For all $j$, $0 \le j \le i+1$, while $k-2-i+j < k$, it holds that
Lemma 1 means that the solutions of equation (39) can be expressed with only one parameter. Note that it will be useful for the study of the equilibrium, especially in calculations to obtain asymptotics. This concludes the third step of the proof.

**STEP 4: A MONOTONICITY ARGUMENT TO CONCLUDE SIMPLE RESERVATION CASE.**
After equations (39) and (40), let us focus now on equation (41). The idea is to prove that the mean number in a tandem of two queues with total capacity K is a strictly increasing function of both arrival-to-service rates. It generalises the monotonicity argument in the similar proof in [10, Section 3.1]. Then using Lemma 1, it concludes to the existence and uniqueness of the equilibrium point.

**Lemma 2 (Monotonicity).** The average number $E(R+V)$ of vehicles and reserved spaces per station, where $(R, V)$ is a random variable with distribution $\pi(\rho_2, \rho_1)$, is a strictly increasing function of both $\rho_2$ and $\rho_1$.

*Proof.* Let $(\rho_1, \rho_2) \mapsto E(R+V)$ be denoted by $g_K$. It is sufficient to prove that, for all $(\rho_1, \rho_2) \in \Gamma$,

$$ \frac{\partial g_K}{\partial \rho_1}(\rho_1, \rho_2) > 0 \text{ and } \frac{\partial g_K}{\partial \rho_2}(\rho_1, \rho_2) > 0 $$

by induction on K.

By a change of indexes, $g_K$ can be rewritten

$$ g_K = \frac{\sum_{i+j \le K} (i+j) \frac{\rho_1^i}{i!} \rho_2^j}{\sum_{i+j \le K} \frac{\rho_1^i}{i!} \rho_2^j} = \frac{\sum_{k=0}^{K} k p_k}{\sum_{k=0}^{K} p_k} $$

where, by definition,

$$ p_k = \sum_{j=0}^{k} \rho_2^j \frac{\rho_1^{k-j}}{(k-j)!} $$

Define also, for $(k, l) \in \mathbb{N}^2$,

$$ r_{l,k} = \frac{p_l}{p_k}. $$

Let $k > 0$ be fixed. We first show that $r_{k,k-1}$ is an increasing function of both $\rho_1$ and $\rho_2$. Indeed

$$ r_{k,k-1} = \frac{p_k}{p_{k-1}} = \frac{\rho_2 p_{k-1} + \rho_1^k / k!}{p_{k-1}} = \rho_2 + \frac{\rho_1^k}{k!} \frac{1}{p_{k-1}}, $$

thus

$$ (47) \qquad \frac{\partial r_{k,k-1}}{\partial \rho_2} = 1 - \frac{\rho_1^k}{k!} \frac{\partial p_{k-1}}{\partial \rho_2} \frac{1}{p_{k-1}^2}. $$

But, by definition of $p_k$,

$$ (48) \qquad \rho_1^k \frac{\partial p_{k-1}}{\partial \rho_2} = \sum_{j=1}^{k-1} j \rho_2^{j-1} \frac{\rho_1^{2k-j-1}}{(k-j-1)!} = \sum_{j=0}^{k-2} (j+1) \rho_2^j \frac{\rho_1^{2k-j-2}}{(k-j-2)!}. $$

And, using that all terms of the sum in the following equation are positive,

$$ (49) \qquad k! p_{k-1}^2 = k! \sum_{1 \le u, v \le k} \rho_2^{u+v-2} \frac{\rho_1^{2k-u-v}}{(k-u)!(k-v)!} \\ > \sum_{i=0}^{k-2} \rho_2^i \rho_1^{2k-i-2} \sum_{j=1}^{i+1} \frac{k!}{(k-2-i+j)!(k-j)!}. $$

For all $j, 0 \le j \le i+1$, while $k-2-i+j < k$, it holds that
$$ \frac{k!}{(k-2-i+j)!(k-j)!} > \frac{1}{(k-i-2)!} $$

and then

$$ \sum_{j=1}^{i+1} \frac{k!}{(k-2-i+j)!(k-j)!} > \frac{i+1}{(k-i-2)!} $$

Plugging in equation (49) and comparing with equation (48), it gives

$$ k!p_{k-1}^2 > \rho_1^k \frac{\partial p_{k-1}}{\partial \rho_2}. $$

Therefore, using (47), it allows to conclude that

$$ \frac{\partial r_{k,k-1}}{\partial \rho_2} > 0. $$

Moreover,

$$ (50) \qquad \frac{\partial r_{k,k-1}}{\partial \rho_1} = \frac{\rho_1^{k-1}}{(k-1)!p_{k-1}} \left( 1 - \frac{\rho_1}{k p_{k-1}} \frac{\partial p_{k-1}}{\partial \rho_1} \right) > 0. $$

because it is easily checked that $\partial p_{k-1}/\partial \rho_1 = p_{k-2}$ and then $k p_{k-1} > \rho_1 p_{k-2}$.

Therefore, if $l > k$, $r_{l,k} = \prod_{i=k+1}^{l} r_{i,i-1}$ is an increasing function of $\rho_2$ and $\rho_1$. This gives that $u_K$ defined by

$$ u_K = \frac{p_K}{\sum_{k=0}^{K} p_k} = \frac{1}{\sum_{k=0}^{K} r_{k,K}} $$

is non decreasing in $x \in \{\rho_1, \rho_2\}$, because $r_{k,K} = 1/r_{K,k}$ is non increasing in $x$.

Note that $g_0$ is constant with $x$ and that

$$ g_K = (1 - u_K)g_{K-1} + Ku_K, $$

which yields that

$$ \frac{\partial g_K}{\partial x} = (K - g_{K-1}) \frac{\partial u_K}{\partial x} + (1 - u_K) \frac{\partial g_{K-1}}{\partial x}. $$

Since $K - g_{K-1} > 0$, $u_K < 1$ and $\partial u_K/\partial x > 0$, by induction we can conclude that $\partial g_K/\partial x > 0$ for all $K \ge 1$. It ends the proof. $\square$

By step 2, it remains to prove that, for any $s > 0$, there exists a unique $(\rho_1, \rho_2)$ solution of both equations (39) and (41). By Lemma 1, equation (39) can be rewritten $\rho_1 = \psi_{\rho_2}(\rho_2)$ with $\rho_2 \ge 0$. Then equation (41) becomes

$$ (51) \qquad s = g_K(\psi_{\rho_2}(\rho_2), \rho_2). $$

Let us denote the right-hand side of equation (51) by function $s_K$ defined on $[0, +\infty[$ which maps $\rho_2$ to $g_K(\psi_{\rho_2}(\rho_2), \rho_2) = E(R+V)$. It is sufficient to prove that $s_K$ is strictly increasing. It is true, using that

$$ s'_{K} = \frac{\partial g_{K}}{\partial \rho_{2}} + \frac{\partial g_{K}}{\partial \rho_{1}} \psi'_{\rho_{2}} $$

as equation (51) holds with $\psi_{\rho_2}$ at $\rho_2$ but also in a neighborhood of $\rho_2$ and using both Lemmas 1 and 2. To conclude the proof, one can easily check that $s_K$ covers the whole interval $[0, +\infty[$. Indeed if $\rho_2 = 0$, the only $\rho_1$ solution of $f(\rho_1, \rho_2) = 0$ is 0 and $E(R+V) = 0$. Similarly, when $\rho_2$ tends to infinity, $\rho_1$ has to tend to $\lambda/\mu$ to keep $f(\rho_1, \rho_2) = 0$ and $E(R+V)$ tends to $+\infty$. It ends the proof for the single reservation case.
$$ \frac{k!}{(k-2-i+j)!(k-j)!} > \frac{1}{(k-i-2)!} $$

and then

$$ \sum_{j=1}^{i+1} \frac{k!}{(k-2-i+j)!(k-j)!} > \frac{i+1}{(k-i-2)!} $$

Plugging in equation (49) and comparing with equation (48), it gives

$$ k! p_{k-1}^2 > \rho_1^k \frac{\partial p_{k-1}}{\partial \rho_2}. $$

Therefore, using (47), it allows to conclude that

$$ \frac{\partial r_{k,k-1}}{\partial \rho_2} > 0. $$

Moreover,

$$ (50) \qquad \frac{\partial r_{k,k-1}}{\partial \rho_1} = \frac{\rho_1^{k-1}}{(k-1)! p_{k-1}} \left( 1 - \frac{\rho_1}{k p_{k-1}} \frac{\partial p_{k-1}}{\partial \rho_1} \right) > 0. $$

because it is easily checked that $\partial p_{k-1}/\partial \rho_1 = p_{k-2}$ and then $kp_{k-1} > \rho_1 p_{k-2}$.

Therefore, if $l > k$, $r_{l,k} = \prod_{i=k+1}^{l} r_{i,i-1}$ is an increasing function of $\rho_2$ and $\rho_1$. This gives that $u_K$ defined by

$$ u_K = \frac{p_K}{\sum_{k=0}^{K} p_k} = \frac{1}{\sum_{k=0}^{K} r_{k,K}} $$

is non decreasing in $x \in \{\rho_1, \rho_2\}$, because $r_{k,K} = 1/r_{K,k}$ is non increasing in $x$.

Note that $g_0$ is constant with $x$ and that

$$ g_K = (1 - u_K)g_{K-1} + Ku_K, $$

which yields that

$$ \frac{\partial g_K}{\partial x} = (K - g_{K-1}) \frac{\partial u_K}{\partial x} + (1 - u_K) \frac{\partial g_{K-1}}{\partial x}. $$

Since $K - g_{K-1} > 0$, $u_K < 1$ and $\partial u_K/\partial x > 0$, by induction we can conclude that $\partial g_K/\partial x > 0$ for all $K \ge 1$. It ends the proof. $\square$

By step 2, it remains to prove that, for any $s > 0$, there exists a unique $(\rho_1, \rho_2)$ solution of both equations (39) and (41). By Lemma 1, equation (39) can be rewritten $\rho_1 = \psi_{\rho_2}(\rho_2)$ with $\rho_2 \ge 0$. Then equation (41) becomes

$$ (51) \qquad s = g_K(\psi_{\rho_2}(\rho_2), \rho_2). $$

Let us denote the right-hand side of equation (51) by function $s_K$ defined on $[0, +\infty[$ which maps $\rho_2$ to $g_K(\psi_{\rho_2}(\rho_2), \rho_2) = E(R+V)$. It is sufficient to prove that $s_K$ is strictly increasing. It is true, using that

$$ s'_{K} = \frac{\partial g_{K}}{\partial \rho_{2}} + \frac{\partial g_{K}}{\partial \rho_{1}} \psi'_{\rho_{2}} $$

as equation (51) holds with $\psi_{\rho_2}$ at $\rho_2$ but also in a neighborhood of $\rho_2$ and using both Lemmas 1 and 2. To conclude the proof, one can easily check that $s_K$ covers the whole interval $[0, +\infty[$. Indeed if $\rho_2 = 0$, the only $\rho_1$ solution of $f(\rho_1, \rho_2) = 0$ is 0 and $E(R+V) = 0$. Similarly, when $\rho_2$ tends to infinity, $\rho_1$ has to tend to $\lambda/\mu$ to keep $f(\rho_1, \rho_2) = 0$ and $E(R+V)$ tends to $+\infty$. It ends the proof for the single reservation case.
**STEP 5: THE DOUBLE RESERVATION CASE.** Straightforwardly equations (32) and (35) leads to

$$\eta_1 = \eta_2 = \frac{\mu}{\nu} \rho_1.$$

Then the main argument is that equations (33) and (34) can be rewritten as equations (39) and (40) where $\rho_1$ is replaced by $\tilde{\rho}_1 = (1 + 2\mu/\nu)\rho_1$ and $\lambda/\mu$ by $a = \lambda/\mu(1 + 2\mu/\nu)$. Indeed, by straightforward algebra, it holds that $\pi_{0V}(\rho) = \pi_{.,0}(\tilde{\rho}_1, \rho_2)$ and $\pi_S(\rho) = \pi_S(\tilde{\rho}_1, \rho_2)$.

Unfortunately, (36) is not rewritten as (41) with the previous change of variables $\tilde{\rho}_1 = (1 + 2\mu/\nu)\rho_1$ but, with careful calculations, as

$$(52) \qquad s = \sum_{i+j \le K} \left( \frac{1 + \mu/\nu}{1 + 2\mu/\nu} i + j \right) \pi_{i,j}(\tilde{\rho}_1, \rho_2).$$

To complete the proof for any $\nu > 0$, it amounts to proving that, for $a > 0$ and $K \in \mathbb{N}$,

$$S(x, y) = \frac{1}{Z(x, y)} \sum_{i+j \le K} (i+2j) \frac{x^i}{i!} y^j$$

is a strictly increasing function of $x \in [0, a]$ under $f(x, y) = 0$, where

$$f(x, y) = (a-x)Z(x, y) - a \sum_{i=0}^{K} \frac{x^i}{i!}, \quad Z(x, y) = \sum_{i+j \le K} \frac{x^i}{i!} y^j.$$

Indeed, the ratio in the right-hand side of equation (52) is a non increasing function from $(0, +\infty)$ to $(1/2, 1)$. By linearity of differentiability and multiplication by a scalar, it suffices to prove that the right-hand side of equation (52) is non-increasing as a function of $\tilde{\rho}_1$, for $\rho_2 = \phi(\rho_1)$ defined in Lemma 1 for the two values $1/2$ and $1$ of this ratio. For value 1, it is exactly Lemma 2. It remains to prove it for value $1/2$, which has been asserted previously.

Nevertheless, due to simple reservation case (steps 2 to 4), by continuity, the proof is complete for $\nu$ enough large.

## ACKNOWLEDGEMENT

The authors would like to thank Cédric Bourdais for his contribution to the analysis of this model during his internship at INRIA.

## REFERENCES

[1] Reza Aghajani, Philippe Robert, and Wen Sun. A large scale analysis of unreliable stochastic networks. Ann. Appl. Probab., 28(2):851–887, 2018.

[2] François Baccelli, Sergey Foss, and Vsevolod Shneer. Migration-Contagion Processes. arXiv preprint, 2022.

[3] Patrick Billingsley. Convergence of probability measures. Wiley Series in Probability and Statistics: Probability and Statistics. John Wiley & Sons Inc., New York, second edition, 1999. A Wiley-Interscience Publication.

[4] Chiara Boldrini, Raffaele Bruno, and Marco Conti. Characterising demand and usage patterns in a large station-based car sharing system. In 2016 IEEE Conference on Computer Communications Workshops (INFOCOM WKSHPS), pages 572–577. IEEE, 2016.

[5] Cédric Bourdais and Christine Fricker. A stochastic model for car-sharing systems. arXiv preprint arXiv:1504.03844, 2015.

[6] Cedric Bourdais, Christine Fricker and Hanene Mohamed. A mean field analysis of a stochastic model for reservation in car-sharing systems. ACM SIGMETRICS Performance Evaluation Review, 48(2):18–20, 2020.

[7] Burak Boyacı, Konstantinos G Zografos, and Nikolas Geroliminis. An optimization framework for the development of efficient one-way car-sharing systems. European Journal of Operational Research, 240(3):718–733, 2015.

[8] Yuntian Deng, Abhishek Gupta, and Ness B Shroff. Fleet sizing and charger allocation in electric vehicle sharing systems. arXiv preprint arXiv:2007.06687, 2020.

[9] Christine Fricker, Nicolas Gast, and Hanene Mohamed. Mean field analysis for inhomogeneous bike sharing systems. Discrete Mathematics & Theoretical Computer Science, 2012.
**STEP 5: THE DOUBLE RESERVATION CASE.** Straightforwardly equations (32) and (35) leads to

$$ \eta_1 = \eta_2 = \frac{\mu}{\nu} \rho_1. $$

Then the main argument is that equations (33) and (34) can be rewritten as equations (39) and (40) where $\rho_1$ is replaced by $\tilde{\rho}_1 = (1 + 2\mu/\nu)\rho_1$ and $\lambda/\mu$ by $a = \lambda/\mu(1 + 2\mu/\nu)$. Indeed, by straightforward algebra, it holds that $\pi_{0V}(\rho) = \pi_{.,0}(\tilde{\rho}_1, \rho_2)$ and $\pi_S(\rho) = \pi_S(\tilde{\rho}_1, \rho_2)$.

Unfortunately, (36) is not rewritten as (41) with the previous change of variables $\tilde{\rho}_1 = (1 + 2\mu/\nu)\rho_1$ but, with careful calculations, as

$$ (52) \qquad s = \sum_{i+j \le K} \left( \frac{1 + \mu/\nu}{1 + 2\mu/\nu} i + j \right) \pi_{i,j}(\tilde{\rho}_1, \rho_2). $$

To complete the proof for any $\nu > 0$, it amounts to proving that, for $a > 0$ and $K \in \mathbb{N}$,

$$ S(x, y) = \frac{1}{Z(x, y)} \sum_{i+j \le K} (i+2j) \frac{x^i}{i!} y^j $$

is a strictly increasing function of $x \in [0, a]$ under $f(x, y) = 0$, where

$$ f(x, y) = (a - x)Z(x, y) - a \sum_{i=0}^{K} \frac{x^i}{i!}, \quad Z(x, y) = \sum_{i+j \le K} \frac{x^i}{i!} y^j. $$

Indeed, the ratio in the right-hand side of equation (52) is a non increasing function from $(0, +\infty)$ to $(1/2, 1)$. By linearity of differentiability and multiplication by a scalar, it suffices to prove that the right-hand side of equation (52) is non-increasing as a function of $\tilde{\rho}_1$, for $\rho_2 = \phi(\rho_1)$ defined in Lemma 1 for the two values $1/2$ and $1$ of this ratio. For value 1, it is exactly Lemma 2. It remains to prove it for value $1/2$, which has been asserted previously.

Nevertheless, due to simple reservation case (steps 2 to 4), by continuity, the proof is complete for $\nu$ enough large.

## ACKNOWLEDGEMENT

The authors would like to thank Cédric Bourdais for his contribution to the analysis of this model during his internship at INRIA.

## REFERENCES

[1] Reza Aghajani, Philippe Robert, and Wen Sun. A large scale analysis of unreliable stochastic networks. *Ann. Appl. Probab.*, 28(2):851–887, 2018.

[2] François Baccelli, Sergey Foss, and Vsevolod Shneer. Migration-Contagion Processes. *arXiv preprint*, 2022.

[3] Patrick Billingsley. *Convergence of probability measures*. Wiley Series in Probability and Statistics: Probability and Statistics. John Wiley & Sons Inc., New York, second edition, 1999. A Wiley-Interscience Publication.

[4] Chiara Boldrini, Raffaele Bruno, and Marco Conti. Characterising demand and usage patterns in a large station-based car sharing system. In *2016 IEEE Conference on Computer Communications Workshops (INFOCOM WKSHPS)*, pages 572–577. IEEE, 2016.

[5] Cédric Bourdais and Christine Fricker. A stochastic model for car-sharing systems. *arXiv preprint arXiv:1504.03844*, 2015.

[6] Cedric Bourdais, Christine Fricker and Hanene Mohamed. A mean field analysis of a stochastic model for reservation in car-sharing systems. *ACM SIGMETRICS Performance Evaluation Review*, 48(2):18–20, 2020.

[7] Burak Boyacı, Konstantinos G Zografos, and Nikolas Geroliminis. An optimization framework for the development of efficient one-way car-sharing systems. *European Journal of Operational Research*, 240(3):718–733, 2015.

[8] Yuntian Deng, Abhishek Gupta, and Ness B Shroff. Fleet sizing and charger allocation in electric vehicle sharing systems. *arXiv preprint arXiv:2007.06687*, 2020.

[9] Christine Fricker, Nicolas Gast, and Hanene Mohamed. Mean field analysis for inhomogeneous bike sharing systems. *Discrete Mathematics & Theoretical Computer Science*, 2012.
[10] Christine Fricker and Nicolas Gast. Incentives and redistribution in homogeneous bike-sharing systems with stations of finite capacity. *Euro journal on transportation and logistics*, 5(3):261–291, 2016.

[11] Christine Fricker and Danielle Tibi. Equivalence of ensembles for large vehicle-sharing models. *The Annals of Applied Probability*, 27(2): 883–916, 2017.

[12] Christine Fricker, Hanene Mohamed, Teodora Popescu, and Martin Trépanier. Stochastic modelling of free-floating car-sharing systems. *CIRRELT Research Report*, 2021.

[13] Gibbens, R. J., Hunt, P. J., and Kelly, F. P. Bistability in communication networks. *Disorder in physical systems*, 113-128, 1990.

[14] Carl Graham and Sylvie Méléard. Propagation of chaos for a fully connected loss network with alternate routing. *Stochastic processes and their applications*, 44(1):159-180, 1993.

[15] Mor Kaspi, Tal Raviv, and Michal Tzur. Parking reservation policies in one-way vehicle sharing systems. *Transportation Research Part B: Methodological*, 62:35–50, 2014.

[16] Mor Kaspi, Tal Raviv, Michal Tzur, and Hila Galili. Regulating vehicle sharing systems through parking reservation policies: Analysis and performance bounds. *European Journal of Operational Research*, 251(3):969–987, 2016.

[17] Davit Martirosyan and Philippe Robert. The equilibrium states of large networks of Erlang queues. *Advances in Applied Probability*, 52(2): 617-654, 2020.

[18] Mehdi Nourinejad, Sirui Zhu, Sina Bahrami, and Matthew J Roorda. Vehicle relocation and staff rebalancing in one-way carsharing systems. *Transportation Research Part E: Logistics and Transportation Review*, 81:98–113, 2015.

[19] Rapport Autolib 2013. https://autolibmetropole.fr/wp-content/uploads/2023/04/Rapport-dactivite-2013.pdf.

[20] Philippe Robert. *Stochastic networks and queues*, volume 52. Springer Science & Business Media, 2013.

[21] Alain-Sol Sznitman. Topics in propagation of chaos. In *Ecole d'été de probabilités de Saint-Flour XIX-1989*, pages 165–251. Springer, 1991.

[22] Autolib Wikipedia. https://en.wikipedia.org/wiki/Autolib%27

(C. FRICKER) DI-ENS, CNRS, INRIA DE PARIS, PARIS, FRANCE
*Email address: christine.fricker@inria.fr*

(H. MOHAMED) MODAL'X, UMR CNRS 9023, UPL, UNIVERSITÉ PARIS NANTERRE
*Email address: hanene.mohamed@parisnanterre.fr*
[10] Christine Fricker and Nicolas Gast. Incentives and redistribution in homogeneous bike-sharing systems with stations of finite capacity. *Euro journal on transportation and logistics*, 5(3):261–291, 2016.

[11] Christine Fricker and Danielle Tibi. Equivalence of ensembles for large vehicle-sharing models. *The Annals of Applied Probability*, 27(2): 883–916, 2017.

[12] Christine Fricker, Hanene Mohamed, Teodora Popescu, and Martin Trépanier. Stochastic modelling of free-floating car-sharing systems. *CIRRELT Research Report*, 2021.

[13] Gibbens, R. J., Hunt, P. J., and Kelly, F. P. Bistability in communication networks. *Disorder in physical systems*, 113-128, 1990.

[14] Carl Graham and Sylvie Méléard. Propagation of chaos for a fully connected loss network with alternate routing. *Stochastic processes and their applications*, 44(1):159-180, 1993.

[15] Mor Kaspi, Tal Raviv, and Michal Tzur. Parking reservation policies in one-way vehicle sharing systems. *Transportation Research Part B: Methodological*, 62:35–50, 2014.

[16] Mor Kaspi, Tal Raviv, Michal Tzur, and Hila Galili. Regulating vehicle sharing systems through parking reservation policies: Analysis and performance bounds. *European Journal of Operational Research*, 251(3):969–987, 2016.

[17] Davit Martirosyan and Philippe Robert. The equilibrium states of large networks of Erlang queues. *Advances in Applied Probability*, 52(2): 617-654, 2020.

[18] Mehdi Nourinejad, Sirui Zhu, Sina Bahrami, and Matthew J Roorda. Vehicle relocation and staff rebalancing in one-way carsharing systems. *Transportation Research Part E: Logistics and Transportation Review*, 81:98–113, 2015.

[19] Rapport Autolib 2013. https://autolibmetropole.fr/wp-content/uploads/2023/04/Rapport-dactivite-2013.pdf.

[20] Philippe Robert. *Stochastic networks and queues*, volume 52. Springer Science & Business Media, 2013.

[21] Alain-Sol Sznitman. Topics in propagation of chaos. In *Ecole d'été de probabilités de Saint-Flour XIX-1989*, pages 165–251. Springer, 1991.

[22] Autolib Wikipedia. https://en.wikipedia.org/wiki/Autolib%27

(C. FRICKER) DI-ENS, CNRS, INRIA DE PARIS, PARIS, FRANCE
*Email address: christine.fricker@inria.fr*

(H. MOHAMED) MODAL'X, UMR CNRS 9023, UPL, UNIVERSITÉ PARIS NANTERRE
*Email address: hanene.mohamed@parisnanterre.fr*