2.1. Signal representation

A regular DWT is defined by two filters: a wavelet filter $\{h_\ell\}_{\ell=0}^{L-1}$ and a scaling filter $\{g_\ell\}_{\ell=0}^{L-1}$ , where L is the number of terms of the filters. In this paper, we use the Daubechies D4 wavelet filter, with $L=4$ terms, given by:

\begin{align*}h_0=\frac{1-\sqrt{3}}{4\sqrt{2}}, \quad h_1=\frac{-3+\sqrt{3}}{4\sqrt{2}}, \quad h_2=\frac{3+\sqrt{3}}{4\sqrt{2}}, \quad h_3=\frac{-1-\sqrt{3}}{4\sqrt{2}}. \end{align*}

The corresponding scaling filter is obtained by

\begin{align*}g_\ell = (\!-1)^{\ell+1}\,h_{L-1-\ell},\qquad \ell=0,1,2,3, \end{align*}

that is,

\begin{align*}g_0=-h_3,\quad g_1=h_2,\quad g_2=-h_1,\quad g_3=h_0. \end{align*}

To see how these filters are utilized for signal representation, consider a telematics signal for a trip i:

\begin{align*}\mathbf{X}_i=(X_{i,0},\ldots,X_{i,T_i-1}). \end{align*}

At the first wavelet level, the MODWT wavelet and scaling filters are applied as local weighted sums of nearby signal values:

(1) \begin{equation}\widetilde W_{i,1,t}=\sum_{\ell=0}^{L-1} \frac{h_\ell}{\sqrt{2}}\,X_{i,(t-\ell)\,\textrm{mod}\, T_i},\qquad\widetilde V_{i,1,t}=\sum_{\ell=0}^{L-1} \frac{g_\ell}{\sqrt{2}}\,X_{i,(t-\ell)\,\textrm{mod}\, T_i},\end{equation}

where the filters $h_{\ell}, g_{\ell}$ are rescaled by $1/\sqrt{2}$ to ensure the variance contributions across levels are comparable. Without scaling, the contribution would be inflated simply because the MODWT produces coefficients at every time point. Here, the wavelet coefficient $\widetilde W_{i,1,t}$ measures a localized change in the signal around time t, while the scaling coefficient $\widetilde V_{i,1,t}$ captures the corresponding local trend. The modulo operation implements circular indexing so that the filter can be applied at every time point of the trip. Thus, the filters $h_\ell$ and $g_\ell$ serve as local scanning weights applied to the signal.

A key feature of the wavelet transform is its multi-resolution decomposition: it represents a time series using components that capture changes over different time scales. We index these scales by a level j, which corresponds to changes in the original time series on a scale of $2^{j-1}$ . At finer levels (small j), the transform captures short-lived, localized changes, while at coarser levels (large j), it captures smoother changes that persist over longer durations.

To account for the multi-scale view, we denote by $\{h_{j,\ell}\}_{\ell=0}^{L_j-1}$ and $\{g_{j,\ell}\}_{\ell=0}^{L_j-1}$ the filters at level j obtained from the DWT filters $\{h_\ell\}$ and $\{g_\ell\}$ by dyadic dilation (i.e., by inserting $2^{j-1}-1$ zeros between successive terms). For the MODWT, we denote the multi-scale filters as

(2) \begin{equation} \widetilde{h}_{j,\ell} = \frac{h_{j,\ell}}{\sqrt{2}},\qquad \widetilde{g}_{j,\ell} = \frac{g_{j,\ell}}{\sqrt{2}},\qquad \ell=0,\ldots,L_j-1,\end{equation}

where the filter width at level j is

\begin{equation*} L_j = (2^j-1)(L-1)+1,\end{equation*}

that is, the number of lags of the original signal references when computing MODWT coefficients at level j. The level-1 filtering formula in (1) generalizes naturally to multiple wavelet levels. At level j, the MODWT uses the rescaled filters $\{\widetilde h_{j,\ell}\}_{\ell=0}^{L_j-1}$ and $\{\widetilde g_{j,\ell}\}_{\ell=0}^{L_j-1}$ defined in (2). The corresponding wavelet and scaling coefficients are

(3) \begin{equation}\widetilde{W}_{i,j,t} \, :\!= \, \sum_{\ell=0}^{L_j-1}\widetilde{h}_{j,\ell}\,X_{i,(t-\ell)\,\textrm{mod}\, T_i},\qquad\widetilde{V}_{i,j,t} \, :\!= \, \sum_{\ell=0}^{L_j-1}\widetilde{g}_{j,\ell}\,X_{i,(t-\ell)\,\textrm{mod}\, T_i},\qquad t=0,\ldots,T_i-1.\end{equation}

The same circular indexing convention allows each level-j filter to scan the entire trip without downsampling.

For any integer $J\ge 1$ , the vectors have the following interpretations:

• • a vector $\widetilde{\mathbf{W}}_{i,j}$ contains the MODWT wavelet coefficients associated with changes in $\mathbf{X}_i$ on a scale of $2^{j-1}$ for $j = 1, \dots, J$ ,

• • and a vector $\widetilde{\mathbf{V}}_{i,J}$ contains the MODWT scaling coefficients associated with changes in $\mathbf{X}_i$ on a scale $2^{J}$ or higher,

where $\widetilde{\mathbf{W}}_{i,j} = (\widetilde{W}_{i,j,t})_{t=0}^{T_i-1}$ and $\widetilde{\mathbf{V}}_{i,J} = (\widetilde{V}_{i,J,t})_{t=0}^{T_i-1}$ (Percival and Walden Reference Percival and Walden2013) and J is the chosen maximum level of MODWT decomposition. Since the filter width $L_j$ increases with j, the coefficients at finer levels represent localized features and coarser levels represent milder features with longer duration. This multi-resolution structure across levels can capture both brief and sustained patterns, all of which are informative driving patterns.

The wavelet filters and coefficients in (3) provide a complete representation of the original time series $\mathbf{X}_i$ . For each level j, we extend the filters $\{\widetilde{h}_{j,\ell}\}_{\ell=0}^{L_j-1}$ and $\{\widetilde{g}_{j,\ell}\}_{\ell=0}^{L_j-1}$ to length $T_i$ by setting $\widetilde{h}_{j,\ell}=0$ and $\widetilde{g}_{j,\ell}=0$ for $\ell=L_j,\ldots,T_i-1$ . Then, we represent $\mathbf{X}_i$ using

\begin{equation*} X_{i,t} \;=\; \sum_{j=1}^{J}\ \sum_{\ell=0}^{T_i-1}\widetilde{h}_{j,\ell}\,\widetilde{W}_{i,j,(t+\ell)\,\textrm{mod}\, T_i} \;+\; \sum_{\ell=0}^{T_i-1}\widetilde{g}_{J,\ell}\,\widetilde{V}_{i,J,(t+\ell)\,\textrm{mod}\, T_i},\end{equation*}

for a given $J\ge1$ . Under the circular indexing in (3), this representation holds for all $t = 0, \dots, T_i-1$ , so $\mathbf{X}_i$ can be reconstructed from $\{\widetilde{\mathbf{W}}_{i,j}\}_{j=1}^{J}$ and $\widetilde{\mathbf{V}}_{i,J}$ . This wavelet transform is particularly useful in UBI settings, where we need to focus on signal-based analysis rather than sparse claim counts, because this representation allows us to map each trip to a coefficient series that retains the full signal information and serves as the foundation for our proposed framework.

2.2. Energy decomposition and variance contributions

As the level j increases, coefficients beyond a certain level tend to represent only broad baseline movement, which is around zero, rather than meaningful risk-relevant patterns in a telematics signal. For this reason, choosing the maximum level of the decomposition J is a natural question to select a meaningful set of levels. We therefore quantify the variance contribution of each level j to trip i using the MODWT energy decomposition. For a fixed integer $J\ge1$ , Percival and Walden (Reference Percival and Walden2013) proved that the MODWT yields the energy decomposition

\begin{equation*}\|\mathbf{X}_i\|^2 \;=\; \sum_{j=1}^{J}\|\widetilde{\mathbf{W}}_{i,j}\|^2 \;+\; \|\widetilde{\mathbf{V}}_{i,J}\|^2,\end{equation*}

because we have

\begin{align*} \|\widetilde{\mathbf{W}}_{i,j}\|^2 + \|\widetilde{\mathbf{V}}_{i,j}\|^2 = \|\widetilde{\mathbf{V}}_{i,j-1}\|^2, \qquad j=1, 2, \dots \end{align*}

with the definition $\widetilde{\mathbf{V}}_{i,0}\, :\!= \,\mathbf{X}_i$ and $\|\cdot\|$ denotes the $\ell_2$ norm.

To quantify contributions for each level, we define a sample variance of $\mathbf{X}_i$ by

\begin{align*}\widehat{\sigma}^2_{\mathbf{X}_i} \;\, :\!= \,\; \frac{1}{T_i}\|\mathbf{X}_i\|^2-\overline{X}_i^{\,2}, \end{align*}

where $\overline{X}_i \, :\!= \, T_i^{-1}\sum_{t=0}^{T_i-1} X_{it}$ . For each level j, define the energy component

\begin{align*}\mathcal{V}_i(j) \;\, :\!= \,\; \frac{1}{T_i}\bigl\|\widetilde{\mathbf{W}}_{i,j}\bigr\|^2, \end{align*}

yielding the variance contribution given by

\begin{equation*}\rho_{i,j} \;\, :\!= \,\; \frac{\mathcal{V}_i(j)}{\widehat{\sigma}^2_{\mathbf{X}_i}}, \qquad j=1,\ldots,J.\end{equation*}

where $\rho_{i,j}\in[0,1]$ and the cumulative sum typically increases as J increases.

This variance contribution $\rho_{i,j}$ will be used for a practical choice of the decomposition depth J in Section 5. In particular, since the incremental contribution $\rho_{i,j}$ typically decreases as j increases, we can identify a cutoff level beyond which additional scales contribute negligibly to the overall variance and are therefore treated as practically insignificant for the representation.

This choice is analogous to the use of explained variance in principal component analysis: each MODWT level contributes a portion of the signal variance, and the cumulative contribution across retained levels provides a practical criterion for selecting the representation depth. The difference is that MODWT levels are ordered by time scale rather than by eigenvalue magnitude, so retaining levels up to J corresponds to retaining localized patterns up to a chosen temporal scale.

2.3. Across-level aggregation

As Weidner et al. (Reference Weidner, Transchel and Weidner2016) emphasize, risk-relevant patterns in telematics signals do not occur on a single time scale within a trip. We therefore aggregate MODWT wavelet coefficients $\widetilde{\mathbf{W}}_{i,j}$ across levels into a single time series. This across-level aggregation summarizes multi-scale patterns and enables a per-exposure behavioral risk index normalized for the trip length via an exposure offset, as will be described in Section 4.

For a set of levels $\mathcal{J}\subset\{1,\ldots,J\}$ , we define the aggregated wavelet coefficient for a trip i at each time point t as

\begin{equation*} C_{i,t} \;\, :\!= \,\; \mathfrak{g} \Bigl(\{\widetilde{W}_{i,j,t}\}_{j\in\mathcal{J}}\Bigr), \qquad t=0,\ldots,T_i-1,\end{equation*}

where $\{\widetilde{W}_{i,j,t}\}_{j\in\mathcal{J}}$ is the collection of the MODWT coefficients across levels at time point t and $\mathfrak{g}: \mathbb{R}^{|\mathcal{J}|}\to\mathbb{R}$ is the aggregation function. Possible choices of $\mathfrak{g}({\cdot})$ include single-level selection, average or weighted pooling, and maximum aggregation across levels. In this paper, we use the maximum rule

\begin{align*}C_{i,t}=\widetilde{W}_{i,j^*(i,t),t}, \end{align*}

where $j^*(i,t)=\arg\max_{j\in\mathcal{J}}|\widetilde{W}_{i,j,t}|$ . This choice naturally preserves the direction of the coefficient, which is needed to distinguish left-tail and right-tail behavior, while retaining the largest localized response across the retained MODWT levels. A single-level representation may miss behaviors that appear at the other time scales, whereas average or weighted pooling can dilute large localized coefficients by combining them with milder coefficients from other levels. The maximum choice therefore provides a better representation that preserves both the magnitude and direction of the most pronounced localized feature at each time point.

This aggregation produces a single time series of the aggregated wavelet coefficients for trip i

(4) \begin{equation} \mathbf{C}_i=(C_{i,0},\ldots,C_{i,T_i-1}),\end{equation}

which summarizes $\{\widetilde{W}_{i,j,t}\}_{j \in \mathcal{J}}$ across J levels at each time point. Regardless of the choice of $\mathfrak{g}({\cdot})$ , $\mathbf{C}_i$ is designed to aggregate driving behaviors that may appear at different temporal scales into a time series. In Sections 3 and 4, $\mathbf{C}_i$ will be an ingredient to build a portfolio distribution that is an input for our proposed framework.

3.1. Portfolio sample construction: hierarchical sampling and thinning

We define a portfolio distribution of aggregated wavelet coefficients (4) as input to the Gaussian–Uniform severity model. The estimation is affected by how coefficients are pooled across drivers and trips. One consideration here is representativeness: if all coefficients are naively pooled, drivers with more trips, longer trips, or over-represented driving styles can contribute disproportionately to the portfolio baseline. For example, if an aggressive driver has a large number of trips, aggressive behavior is over-represented in the portfolio, making the actually extreme patterns appear typical relative to the fitted baseline. The other consideration is temporal dependence because the aggregated coefficients within a trip form a time series. A naive pooling for all coefficients can violate the i.i.d. sampling assumption underlying likelihood fitting. We therefore construct the portfolio distribution using a sampling scheme that addresses both representativeness across drivers and trips and temporal dependence within trips.

Specifically, we propose a hierarchical sampling stage together with thinning within the trip. The hierarchical sampling stage is designed to improve representativeness: a driver is selected uniformly, and then a trip from that driver is selected uniformly. This gives drivers and trips comparable opportunity to contribute to the portfolio baseline, rather than weighting the baseline only by the number or length of available trips. In larger naturalistic portfolios, this helps prevent the baseline from being dominated by over-represented driver groups, routes, or driving styles. In the controlled empirical application in Section 5, we retain all available trips because the data are collected under a balanced experimental design with prespecified routes and instructed driving states.

The thinning stage then addresses the temporal dependence by retaining coefficients within each selected or available trip at time points that are spaced far enough apart so that the empirical autocorrelation is reduced below a chosen threshold. However, if the thinning lag is chosen separately for each trip, the exposure may reflect the signal complexity of that trip, rather than the observed amount of measurements only. To avoid this issue, we propose a common global thinning lag rather than using a different thinning lag for each trip. For each trip i, we first compute the sample autocorrelation function (ACF) of the aggregated coefficient series $\mathbf{C}_i$ and find the smallest lag $L_i$ such that the absolute autocorrelations for $K_{\mathrm{acf}}$ consecutive lags remain below a prespecified threshold. We then define the common thinning lag as $L_{\mathrm{global}} \, :\!= \, \max_i L_i.$ This conservative thinning lag ensures that all trips satisfy the ACF threshold criterion, because $L_{\mathrm{global}}$ is at least as large as every trip-specific thinning lag. This avoids making each trip’s exposure depend on its own signal complexity. Given $L_{\mathrm{global}}$ , each trip is thinned by retaining coefficients separated by this common lag. Specifically, for trip i, we randomly choose a starting offset $s_i\in\{0,1,\ldots,L_{\mathrm{global}}-1\}$ and define the retained time point set as

\begin{align*}\mathcal{T}_i=\left\{s_i+qL_{\mathrm{global}}\,:\,q=0,1,\ldots,\left\lfloor\frac{T_i-1-s_i}{L_{\mathrm{global}}}\right\rfloor\right\}. \end{align*}

Pooling the retained coefficients yields the portfolio sample

\begin{align*}\mathbf{C}\, :\!= \,\{C_r\}_{r=1}^n=\bigcup_i \{C_{i,t}\,:\,t\in\mathcal{T}_i\},\qquad n=\sum_i |\mathcal{T}_i|, \end{align*}

where r is a re-indexing of the retained (i, t) pairs. We treat the thinned portfolio sample $\mathbf{C}$ as an approximately i.i.d. input when fitting the Gaussian–Uniform severity model in the following sections. This approximate i.i.d. assumption is used as a working assumption for estimating the severity model. Although retained coefficients are still not exactly independent observations even after thinning, the thinning stage with a common lag provides a practical approximation to an i.i.d. portfolio sample while keeping exposure comparable across trips. This thinned portfolio sample $\mathbf{C}$ will be modeled by the Gaussian–Uniform mixture model in the following section.

3.2. Gaussian–Uniform severity model, parameter constraints, and likelihood function

3.2.2. Parameter constraints

Several constraints are imposed so that the fitted mixture is both numerically well behaved and interpretable as a severity model. First, we impose a tail-separation constraint to preserve the intended bulk–tail roles of the Gaussian–Uniform mixture (Coretto and Hennig Reference Coretto and Hennig2011; Coretto Reference Coretto2022). The Gaussian components represent the central, typical part of the portfolio distribution, while the Uniform components are restricted to tail regions. After ordering the Gaussian means as $\mu_1\le\cdots\le\mu_G$ , we require

\begin{align*}u^-_1 \le \mu_1-\delta\sigma_1,\qquad\mu_G+\delta\sigma_G \le u^+_1,\qquad \delta\gt0. \end{align*}

This keeps the closest Uniform tail endpoints sufficiently separated from the Gaussian core, preventing ordinary observations near the core from being absorbed by Uniform components.

Second, the Uniform supports are constrained to be ordered and non-overlapping within each tail. The left-tail and right-tail supports are consecutive intervals,

\begin{align*}\theta^-_{m'}=(u^-_{m'+1},u^-_{m'}),\qquad m'=1,\ldots,M^-, \qquad \text{and} \qquad\theta^+_{m}=(u^+_{m},u^+_{m+1}),\qquad m=1,\ldots,M^+. \end{align*}

This ordering allows the fitted intervals to be interpreted as shallow-to-deep severity layers.

Third, the fitted tail-layer probabilities are constrained to be compatible with the severity ordering. With layers ordered from the Gaussian core toward the deeper tails, we impose

\begin{align*}\pi^-_1\ge \pi^-_2\ge\cdots\ge \pi^-_{M^-},\qquad\pi^+_1\ge \pi^+_2\ge\cdots\ge \pi^+_{M^+}. \end{align*}

This ensures that deeper tail layers correspond to rarer portfolio behavior than shallower layers, which is important because the severity weights used later are constructed from the fitted layer probabilities.

Finally, component scales are bounded away from zero to avoid likelihood degeneracy in constrained mixture estimation (Tanaka and Takemura Reference Tanaka and Takemura2006). Denote the scales as

\begin{align*}v_j = \begin{cases}(u^-_{m'} - u^-_{m'+1})/\sqrt{12} & \text{if } j \text{ indexes } m' = 1, \dots, M^-, \\[6pt]\sigma_g & \text{if } j \text{ indexes } g = 1,\dots,G, \\[6pt](u^+_{m+1} - u^+_{m})/\sqrt{12} & \text{if } j \text{ indexes } m = 1, \dots, M^+,\end{cases} \end{align*}

For a fixed $d\in(0,1)$ , we require $\exp(\!-n^d)\le v_j$ for all j, where n is the size of the portfolio sample. This constraint prevents Gaussian standard deviations and Uniform interval widths from collapsing to zero. More detailed information on these constraints is given in the Online Supplementary Material.

3.2.3. Likelihood function

The resulting finite mixture density is written as

\begin{equation*} f(\!\bullet\,;\,\eta) \;\, :\!= \,\; \sum_{m'=1}^{M^-}\pi^-_{m'}\,u\!\left(\bullet\,;\,\theta^-_{m'}\right) \;+\; \sum_{g=1}^{G}\pi_g\,\phi\!\left(\bullet\,;\,\theta_g\right) \;+\; \sum_{m=1}^{M^+}\pi^+_{m}\,u\!\left(\bullet\,;\,\theta^+_{m}\right). \end{equation*}

where $\eta=(\boldsymbol{\pi},\boldsymbol{\Theta})$ is the vector of model parameters with

\begin{align*}\boldsymbol{\pi}=\left(\{\pi^-_{m'}\}_{m'=1}^{M^-},\{\pi_g\}_{g=1}^{G},\{\pi^+_{m}\}_{m=1}^{M^+}\right) \quad \text{and} \quad\boldsymbol{\Theta}=\left(\{\theta^-_{m'}\}_{m'=1}^{M^-},\{\theta_g\}_{g=1}^{G},\{\theta^+_{m}\}_{m=1}^{M^+}\right). \end{align*}

Given the portfolio sample $\mathbf{C}$ , the observed-data log-likelihood is

\begin{align*}l_n(\eta)=\sum_{r=1}^n \log f(C_r\,;\,\eta). \end{align*}

This likelihood serves as an objective function used to estimate the Gaussian core, the ordered Uniform tail layers, and the corresponding membership probabilities. For notational convenience, we combine and re-index those parameters, where $M=M^-+M^+$ , by

\begin{align*}\pi_m =\begin{cases}\pi^-_{m}, & m=1,\ldots,M^-,\\[2pt]\pi^+_{m-M^-}, & m=M^-+1,\ldots,M,\end{cases}\qquad\Theta_m =\begin{cases}\theta^-_{m}=(u^-_{m+1},u^-_{m}), & m=1,\ldots,M^-,\\[2pt]\theta^+_{m-M^-}=(u^+_{m-M^-},u^+_{m-M^-+1}), & m=M^-+1,\ldots,M.\end{cases} \end{align*}

The Gaussian–Uniform mixture converts the continuous portfolio sample $\mathbf{C}$ into a common set of portfolio-anchored severity layers. The Gaussian components define the central region of typical behavior, while the fitted Uniform supports $\widehat{\Theta}_1,\ldots,\widehat{\Theta}_M$ partition the left and right tails into abnormal behavioral regions. For a trip i, a retained coefficient contributes to layer m when

\begin{align*}C_{i,t}\in\widehat{\Theta}_m,\qquad m=1,\ldots,M. \end{align*}

Thus, the fitted supports determine where tail-layer counts are observed. Here, severity is interpreted as portfolio rarity of a telematics pattern, not as insurance claim size. The fitted Uniform mixture probabilities $\widehat{\pi}_m$ quantify this rarity: smaller $\widehat{\pi}_m$ corresponds to rarer tail behavior and therefore to a more severe behavioral deviation. These probabilities are later converted into severity weights so that the risk index combines how often a trip enters each fitted tail layer with how rare that layer is in the portfolio.

This construction connects the severity model to the behavioral risk index developed in Section 4. For each trip i, the fitted severity layers allow us to summarize the retained coefficient series $\{C_{i,t}\,:\,t\in\mathcal{T}_i\}$ by counting how many retained coefficients fall into each fitted severity layer. These layer-wise counts describe how often trip i enters each tail region, while the fitted probabilities $\widehat{\pi}_m$ describe how rare those regions are in the portfolio. The risk index then combines these two pieces of information: the frequency with which a trip visits the fitted tail layers and the severity implied by the rarity of those layers.

3.3. Parameter estimation by MU-MEMR

Uniform endpoints make the mixture likelihood discontinuous, so Coretto (Reference Coretto2022) proposed MEMR to make fitting feasible for a model with a single Uniform component. In this paper, we propose a MU-MEMR algorithm to extend the MEMR to multiple Uniform components. The purpose of MU-MEMR is to provide a feasible likelihood-based fitting procedure for a Gaussian–Uniform mixture with multiple ordered Uniform tail components. It does so by replacing the continuous endpoint-estimation problem with a finite, constrained search over candidate tail partitions, while using EM updates to the parameters and membership probabilities of components within each candidate configuration. This extension is central to our framework because it enables the fitted portfolio distribution to contain multiple thresholds for the severity of telematics patterns.

The MU-MEMR estimates the parameter vector $\eta=(\boldsymbol{\pi},\boldsymbol{\Theta})$ under a finite constrained likelihood search. The Gaussian core is initialized by trimmed k-means (Cuesta-Albertos et al. Reference Cuesta-Albertos, Gordaliza and Matrán1997) so that extreme tail observations do not dominate the central fit. Then, for a fixed $(M^-,G,M^+)$ , we use the trimmed observations to construct finite candidate grids for the left- and right-tail endpoints,

\begin{align*}\mathcal U^-=\{b^-_{(1)},\ldots,b^-_{(q)}\},\qquad\mathcal U^+=\{b^+_{(1)},\ldots,b^+_{(p)}\}, \end{align*}

where $\mathcal U^-$ and $\mathcal U^+$ are finite grids of candidate cut points for the left- and right-tail Uniform endpoints. They are constructed from the trimmed tail observations by dividing the empirical left- and right-tail ranges into q and p ordered grid locations, respectively, and then snapping each grid location to a nearby observed coefficient. Thus, q and p control the resolution of the tail-endpoint search: larger values provide a finer approximation to the continuous endpoint space, at the cost of more candidate partitions. For a given $(M^-,G,M^+)$ , one candidate endpoint configuration is obtained by selecting

\begin{align*}1\le \Delta^-_1\lt\cdots\lt\Delta^-_{M^-}\le q,\qquad1\le \Delta^+_1\lt\cdots\lt\Delta^+_{M^+}\le p, \end{align*}

which corresponds to choosing one of the $\binom{q}{M^-}$ left-tail endpoint combinations and one of the $\binom{p}{M^+}$ right-tail endpoint combinations. The selected grid points then determine the Uniform endpoints,

\begin{align*}u^-_{M^-+1-m'} = b^-_{(\Delta^-_{m'})},\quad m'=1,\ldots,M^-,\qquad u^+_{m}=b^+_{(\Delta^+_{m})},\quad m=1,\ldots,M^+, \end{align*}

where the outer fixed tail endpoints are fixed at the tail extremes $u^-_{M^-+1}=c^-_{\min}$ and $u^+_{M^++1}=c^+_{\max}.$

For each candidate endpoint configuration, the Uniform supports are held fixed, and the remaining Gaussian parameters and membership probabilities are estimated. More specifically, for each retained coefficient $C_r$ , the E-step computes the posterior membership probability of component k,

\begin{align*}\tau_k(C_r\,;\,\eta^{(h)})=\frac{\pi_k^{(h)} f_k(C_r\,;\,\theta_k^{(h)})}{\sum_{\ell}\pi_\ell^{(h)} f_\ell(C_r\,;\,\theta_\ell^{(h)})}, \end{align*}

where h indexes the iteration and k ranges over the $(M^-+G+M^+)$ components in the mixture model. In the M-step, these posterior memberships are used to update the component probabilities and the Gaussian parameters, while the candidate Uniform endpoints remain fixed for that run.

Thus, MU-MEMR is a constrained likelihood-based implementation of the Gaussian–Uniform mixture model, building on the MEMR idea of Coretto (Reference Coretto2022) and extending it to multiple ordered Uniform tail layers. Although this discrete candidate search procedure does not guarantee a global optimum over the continuous space of all possible Uniform endpoints, we can approximate the continuous space when the resolution q and p increase. The full algorithm is provided in the Online Supplementary Material.

4.1. Trip-level behavioral risk index

In Section 2, we introduced the series of aggregated wavelet coefficients for trip i as $\mathbf{C}_i=\{C_{i,t}\}^{T_i-1}_{t=0}$ , where $t=0,\ldots,T_i-1$ indexes time points and $\mathcal{T}_i$ is the set of retained coefficients obtained by the proposed sampling scheme. In this section, we set $E_i\, :\!= \,|\mathcal{T}_i|$ to interpret as an exposure measure. For each individual trip i, we want to examine the frequency profile for each severity layer, so we define the trip-level multi-layer tail counts (MLTCs) by

\begin{equation*} N_{im} \;\, :\!= \,\; \sum_{t=0}^{E_i-1}\mathbf{1}\{\,C_{i,t}\in\Theta_m\,\}, \qquad m=1,\ldots,M,\end{equation*}

where $N_{im}$ counts the number of coefficients that fall in layer $\Theta_m$ . The vector $\mathbf{N}_i=(N_{i1},\ldots,N_{iM})$ summarizes the tail profile across severity layers for the trip i. Because the layers ${\Theta_m}$ are disjoint, $\mathbf{N}_i$ captures not only how often tail patterns occur but also how the tail counts are distributed across depths.

A Poisson model with exposure offset $E_i$ provides a simple way to model MLTCs and to compare trips of different lengths through per-exposure intensities $\lambda_{im}$ . Conditional on the per-exposure intensities $\lambda_{im}\gt0$ , we therefore propose to model the MLTC by

\begin{align*}N_{im}\mid \lambda_{im},E_i \;\sim\; \mathrm{Poisson}(E_i\,\lambda_{im}),\qquad m=1,\ldots,M. \end{align*}

To update the risk profile using a closed form, we propose to use independent Gamma priors on intensities $\lambda_{im}$ for each severity layer m as

\begin{align*}\lambda_{im}\sim \mathrm{Gamma}(\alpha_{0m},\beta_{0m}),\qquad m=1,\ldots,M, \end{align*}

which yields conjugate posteriors and supports stable trip- and sequential driver-level risk evaluation. By Gamma–Poisson conjugacy, the intensity is updated with the information $(E_i, \mathbf{N}_i)$ of trip i:

\begin{align*}\lambda_{im}\mid E_i,N_{im}\;\sim\;\mathrm{Gamma}(\alpha_{0m}+N_{im},\;\beta_{0m}+E_i),\qquad m=1,\ldots,M. \end{align*}

The layer-wise hyperparameters $(\alpha_{0m},\beta_{0m})$ are initialized by a winsorized-moment empirical Bayes step (Phipson et al. Reference Phipson, Lee, Majewski, Alexander and Smyth2016) to reduce the influence of extreme rates on the prior fit. For each layer m, we compute per-exposure rates $N_{im}/E_i$ across trips and replace the values outside the $[\Omega, 1-\Omega]$ range with the appropriate quantile cutoff $q_{\Omega,m}$ or $q_{1-\Omega,m}$ . We then initialize the hyperparameters $(\alpha_{0m}, \beta_{0m})$ by moment matching using the sample mean and variance of the winsorized rates.

The Poisson specification is used as a convenient working frequency model. Because MLTCs summarize visits of MODWT coefficients to severity layers rather than fully independent observations, residual dependence after thinning may induce overdispersion. Section 5.5 therefore compares the Poisson–Gamma model with a Negative Binomial–Beta alternative. Since the main classification and risk index conclusions are not materially changed, we use the simpler Poisson–Gamma model as an interpretable baseline frequency specification.

The MLTCs are penalized based on their portfolio-level rarity. For the penalty, we define severity weights by

\begin{equation*}w_m\;\, :\!= \,\;\frac{(1/\pi_m)^{\gamma}}{\sum_{l=1}^M (1/\pi_l)^{\gamma}},\qquad m=1,\ldots,M,\end{equation*}

where the membership probabilities $\{\pi_m\}_{m=1}^M$ are estimated by the severity model and $\gamma \ge 0$ is selected over a fixed grid using leave-one-driver-out classification validation in the empirical validation in Section 5. In general, when a validation dataset is not available for choosing $\gamma$ , we suggest choosing a value which is able to balance between information contained in the shallower and deeper layers, for example, $\gamma = 1$ . The parameter $\gamma$ controls the dispersion of severity weights over the layers: smaller $\gamma$ spreads severity weights more uniformly across layers, while larger $\gamma$ puts more weight on deeper layers.

For the trip i, we now propose the trip-level behavioral risk index that summarizes severity-weighted layer-wise risk rates per exposure:

\begin{align*}S_i \;\, :\!= \,\; \sum_{m=1}^M w_m\,\lambda_{im}. \end{align*}

This index aggregates the per-exposure intensities $\{\lambda_{im}\}_{m=1}^M$ of an individual trip, while using the portfolio-level severity weights $\{w_m\}_{m=1}^M$ that commonly penalize the layers across trips. Since we assume conditional independence of $N_{im}$ given $(\lambda_{im}, E_i)$ and independent Gamma priors across layers m, we can estimate $S_i$ by a sum of layer-wise posterior means:

(5) \begin{equation}\widehat{S}_i\;=\;\mathbb{E}[S_i\mid E_i,\mathbf{N}_i]\;=\;\sum_{m=1}^M w_m\, \frac{\alpha_{0m} + N_{im}}{\beta_{0m} + E_i},\end{equation}

where $\dfrac{\alpha_{0m}+N_{im}}{\beta_{0m}+E_i} = \mathbb{E}[\lambda_{im}\mid E_i,\mathbf{N}_i]$ . The estimated index $\widehat{S}_i$ admits the layer-wise decomposition, so each term $w_m\,\mathbb{E}[\lambda_{im}\mid E_i,\mathbf{N}_i]$ quantifies the contribution of severity layer m to the overall risk. Consequently, $\widehat{S}_i$ jointly combines the frequency information in the MLTC profile $\mathbf{N}_i$ with the portfolio-level severity weights $\{w_m\}^M_{m=1}$ on a common per-exposure scale through $E_i$ , yielding a frequency–severity joint trip-level behavioral risk index.

4.2. Driver-level behavioral risk index

Once the trip-level index is defined, we extend this to driver-level through sequential updating of $\lambda_{im}$ over a driver’s trips. Suppose trips of driver d are observed in time order $i=1,\ldots,k$ , and trip i provides exposure $E_i$ and MLTC $\mathbf{N}_{i} = (N_{i1}, \dots, N_{iM})$ . For each layer m, let $\lambda_{dm}$ denote the driver’s per-exposure intensity of tail layer m. After k trips, conjugacy gives

\begin{align*}\lambda_{dm} \mid \mathcal{H}_d^{(k)} \sim \mathrm{Gamma}(\alpha_{dm}^{(k)},\beta_{dm}^{(k)}), \qquad \mathcal{H}_d^{(k)}=\{(E_i,\mathbf{N}_i) \, : \, i=1,\ldots,k\}, \end{align*}

with sequential updates

\begin{align*}\alpha_{dm}^{(k)}=\alpha_{dm}^{(k-1)}+N_{km}, \qquad \beta_{dm}^{(k)}=\beta_{dm}^{(k-1)}+E_k, \end{align*}

where $\alpha_{dm}^{(0)}=\alpha_{0m}$ and $\ \beta_{dm}^{(0)}=\beta_{0m}$ initialized in the same way for the trip-level index.

To estimate a driver’s long-term risk profile, we summarize each layer intensity $\lambda_{dm}$ after k trips by its posterior mean. Under the Poisson–Gamma model, this model choice yields a closed-form update that can be implemented online as each new trip arrives. Moreover, the posterior mean admits a credibility representation because it blends the accumulated history $\widehat{\lambda}^{(k-1)}_{dm}$ with the new trip’s empirical rate $N_{km}/E_k$ through exposure-based weights. The posterior mean intensity after k trips is

\begin{equation*}\widehat{\lambda}_{dm}^{(k)} \;\, :\!= \,\; \mathbb{E}\!\left[\lambda_{dm}\mid \mathcal{H}_d^{(k)}\right] \;=\; \frac{\alpha_{dm}^{(k)}}{\beta_{dm}^{(k)}} \;=\; \underbrace{\frac{\beta_{dm}^{(k-1)}}{\beta_{dm}^{(k-1)}+E_k}}_{\text{weight on past}}\,\widehat{\lambda}_{dm}^{(k-1)} \;+\; \underbrace{\frac{E_k}{\beta_{dm}^{(k-1)}+E_k}}_{\text{weight on new}} \,\frac{N_{km}}{E_k},\end{equation*}

where the weights depend on accumulated precision and exposure: larger $E_k$ increases the influence of the new trip, whereas small $E_k$ yields a more conservative update driven by the past through $\beta_{dm}^{(k-1)}$ .

Given the updated layer intensities, we estimate the driver-level behavioral risk index after k trips:

\begin{equation*}\widehat{S}_d^{\,(k)} \;\, :\!= \,\; \sum_{m=1}^M w_m\,\widehat{\lambda}_{dm}^{(k)}= \sum_{m=1}^M w_m\,\frac{\alpha_{dm}^{(k)}} {\beta_{dm}^{(k)}}.\end{equation*}

This index provides a time-evolving summary of the driver’s overall risk profile, aggregating evidence accumulated across trips. The trip-level index is recovered as the special case $k=1$ , then $\alpha^{(1)}_{dm} = \alpha_{0m} + N_{1m}$ and $\beta^{(1)}_{dm} = \beta_{0m} + E_{1}$ , so $\hat{S}^{(1)}_d$ matches (5).

In Section 5, we will use these outputs: MLTC profiles $\mathbf{N}_i$ , trip-level indices $\widehat{S}_i$ , and sequentially updated driver-level indices $\widehat{S}^{(k)}_d$ , to examine how different behavioral states exhibit distinct profile in terms of both frequency and severity profiles and to evaluate how well the proposed behavioral risk index distinguishes risky from normal trips in classification.

5.1. Portfolio construction and model setup

We select the MODWT decomposition depth J using the variance contribution ${\rho_{i,j}}$ defined in Section 2 so that we retain as much risk-relevant information as possible while keeping the representation computationally efficient. We begin from a conservative reference trip by Driver 4, the aggressive trip on motorway, because it has the longest trip length and therefore admits the largest possible decomposition depth under the conventional rule for a full depth, $J_{\text{full}}=\lfloor \log_2(T_i)\rfloor$ . The cumulative variance of this trip is $99.2\%$ at $J_{\text{full}}=13$ , while truncation at $J=6$ and $J=7$ retains $87.8\%\ (\!=\sum^6_{j=1}\rho_{i,j})$ and $94.6\%\ (\!=\sum^7_{j=1}\rho_{i,j})$ of the sample variance $\mathbf{X}_i$ , respectively. Figure 3 shows a sharp drop in the incremental variance contribution when moving from level 6 to level 7. We therefore select $J=6$ as a conservative and efficient choice.

Following the sampling scheme proposed in Section 3, we construct the portfolio sample $\mathbf{C}$ from the aggregated coefficient series $\mathbf{C}_i=\{C_{it}\}_{t=0}^{T_i-1}$ . For each trip, $\mathbf{C}_i$ is obtained by applying the maximum aggregation rule over the retained MODWT levels $J=\{1,\ldots,6\}$ , as specified in Section 2.

At the hierarchical sampling stage, we select all 40 available trips, as discussed in Section 3.1. At the thinning stage, we set the ACF threshold as $0.1$ and require the threshold condition to hold for $K_{\mathrm{acf}}=3$ consecutive lags. For each trip i, we compute the candidate lag $L_i$ from the sample ACF of $\mathbf{C}_i$ and set $L_{\mathrm{global}}=\max_i L_i$ for thinning. The dataset gives $L_{\mathrm{global}}=15$ , so each trip is thinned by retaining coefficients every 15 time points. Using this common lag yields 20,759 retained coefficients for the portfolio sample $\mathbf{C}$ from the original 311,381 observations. We then define exposure as $E_i\, :\!= \,|\mathcal{T}_i|$ and compute the MLTC profiles and risk indices on this sample.

We now fit the severity model to the portfolio sample $\mathbf{C}$ constructed using the maximum aggregation rule, $J=6$ , and $L_{\text{global}} = 15.$ For the model tuning parameters, we set $(d,\alpha,\delta,q,p)=(0.5,\,0.05,\,1.96,\,12,\,10)$ and fit over $G \in \{1,2\}$ , $M^- \in \{1, 2, \dots,8\}$ , and $M^+ \in \{1, 2, \dots, 6\}$ . We set $d=0.5$ as recommended by Tanaka and Takemura (Reference Tanaka and Takemura2006) and $\alpha=0.05$ to use the outer 5% of the portfolio to form the base grid for candidate search for Uniform parameters. The grid (q, p) controls the tail partition resolution for the base grid, and $(M^-, M^+)$ controls the number of severity layers on each tail. We use finer resolution on the left tail by setting $q=12$ and allowing more layers $M^- \in \{1, 2, \dots,8\}$ because the portfolio sample has a longer left tail, while we set $p=10$ with $M^+ \in \{1, 2, \dots,6\}$ on the right tail.

5.2. Fitted severity model and portfolio-level severity weights

Under the model setup described above, we fit the proposed MU-MEMR severity model over $(M^-,G,M^+)$ and consider AIC/BIC/log-likelihood as the model selection criteria. In what follows, we report results under $(M^-,G,M^+)=(4,2,4)$ selected by BIC. Compared with AIC and log-likelihood, BIC discourages overly complex mixtures and reduces the chance of creating extremely sparse deep layers. This helps prevent the risk index from depending too heavily on the deepest severity layers only.

Using the estimated membership probabilities of the fitted severity layers, we define the severity weights $\mathbf{w}(\gamma)=(w_1(\gamma),\ldots,w_M(\gamma))$ . The tuning parameter $\gamma$ controls the dispersion of severity weights across layers and therefore directly affects risk assessment: if $\gamma$ is too small, the weights are nearly uniform and the severity component has little effect; if $\gamma$ is too large, the risk index may be driven by only a few observations in the most extreme layers, reducing the benefit of the multi-layer construction. To obtain the best severity scale for a robust risk assessment, we select $\gamma$ by maximizing the discrimination accuracy between normal and risky trips in the binary classification validation. The value is selected as $\gamma^\star=0.9$ , and thus we apply these severity weights $\mathbf{w}=\mathbf{w}(\gamma^\star)$ consistently to all trips as a common portfolio-level severity scale. This portfolio-anchored model makes the proposed index comparable across trips. More detailed information is described in the Online Supplementary Material.

The resulting fitted severity model, shown in Figure 4 and Table 3, illustrates how the Gaussian–Uniform mixture separates typical driving dynamics from rare tail deviations. The bulk–tail separation in Figure 4 shows that the Gaussian components stably capture the central portfolio baseline, while the Uniform components absorb abnormal tail behavior. Although the empirical distribution has a longer left tail, the two Gaussian components remain nearly symmetric around zero with similar variances, as reported in Table 3. This indicates that tail observations do not inflate the Gaussian scale because the Uniform components absorb the contaminated tail region. This separation is important for risk assessment: without Uniform absorption, the fitted baseline could treat the tail region as ordinary behavior and become undesirably generous for tail patterns. The fitted severity model therefore provides a stable portfolio-level reference for typical behavior and a set of ordered tail layers for quantifying signal severity.

Table 3.

Estimated parameters of the Gaussian–Uniform severity model based on BIC: Gaussian bulk components with their means $\mu_g$ , scales $\sigma_g$ , and membership probabilities $\pi_g$ , and Uniform tail layers with their intervals, layer probabilities $\pi_m$ , and severity weights $w_m$ . The severity weights are computed using the fixed value $\gamma^\star=0.9$ .

Table 3. Long description

A table titled Estimated parameters of the GaussianUniform severity model based on BIC. The table is divided into two main sections: Gaussian component and Uniform component. The Gaussian component section has two columns labeled g = 1 and g = 2, each with three rows for mean, scale, and membership probability. The Uniform component section has intervals, layer probabilities, and severity weights for left-tail and right-tail layers. The Gaussian component section shows values for means, scales, and membership probabilities for g = 1 and g = 2. The Uniform component section lists intervals, layer probabilities, and severity weights for left-tail layers with m’ values from 4 to 1 and right-tail layers with m values from 1 to 4.

Figure 4.

Fitted Gaussian–Uniform mixture under BIC for the portfolio sample $\boldsymbol{C}$ of aggregated wavelet coefficients: the Gaussian bulk captures typical driving dynamics, while multi-layer Uniform tails define ordered severity layers; rug marks along the axis show the assignment of individual coefficients.

5.3. Trip- and driver-level behavioral risk indices

We estimate the proposed trip- and driver-level behavioral risk indices by combining the Poisson–Gamma frequency model in Section 4 with the portfolio-level severity layers and severity weights from the fitted Gaussian–Uniform mixture model in Section 3.

Table 4.

Numerical result on secondary road. Behavioral Risk Index: Trip (trip-level) and Updated (driver-level) with exposure given by the number of retained coefficients $E_i \, :\!= \,|\mathcal{T}_i|$ and MLTC counts in the left tail (Lv4- $\sim$ Lv1-) and right tail (Lv1+ $\sim$ Lv4+) under the layers $(M^-{=}4,\,M^+{=}4)$ . The severity weights are computed using the fixed value $\gamma^\star=0.9$ . Zeros in the MLTC columns are left blank.

Table 4. Long description

The table presents numerical results on secondary roads, focusing on behavioral risk indices and MLTC counts for different drivers and trip conditions. It includes trip-level and driver-level behavioral risk indices, exposure values, and MLTC counts in the left and right tails under various layers. The table has 6 rows and 11 columns. Column headers are Exposure, Trip, Updated, Lv4, Lv3, Lv2, Lv1, Lv1+, Lv2+, Lv3+, and Lv4+. Row labels include Driver Label and specific conditions like normal1, aggressive, normal2, and drowsy. Row 1: D1, normal1, 411, 3.47, 3.47, 1, , , , , , . Row 1: D1, aggressive, 354, 15.64, 9.49, 4, 5, 9, 2, 2, 1. Row 1: D1, normal2, 435, 5.42, 7.75, 1, 1, 3, 3, , . Row 1: D1, drowsy, 335, 14.89, 9.75, 3, 3, 7, 4, 3, 2. Row 2: D2, normal1, 434, 4.32, 4.32, 2, 1, 1, , , . Row 2: D2, aggressive, 411, 34.84, 20.83, 4, 17, 27, 30, 9, 4. Row 2: D2, normal2, 458, 4.31, 15.01, 4, 3, 1, , , . Row 2: D2, drowsy, 434, 20.52, 17.52, 3, 5, 9, 8, 4, 4, 1. Row 3: D3, normal1, 466, 3.09, 3.09, 1, 1, , , , . Row 3: D3, aggressive, 452, 10.44, 6.21, 3, 14, 8, 4, 1, , . Row 3: D3, normal2, 446, 4.28, 4.99, 7, 2, , , , . Row 3: D3, drowsy, 442, 7.98, 5.51, 1, 3, 1, 3, , . Row 4: D4, normal1, 486, 3.14, 3.14, 2, 1, , , , . Row 4: D4, aggressive, 441, 11.09, 6.46, 1, 2, 13, 14, 4, , . Row 4: D4, normal2, 460, 3.69, 4.92, 5, 1, , , , . Row 4: D4, drowsy, 458, 6.41, 4.95, 1, 1, 1, 3, 1, 1. Row 5: D5, normal1, 467, 3.97, 3.97, 5, 1, 1, , , . Row 5: D5, aggressive, 282, 41.59, 19.93, 1, 3, 4, 5, 3, 3, 4, 1. Row 5: D5, normal2, 462, 6.27, 14.72, 1, 1, 6, 6, , . Row 5: D5, drowsy, 467, 6.43, 12.27, 1, 1, 4, 2, 1, , . Row 6: D6, normal, 519, 6.75, 6.75, 10, 8, 1, 1, , , . Row 6: D6, drowsy, 480, 11.20, 8.99, 3, 11, 4, 6, 1, , .

The frequency and severity models are fitted using the full set of 40 trips pooled across both secondary and motorway roads, so the estimated portfolio baseline and severity scale reflect the entire dataset. For exposition and interpretability, we focus on the secondary road when presenting MLTC profiles and risk indices. This focus can isolate behavioral differences and enable characterization of the resulting frequency–severity profiles across the three driving behavioral states by examining the trips on the same route and speed regime. Moreover, D1–D5 have four trips on secondary road for each (vs. three motorway trips), which provides a richer sequence within drivers and allows the sequential updating mechanism to be illustrated more clearly: normal trips tend to moderate the driver-level index, while risky trips tend to increase it.

From the estimated MLTC profiles and risk indices, we demonstrate both (i) how tail depth changes ranking of trip-level indices and (ii) how sequential updating aggregates evidence into a driver-level risk profile over time. Table 4 illustrates how the multi-layer tail structure shapes risk evaluation by accounting for which severity layers are reached, not only how many tail counts occur in total. Across drivers, normal trips are characterized by small MLTC totals concentrated in shallow layers (typically Lv1 $\pm$ or Lv2 $\pm$ ), whereas the risky trips (i.e., aggressive and drowsy) tend to exhibit either higher MLTC totals and/or a broader spread of tail mass into deeper layers. For risky trips, aggressive trips typically distribute MLTC mass more heavily across multiple layers, whereas drowsy trips more often display a lighter overall mass but still reach deep layers.

This layer-wise structure characterizes interpretable frequency–severity profiles across different behavioral driving states. Here, a trip’s frequency profile refers to the overall MLTC mass, while its severity profile describes how that tail mass is distributed across layers. The normal trips tend to be low-frequency/low-severity because most of their MLTC mass is concentrated in the shallow layers. The aggressive trips tend to be high-frequency/high-severity because they show many tail counts that spread across layers and commonly reach deeper tail layers (e.g., Lv3 $\pm$ or beyond). The drowsy trips often exhibit a distinct low-frequency/high-severity profile because their counts are often similar in shallow layers compared with normal, but the observed counts are more likely to appear in deeper layers. This is exactly where the proposed framework is most informative because drowsy driving can look “quiet” under the frequency-only perspective, while its high-severity profile reveals itself in the tails.

Once the frequency model estimates layer-wise risk rates $\lambda_{im}$ for a trip, we compute the trip-level behavioral risk index (Trip) by averaging these rates with the portfolio-level severity weights, that is, a severity-weighted average of risk rates. Table 4 shows that this construction yields a clear within-driver separation between normal and risky trips. For each driver, every risky trip has a larger trip-level index than the driver’s normal trips. That is, it consistently identifies the instructed risky trips as riskier than the corresponding normal trips within each driver. This is a meaningful property in the controlled setting because drivers and vehicles may have different baseline driving styles, creating heterogeneity across drivers. Even with this heterogeneity, the proposed index can recover the behavioral driving state within each driver by assigning higher trip-level indices to the risky trips than to the corresponding normal trips.

The MLTC profiles explain why this separation is obtained. For Drivers 1–5, the aggressive trip has the largest trip-level index within the driver, reflecting heavier and broader tail behavior. For instance, Driver 3’s aggressive trip has the largest index within Driver 3, $10.44$ , with tail mass spread across both left- and right-tail layers, including Lv2− $=3$ , Lv1− $=14$ , Lv1+ $=8$ , Lv2+ $=4$ , and Lv3+ $=1$ . At the same time, the MLTC profile also demonstrates that the proposed index is not merely a frequency-based measure. For Driver 3, normal2 has more total MLTCs than drowsy (9 vs. 8) and both have similar exposures, but drowsy receives a larger index ( $7.98$ vs. $4.28$ ) because its tail counts reach deeper severity layers, including Lv3− $=1$ and Lv3+ $=3$ , whereas normal2 is concentrated only in shallower layers (Lv1- $=7$ and Lv1+ $=2$ ). Similarly, this interpretation works even for the same label, Driver 2’s drowsy trip receives a much larger index than Driver 5’s drowsy trip ( $20.52$ vs. $6.43$ ), because Driver 2’s drowsy trip has substantially heavier and deeper tail-layer counts, including Lv3- $=3$ , Lv2- $=5$ , Lv1- $=9$ , Lv1+ $=8$ , Lv2+ $=4$ , Lv3+ $=4$ , and Lv4+ $=1$ . These examples show that the index provides an interpretable ordering not only across behavioral states but also within the same state by distinguishing relatively mild from more severe tail behavior.

The MLTC breakdown also reveals which side of the tail drives risk, linking back to the physical interpretation of longitudinal acceleration. For example, Driver 1’s aggressive has a stronger right-tail profile than left-tail profile: its left-tail counts are Lv2− $=4$ and Lv1− $=5$ (total 9), whereas its right-tail counts are Lv1+ $=9$ , Lv2+ $=2$ , Lv3+ $=2$ , and Lv4+ $=1$ (total 14). This indicates that abnormal acceleration patterns contribute more to this trip’s behavioral risk than deceleration patterns. Thus, the MLTC profile helps identify which type of maneuver is driving the trip-level risk.

As a long-term risk measure, it is essential to assess a driver’s overall risk profile. We compute a driver-level (Updated) behavioral risk index via sequential updating to quantify how a driver’s overall profile evolves as trips accumulate. We illustrate this sequential updating dynamic by assuming an illustrative order of trips for each driver. Specifically, we arrange the four secondary road trips as normal1 $\rightarrow$ aggressive $\rightarrow$ normal2 $\rightarrow$ drowsy because this alternating normal–risky pattern can demonstrate how the driver-level profile reacts as different risk evidence is incorporated over time.

Under this illustrative order, Table 4 shows that the updated index typically increases when a risky trip is incorporated and is moderated when a normal trip is incorporated. For example, Driver 1 changes from $3.47$ after normal1 to $9.49$ after aggressive, then to $7.75$ after normal2, and finally to $9.75$ after drowsy. Driver 2 similarly changes from $4.32$ to $20.83$ , then to $15.01$ , and finally to $17.52$ . These sequences show the credibility-style nature of the update: risky trips tend to increase the driver-level profile, while normal trips can moderate it by adding lower risk evidence.

To understand the updating mechanism at the layer level, consider Driver 5 around the transition normal2 $\rightarrow$ drowsy. Let $\widehat{\lambda}_{D5,m}^{(k)}=\mathbb{E}\!\left(\lambda_{D5,m}\mid \mathcal{H}_{D5}^{(k)}\right)$ be Driver 5’s posterior mean intensity for severity layer m after the first k trips. Here, $k=3$ corresponds to the history $\{\texttt{normal1},\texttt{aggressive},\texttt{normal2}\}$ , and $k=4$ additionally includes the drowsy trip. The updated index decreases from $14.72$ to $12.27$ because the fourth trip reduces several posterior layer intensities that carry high severity weights. In units of $10^{-4}$ , the deepest left-tail intensity decreases from $\widehat{\lambda}_{D5,\mathrm{Lv4-}}^{(3)}=8.25$ to $\widehat{\lambda}_{D5,\mathrm{Lv4-}}^{(4)}=5.96$ , and the deepest right-tail intensity decreases from $\widehat{\lambda}_{D5,\mathrm{Lv4+}}^{(3)}=7.24$ to $\widehat{\lambda}_{D5,\mathrm{Lv4+}}^{(4)}=5.65$ . The Lv3+ intensity also decreases from $31.92$ to $29.31$ . Although Lv2+ increases slightly from $34.31$ to $36.42$ , this layer has a smaller severity weight than Lv4− and Lv4+, so the increase does not offset the reductions in the deeper, higher–weight layers. Consequently, the severity-weighted sum $\widehat{S}_{D5}^{(k)}=\sum_m w_m\widehat{\lambda}_{D5,m}^{(k)}$ moves downward after the drowsy trip. As more trips are accumulated, the driver-level index provides a more stable summary of the driver’s long-term behavioral risk profile.

5.4. Binary classification validation: normal versus risky

In this subsection, we demonstrate how effective the proposed severity weighting is at separating instructed behavioral states through a binary trip-level classification task (normal vs. risky). We use all 40 trips from the controlled datasets and compare three MLTC-based inputs under the same Poisson–Gamma classifier: Model (A) Threshold-count frequency (single layer) uses the scalar total MLTC count $N_{i1}^{\mathrm{total}}=\sum_{m=1}^M N_{im},$ with $M=1$ and $w_1\equiv 1$ . This model removes the layer-wise structure and severity weighting and therefore serves as a simple threshold-count benchmark that aggregates the counts of harsh braking or acceleration. Model (B) Layer-wise frequency (unweighted) uses the full MLTC vector $\mathbf{N}_i=(N_{i1},\ldots,N_{iM})$ with uniform weights $w_m\equiv 1/M$ , preserving the multi-layer structure without severity weighting. Model (C) Severity-weighted frequency (proposed) uses the same MLTC vector but applies rarity-based severity weights $w_m(\gamma)$ from the fitted severity model. The step (A) $\to$ (B) isolates the gain from moving beyond a single threshold-count benchmark to a multi-layer threshold representation, and the step (B) $\to$ (C) isolates the gain from incorporating portfolio-level severity weights into risk evaluation; see the Online Supplementary Material for algorithmic details.

We evaluate the three models using leave-one-driver-out cross-validation (LODO CV), where the test fold contains all trips from one held-out driver. This validation scheme is used to assess whether the proposed framework generalizes to an unseen driver by preventing driver leakage across trips from the same driver. The severity model is re-estimated by MU-MEMR within each training fold. In each LODO fold, the held-out driver is excluded before constructing the portfolio sample and before fitting the Gaussian–Uniform mixture. The fitted severity layers and membership probabilities are then applied to the held-out driver’s trips only for evaluation. In this setting, the portfolio baseline, severity layers, membership probabilities, and severity weights are learned from the training drivers only. The out-of-fold predictions from the six held-out drivers are then pooled to construct Figure 5 and to compute the reported balanced accuracy (BA).

The classification result gives the best $\gamma$ as $\gamma^\star=0.9$ over $\{0.1,0.2,\ldots,2.0\}$ under the BA. Figure 5 summarizes the performance of the three models: Model A uses only the total number of MLTCs and achieves BA $76.7\%$ ; Model B preserves the layer-wise MLTC structure but does not use severity weights and achieves BA $76.0\%$ ; Model C incorporates portfolio-level severity weights and achieves BA $92.7\%$ in these controlled datasets. The comparison between Models A and B shows that retaining the layer structure alone is not sufficient, whereas the improvement from Model B to Model C shows that severity weighting is crucial for improving discrimination beyond total frequency counting.

Figure 5.

Binary classification (normal vs. risky) using MLTC-based inputs under leave-one-driver-out cross validation (LODO CV). The figure shows the out-of-fold distributions of $P(\texttt{risky})$ for Total frequency, Layer-wise frequency, and Severity-weighted frequency. For Model C, the severity weights use the best fixed value $\gamma^\star=0.9$ from the $\gamma$ grid under the Poisson–Gamma classifier. Dashed horizontal lines indicate the model-specific decision thresholds selected from the training folds.

Specifically, Figure 5 compares the out-of-fold distributions of $P(\texttt{risky})$ for all 40 trips and highlights both the between-class separation and the within-class concentration. Viewed from these two perspectives, Model A (BA $=76.7\%$ ) shows weak separation on both fronts. Its predicted probabilities are highly polarized, but the polarization is not well aligned with the true labels: risky trips appear both near 0 and near 1, while normal trips are mostly near 0 but still include observations above the decision threshold. As a result, the two classes overlap at both low and high probability levels, and the risky class exhibits poor within-class concentration. Model B (BA $=76.0\%$ ) improves the representation by retaining the layer-wise MLTC structure, and this yields a meaningful within-class concentration for normal trips. However, the risky trips remain widely dispersed from low to high $P(\texttt{risky})$ , so the between-class overlap is still substantial. This explains why Model B does not improve the overall discrimination relative to Model A despite using richer input information. Model C (BA $=92.7\%$ ) shows the clearest improvement in both between-class separation and within-class concentration. After severity weighting is introduced, the normal trips form a relatively compact band below the threshold, whereas the risky trips are concentrated above the threshold, with many lying in the upper-probability region and only limited overlap near the boundary. Overall, the distributional view reinforces the BA results by showing that severity weighting is the key step that turns MLTC layer information into a more stable and more discriminative indicator for identifying the behavioral risk.

5.5. Robustness and sensitivity analysis

To assess the stability of the proposed framework, we examine how the index behaves under three implementation and modeling choices. First, because the thinning stage retains coefficients from a selected starting phase, changing the offset can alter the thinned portfolio sample and fitted severity model; we therefore examine thinning-phase sensitivity by refitting the severity model over all possible offsets. Second, because $\gamma$ controls the dispersion of severity weights across fitted tail layers, we examine the sensitivity of the proposed index to $\gamma$ . Third, because residual dependence after thinning may induce overdispersion in the MLTCs, we assess robustness to the frequency model specification by comparing a Negative Binomial–Beta alternative. Figure 6 summarizes these sensitivity and robustness checks.

Figure 6.

Robustness and sensitivity diagnostics. Left: Model C BA over $\gamma\in\{0.1,\ldots,2.0\}$ under the Poisson–Gamma and Negative Binomial–Beta classifiers; the dashed line is the unweighted Model B baseline (BA $=76.0\%$ ). Right: Spearman’s correlation of trip-level index rankings across thinning offsets relative to offset 0 under $L_{\mathrm{global}}=15$ .

We first assess the sensitivity of the severity model and resulting indices when the thinning phase changes. Because the common lag is $L_{\mathrm{global}}=15$ , there are 15 deterministic starting offsets, $s=0,\ldots,14$ . For each offset, we reconstruct the thinned portfolio sample, refit the Gaussian–Uniform severity model, and recompute the trip-level indices using the fixed value $\gamma^\star=0.9$ . The fitted layer complexity changes only mildly across offsets: the selected numbers of layers vary between $M^-\in\{4,5\}$ and $M^+\in\{4,5\}$ , corresponding to 8–9 total severity layers. More importantly, the downstream trip rankings are highly stable: using the offset-0 ranking as reference, the minimum Spearman rank correlation is $0.8972$ , and the average is $0.9468$ over the non-reference offsets. Thus, this empirically indicates that the relative ordering of the trip-level indices is not an artifact of a particular thinning phase.

We then assess whether the classification gain of Model C is sensitive to $\gamma$ or to the assumed frequency model. For each fixed $\gamma\in\{0.1,0.2,\ldots,2.0\}$ , we repeat the LODO CV under two classifiers. Figure 6 (left panel) shows that Model C remains at least as accurate as the Model B baseline for 17 of the 20 values under both the Poisson–Gamma and Negative Binomial–Beta classifiers. Under the Poisson–Gamma classifier, Model C reaches its maximum BA of $92.71\%$ at $\gamma=0.9$ . Under the Negative Binomial–Beta classifier, the maximum BA is $90.54\%$ at $\gamma=0.7,0.8,$ and $0.9$ . Hence, the discrimination improvement from severity weighting is not confined to a single tuning value and persists under an alternative frequency model that explicitly allows overdispersed counts.

The shape of the $\gamma$ curve is also informative. When $\gamma$ is very small, especially at $\gamma=0.1$ , the severity weights are nearly uniform, so the severity component is effectively muted and Model C essentially coincides with the unweighted Model B benchmark. As $\gamma$ increases into the intermediate range, rarity-based weighting becomes strong enough to make severe tail-layer behavior count more heavily without collapsing the signal to only the deepest layers; this is precisely the range in which the class distributions are most cleanly separated. By contrast, for very large values of $\gamma$ , the weights become too concentrated on only the deepest layers, and performance drops below the Model B baseline at $\gamma=1.8,1.9,$ and $2.0$ because shallower but still informative layers are underused. The selected value $\gamma^\star=0.9$ balances the use of the full layer-wise severity structure against the risk of over-concentration.

The comparison between the Poisson–Gamma and Negative Binomial–Beta classifiers leads to the same qualitative conclusion. Although the Negative Binomial–Beta specification allows overdispersion, it does not improve on the Poisson–Gamma maximum in this controlled dataset, and both frequency models produce the same broad sensitivity pattern across $\gamma$ . This suggests that the principal source of the gain is the portfolio-level severity weighting itself, not the particular distributional assumption used for the frequency model. In other words, the improvement comes from converting layer-wise MLTC counts into a severity-aware representation that better separates normal and risky driving while keeping predictions more concentrated within each class.

Taken together, these checks show that the proposed behavioral risk index is not driven by a particular thinning offset, an isolated value of $\gamma$ , or the Poisson–Gamma frequency model specification. The thinning-phase analysis supports the stability of the trip-level ranking, while the $\gamma$ and frequency model analyses show that severity weighting provides a robust improvement over the unweighted layer-wise benchmark. These findings strengthen the interpretation that the empirical performance of Model C reflects persistent portfolio-level severity information captured by the MLTC representation.