2.1. Microscopic traffic networks

We consider a road network that is typically embedded into a two-dimensional area $A\subseteq\mathbb{R}^2$ . The network may consist of roads, junctions, roundabouts, intersections, highways, etc. on which vehicles move. The collection of all vehicles in the network is denoted by $\mathcal{M}$ . Each vehicle $i\in\mathcal{M}$ is assigned an origin–destination pair $(O^i,D^i)\in A$ .

We consider a fixed time horizon $T>0$ . Vehicles move over time from their origin to their destination on a (potentially changing) path. We denote by $x^i(t)$ the position of vehicle i at time $t\in[0,T]$ , by $v^i(t)=\frac{\textrm{d}}{\textrm{d}t}x^i(t)$ their velocity, and by $a^i(t)=\frac{\textrm{d}}{\textrm{d}t}v^i(t)$ their acceleration. We make the implicit assumption that vehicles are located in $O^i$ until some release time and remain in $D^i$ once reached. Thus, we let $\mathcal{M}(t)=\left\{i\in\mathcal{M}\colon x^i(t)\notin\{O^i,D^i\}\right\}$ be those vehicles which are currently inside the network, that is, they have left their origin but not reached their destination, yet. If we only consider vehicles which belong to a certain group of vehicles (also called a fleet) $\Phi\subseteq\mathcal{M}$ , we will write $\mathcal{M}^\Phi(t)$ . This could, for example, be a fleet of vehicles with the same driving characteristics.

At the core of many microscopic traffic models are car-following models. Prominent examples include the Intelligent Driver Model (Treiber et al., Reference Treiber, Hennecke and Helbing2000), the Optimal Velocity Model (Bando et al., Reference Bando, Hasebe, Nakayama, Shibata and Sugiyama1994 and 1995), and the Krauß model (Krauß, Reference Krauß1998). Car-following models determine the acceleration behavior of an individual vehicle i along its path on the basis of information on the positions and velocities of the vehicles, typically in a neighborhood of i, and the properties of the system. Often, only the preceding vehicle on the road is relevant, and the acceleration of i is constructed such that vehicles move forward while maintaining a minimal distance. Through specific choices more complex traffic scenarios (e.g., intersections, overtaking) can still be represented in such a manner. Mathematically, car-following models correspond to systems of coupled ordinary differential equations.

Traffic State. We denote by $\gamma(t)=\left(x^i(t),v^i(t),a^i(t)\right)_{i\in\mathcal{M}}$ the state of the traffic system at time t. It records the position, velocity, and acceleration of any vehicle. The evolution of the traffic system over time is depicted by the (high-dimensional) trajectory $t\mapsto \gamma(t)$ .

Macroscopic traffic statistics aggregate these microscopic data. Typical examples include traffic flow (number of vehicles that pass a certain point per time unit), traffic density (number of vehicles per length unit), and average speed. These measures quantify the performance of traffic systems.

Local Traffic Conditions. In order to model the occurrence of accidents depending on local traffic conditions, we partition A in regions. More precisely, we partition A into a finite number of disjoint sets $A_r\subseteq A$ such that $ A = \bigcup_{r=1}^R A_r$ and $R \in\mathbb{N}$ . We call the elements $A_r$ of the partition a traffic module.

We let $\mathcal{M}_r(t)=\left\{i\in\mathcal{M}(t)\colon x^i(t)\in A_r\right\} \subseteq \mathcal{M}(t)$ denote those vehicles that are in $A_r$ at time t (with $\mathcal{M}^\Phi_r(t)$ defined in analogy). The local traffic state of the module is $\gamma_r(t) = ( x^i(t), v^i(t), a^i(t))_{i \in \mathcal{M}_r(t)}$ . Key local traffic characteristics (density, flow, speed, etc.) can then be expressed as functions of $\gamma_r(t)$ and its evolution over small time windows.

2.2. Microscopic traffic model with accident losses

So far, the evolution of the traffic system is a deterministic function of time. The advantage of microscopic models is that they enable a detailed simulation of traffic systems. The driving behavior of the vehicles can be varied, likewise their number and paths, road conditions, etc., in order to generate many different scenarios. Such models provide a detailed picture, similar to digital twins of reality, and can be used to analyze potential future traffic systems or to understand the impact of new technologies.

We consider a finite collection of different traffic scenarios $\gamma^k \,:\!=\, (\gamma^k(t))_{t\in [0,T]}$ with $k \in \{1,2, \dots, K \} $ for a short time horizon $T>0$ . The aim is to analyze characteristics of traffic over the long time horizon N T, for example, 1 year, for some large $N\in \mathbb{N}$ ; this is modeled by a finite sequence of traffic scenarios $(k_1, k_2, \dots, k_N) \in \{1,2, \dots, K \}^N $ . The N subintervals of length T are called time buckets. We will be interested in quantities aggregated or averaged over the whole time horizon NT. Examples include the average traffic flow, the total number of accidents, the aggregate losses due to accidents, etc. These quantities depend not only on the order of the traffic scenarios during this time period but also on their number of occurrences.

We denote by $\mu^k$ the number of occurrences of scenario k divided by N, that is, the relative frequency of this traffic scenario over the considered time horizon N T. The vector $\mu= (\mu^1, \mu^2, \dots, \mu^K)^\top$ lies in the simplex $\Delta^{K-1} = \{x \in \mathbb{R}^K_+: \sum_{k=1}^K x^k =1 \}$ . We assume that $\mu$ is not deterministic, but a random variable. This is to account for the fact that the relative frequencies of traffic scenarios fluctuate over different years due to varying weather conditions, random changes in traffic demand, or other factors. From a mathematical point of view, this construction leads to a mixture model with exogenous factor $\mu$ .

Accident Occurences. We now introduce our traffic accident model that will permit an analysis of aggregate losses and corresponding insurance contracts. The likelihood of the occurrence of an accident is modeled as a function of the traffic scenario. We consider two specifications: a binomial and a Poisson model.

1. (i) Binomial Model. Accidents are rare events. For a given traffic scenario k, we assume that the probability $p^k$ of an accident is close to zero. This probability may, of course, depend on the evolution of the traffic scenario, that is, on the path $t\mapsto \gamma^k(t)$ , and we will discuss concrete specifications later. Given a realization of $\mu$ , accidents are assumed to be independent across time buckets. This implies that traffic scenario k occurs for $N \mu^k$ time buckets corresponding to a duration of $N T \mu^k$ , and the number of accidents $C^k$ during this period has a conditional binomial distribution with parameters $p^k$ and $N \mu^k$ :

   \begin{equation*}C^k \mid \mu \sim {\rm Bin} \left(p^k, N \mu^k \right).\end{equation*}

2. (ii) Poisson Model. An alternative model assumes that accidents occur at random times with a distribution governed by an intensity $\lambda^k/T$ that depends on the traffic scenario k. More specifically, the number of accidents $C^k$ during the period governed by scenario k of duration $N T \mu^k$ is conditionally Poisson distributed with parameter $\lambda^k N \mu^k$ :

   \begin{equation*}C^k \mid \mu \sim {\rm Poiss} \left(\lambda^k N \mu^k\right).\end{equation*}

Accident Losses. Loss sizes conditional on the occurrence of accidents are assumed to be independent across traffic scenarios and across time buckets. We assume that the conditional loss distribution with conditional distribution function $F^k$ depends only on the traffic scenario k. We will discuss examples below. Random total losses over the considered time horizon NT are equal to

\begin{equation*}L\; = \; \sum_{k=1}^K \; \sum_{c=1}^{C^k} \; X^k_c \end{equation*}

where the random variables $X^k_c$ , $k=1,2, \dots, K$ , $c\in \mathbb{N}$ , are independent and $X^k_c \sim F^k$ , $c\in\mathbb{N}$ , for any k.

Concrete Specifications. Microscopic traffic models are experimental environments that allow to simulate the behavior of systems where no real data are yet available. Traffic planning can be supported by such models, and the impact of new technologies can be tested in a counterfactual analysis. Here, we specify the general principles how accident occurrences and losses can be based on microscopic traffic models. An implementation will in this paper be based on SUMO, see Section 2.3, but could also rely on any other suitable traffic model.

Initially, K traffic scenarios need to be selected as a basis for the model. While running any deterministic traffic scenario k over the time window [0, T], information can be extracted about the traffic states $\gamma_r^k$ in each module $r\in \{1,2, \dots, R\}$ . In SUMO typically not complete data on the whole paths are extracted, but only selected information at loop detectors in the network that are part of the implementation.

In reality, the likelihood of accidents typically increases with higher traffic density and higher velocities, ceteris paribus. Also the distribution of losses is influenced by quantities of this type. Examples are described in Section 4. This allows a computation of $p^k$ and $\lambda^k$ as a function of the data. Using the data associated with the modules, we may specify probabilities and intensities for the modules such that $p^k= \sum_{r=1}^R p^k_r$ and $\lambda^k = \sum_{r=1}^R \lambda^k_r$ . In the binomial model, $p^k_r/p^k$ is the conditional probability that the accident is in module r given that an accident occurs. In the Poisson model, the intensities $\lambda^k_r$ , $r=1,2, \dots, R$ , determine the accident times for each module. The resulting sequence of random times in the whole traffic system possesses the intensity $\lambda^k$ . Conversely, if we first simulate random times with intensity $\lambda^k$ and then randomly choose a corresponding module with probability $\lambda^k_r/\lambda^k$ in a second step, the random times associated with each module r possess intensity $\lambda^k_r$ . Both procedures produce a number of accidents $C^k$ that occur during the period governed by scenario k.

The distributions of losses given a single random event will be chosen as follows. For each traffic scenario, we consider a collection of distribution functions $(F^{k,\psi})_{\psi\in \Psi}$ where $\psi$ corresponds to data that may be extracted from the traffic simulation. In order to do so, we uniformly simulate a random time in [0, T] and extract at this time the data from scenario k that determine $\psi$ . The resulting distribution $F^k$ is a mixture of the distributions $(F^{k,\psi})_{\psi\in \Psi}$ . The mixing distribution is derived on the basis of the traffic data of scenario k that are generated from our microscopic traffic model.

Insurance Contracts and Statistical Functionals. The microscopic traffic model with accidents will be the basis for the generation of aggregate losses. We study statistical functionals and insurance contracts. The analysis will be based on Monte Carlo simulations, but we also compare approximation techniques that we describe in Section 3.

The focus is on functions of aggregate losses L. Letting $h\colon \mathbb{R} \to \mathbb{R} $ be an increasing function, we analyze h(L). In particular, we investigate the following functions corresponding to three types of insurance coverage: $h(x)=x$ (full coverage), $h(x)=\max\!(x-\theta,0)$ , $\theta\geq 0$ (constant deductible), $h (x) =\min(x,\theta)$ , $\theta\geq 0$ (stop loss). In each case, we evaluate various statistical functionals:

1. (i) Expectation. $\mathbb{E}(h(L))$ ,

2. (ii) Variance. $\textrm{Var}(h(L))=\mathbb{E}((h(L))^2))-\mathbb{E}(h(L))^2$ ,

3. (iii) Skewness. $\varsigma_{h(L)}=\frac{\mathbb{E}\left[(h(L)-\mathbb{E}(h(L)))^3\right]}{(\textrm{Var}(h(L)))^{3/2}}$ ,

4. (iv) Value-at-Risk. $\textrm{VaR}_{p}(h(L))=\inf \{x\in \mathbb{R}\colon P(h(L)\leq x)\geq p\}$ ,

5. (v) Expected Shortfall. $\textrm{ES}_{p}(h(L))=\frac{1}{1-p}\int_p^1 \textrm{VaR}_q(h(L))\,\textrm{d}q$ .

These functionals allow also the computation of insurance premiums on the basis of premium principles such as the expectation principle, the variance principle, or the standard deviation principle.

2.3. Traffic scenarios in SUMO

3.1. Monte Carlo methods

The Monte Carlo simulation of L requires sampling the number of accidents $C^k$ in either the binomial or Poisson model and sampling the independent conditional losses $X^k_c$ , $c\in \mathbb{N}$ , for each traffic scenario $k=1,2, \dots, K$ from the corresponding distribution $F^k$ . These tasks can be performed separately but require both a prior evaluation of the microscopic traffic model.

1. (i) Prior Evaluation of Traffic Model. For each traffic scenario k, a single run delivers data that are the basis for a computation of, respectively, the accidents probabilities $p^k$ and intensities $\lambda^k$ as well as the corresponding values $p^k_r$ and $\lambda^k_r$ on the level of the modules $r=1,2, \dots, R$ .

2. (ii) Number of Accidents. Sampling from $\mu$ and using the results of the prior evaluation allow to sample the number of accidents $C^k$ for each traffic scenario k in both the binomial and Poisson model.

3. (iii) Conditional Accident Losses. Based on the precomputed values of the accident probabilities and the accident intensities, respectively, we simulate for each traffic scenario k the random locations and times of accidents. These data can be stored. For these locations and times, traffic data $\psi$ are extracted from an additional single run of traffic scenario k. Given $\psi$ , the losses are generated according to conditional loss distributions $F^{k, \psi}$ as described in Section 2.2. Details of the implementation are explained in Section 4.

These Monte Carlo methods can be flexibly applied to all considered functionals. In the special cases of expectation and variance, Wald’s equation can simplify the computation, since the losses L are given in the form of a collective model.

3.2. Gaussian approximation

Another way to compute $\mathbb{E}(h(L))$ for the considered types of insurance coverage is a Gaussian approximation, possibly improved by a correction term. A Gaussian approximation can easily be motivated within the binomial model. For N sufficiently large, the distribution of L given $\mu$ is approximately normal, implying that L is a mean-variance mixture of Gaussian distributions. This is an important structural insight from this approximation.

In this section, we condition on $\mu$ , that is, suppose that $\mu$ is fixed and given. The general results for $\mu$ random are then a corollary by considering suitable mixtures according to the distribution of $\mu$ . The total random losses in the binomial model can be rewritten as:

\begin{equation*}L\; = \; \sum_{k=1}^K \; \sum_{c=1}^{N\mu^k} \; \textbf{1}^k_c \cdot X^k_c \end{equation*}

where the random variables $\textbf{1}^k_c$ , $X^k_c$ , $k=1,2, \dots, K$ , $c\in \mathbb{N}$ , are independent, $X^k_c \sim F^k$ , and $\textbf{1}^k_c$ are Bernoulli random variables taking the value 1 with probability $p^k$ and the value 0 otherwise, $c\in\mathbb{N}$ , for any k. Setting

\begin{equation*}Y^k_c \,:\!=\, \textbf{1}^k_c \cdot X^k_c, \quad m^k \,:\!=\, \mathbb{E}\left(Y^k_c\right), \quad {\left(\sigma^k\right)}^2 \,:\!=\, \mathbb{E}\left(\left[Y^k_c - m^k \right]^2\right), \quad {\left(\zeta^k\right)}^3 \,:\!=\, \; \; \mathbb{E}\left(\left[Y^k_c - m^k \right]^3\right),\end{equation*}

$k=1,2, \dots, K$ , $c\in \mathbb{N}$ , a classical normal approximation of L is $\sum_{k=1}^K N\mu^k m^k + Z $ with

\begin{equation*}Z \quad \sim \quad \mathcal{N} \left( 0, \; \sum_{k=1}^K N\mu^k {\left(\sigma^k\right)}^2 \right).\end{equation*}

We focus on three types of insurance coverage, $h(x)=x$ (full coverage), $h(x)=\max(x-\theta,0)$ , $\theta\geq 0$ (constant deductible), $h=\min(x,\theta)$ , $\theta\geq 0$ (stop loss), and obtain an approximation $\mathbb{E}(h(Z))$ of $\mathbb{E}(h(L))$ in each of these cases.

On the basis of Stein’s method (see Chen et al., Reference Chen, Goldstein and Shao2011 and Ross, Reference Ross2011 for an overview), El Karoui and Jiao (Reference El Karoui and Jiao2009) suggest correction terms in order to improve the approximation, that is, the approximation $\mathbb{E}\left(h\left(\sum_{k=1}^K N\mu^k m^k + Z\right)\right)$ is replaced by the corrected approximation:

\begin{equation*}\mathbb{E}\left(h\left(\sum_{k=1}^K N\mu^k m^k + Z\right)\right)\; + \; C_h. \end{equation*}

The correction $C_h$ depends on the degree of smoothness of the derivatives of the function h and thus differs (see Theorem 3.1 and Proposition 3.6 in El Karoui and Jiao, Reference El Karoui and Jiao2009) for the three types of coverage. We define

\begin{equation*}d_1 \; = \; \sum_{k=1}^K N\mu^k m^k, \quad d_2 \;= \; \sum_{k=1}^K N\mu^k {\left(\sigma^k\right)}^2, \quad d_3 \; = \; \sum_{k=1}^K N\mu^k \left(\zeta^k\right)^3, \quad \tilde h (x) \; = \; h\left( x + d_1 \right) .\end{equation*}

We obtain the following correction terms:

1. (i) Full Coverage. In the case of full coverage, the correction term of El Karoui and Jiao (Reference El Karoui and Jiao2009) disappears. In general, if h is some Lipschitz function with bounded third derivative, the correction term equals $C_h \; = \; \dfrac{d_3}{2 d_2^2} \cdot \mathbb{E}\left( \left\{ \dfrac{Z^2}{3 d_2} - 1 \right\} Z \tilde h(Z) \right).$

2. (ii) Constant Deductible. $C_h \; = \; \dfrac{ (\theta - d_1)d_3 }{6 d_2} \cdot \dfrac{1}{\sqrt{2 \cdot \pi \cdot d_2}} \cdot \exp \left\{ - \dfrac{(\theta - d_1)^2}{2 \cdot d_2}\right\}$

3. (iii) Stop Loss. A stop loss $x\mapsto\min\!(x,\theta)$ can be written as the difference between full coverage $x\mapsto x$ and a constant deductible $x\mapsto \max(x-\theta,0)$ . This implies that the correction term for a constant deductible appears with a negative sign in this case.

The advantage of the (corrected) Gaussian approximation in comparison with pure Monte Carlo is that, once the numbers $m^k$ , ${\left(\sigma^k\right)}^2$ , and $\left(\zeta^k\right)^3$ have been computed for each traffic scenario $k=1,2, \dots, K$ , no further data need to be stored or sampled in order to compute $\mathbb{E}(h(L))$ . The approximate representation of the distribution of L as a mean-variance mixture is a considerable simplification.

4.1. SUMO scenario and accident data

Wildau is a small German city of approximately 10,000 inhabitants, located around 30 km southeast of the capital Berlin. A SUMO model of the city was developed within a study project by the Technical University of Applied Sciences Wildau and is publicly available (see github.com/DLR-TS/sumo-scenarios/tree/main/Wildau).

SUMO Scenario. The implemented road network is visualized in Figure 1. It is specified using 646 nodes connected by 1426 edges. The city itself is crossed by the railway; the tracks are represented by the gray line. Vehicles in the present scenario are calibrated from real traffic counts provided by local authorities, see also Behrisch and Hartwig (Reference Behrisch and Hartwig2022). Such numbers of vehicles per time unit for certain positions are then assigned to different routes through the city. The original scenario has a duration of $7010\,\textrm{s}$ . Empty in the beginning, vehicles enter the system with a peak of approximately 240 vehicles that drive simultaneously. In total, 2502 vehicles are generated.

In the following section, we describe in detail how we adjust this SUMO scenario to obtain suitable ingredients for our case studies. This yields a collection of traffic scenarios that allow us to compare the effects of different driving characteristics and fleet sizes on the total loss and related insurance premiums.

Varying Traffic Conditions. For different collections of model parameters representing different traffic systems we generate adjusted SUMO scenarios. For each choice of parameters, we proceed as follows to produce $K=100$ traffic scenarios of length $T=60\,\textrm{s}$ . Traffic scenarios $k=1,2, \dots, 50$ correspond to selected time intervals from the SUMO scenario. Traffic scenarios $k= 51, 52, \dots, 100$ represent higher traffic volumes. They are generated by replacing the original route file by a route file that consists of two copies of the original route file. This simple procedure generates a larger amount of vehicles along the original paths. The traffic scenarios are again selected time intervals from the corresponding SUMO scenario.

To represent the full year, we set $N=365\cdot 24\cdot 60=525,600$ . We need to specify the random vector $\mu=(\mu^1,\dots,\mu^K)^\top$ describing the number of occurrences of the individual traffic scenarios divided by N. For the purpose of illustration, we specifically assume that two probability measures $\nu_g$ , $\nu_b$ are given on $\{1,2, \dots, K \}$ which approximately correspond to the relative frequencies of traffic scenarios in two prototypical years $y=g,b$ . In addition, we suppose that the type y of the current year is random where both values g and b have probability $1/2$ . Given y, we generate $\mu=(\mu^1,\dots,\mu^K)^\top$ from a multinomial distribution corresponding to $\nu_y$ . That is, for all time buckets $n=1,2, \dots, N$ , a traffic scenario k is chosen independently from the distribution $\nu_y$ on $\{1,2, \dots, K\}$ . Dividing the number of occurrences of a scenario k by N, one obtains its random relative frequency $\mu^k$ for any $k=1,2, \dots, K$ . In our case study, the distribution $\nu_g$ corresponds to lower traffic densities on average, while $\nu_b$ is associated with higher traffic densities, that is, we set

\begin{equation*}\nu_g \,:\!=\, \left\{\begin{array}{l@{\quad}l}1/75, & k=1,2, \dots, 50,\\[5pt]1/150, & k =51,52, \dots, 100\end{array}\right., \quad \nu_b \,:\!=\, \left\{\begin{array}{ll}1/150, & k=1,2, \dots, 50,\\[5pt]1/75, & k =51,52, \dots, 100.\end{array}\right.\end{equation*}

Accident Data. The German Accident Atlas provided by Statistical Offices of the Federation and the Länder (2022) depicts the locations of all police-reported accidents involving personal damage that occurred within 1 year. In 2020, within the modeled area of Wildau (approximately) 48 accidents were registered. There are also aggregate statistics for Germany for all police-reported accidents. In 2020, approximately $11.8\,{\%}$ of all road accidents involved personal damage (see Statistisches Bundesamt, 2023). We use an estimate of $\bar{c}_{\textrm{year}}=48/11.8{\%}\approx 407$ accidents for calibration purposes.

4.2. Model specification

Our goal is to analyze accident losses for a fleet $\Phi$ over the time horizon of 1 year.

Fleet Definition. In the Wildau scenario, vehicles are defined using repeated flows from origins to destinations. Passenger cars (next to trucks and the train) are defined via 90 different flows. A flow generates vehicles of a given type at a given position. From this position, they navigate to a specified destination.

In SUMO, a passenger car represents a certain vehicle type. Initially, all passenger cars belong to the same vehicle type and, consequently, have the same driving characteristics. To introduce a fleet $\Phi$ of vehicles whose driving characteristics we can vary, we define a new vehicle type $\Phi$ and construct corresponding new SUMO scenarios. Fixing a fraction $\rho^\Phi\in[0,1]$ of vehicles belonging to $\Phi$ , we retain approximately $1-\rho^\Phi$ of the existing flow definitions and modify $\rho^\Phi$ of the flow definitions suitably in order to model the fleet. In our case studies, we consider $\rho^\Phi = 10 {\%},\, 50{\%},\, 90 {\%}$ .

Driving Configuration. Vehicles in a fleet $\Phi$ are of the same type. Various characteristics can be varied in SUMO; we focus on maximal speed $v_\textrm{max}$ , maximal acceleration $a_\textrm{max}>0$ , and time headway $\zeta>0$ . The time headway is the distance which is kept to the preceding vehicle measured in time, that is, a velocity-weighted safety distance. We refer to a fixed selection of driving characteristics as a driving configuration.

In our case studies, we will vary the driving configuration for all vehicles in the fleet $\Phi$ and keep all other vehicles as originally introduced; we use the implementation of an Intelligent Driver Model without any speed deviation that does not include further random effects. A driving configuration of vehicles in fleet $\Phi$ is denoted by $\xi=(v_\textrm{max},a_\textrm{max},\zeta)$ . Specifically, we consider

\begin{align*}\xi^{1a}&=\left(5\,\textrm{m/s},\, 0.8\,\textrm{m/s}^2,\, 3.0\,\textrm{s}\right),\,\xi^{2a}=\left(10\,\textrm{m/s},\,0.8\,\textrm{m/s}^2,\, 2.0\,\textrm{s}\right),\, \xi^{3a}=\left(15\,\textrm{m/s},\,0.8\,\textrm{m/s}^2,\,1.0\,\textrm{s}\right), \\\xi^{1b}&=\left(5\,\textrm{m/s},\,2.6\,\textrm{m/s}^2,\,3.0\,\textrm{s}\right),\,\xi^{2b}=\left(10\,\textrm{m/s}, 2.6\,\textrm{m/s}^2,\, 2.0\,\textrm{s}\right),\, \xi^{3b}=\left(15\,\mathrm{m/s},\, 2.6\,\textrm{m/s}^2,\, 1.0\,\textrm{s}\right).\end{align*}

The configurations 1,2,3 increase in terms of “aggressiveness” from driving slowly with a large headway to fast with a small headway – with two options a and b for the maximal acceleration.

Accident Occurrence. The best estimate for the total number of accidents in Wildau is $\bar{c}_\textrm{year} \approx 407$ . From this, we derive a uniform and a nonuniform accident occurrence model. In both cases, we specify accident probabilities $p^{1,k}$ and $p^{2,k}$ for the binomial model as well as accident intensities $\lambda^{1,k}$ and $\lambda^{2,k}$ for the Poisson model.

1. (i) Uniform Accident Occurrence. Assuming that accidents occur uniformly over the year, we obtain a probability per time bucket of an accident in the system of $p^1 =\bar{c}_\textrm{year}/N \approx 7.7 \times 10^{-4}$ . This is the accident probability that we allocate to each traffic scenario k. We also suppose that accidents occur uniformly across all vehicles in the system. This implies that the probability that any accident occurs in scenario k within the fleet $\Phi$ is

   \begin{equation*} p^{\Phi, 1,k} \; = \; \rho^\Phi\cdot p^1, \quad k\in\{1,\dots,K\}. \end{equation*}
   This probability is used in the binomial model. For the Poisson model, we set $\lambda^{\Phi,1,k}=p^{\Phi, 1,k}$ , $k=1,2, \dots, K$ , since the intensity approximately equals the probability of an accident per time bucket.

   In the case of uniform accident occurrence, we do not consider any spatial variations of the likelihood of accidents due to different traffic conditions. This means that we do not distinguish any modules, that is, we set $R=1$ .

2. (ii) Non-Uniform Accident Occurrence. In reality, the likelihood of accidents depends on external factors such as weather and local traffic conditions, for example, the velocity of vehicles and traffic density. The quantities vary spatially and over time.

   From SUMO runs, we obtain for each traffic scenario $k=1,2, \dots, K$ and each module $r=1,2, \dots, R$ pairs $(d_r^k,\bar{v}_r^k)$ on the average density and velocity. These statistics can be computed in SUMO, for example, from data that are obtained at induction loop detectors which are placed within the modules; as a proxy for density, we extract the occupancy of the loop detector, that is, the fraction of time which it is occupied by a vehicle.

   For $r=1,\dots,R$ , we choose benchmark values $d^*_r$ and $\bar{v}^*_r$ for the density and velocity and specify occurrence probabilities and intensities that vary spatially and over time:

   \begin{equation*} \lambda^{\Phi, 2,k}_r \; = \; p^{\Phi, 2,k}_r \,:\!=\, \frac{p^{\Phi, 1,k}}{R} \cdot \frac{\bar{v}_r^k}{\bar{v}_r^*}\cdot\frac{d_r^k}{d_r^*}\cdot e^{-(\zeta^\Phi-1)}, \quad k\in\{1,\dots,K\},\quad r\in\{1,\dots,R\}.\end{equation*}
   The last term refers to deviations of the time headway from SUMO’s default value of $1.0\,\textrm{s}$ : a larger time headway is associated with less risky driving. We set $p^{\Phi, 2,k}=\sum_{r=1}^R p^{\Phi, 2,k}_r$ and $\lambda^{\Phi, 2,k}=\sum_{r=1}^R \lambda^{\Phi, 2,k}_r$ .

   In our case studies, we will consider a grid of $R=4$ modules and compute $d_r^k$ and $\bar{v}_r^k$ as averages over measurements from 10 induction loop detectors that are placed in each module (see also Figure 2). We use the scenario averages $d^*_r=\frac{1}{K}\sum_{k=1}^K d_r^k$ and $\bar{v}_r^*=\frac{1}{K}\sum_{k=1}^K \bar{v}_r^k$ . If $\sum_{k=1}^K \mathbb{E}(\mu^k) \bar{v}_r^k = \bar{v}_r^* $ , $\sum_{k=1}^K \mathbb{E}(\mu^k) d_r^k = d_r^*$ , and $\zeta^\Phi =1$ , then we essentially recover on average the case of uniform accident occurrence.

Figure 2. Partition of Wildau and placement of induction loop detectors.

Accident Losses. The distributions $F^k$ of accident losses associated with a traffic scenario $k = 1,2, \dots, K$ are constructed on the basis of traffic data that are extracted from the SUMO runs. The general procedure was described in Section 2.2; here, we explain the specific implementation that we use in our case studies.

The likelihood of accident occurrence was discussed in the previous section. $F^k$ is the conditional distribution of an accident loss in traffic scenario k if an accident occurs. Time in traffic scenario k is enumerated by $t\in [0,T]$ , and we assume that the time $\tau$ of the accident conditional on its occurrence is uniformly distributed on [0, T], that is, $\tau\; \sim \; \textrm{Unif}[0,T].$ We choose a module $\mathcal{R}$ at random in which the accident occurs and assume, respectively, that

\begin{equation*}P(\mathcal{R}=r) \; = \; \frac{p^{\Phi, \cdot, k}_r}{p^{\Phi, \cdot, k}}, \quad\quad P(\mathcal{R}=r)=\frac{\lambda^{\Phi, \cdot, k}_r}{\lambda^{\Phi, \cdot, k}}, \quad \quad\quad r\in\{1,\dots,R\}. \end{equation*}

These ratios depend in the case of nonuniform accident occurrence on the specific fleet, since the properties of the fleet alter the route file that is used to generate the SUMO scenarios; this is true, although the multiplicative terms $\rho^\Phi$ appear in both the numerator and denominator and cancel out.

In the chosen module, we pick one or more vehicles at random and extract from the traffic scenario data for these vehicles. In our concrete implementation, we simply choose at time $\tau$ a single vehicle I uniformly at random in module $\mathcal{R}$ , that is, its conditional distribution is

\begin{equation*}I\mid \tau,\mathcal{R} \; \sim \; \textrm{Unif}\!\left(\mathcal{M}_\mathcal{R}^\Phi(\tau)\right). \end{equation*}

For the purpose of illustrating our approach, the only data we extract are the velocities $v^I$ of the randomly chosen vehicles that are involved in accidents. We set $\psi= v^I$ and assume that the conditional loss distribution $F^{k,\psi} $ is known (we assume that $F^{k,0}$ corresponds to a Dirac measure in 0; if $\mathcal{M}_\mathcal{R}^\Phi(\tau)=\emptyset$ , we set $\psi=0$ , resulting in 0 losses). If we denote the distribution of $\psi$ by $\mathcal{L}^{k}$ , we obtain the distribution $F^k$ as a mixture:

\begin{equation*} F^k \; = \; \int F^{k, \psi} \;d \mathcal{L}^{k} .\end{equation*}

In our case studies, we will assume that $F^{k, \psi} = F^\psi$ for all k; however, the mixing distribution $\mathcal{L}^{k}$ will depended on the traffic scenario k. We consider the following examples for $F^\psi$ :

1. (i) Gamma Distribution. We define distributions with varying levels of dispersion. A measure for the dispersion of a random variable X is the coefficient of variation defined by $c_v=\sqrt{\textrm{Var}(X)}/\mathbb{E}(X)$ . For $c_v\in\{1/2,\,1,\,2\}$ , we choose $F^\psi \; = \; \Gamma \left( \frac{1}{c_v^2}, \; \frac{1}{c_v^2 \psi^2} \right).$ The expectation of this distribution is $\psi^2$ and increases quadratically with $\psi$ , the velocity of the vehicle involved in an accident; this is consistent with the fact that losses scale with kinetic energy. The variance of the distribution $F^\psi$ equals $c_v^2\psi^4$ ; hence, the coefficient of variation is indeed $c_v$ .

2. (ii) Log-Normal Distribution. We consider log-normal distributions with expectation $\psi^2$ and variance $c_v^2 \psi^4$ , implying that the coefficient of variation is again $c_v\in\{1/2,\,1,\,2\}$ . This log-normal distribution is obtained as the distribution of $\exp\!(Z)$ for a normal random variable Z with expectation $\ln(\psi^2/\sqrt{1+c_v^2})$ and variance $\ln\left(1+c_v^2\right)$ , that is,

   \begin{equation*}F^\psi \; = \; \mathcal{LN}\left( \ln \left( \frac{\psi^2}{\sqrt{1+c_v^2}} \right), \; \ln(1+c_v^2) \right) .\end{equation*}

4.3. Case studies

4.3.1. Overview

We illustrate our modeling approach in case studies on multiple levels. A selection of case studies is discussed in detail in Sections 4.3.2 & 4.3.3. All numerical results for the following choices are documented online in tables in Appendix F:

1. (i) Fleet Models. We analyze six driving configurations with three different fleet proportions.

2. (ii) Accident Occurrence. In our traffic system, accidents occur uniformly or nonuniformly in space. Their number is given by a binomial or a Poisson model with parameters depending on traffic condition.

3. (iii) Accident Losses. We study two parametric families of loss distributions with three different choices for the coefficient of variation.

4. (iv) Insurance Design. Insurance losses are a function of the total losses; we distinguish three contract designs.

Denoting aggregate losses by L, we evaluate for each type of insurance coverage h the resulting insurance losses h(L) in terms of their expectation, variance, and skewness, and the monetary risk measures Value-at-Risk and Average Value-at-Risk, also called Expected Shortfall. To analyze the distributions in detail, we provide qq-plots and estimates of cumulative distribution functions (CDFs) and densities. The main tool to access the random variable h(L) is Monte Carlo sampling; we provide a pseudocode how we obtain samples of L in Algorithm 1 in Appendix E. In Section 4.3.3, we compare this approach to the normal mean-variance mixture approximation introduced in Section 3.2.

Figure 3. Impact of fleet size and driving configuration on the total loss L (solid lines represent normalized values (left y-axis) and dashed lines unnormalized ones (right y-axis)).

This paper explores the analysis and management of risks that occur in vehicle fleets in traffic systems. We distinguish to perspectives:

1. (i) The Engineering Perspective. In these case studies, we fix the accident occurrence and accident loss distributions and vary the fleet models, that is, driving configurations and fleet proportions. We focus on $\mathbb{E}(L)$ , $\textrm{Var}(L)$ , and complement these with analyses of the performance of traffic system.

2. (ii) The Actuarial Perspective. In these case studies, we fix the fleet model and vary accident occurrence and accident loss distributions as well as the insurance design. We study the distribution of L and the insurance prices $\mathbb{E}(h(L))$ .

4.3.2. The engineering perspective

Our micro-modeling approach allows us to study the effects of different traffic-related controls on total losses L; we investigate the effects of fleet size and traffic configuration. Throughout this section, we consider nonuniform accident occurrence in the binomial model with Gamma-distributed accident losses and a coefficient of variation $c_v=1$ .

Losses. We evaluate expected loss $\mathbb{E}(L)$ and standard deviation $\textrm{std}(L)$ for different fleet models. To compare losses for different fleet sizes, we normalize losses per 100 expected insured vehicles: for each traffic scenario k, the number of insured vehicles is the number of vehicles belonging to $\Phi$ as given in the underlying route file. We refer to Appendix F for more details. The model specification was explained in Section 4.2 which includes in particular a description of the driving configurations. The results are documented in Figure 3. The solid lines are the normalized quantities.

In Figure 3(a), we see that increasing the aggressiveness of driving increases both the total and normalized expected loss. An impact of the maximal acceleration on losses is only substantial for the most aggressive driving configurations $\xi^{3\cdot}$ . Increasing the fleet size increases the expected loss which is primarily due to the fact that we count losses only within the fleet and a higher volume is associated with higher losses. More interesting is the normalized case: apparently higher speeds also increase the normalized losses.

In Figure 3(b), we have the corresponding standard deviations. Increasing the aggressiveness of driving increases the standard deviation of the total loss. The standard deviations of the normalized losses are decreasing in the fleet size. The main reason is that fluctuations normalized for a fixed volume are larger for smaller pools than for larger pools; a rational for this is provided by the law of large numbers and the central limit theorem.

The frequency and severity of accidents are, of course, increasing in the aggressiveness of driving. To demonstrate this, we evaluate the expectation and the standard deviation of the average accident frequency $\sum_{k=1}^K\mu^kp^k$ (both normalized and unnormalized) and the average accident severity $\sum_{k=1}^K \mu^k X_1^k$ , as displayed in Figures 4 and 5. A larger fleet increases frequency and, in aggressive scenarios, also the expected accident severity. This is, of course, due to the specific choice of the driving behavior of the considered fleets in comparison with the driving behavior of the remaining vehicles and does not necessarily hold for all traffic systems in general. In this section, we discussed losses for selected special cases. A comprehensive set of tables for all other cases and statistical functionals is provided online in Appendix F.

Figure 4. Impact of fleet size and driving configuration on expectation and standard deviation of average accident frequency (solid lines represent normalized values (left y-axis) and dashed lines unnormalized ones (right y-axis)).

Changing the characteristics of the fleet not only affects the losses. At the same time, this has an impact on the performance of the traffic system. We present further results in Appendix C.

4.3.3. The actuarial perspective

From an actuarial perspective, it is relevant to understand the risk that corresponds to the insurance losses. This requires a more detailed analysis of the probability distributions. To do this, we pick a particular fleet model and use probabilistic techniques to evaluate the distribution of h(L). From now on, we consider $\rho^\Phi=0.5$ with driving configuration $\xi^{2a}$ and nonuniform accident occurrence.

Distributional Analysis of Losses. We start our investigations with the total losses L. Table 1 shows the evaluation of statistical functionals for different accident losses in the case of the binomial model. These numbers quantify the risk entailed in the total losses. Both the distributional family and the chosen coefficient of variation for the accident loss model have a substantial effect on the risk.

Table 1. Statistical Functionals of L for $\rho^\Phi=0.5$ and $\zeta^{2a}$ .

The statistical functionals of the total loss are approximated using 10,000 independent samples of L.

Figure 5. Impact of fleet size and driving configuration on expectation and standard deviation of average accident severity.

A visual impression of the distributions is provided in Figure 6. We compare different accident loss models while fixing the coefficient of variation $c_v=2$ . We plot the empirical distribution functions as estimates of the CDF and a kernel density estimate of the corresponding densities. Moreover, Figure 6(c) shows qq-plots for the quantiles of standardized values of the losses against quantiles of a standard normal distribution (samples are standardized by subtracting their sample average and dividing by their sample standard deviation).

Figure 6. Distribution of the total loss for fixed coefficient of variation $c_v=2$ .

We find that the binomial and the Poisson model do not differ too much. Yet, log-normal accident losses produce heavier tails than the corresponding Gamma losses. The qq-plots reveal that the right tails are heavier compared to a normal distribution, while the left tails are lighter. The latter observation simply relates to the fact that the original losses are nonnegative, while the normal distribution takes values on the whole real line.

In Figure 7, we analyze the impact of the coefficient of variation while fixing the log-normal distribution for the accident losses. We see again that binomial and Poisson model do not differ substantially. However, the effect of the coefficient of variation is clearly visible: increasing $c_v$ produces heavier right tails. Introducing dispersion to accident losses substantially changes the distribution of the total losses.

Figure 7. Distribution of the total loss for log-normal accident losses and varying coefficient of variation.

Comparison of Losses and Normal Mean-Variance Mixture Approximation. In Section 3.2, we suggested a normal mean-variance mixture approximation for the total loss. To study the quality of this approximation, we generate 10,000 samples from the approximation; in the following, we focus on the case of Gamma-distributed accident losses with coefficient of variation $c_v=1$ .

To sample from the approximation, we rely on the following computations of $m^k$ and $(\sigma^k)^2$ . Using $\mathbb{E}(\textbf{1}_c^k)=p^k$ and $\mathbb{E}(X_c^k) = \mathbb{E}(\mathbb{E}(X_c^k\mid \psi)) = \mathbb{E}(\psi^2) = \int \psi^2 d \mathcal{L}^k $ , we obtain for the Gamma losses $m^k =p^k\cdot \int \psi^2 d \mathcal{L}^k $ and $(\sigma^k)^2 = p^k\cdot c_v^4\cdot \int \psi^4 d \mathcal{L}^k \cdot \left(1+\frac{1}{c_v^2}\right)\frac{1}{c_v^2}-(m^k)^2$ .

Remark 4.3. For the notation, we refer to Sections 3.2 & 4.2. The computations are valid for any coefficient of variation $c_v$ , but we use only $c_v=1$ in the numerical case study.

The involved moments of $\psi$ are approximated using 10,000 samples from the traffic simulation, for each $k=1,\dots,K$ . A sample from the normal mean-variance mixture approximation is generated by, first, sampling $\mu$ and, second (conditional on $\mu$ ), sampling the normal random variable $\sum_{k=1}^K N\mu^k m^k + Z $ with $Z \sim \mathcal{N} \left( 0, \; \sum_{k=1}^K N\mu^k {\left(\sigma^k\right)}^2 \right).$

Figure 8 shows the qq-plot comparing quantiles of the crude Monte Carlo simulation with quantiles of the approximation. This demonstrates the quality of our suggested approximation. It is almost exact between the 5% and 95% quantile as the values lie on the half-line. It is still very good for the 1–5% and 95–99% quantiles and is only less accurate in the extreme tails where also in the Monte Carlo simulation only few data points are available. These analyses, on the one hand, confirm the postulated structural model insight. On the other hand, they also validate the implementation of our crude Monte Carlo sampling.

Figure 8. QQ-Plot of 10,000 Monte Carlo samples (y-axis) versus 10,000 samples of mixture approximation (x-axis) (the red dots mark the $1\%,\,5\%,\,10\%,\,90\%,\,95\%,\,99\%$ quantiles).

Pricing and Evaluation Methods. To conclude our case studies, we study prices for various insurance contracts. We compare $\mathbb{E}(L)$ (full coverage), $\mathbb{E}(\!\max\!(L-\theta,0))$ (constant deductible), and $\mathbb{E}(\min(L,\theta))$ (stop loss) for different values of $\theta$ . The results are given in Figure 9. We obtain the typical hockey stick profiles satisfying the parity $\mathbb{E}(L)=\mathbb{E}(\!\max\!(L-\theta,0))+\mathbb{E}(\min(L,\theta))$ . We note that other insurance contracts can easily be represented in our framework; also deductibles per accident can be implemented by changing the accident loss distributions accordingly.

Figure 9. Insurance prices.

Besides Monte Carlo methods, our normal mean-variance mixture approximation, and the correction suggested in Section 3.2 provide alternative techniques to compute the prices $\mathbb{E}(h(L))$ for the different types of coverage h. To compute the correction term $C_h$ , we need to compute $(\zeta^k)^3$ (for the notation we refer to Sections 3.2 & 4.2). For Gamma-distributed losses, we obtain

\begin{equation*}(\zeta^k)^3=p^k\cdot c_v^6\cdot \int \psi^6 d \mathcal{L}^k \cdot \left(2+\frac{1}{c_v^2}\right)\left(1+\frac{1}{c_v^2}\right)\frac{1}{c_v^2} - 3m^k(\sigma^k)^2-(m^k)^3.\end{equation*}

Since $\mathbb{E}(h(L))=\mathbb{E}(\mathbb{E}(h(L)\mid\mu))$ , we may generate samples of $\mu$ and evaluate the corresponding conditional expectations $\mathbb{E}(h(L)\mid\mu)$ ; in the normal mean-variance mixture approximation, these are expectations of functions of normally distributed random variables. We compute these expectations numerically as integrals with respect to Lebesgue measure using a normal density.

We compare the estimation errors of the different approaches in Figure 10. We produce 100,000 samples of L to approximate the “true” value of $\mathbb{E}(\!\max\!(L-\theta,0))$ of coverage with constant deductible, for different values of $\theta$ . Independently, we generate 10,000 samples and consider Monte Carlo approximations based on all 10,000 samples and based only on the first 1000 samples. We also study the normal mean-variance mixture approximation with and without correction using 1000 samples of $\mu$ .

Figure 10. Comparison of estimation errors.

While the absolute error generally decreases in the deductible $\theta$ for all methods (apart from some local effects), the relative error increases for larger values of $\theta$ which are associated with a lower price of the contract. At the same time, we observe that, in terms of the relative error, the normal mean-variance mixture approximation produces reasonable estimation results (compared to 1000 samples) for moderate values of $\theta$ . This is in line with our previous observations on the quality of our normal mean-variance mixture approximation which becomes worse in the extreme tail of L. Interestingly, the estimation error can largely be reduced using the correction; with the correction, the estimation becomes quite good even for large values of $\theta$ .