* [1] [author] Asmussen, SørenS. (2003). Applied probability and queues, second ed. Applications of Mathematics (New York) 51\. Springer-Verlag, New York. Stochastic Modelling and Applied Probability. MR1978607 (2004f:60001)
  * [2] [author] Baccelli, F.F. and Brémaud, P.P. (2003). Elements of queueing theory: Palm martingale calculus and stochastic recurrences, 2nd ed. Applications of Mathematics 26\. Springer, Berlin.
  * [3] [author] Barbour, A. D.A. D. (1990). Stein’s method for diffusion approximations. Probab. Theory and Related Fields 84 297-322. 10.1007/BF01197887
  * [4] [author] Barbour, ADA., Ross, NathanN. and Zheng, GuangquG. (2023). Stein’s method, Gaussian processes and Palm measures, with applications to queueing. The Annals of Applied Probability 33 3835–3871.
  * [5] [author] Barbour, A. D.A. D. (1988). Stein’s Method and Poisson Process Convergence. Journal of Appl. Probab. 25 175-184.
  * [6] [author] Bertsimas, D.D. and Gamarnik, D.D. (2022). Queueing theory: classical and modern methods, 1st ed. Dynamic Ideas LLC, Belmont, MA.
  * [7] [author] Besançon, EustacheE., Decreusefond, LaurentL. and Moyal, PascalP. (2020). Stein’s method for diffusive limits of queueing processes. Queueing Systems 95 173–201.
  * [8] [author] Boon, MAAM., Janssen, AJEMA. and van Leeuwaarden, JSHJ. (2023). Heavy-traffic single-server queues and the transform method. Indagationes Mathematicae 34 1014–1037.
  * [9] [author] Braverman, AntonA. (2022). The Prelimit Generator Comparison Approach of Stein’s Method. Stochastic Systems 12 181-204. 10.1287/stsy.2021.0085
  * [10] [author] Braverman, AntonA., Dai, J. G.J. G. and Miyazawa, MasakiyoM. (2017). Heavy traffic approximation for the stationary distribution of a generalized Jackson network: the BAR approach. Stochastic Systems 7 143–196.
  * [11] [author] Braverman, AntonA., Dai, JGJ. and Miyazawa, MasakiyoM. (2024). The BAR approach for multiclass queueing networks with SBP service policies. Stochastic Systems.
  * [12] [author] Braverman, AntonA. and Dai, J. G.J. G. (2017). Stein’s method for steady-state diffusion approximations of ${M}/\mathit{Ph}/n+{M}$ systems. Ann. of Appl. Probab. 27 550-581. 10.1214/16-AAP1211
  * [13] [author] Braverman, A.A., Dai, J. G.J. G. and Feng, J.J. (2016). Stein’s method for steady-state diffusion approximations: An introduction through the Erlang-A and Erlang-C models. Stoch. Syst. 6 301-366.
  * [14] [author] Brown, Timothy C.T. C. and Xia, AihuaA. (2001). Stein’s method and birth-death processes. Ann. Probab. 29 1373–1403. 10.1214/aop/1015345606
  * [15] [author] Chen, YanY. and Whitt, WardW. (2020). Extremal models for the GI/ GI/ K waiting-time tail-probability decay rate. Operations Research Letters 48 770–776.
  * [16] [author] Chen, YanY. and Whitt, WardW. (2021). Extremal GI/GI/1 queues given two moments: exploiting Tchebycheff systems. Queueing Systems 97 101–124.
  * [17] [author] Chen, YanY. and Whitt, WardW. (2022). Correction to: Extremal GI/GI/1 queues given two moments: exploiting Tchebycheff systems. Queueing Systems 102 553–556.
  * [18] [author] Chen, YanY. and Whitt, WardW. (2022). SET-VALUED PERFORMANCE APPROXIMATIONS FOR THE QUEUE GIVEN PARTIAL INFORMATION. Probability in the Engineering and Informational Sciences 36 378–400.
  * [19] [author] Dai, J. G.J. G. and Meyn, Sean P.S. P. (1995). Stability and convergence of moments for multiclass queueing networks via fluid limit models. IEEE Transactions on Automatic Control 40 1889–1904.
  * [20] [author] Daley, DJD., Kreinin, A YaA. Y. and Trengove, CDC. (1992). Inequalities concerning the waiting time in single-server queues: A survey. Oxford Statistical Science Series 1 177–177.
  * [21] [author] Gast, NicolasN. (2017). Expected values estimated via mean-field approximation are 1/n-accurate. Proceedings of the ACM on Measurement and Analysis of Computing Systems 1 1–26.
  * [22] [author] Gaunt, Robert E.R. E. and Walton, NeilN. (2020). Stein’s method for the single server queue in heavy traffic. Statistics & Probability Letters 156 108566\. https://doi.org/10.1016/j.spl.2019.108566
  * [23] [author] Götze, F.F. (1991). On the rate of convergence in the multivariate CLT. Ann. Probab. 19 724–739. 10.1214/aop/1176990448
  * [24] [author] Gurvich, ItaiI. (2014). Diffusion models and steady-state approximations for exponentially ergodic Markovian queues. Ann. Appl. Probab. 24 2527–2559. 10.1214/13-AAP984
  * [25] [author] Gurvich, ItaiI., Huang, JunfeiJ. and Mandelbaum, AvishaiA. (2014). Excursion-Based Universal Approximations for the Erlang-A Queue in Steady-State. Math. Oper. Res. 39 325-373. 10.1287/moor.2013.0606
  * [26] [author] Gut, AllanA. (1974). On the moments and limit distributions of some first passage times. The Annals of Probability 277–308.
  * [27] [author] Harrison, J. MichaelJ. M. and Reiman, Martin I.M. I. (1981). Reflected Brownian motion on an orthant. Ann. Probab. 9 302–308. 606992
  * [28] [author] Huang, JunfeiJ. and Gurvich, ItaiI. (2018). Beyond Heavy-Traffic Regimes: Universal Bounds and Controls for the Single-Server Queue. Oper. Res. 66 1168-1188. 10.1287/opre.2017.1715
  * [29] [author] Kingman, J. F. C.J. F. C. (1961). The single server queue in heavy traffic. Mathematical Proceedings of the Cambridge Philosophical Society 57 902-904. 10.1017/S0305004100036094
  * [30] [author] Kingman, J. F. C.J. F. C. (1962). On queues in heavy traffic. J. Roy. Statist. Soc. Ser. B 24 383–392. 0148146
  * [31] [author] Li, JingwenJ. and Ou, JihongJ. (1995). Characterizing the idle-period distribution of GI/G/1 queues. Journal of applied probability 32 247–255.
  * [32] [author] Loulou, RichardR. (1978). An explicit upper bound for the mean busy period in a $GI/G/1$ queue. Journal of Applied Probability 15 452–455. 10.2307/3213419
  * [33] [author] Miyazawa, MasakiyoM. (1994). Rate conservation laws: a survey. Queueing Systems 15 1-58.
  * [34] [author] Miyazawa, MasakiyoM. (2015). Diffusion approximation for stationary analysis of queues and their networks: a review. J. Oper. Res. Soc. Japan 58 104–148. 10.15807/jorsj.58.104 3363214
  * [35] [author] Miyazawa, MasakiyoM. (2017). A unified approach for large queue asymptotics in a heterogeneous multiserver queue. Adv. in Appl. Probab. 49 182–220. 10.1017/apr.2016.84 3631221
  * [36] [author] Peköz, Erol A.E. A. and Röllin, AdrianA. (2011). New rates for exponential approximation and the theorems of Rényi and Yaglom. The Annals of Probability 39 587 – 608. 10.1214/10-AOP559
  * [37] [author] Ross, NathanN. (2011). Fundamentals of Stein’s method. Probab. Surv. 8 210–293. 10.1214/11-PS182 2861132 (2012k:60079)
  * [38] [author] Stein, CharlesC. (1972). A bound for the error in the normal approximation to the distribution of a sum of dependent random variables. In Proceedings of the Sixth Berkeley Symposium on Mathematical Statistics and Probability, Volume 2: Probability Theory 583–602. University of California Press, Berkeley, Calif.
  * [39] [author] Stolyar, Alexander L.A. L. (2015). Tightness of stationary distributions of a flexible-server system in the Halfin-Whitt asymptotic regime. Stoch. Syst. 5 239–267. 10.1214/14-SSY139
  * [40] [author] Wolff, Ronald WR. W. and Wang, Chia-LiC.-L. (2003). Idle period approximations and bounds for the GI/G/1 queue. Advances in Applied Probability 35 773–792.
  * [41] [author] Ying, LeiL. (2016). On the Approximation Error of Mean-Field Models. In Proceedings of the 2016 ACM SIGMETRICS International Conference on Measurement and Modeling of Computer Science 285–297. ACM.
  * [42] [author] Ying, LeiL. (2017). Stein’s Method for Mean Field Approximations in Light and Heavy Traffic Regimes. Proc. ACM Meas. Anal. Comput. Syst. 1 12:1–12:27. 10.1145/3084449
  * [43] [author] Zhou, XingyuX. and Shroff, NessN. (2020). A Note on Stein’s Method for Heavy-Traffic Analysis. arXiv preprint arXiv:2003.06454.


# Irregular proton injection to high energies at interplanetary shocks

Domenico Trotta The Blackett Laboratory, Department of Physics, Imperial
College London, London SW7 2AZ, UK  Timothy S. Horbury The Blackett
Laboratory, Department of Physics, Imperial College London, London SW7 2AZ, UK
David Lario Heliophysics Science Division, NASA Goddard Space Flight Center,
Greenbelt, MD 20771, USA  Rami Vainio Department of Physics and Astronomy,
University of Turku, Finland  Nina Dresing Department of Physics and
Astronomy, University of Turku, Finland  Andrew Dimmock Swedish Institute of
Space Physics, Uppsala, Sweden  Joe Giacalone Lunar and Planetary Laboratory,
University of Arizona, Tucson, USA  Heli Hietala School of Physics and
Astronomy, Queen Mary University of London, London E1 4NS, UK  Robert F.
Wimmer-Schweingruber Institute of Experimental and Applied Physics, Kiel
University, 24118 Kiel, Germany  Lars Berger Institute of Experimental and
Applied Physics, Kiel University, 24118 Kiel, Germany  Liu Yang Institute of
Experimental and Applied Physics, Kiel University, 24118 Kiel, Germany

(Received XXX; Revised YYY; Accepted ZZZ)

###### Abstract

How thermal particles are accelerated to suprathermal energies is an unsolved
issue, crucial for many astrophysical systems. We report novel observations of
irregular, dispersive enhancements of the suprathermal particle population
upstream of a high-Mach number interplanetary shock. We interpret the observed
behavior as irregular “injections” of suprathermal particles resulting from
shock front irregularities. Our findings, directly compared to self-consistent
simulation results, provide important insights for the study of remote
astrophysical systems where shock structuring is often neglected.

Acceleration of particles — plasmas – shock waves — Sun: heliosphere — Sun:
solar wind

††journal: ApJL

revtex4-1.clsrevtex4-2.cls

##  1 Introduction

Collisionless shock waves are fundamental sources of energetic particles,
which are ubiquitously present in our universe and pivotal to explain many of
its features, such as the non-thermal radiation emission common to many
astrophysical sources, as revealed by decades of remote and direct
observations (Reames, 1999; Amato & Blasi, 2018). Particle acceleration to
suprathermal energies from thermal plasma, less understood than particle
acceleration starting from an already energised population, remains a puzzle,
and has been object of extensive theoretical and numerical investigations
(Drury, 1983; Caprioli & Spitkovsky, 2014; Trotta et al., 2021).

Shocks in the heliosphere, unique as directly accessible by spacecraft
(Richter et al., 1985), provide the missing link to remote observations of
astrophysical systems. Direct observations of the Earth’s bow shock using
single and multi-spacecraft approaches (e.g., Johlander et al., 2016) reveal a
complex scenario of energy conversion and particle acceleration at the shock
transition (Amano et al., 2020; Schwartz et al., 2022). The emerging picture,
well supported by theory and modelling, is that small scale irregularities in
the spatial and temporal evolution of the shock environment (Greensadt et al.,
1980; Matsumoto et al., 2015) are fundamental for efficient ion injection to
high energies (Dimmock et al., 2019). This idea of irregular particle
injection has been investigated in the past for the Earth’s bow shock
(Madanian et al., 2021) and in numerical simulations (Guo & Giacalone, 2013),
thus suggesting that particle behaviour at shocks is much more complex than
what is expected neglecting space-time irregularities, as suggested by early
theoretical and numerical works (Decker, 1990; Ao et al., 2008; Lu et al.,
2009).

Such a complex picture is not as well observed and understood for shocks
beyond the Earth’s bow shock. In particular, shock structuring at
Interplanetary (IP) shocks, generated as a consequence of phenomena such as
Coronal Mass Ejections (CMEs, Gosling et al., 1974) and its role in particle
acceleration remains elusive (Blanco-Cano et al., 2016; Kajdič et al., 2019).
IP shocks are generally weaker and have larger radii of curvature with respect
to Earth’s bow shock, allowing for direct observations of collisionless shocks
in profoundly different regimes (e.g., Kilpua et al., 2015; Yang et al.,
2020), and are more relevant to astrophysical environments such as galaxy
cluster shocks, where shock irregularities are not resolved, but they are
likely to play a crucial role in efficient particle acceleration (Brunetti &
Jones, 2014). Therefore, the study of particle injection at IP shocks is
fundamental to test our current understanding built on Earth’s bow shock, as
well for addressing shocks at objects currently beyond reach. This paper
demonstrates that, in order to address the suprathermal particle production
upstream of supercritical collisionless shocks, the inherent variability of
the injection process in both time and space must be taken into account.

The Solar Orbiter mission (SolO, Müller et al., 2020) probes the inner
heliosphere with unprecedented levels of time-energy resolution for energetic
particles, thus opening a new observational window for particle acceleration.
In this work, we study the acceleration of low-energy ($\sim$ 1 keV) particles
to supra-thermal energies ($\sim$ 50 keV) at a strong IP shock observed by
SolO at heliocentric distance of about 0.8 AU on 2021 October 30th at 22:02:07
UT. We use the SupraThermal Electrons and Protons sensor (STEP) of the
Energetic Particle Detector (EPD) suite (Rodríguez-Pacheco et al., 2020),
measuring particles in the 6 - 60 keV energy range (close to the injection
range), at the very high time resolution of 1 s, close to suprathermal
particle gyroscales. Our work exploits such novel, previously unavailable
datasets for suprathermal particles upstream of IP shocks. We resolve upstream
enhancements in the suprathermal particle population with dispersive velocity
signatures, and link them to irregular proton injection along the shock front.
Our findings are corroborated by kinetic simulations showing similar irregular
proton energization upstream close to the shock, thus elucidating the
mechanisms responsible for this behaviour. This letter is organised as
follows: results are presented in Section 2. SolO observations are shown and
discussed in Section 2.1, while modelling results are reported in 2.2. The
conclusions are in Section 3.

##  2 Results

###  2.1 Solar Orbiter Observations

Figure 1: Event overview. (a) EPD-Electron Proton Telescope (EPT) particle
flux (sunward aperture). (b) EPD-STEP particle flux (magnet channel averaged
over the entire field of view). (c) Pitch angle distributions for ions with an
energy of 0.011 - 0.019 MeV in the spacecraft frame. (d) Time profile of the
STEP energy flux in the 0.012 - 0.015 MeV energy channel at full resolution
(blue), and time-averaged using a 1 minute window. (e) SWA-PAS ion energy flux
(Owen et al., 2020). (f) SWA-PAS proton density. (g) MAG burst magnetic field
data in RTN coordinates (Horbury et al., 2020). The magenta line marks the
shock crossing, and the black rectangle highlights the dispersive energetic
particle enhancements observed by STEP. Differential fluxes are in
$\mathrm{E^{2}\cdot cm^{-2}s^{-1}sr^{-1}MeV}$ for the EPD instruments and
$\mathrm{cm^{-2}s^{-1}eV}$ for PAS.

Fig. 1 shows a 30 minute overview across the shock transition. Panels (a)-(b)
reveal the presence of shock accelerated particles at energies of up to 100
keV, while particle fluxes at higher energies do not respond to the shock
passage. At these high energies the fluxes were enhanced following a large
Solar Energetic Particle (SEP) event (see Klein et al., 2022).

The most striking feature of the period prior to the shock arrival at SolO is
the irregular energetic particle enhancements particularly evident at 10 - 30
keV energies (Fig. 1 (b), black box), found in the time interval $~{}\sim$ 15
minutes before the shock crossing, corresponding to $2\times 10^{5}$ km or
2500 ion inertial lengths, $d_{i}$. These particle enhancements have the novel
feature of being dispersive in energy and are the focus of this work. The
typical timescales at which the irregularities are observed are of 10-20
seconds, corresponding to spatial scales of about 50 $d_{i}$. Such signatures
were previously inaccessible to observations, as shown in Fig. 1 (c), where
the time profile of ion differential flux in the 0.012 - 0.015 MeV channel,
rising exponentially up to the shock (Giacalone, 2012), is shown at full
resolution (blue) and averaged using a $\sim$ 1 minute window, typical of
previous IP shock measurements. Fig. 1(d) shows pitch angle intensities for
0.011 – 0.019 MeV ions (i.e., energies at which the irregular enhancements are
observed). Pitch angles are computed in the plasma rest frame assuming that
all ions are protons, and performing a Compton-Getting correction (Compton &
Getting, 1935a), thereby combining magnetic field data from the magnetometer
(MAG, Horbury et al., 2020), and solar wind plasma data from the Proton and
Alpha particle Sensor (PAS) on the Solar Wind Analyser (SWA) instrument suite
(Owen et al., 2020), and particle data from EPD/STEP (Yang, L. et al., 2023).
For the interval studied, low pitch angles are in the 30∘ field of view of
STEP, relevant for shock reflected particles. The irregular enhancements of
energetic particles are field aligned, as is evident for the strongest signal
close to the shock transition. The flux enhancement visible in PAS (Fig. 1(e))
at lower energies starting immediately before the shock (22:00 UT) also
reveals a field-aligned population. The study of the PAS low-energy population
and the behaviour very close to the shock transition is object of another
investigation (Dimmock et al., 2023).

Figure 2: Left: Spectrogram of the irregular signal in seconds from shock vs
$1/v$ axes, with the velocity dispersion shown by the solid magenta line
(top). Time series showing the local $\theta_{Bn}(t)$ angle. The red and grey
dashed lines represent the average $\theta_{Bn}$ and a 90∘ angle, respectively
(bottom). Right: Cartoon showing the corrugated shock front with local shock
normal, trajectory of a reflected particle and the Solar Orbiter trajectory
(SolO model: esa.com).

The magnetic field reveals a wave foreshock $\sim$ 2 minutes upstream of the
shock, in conjunction with a population of low-energy ($\sim$ 4 keV) reflected
particles seen by SWA/PAS, visible as the light blue enhancement in Fig. 1(e)
around 22:00 UT. Interestingly, the magnetic field is quieter where signals of
irregular injection are found, indicating that efficient particle scattering
may be reduced in this region (Lario et al., 2022). In this “quiet” shock
upstream, we found two structures compatible with shocklets in the process of
steepening ($\sim$21:57 UT), very rarely observed at IP shocks (Wilson et al.,
2009; Trotta et al., 2023a).

The shock parameters were estimated using upstream/downstream averaging
windows varied systematically between 1 and 8 minutes (Trotta et al., 2022a).
The shock was oblique, with a normal angle $\theta_{Bn}=44\pm 1.5^{\circ}$
(obtained with the Mixed Mode 3 technique (MX3 Paschmann & Schwartz, 2000),
compatible with MX1,2 and Magnetic Coplanarity). The shock speed in the
spacecraft frame and along the shock normal is $\mathrm{V_{shock}=400\pm
5\,km/s}$. The shock Alfvénic and fast magnetosonic Mach numbers are $\rm
M_{A}\sim 7.6$ and $\rm M_{fms}\sim 4.6$, respectively. Thus, the event
provides us with the opportunity to study a shock with particularly high Mach
number in comparison with other IP shocks, while the shock speed is moderate
with respect to typical IP shocks (Kilpua et al., 2015). The shock is
supercritical, and therefore expected to have a corrugated, rippled front
(Trotta & Burgess, 2019; Kajdic et al., 2021). The presence of reflected
particles, enhanced wave activity in close proximity (1 minute) to the shock
transition and upstream shocklets in the process of steepening is consistent
with the local shock parameters (Blanco-Cano et al., 2016).

To further elucidate the dispersive nature of the suprathermal particles, we
show the STEP energy spectrogram in $1/v$ vs $t$ space (Fig. 2). Here,
particle speeds are referred to the center of the relative energy bin and
computed in the spacecraft rest frame, assuming that all particles detected
are protons (see Wimmer-Schweingruber et al., 2021, for further details).
During the period of irregular particle enhancements, we also combined
magnetic field and plasma data to compute the particle pitch angles in the
solar wind frame (Compton & Getting, 1935b), revealing that the particles
detected by STEP are closely aligned with the field (not shown here).
Interestingly, by visual inspection, it can be seen that these dispersive
signals are shallower going far upstream, consistent with the fact that they
are injected from more distant regions of the shock.

The dispersive flux enhancements are associated with irregular acceleration of
protons along the shock front. Indeed, due to their dispersive nature, the
particles detected by STEP cannot be continuously produced at the shock and
propagated upstream, but they must come from a source that is only temporarily
magnetically connected to the spacecraft due to time and/or space
irregularities. Then, the fastest particles produced at the irregular source
are detected first by the spacecraft, followed by the slower ones, yielding
the observed dispersive behaviour. Given the short timescales at which
energetic particle enhancements are observed with respect to the shock and the
quiet behaviour of upstream magnetic field in the 10 minutes upstream of the
shock, we assume that particles do not undergo significant scattering from
their (irregular) production to the detection at SolO. It is then natural to
investigate the connection with the shock. The bottom-left panel of Fig. 2
shows the local
$\theta_{Bn}(t)\equiv\cos^{-1}\left(\mathbf{B}(t)\cdot\hat{\mathbf{n}}_{\rm
shock}/|\mathbf{B}(t)|\right)$ changing significantly when the dispersive
signals are observed, indicating that the spacecraft was indeed connected to
different portions of the (corrugated) shock front, which in turn is expected
to respond rapidly to upstream changes, as recent simulation work elucidated
(e.g., Trotta et al., 2023b). Note that, given the single-spacecraft nature of
the observations, the average shock normal computed with MX3 for both local
and average $\theta_{Bn}$estimation was used.

To further support this idea, similarly to Velocity Dispersion Analyses (VDA)
used to determine the injection time of SEP events (e.g., Lintunen & Vainio,
R., 2004; Dresing et al., 2023), we chose the clearest dispersive signal
($\sim$ 100 seconds upstream of the shock) and we superimpose the following
relation (indicated by the magenta line in Fig. 2):

$t_{\mathrm{O}}(v)=t_{i}+\frac{s}{v},$ (1)

where $t_{\mathrm{O}}$ represents the time at which the flux enhancement is
observed for a certain speed $v$, $t_{i}$ is the time of injection at the
source, and $s$ is the distance travelled by the particles from the source to
the spacecraft. Thus, the argument is that the dispersive signals are due to
accelerated particles produced by different portions of the shock front
temporarily connected with the spacecraft, as sketched in Fig. 2 (right). We
note that, due to the very high energy-time resolution of STEP, it was
possible to perform the VDA on such small ($\sim$seconds) time scales.
Determining $t_{i}$ based on the time when the highest energy particles are
observed ($t_{i}\sim-130s$), the source distance that we obtain through
Equation 1 is $s\approx 4\times 10^{4}$ km ($\sim 500d_{i}$), compatible with
their generation at the approaching shock, for which we would expect
$s\sim\mathrm{V_{shock}\Delta t/sin(\theta_{Bn})}$, where $\mathrm{V_{shock}}$
is the average shock speed, $\Delta t$ is the time delay between the
observation of the dispersive signal and the shock passage. This is also
compatible with the fact that the other dispersive signals observed further
upstream, such as the one before 21:54, about 500 seconds upstream of the
shock (see Fig. 2), show a shallower inclination, though a more precise,
quantitative analysis of this behaviour is complicated by the high noise
levels of the observation, and will be the object of later statistical
investigation employing more shock candidates (Yang, L. et al., 2023).

Figure 3: Top: Simulation snapshot of proton density (colormap). The inset
shows a zoom around the shock transition (grey), and the local shock position
is superimposed, with a colormap correesponding to the local $\theta_{Bn}$.
Bottom: Density map of upstream superathermal protons (colormap) and magnetic
field lines (magenta) computed at the same simulation time as (a). The inset
shows the upstream particle energy spectrum, with the dashed blue lines
indicating the suprathermal energy range considered.

###  2.2 Shock Modelling

Further insights about shock front irregularities are limited by the single-
spacecraft nature of these observations. Therefore, we employ 2.5-dimensional
kinetic simulations, with parameters compatible with the observed ones, to
model the details of the shock transition, where proton injection to
suprathermal energies takes place, relevant to our interpretation of the
dispersive signals and enabling us to see how the shock surface and normal
behave at small scales (see Fig. 2). In the simulations, protons are modelled
as macroparticles and advanced with the Particle-In-Cell (PIC) method, while
the electrons are modelled as a massless, charge-neutralizing fluid (Trotta et
al., 2020).

In the model, distances are normalised to the ion inertial length $d_{i}$,
times to the upstream inverse cyclotron frequency ${\Omega_{ci}}^{-1}$,
velocity to the Alfvén speed $v_{A}$, and the magnetic field and density to
their upstream values $B_{0}$ and $n_{0}$. The shock is launched with the
injection method (Quest, 1985), where an upstream flow speed
$V_{\mathrm{in}}=4.5v_{A}$ was chosen, corresponding to $M_{A}\sim 6$. The
shock nominal $\theta_{Bn}$ is 45∘. The simulation domain is 512 $d_{i}$
$\times$ 512 $d_{i}$, with resolution $\Delta x$ = $\Delta y$ = 0.5 $d_{i}$
and a particle time-step $\Delta t_{pa}$ = 0.01 $\Omega_{ci}^{-1}$. The number
of particles per cell used is always greater than 300. This choice of
parameters is compatible with the local properties of the IP shock as
estimated from the SolO measurements. However, inherent variability routinely
found in the simulations at small scales and in the observations at larger
scales must be considered when comparing numerical and observational results.
We note that these simulations are initialised with a laminar upstream, and
therefore the fluctuations that impact the shock are self-generated (due to
particle reflection and subsequent upstream propagation). An exhaustive
characterization of these self-induced fluctuations is discussed in Kajdic et
al. (2021).

Simulation results are shown in Fig. 3. In the top panel, we present the
proton density for a simulation snapshot where the shock transition is well-
developed, showing the strongly perturbed character of the shock front. In
such an irregular shock transition, particle dynamics become extremely complex
(e.g., Lembege & Savoini, 1992). To further elucidate the irregularities of
the shock front, we computed the shock position in the simulation domain (with
the criterion $B>3B_{0}$, as in Trotta et al. (2023b)) and evaluated the local
$\theta_{Bn}$ along it (Fig. 3(a), inset), showing high variability (see the
sketch in Fig. 2).

In the bottom panel of Fig. 3, we study the self-consistently shock-
accelerated protons. The upstream energy spectrum is shown in the inset, with
a peak at the inflow population energies and a suprathermal tail due to the
accelerated protons. To address particle injection, we analyse the upstream
spatial distribution of such suprathermal protons (Fig. 3(b)) at the energies
highlighted in the inset, which are a factor of 10 larger than the typical
energies of particles in the upstream inflow population, in a similar fashion
as the energy separation between the STEP energies at which the irregular
enhancements are observed ($\sim$ 10 keV) and the Solar wind population
energies measured by PAS ($\sim$ 1 keV) . It can be seen that suprathermal
particles are not distributed uniformly, and their spatial distribution varies
with their locations along the shock front, another indication of irregular
injection. Furthermore, we observed that the length scale of the
irregularities is of 50 $d_{i}$, directly comparable with the irregularities
seen in the STEP fluxes (see Fig. 1). Higher energy particles also show
irregularities.

##  3 Conclusions

We studied irregular particle acceleration from the thermal plasma using novel
SolO observations. Particle injection to high energies is an extremely
important issue for a large collection of astrophysical systems making the
SolO shock on 2021 October 30th an excellent event to tackle this interesting
problem. The capabilities of the SolO EPD suite were exploited to probe the
complex shock front behaviour in the poorly investigated IP shock case. From
this point of view, in-situ observations of irregular particle enhancements
have been used as a tool to address the (remote) structuring of the shock, an
information not available by simply looking at the spacecraft shock crossing
of in one point in space and time. Such an approach is reminiscent to the ones
used to reconstruct the properties of SEP events (Krucker et al., 1999), and
even to the ones looking at the properties of the heliospheric termination
shock with the Interstellar Boundary Explorer mission (IBEX, McComas et al.,
2009), where particles produced at different portions of the shock are used to
understand its dynamics (Zirnstein et al., 2022).

The hybrid kinetic simulations are consistent with this complex scenario of
proton acceleration, with irregularly distributed suprathermal particles along
the shock front, an invaluable tool to elucidate the small-scale behaviour of
this IP shock and of shock transitions in a variety of astrophysical systems.
Our model highlights the very small-scale behaviour of the shock, but neglects
other effects like pre-existing turbulence and interplanetary disturbances
that may be important (Lario & Decker, 2002; Trotta et al., 2022b; Nakanotani
et al., 2022; Trotta et al., 2023b). The direct investigation of shock
acceleration in systems other than the Earth’s bow shock (having a small
radius of curvature and many other properties important for planetary bow
shocks) is important to build a comprehensive understanding of collisionless
shocks energetics. This work significantly strengthens an evolving theory of
collisionless shock acceleration. Combining high resolution energetic particle
data upstream of heliospheric shocks with hybrid simulations, we have shown,
for interplanetary shocks, that the inherent variability of the injection
process in both time and space must be considered to solve the problem of how
suprathermal particle injection occurs in astrophysical systems. The process
analysed here is general, as it does not depend on how shock irregularities
are generated. Indeed, this study is relevant for astrophysical systems where
shock front irregularities cannot be resolved but are likely to play an
important role for particle acceleration from the thermal distribution, such
as galaxy cluster shocks, where efficient particle acceleration, which is
inferred to happen at very large, $\sim$ Mpc scales, remains a puzzle,
particularly in the absence of pre-existing cosmic rays (Botteon et al.,
2020).

This study has received funding from the European Unions Horizon 2020 research
and innovation programme under grant agreement No. 101004159 (SERPENTINE,
www.serpentine-h2020.eu). Part of this work was performed using the DiRAC Data
Intensive service at Leicester, operated by the University of Leicester IT
Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk),
under the project “dp031 Turbulence, Shocks and Dissipation in Space Plasmas”.
N.D. acknowledges the support of the Academy of Finland (SHOCKSEE, grant nr.
346902). H.H. is supported by the Royal Society University Research Fellowship
URF\R1\180671. D.L. acknowledges support from NASA Living With a Star (LWS)
program NNH19ZDA001N-LWS, and the Goddard Space Flight Center Heliophysics
Innovation Fund (HIF) program.

## References

  * Amano et al. (2020) Amano, T., Katou, T., Kitamura, N., et al. 2020, Phys. Rev. Lett., 124, 065101, doi: 10.1103/PhysRevLett.124.065101
  * Amato & Blasi (2018) Amato, E., & Blasi, P. 2018, Advances in Space Research, 62, 2731 , doi: https://doi.org/10.1016/j.asr.2017.04.019
  * Ao et al. (2008) Ao, X., Zank, G. P., Pogorelov, N. V., & Shaikh, D. 2008, Physics of Fluids, 20, 127102, doi: 10.1063/1.3041706
  * Blanco-Cano et al. (2016) Blanco-Cano, X., Kajdič, P., Aguilar-Rodríguez, E., et al. 2016, Journal of Geophysical Research: Space Physics, 121, 992, doi: https://doi.org/10.1002/2015JA021645
  * Botteon et al. (2020) Botteon, Brunetti, G., Ryu, D., & Roh, S. 2020, A&A, 634, A64, doi: 10.1051/0004-6361/201936216
  * Brunetti & Jones (2014) Brunetti, G., & Jones, T. W. 2014, International Journal of Modern Physics D, 23, 1430007, doi: 10.1142/S0218271814300079
  * Caprioli & Spitkovsky (2014) Caprioli, D., & Spitkovsky, A. 2014, The Astrophysical Journal, 783, 91, doi: 10.1088/0004-637X/783/2/91
  * Compton & Getting (1935a) Compton, A. H., & Getting, I. A. 1935a, Phys. Rev., 47, 817, doi: 10.1103/PhysRev.47.817
  * Compton & Getting (1935b) —. 1935b, Phys. Rev., 47, 817, doi: 10.1103/PhysRev.47.817
  * Decker (1990) Decker, R. B. 1990, Journal of Geophysical Research: Space Physics, 95, 11993, doi: https://doi.org/10.1029/JA095iA08p11993
  * Dimmock et al. (2023) Dimmock, Gedalin, M., Lalti, A., et al. 2023, A&A, doi: 10.1051/0004-6361/202347006
  * Dimmock et al. (2019) Dimmock, A. P., Russell, C. T., Sagdeev, R. Z., et al. 2019, Science Advances, 5, eaau9926, doi: 10.1126/sciadv.aau9926
  * Dresing et al. (2023) Dresing, N., Rodríguez-García, L., Jebaraj, I. C., et al. 2023, A&A, accepted, doi: 10.1051/0004-6361/202345938
  * Drury (1983) Drury, L. O. 1983, Reports on Progress in Physics, 46, 973, doi: 10.1088/0034-4885/46/8/002
  * Giacalone (2012) Giacalone, J. 2012, The Astrophysical Journal, 761, 28, doi: 10.1088/0004-637X/761/1/28
  * Gosling et al. (1974) Gosling, J. T., Hildner, E., MacQueen, R. M., et al. 1974, Journal of Geophysical Research, 79, 4581, doi: 10.1029/JA079I031P04581
  * Greensadt et al. (1980) Greensadt, E. W., Russell, C., Gosling, J., et al. 1980, Journal of Geophysical Research: Space Physics, 85, 2124, doi: https://doi.org/10.1029/JA085iA05p02124
  * Guo & Giacalone (2013) Guo, F., & Giacalone, J. 2013, The Astrophysical Journal, 773, 158, doi: 10.1088/0004-637X/773/2/158
  * Horbury et al. (2020) Horbury, T. S., O’Brien, H., Carrasco Blazquez, I., et al. 2020, Astronomy & Astrophysics, 642, A9, doi: 10.1051/0004-6361/201937257
  * Johlander et al. (2016) Johlander, A., Vaivads, A., Khotyaintsev, Y. V., Retinó, A., & Dandouras, I. 2016, The Astrophysical Journal Letters, 817, L4, doi: 10.3847/2041-8205/817/1/L4
  * Kajdic et al. (2021) Kajdic, P., Pfau-Kempf, Y., Turc, L., et al. 2021, Journal of Geophysical Research: Space Physics, 126, e2021JA029283, doi: https://doi.org/10.1029/2021JA029283
  * Kajdič et al. (2019) Kajdič, P., Preisser, L., Blanco-Cano, X., Burgess, D., & Trotta, D. 2019, ApJ, 874, L13, doi: 10.3847/2041-8213/ab0e84
  * Kilpua et al. (2015) Kilpua, E. K., Lumme, E., Andreeova, K., Isavnin, A., & Koskinen, H. E. 2015, Journal of Geophysical Research: Space Physics, 120, 4112, doi: 10.1002/2015JA021138
  * Klein et al. (2022) Klein, K.-L., Musset, S., Vilmer, N., et al. 2022, Astronomy & Astrophysics, 663, A173, doi: 10.1051/0004-6361/202243903
  * Krucker et al. (1999) Krucker, S., Larson, D. E., Lin, R. P., & Thompson, B. J. 1999, The Astrophysical Journal, 519, 864, doi: 10.1086/307415
  * Lario & Decker (2002) Lario, D., & Decker, R. B. 2002, Geophysical Research Letters, 29, 31, doi: https://doi.org/10.1029/2001gl014017
  * Lario et al. (2022) Lario, D., Richardson, I. G., Wilson, L. B., I., et al. 2022, ApJ, 925, 198, doi: 10.3847/1538-4357/ac3c47
  * Lembege & Savoini (1992) Lembege, B., & Savoini, P. 1992, Physics of Fluids B, 4, 3533, doi: 10.1063/1.860361
  * Lintunen & Vainio, R. (2004) Lintunen, & Vainio, R. 2004, A&A, 420, 343, doi: 10.1051/0004-6361:20034247
  * Lu et al. (2009) Lu, Q., Hu, Q., & Zank, G. P. 2009, The Astrophysical Journal, 706, 687, doi: 10.1088/0004-637X/706/1/687
  * Madanian et al. (2021) Madanian, H., Schwartz, S. J., Fuselier, S. A., et al. 2021, The Astrophysical Journal Letters, 915, L19, doi: 10.3847/2041-8213/ac0aee
  * Matsumoto et al. (2015) Matsumoto, Y., Amano, T., Kato, T. N., & Hoshino, M. 2015, Science, 347, 974, doi: 10.1126/science.1260168
  * McComas et al. (2009) McComas, D. J., Allegrini, F., Bochsler, P., et al. 2009, Space Sci. Rev., 146, 11, doi: 10.1007/s11214-009-9499-4
  * Müller et al. (2020) Müller, St. Cyr, O. C., Zouganelis, I., et al. 2020, A&A, 642, A1, doi: 10.1051/0004-6361/202038467
  * Nakanotani et al. (2022) Nakanotani, M., Zank, G. P., & Zhao, L.-L. 2022, The Astrophysical Journal, 926, 109, doi: 10.3847/1538-4357/ac4781
  * Owen et al. (2020) Owen, C. J., Bruno, R., Livi, S., et al. 2020, Astronomy & Astrophysics, 642, doi: 10.1051/0004-6361/201937259
  * Paschmann & Schwartz (2000) Paschmann, G., & Schwartz, S. J. 2000, ESA Special Publication, Vol. 449, ISSI Book on Analysis Methods for Multi-Spacecraft Data, ed. R. A. Harris, 99
  * Quest (1985) Quest, K. B. 1985, Phys. Rev. Lett., 54, 1872, doi: 10.1103/PhysRevLett.54.1872
  * Reames (1999) Reames, D. V. 1999, Space Sci. Rev., 90, 413, doi: 10.1023/A:1005105831781
  * Richter et al. (1985) Richter, A. K., Hsieh, K. C., Luttrell, A. H., Marsch, E., & Schwenn, R. 1985, Review of Interplanetary Shock Phenomena Near and within 1 AU (American Geophysical Union (AGU)), 33–50, doi: https://doi.org/10.1029/GM035p0033
  * Rodríguez-Pacheco et al. (2020) Rodríguez-Pacheco, Wimmer-Schweingruber, R. F., Mason, G. M., et al. 2020, A&A, 642, A7, doi: 10.1051/0004-6361/201935287
  * Schwartz et al. (2022) Schwartz, S. J., Goodrich, K. A., Wilson III, L. B., et al. 2022, Journal of Geophysical Research: Space Physics, 127, doi: https://doi.org/10.1029/2022JA030637
  * Trotta & Burgess (2019) Trotta, D., & Burgess, D. 2019, MNRAS, 482, 1154, doi: 10.1093/mnras/sty2756
  * Trotta et al. (2020) Trotta, D., Burgess, D., Prete, G., Perri, S., & Zimbardo, G. 2020, MNRAS, 491, 580, doi: 10.1093/mnras/stz2760
  * Trotta et al. (2023a) Trotta, D., Hietala, H., Horbury, T., et al. 2023a, Monthly Notices of the Royal Astronomical Society, 520, 437, doi: 10.1093/mnras/stad104
  * Trotta et al. (2021) Trotta, D., Valentini, F., Burgess, D., & Servidio, S. 2021, Proceedings of the National Academy of Sciences, 118, e2026764118, doi: 10.1073/pnas.2026764118
  * Trotta et al. (2022a) Trotta, D., Vuorinen, L., Hietala, H., et al. 2022a, Frontiers in Astronomy and Space Sciences, 9, doi: 10.3389/fspas.2022.1005672
  * Trotta et al. (2022b) Trotta, D., Pecora, F., Settino, A., et al. 2022b, The Astrophysical Journal, 933, 167, doi: 10.3847/1538-4357/ac7798
  * Trotta et al. (2023b) Trotta, D., Pezzi, O., Burgess, D., et al. 2023b, Monthly Notices of the Royal Astronomical Society, 525, 1856, doi: 10.1093/mnras/stad2384
  * Wilson et al. (2009) Wilson, L. B. I., Cattell, C. A., Kellogg, P. J., et al. 2009, Journal of Geophysical Research: Space Physics, 114, doi: https://doi.org/10.1029/2009JA014376
  * Wimmer-Schweingruber et al. (2021) Wimmer-Schweingruber, Janitzek, N. P., Pacheco, D., et al. 2021, A&A, 656, A22, doi: 10.1051/0004-6361/202140940
  * Yang et al. (2020) Yang, L., Berger, L., Wimmer-Schweingruber, R. F., et al. 2020, The Astrophysical Journal Letters, 888, L22, doi: 10.3847/2041-8213/ab629d
  * Yang, L. et al. (2023) Yang, L., Heidrich-Meisner, V., Berger, L., et al. 2023, A&A, 673, A73, doi: 10.1051/0004-6361/202245681
  * Zirnstein et al. (2022) Zirnstein, E. J., Shrestha, B. L., McComas, D. J., et al. 2022, Nature Astronomy, 6, 1398, doi: 10.1038/s41550-022-01798-6


11institutetext: Fabio Baione, 22institutetext: Department of Statistics,
Sapienza University of Rome, Viale Regina Elena, 295 00161 Roma, Tel. +39 06
49255317, 22email<EMAIL_ADDRESS>33institutetext: Davide
Biancalana, 44institutetext: Department of Statistics, Sapienza University of
Rome, Italy 44email<EMAIL_ADDRESS>55institutetext: Paolo De
Angelis, 66institutetext: Department of methods and models for economics
territory and finance, Sapienza University of Rome, Italy 66email:
<EMAIL_ADDRESS>
# An application of Zero-One Inflated Beta regression models for predicting
health insurance reimbursement

Fabio Baione  Davide Biancalana and Paolo De Angelis

###### Abstract

In actuarial practice the dependency between contract limitations
(deductibles, copayments) and health care expenditures are measured by the
application of the Monte Carlo simulation technique. We propose, for the same
goal, an alternative approach based on Generalized Linear Model for Location,
Scale and Shape (GAMLSS). We focus on the estimate of the ratio between the
one-year reimbursement amount (after the effect of limitations) and the one
year expenditure (before the effect of limitations). We suggest a regressive
model to investigate the relation between this response variable and a set of
covariates, such as limitations and other rating factors related to health
risk. In this way a dependency structure between reimbursement and limitations
is provided. The density function of the ratio is a mixture distribution,
indeed it can continuously assume values mass at 0 and 1, in addition to the
probability density within $\left(0,1\right)$. This random variable does not
belong to the exponential family, then an ordinary Generalized Linear Model is
not suitable. GAMLSS introduces a probability structure compliant with the
density of the response variable, in particular zero-one inflated beta density
is assumed. The latter is a mixture between a Bernoulli distribution and a
Beta distribution.

##  1 Introduction

Insurance deductibles are a very important contract boundary which aim is to
limit the abuse of reimbursement requests especially in health insurance
policies. Generally, the data set of such contracts contains only censored
and/or truncated values, but when the data contains the medical invoices, as
well as the reimbursement values, it is possible to investigate the effect of
deductibles in terms of the ratio of the truncated values and the overall
loss, briefly referred as indicated deductible relativity (IDR). A standard
reference for deductible pricing in actuarial science is in BrownLennox . A
practical approach to modeling the deductible rates in a ratemaking process is
the adoption of regression models such as Generalized Linear Models (GLM), see
Lee , FreesLee , FreesLee2  among others. Anyway, GLMs are used to predict
aggregate claims of the truncated response variable using a log deductible
covariate and their use is limited to the class of exponential family
distributions. Our aim is to propose a different approach to measure the
influence of deductibles on health care expenditures by focusing on a
regression model, where IDR is the response variable. Anyway, the probability
function of IDR in our framework belongs to a particular mixture distribution
named Zero-One Inflated Beta (ZOIB) Ospina . As a consequence we suggest to
use a regression model of Generalized Linear Model for Location, Scale and
Shape (GAMLSS) Rigby  type.

##  2 Actuarial Framework

We consider an Health insurance company, whose portfolio is allocated on
$H$-homogeneous sub-portfolios. Let:

  * •

$Y_{i}$ is the random variable (r.v.) expenditure for single episode, i.e.
before the application of deductibles;

  * •

$L_{i}$ be the r.v. reimbursement for single episode, i.e. after the
application of deductibles;

  * •

$\textbf{x}_{i}$: be the row vector of the design matrix, providing
information about policyholder or contract rating factors;

  * •

$f$ the amount of deductible ;

  * •

$M$ the amount of the out of pocket maximum.

Hence, we can state:

$L_{i}=\min\left\\{I_{\left(Y_{i}-f>0\right)}\cdot Y_{i};M\right\\}$ (1)

where $I_{\left(Y_{i}-f>0\right)}$ is the indicator r.v. for the event
$Y_{i}-f>0$. Hence, if the invoice amount is higher than $f$ the health
insurance pays the full expenditure within the limit fixed by $M$ and 0
otherwise. In order to assess the influence of the deductible and the out of
pocket maximum on the loss of the company, we focus on the random variable
$R_{i}=\frac{Y_{i}}{L_{i}}$, that is the proportion of expenditure reimbursed
(Indicated Deductible Relativity IDR).
Letting $h$ index the number of the $H$ risk classes and $\boldmath{D}_{h}$ be
the $g_{h}$-elements set of policyholders belonging to $h$-th risk class, we
can state $\textbf{x}_{i\in\boldmath{D}_{h}}=\textbf{x}_{h},\forall i\in
D_{h}$. Hence, in the following,
$E(R|\textbf{x}_{i\in\boldmath{D}_{h}})=E(R|\textbf{x}_{h})=E(R_{h}),\forall
i\in\boldmath{D}_{h}$ is the mean conditioned to the $\textbf{x}_{h}$ vector
of risk factors.
It is important to note that the r.v. $R_{h}$ can continuously assume values
between zero and one, but $Prob\left(R_{h}=0\right)>0$ and
$Prob\left(R_{h}=1\right)>0$. It means that the density function of $R_{h}$ is
a mixture distribution, because $R_{h}$ is derived from a collection of two
random variables: the first one representing the event $(R_{h}=0\cup R_{h}=1)$
and the second defining the level of proportion conditioned to the event
$\left(0<R_{h}<1\right)$. This mixture belongs to the so called Zero-One
Inflated distributions as the density function of the IDR can continuously
assume values mass at 0 and 1, in addition to the probability density within
$\left(0,1\right)$. To model the density between $\left(0,1\right)$ we
consider a Beta distribution. Then the r.v. $R_{h}$ is assumed to be Zero-and-
One Inflated Beta distributed.

We introduce a mixture between a Bernoulli distribution and a Beta
distribution as proposed in (Ospina and Ferrari Ospina ). Specifically,
assuming that the cumulative distribution function (hereafter cdf) of the
generic random variable $R$ is:

$BEINF\left(r;p_{o},p_{1},a,b\right)=\left(p_{0}+p_{1}\right)\cdot
Ber\left(r;\frac{p_{1}}{p_{0}+p_{1}}\right)+\left(1-p_{0}+p_{1}\right)\cdot
Beta\left(r;a,b\right)$ (2)

where $Ber\left(\cdot;\frac{p_{1}}{p_{0}+p_{1}}\right)$ represents the cdf of
a Bernoulli random variable with parameter $\frac{p_{1}}{p_{0}+p_{1}}$ and
$Beta\left(r;a,b\right)$ is the Beta cdf, whose density function is:

$beta(r;a,b)=\frac{\Gamma(a+b)}{\Gamma(a)\cdot\Gamma(b)}\cdot
r^{a-1}\cdot(1-r)^{b-1}\ ,\;\;\;0<r<1$ (3)

We say that $R$ has a Zero-and-One Inflated Beta distribution (henceforth
BEINF) (i.e $R\sim BEINF\left(p_{o},p_{1},a,b\right)$) with parameters
$p_{0}$, $p_{1}$, $a$ and $b$, if its density function (hereafter
$beinf\left(r;p_{o},p_{1},a,b\right)$) with respect to the measure generated
by the mixture is given by:

$beinf\left(r;p_{o},p_{1},a,b\right)=\left\\{\begin{array}[]{cc}p_{0}&r=0\\\
p_{1}&r=1\\\ \left(1-p_{0}-p_{1}\right)\cdot
f(r;a,b)&\,\,\,\,\,\,otherwise\end{array}\right.$ (4)

It is worth noting that is
$p_{0}=Prob\left(R=0\right)$,$p_{1}=Prob\left(R=1\right)$.

In order to consider a BEINF distribution for each $R_{h}$, we need to specify
the dependence structure between the response variables and the set of
covariates $\textbf{x}_{h}$ by means of a regression model.

##  3 GAMLSS theory

Considering the features of the r.v.s previously introduced a Generalized
Linear Model for Location, Scale and Shape (GAMLSS) seems to be an appropriate
choice (see Rigby and Stasinopoulos Rigby  and Stas ) for our goal. A GAMLSS
model assumes independent observations $r_{i}$, $i=1,\dots,I$ from a random
variable $R_{i}$, with probability density function
$f\left(r_{i}|\boldmath{\theta^{i}}\right)$, conditional on
$\boldmath{\theta^{i}}=\left(\theta_{1,i},\theta_{2,i},\theta_{3,i},\theta_{4,i}\right)=\left(\mu_{i},\sigma_{i},\nu_{i},\tau_{i}\right)$
a vector of four distribution parameters, each of which can be a function to
the explanatory variables. We shall refer to
$\left(\mu_{i},\sigma_{i},\nu_{i},\tau_{i}\right)$ as the distribution
parameters. The first two population distribution parameters $\mu_{i}$ and
$\sigma_{i}$ are usually characterized as location and scale parameters, while
the remaining parameter(s), if any, are characterized as shape parameters,
e.g., skewness and kurtosis parameters, although the model may be applied more
generally to the parameters of any population distribution, and can be
generalized to more than four distribution parameters. Rigby and Stasinopoulos
define the original formulation of a GAMLSS model as follows. Let
$\boldmath{r}^{T}=\left(r_{1},\ldots,r_{I}\right)$ be the $I$ length vector of
the response variable. Also for $k=1,2,3,4$, let $g_{k}(.)$ be known monotonic
link functions relating the distribution parameters to explanatory variables
by:

$g_{k}\left(\boldmath{\theta}_{k}\right)=\eta_{k}=\boldmath{X}_{k}\beta_{k}+\sum_{j=1}^{J_{k}}Z_{jk}\gamma_{jk},\,\,\,\,\,with\,\,\,\,k=1,2,3,4.$
(5)

where $\boldmath{\mu}$, $\boldmath{\sigma}$, $\boldmath{\nu}$,
$\boldmath{\tau}$ and $\boldmath{\eta}_{k}$ are vectors of length $I$,
${\beta_{k}^{T}}=\left(\beta_{1k},\beta_{2k},\ldots,\beta_{m_{k}k}\right)$ is
a parameter vector of length $m_{k}$, $\boldmath{X}_{k}$ is a fixed known
design matrix of order $I\times m_{k}$, $Z_{jk}$ is a fixed known $I\times
q_{jk}$ design matrix and $\gamma_{jk}$ is a $q_{jx}$ dimensional random
variable which is assumed to be distributed as $\gamma_{jk}\sim
N_{q_{jk}}\left(\boldmath{0},\boldmath{G}_{jk}^{-1}\right)$, where
$\boldmath{G}_{jk}^{-1}$ is the (generalized) inverse of a $q_{jk}\times
q_{jk}^{\prime}$ symmetric matrix
$\boldmath{G}_{jk}=\boldmath{G}_{jk}\left(\boldmath{\lambda}_{jk}\right)$,
which may depend on a vector of hyperparameters $\boldmath{\lambda}_{jk}$, and
where if $\boldmath{G}_{jk}$ is singular then $\gamma_{jk}$ is understood to
have an improper prior density function proportional to
$exp\left(-\frac{1}{2}\gamma_{jk}^{T}\boldmath{G}_{jk}\gamma_{jk}\right)$.

The population probability (density) function
$f\left(y|\boldmath{\theta}\right)$ in model (5) is deliberately left general
with no explicit conditional distributional form for the response variable. As
in the previous section we have introduced a density function particularly
useful for our ratemaking purposes, the Zero-One Inflated Beta, then we can
express the following parametrization as proposed in Rpack :

$\boldmath{\mu}=\frac{a}{a+b};\;\;\;\boldmath{\sigma}=\frac{1}{a+b+1};\;\;\;\boldmath{\nu}=\frac{p_{0}}{1-p_{0}-p_{1}};\;\;\;\boldmath{\tau}=\frac{p_{1}}{1-p_{0}-p_{1}}$
(6)

##  4 Some numerical Results

We set our model on a database from an Italian health insurance company
between years 2017 and 2018, the company offers two kind of covers: the first
for surgery and the second for diagnostic. For this two products, we observe a
number of policyholder exposed $I_{s}=63,790$ and $I_{d}=58,994$ respectively.
We apply a GAMLSS on the observed IDR as response variable, where the
covariates are a class of deductibles with 3 levels (1 lowest).

A logit function is used for $\mu$ and $\sigma$, as these parameters must be
included between 0 and 1 (see 5). For $\nu$ and $\tau$ we choose a logarithmic
function, since the lowest value of the AIC is obtained. Random effect are not
considered in this study.
Hence, equations (7), (8), (9) and (10) reports the model assumed for the four
parameters

$\displaystyle logit\left(\boldmath{\mu}\right)$ $\displaystyle=$
$\displaystyle\beta_{1,0}+deductible_{2}\cdot\beta_{1,1}+deductible_{3}\cdot\beta_{1,2}$
(7) $\displaystyle logit\left(\boldmath{\sigma}\right)$ $\displaystyle=$
$\displaystyle\beta_{2,0}$ (8) $\displaystyle log\left(\boldmath{\nu}\right)$
$\displaystyle=$
$\displaystyle\beta_{3,0}+deductible_{2}\cdot\beta_{3,1}+deductible_{3}\cdot\beta_{3,2}$
(9) $\displaystyle log\left(\boldmath{\tau}\right)$ $\displaystyle=$
$\displaystyle\beta_{4,0}+deductible_{2}\cdot\beta_{4,1}+deductible_{3}\cdot\beta_{4,2}$
(10)

In Table 1 and Figure 1 the fitted values are reported and compared with the
observed values.

Table 1: Observed vs Fitted values. Branch | Deductible | Beta | $p_{0}$ | $p_{1}$
---|---|---|---|---
|  | Exposure | Observed | Fitted | Exposure | Observed | Fitted | Exposure | Observed | Fitted
Surgery | Level 1 | 366 | 88.72% | 88.08% | 47,845 | 97.88% | 97.88% | 47,845 | 1.35% | 1.35%
| Level 2 | 167 | 79.66% | 80.30% | 9,566 | 98.24% | 98.24% | 9,566 | 0.01% | 0.01%
| Level 3 | 867 | 76.90% | 76.37% | 6,379 | 86.41% | 86.41% | 6,379 | 0.00% | 0.00%
Diagnostic | Level 1 | 14,358 | 77.04% | 76.58% | 44,245 | 52.67% | 52.67% | 44,245 | 14.88% | 14.88%
| Level 2 | 8,649 | 70.08% | 69.47% | 8,889 | 2.70% | 2.70% | 8,889 | 0.00% | 0.01%
| Level 3 | 4,298 | 64.10% | 64.12% | 5,860 | 26.66% | 26.66% | 5,860 | 0.00% | 0.00%

t] 25pt20pt

(a) $p_{0}$: surgery

(b) $p_{0}$: diagnostic

(c) $p_{1}$: surgery

(d) $p_{1}$: diagnostic

(e) $\mu$: surgery

(f) $\mu$: diagnostic

Figure 1: Goodness of fit of GAMLSS with BEINF distribution.

As observable, the GAMLSS shows a very good quality of fitting as the
probability of a 0 and 1 are perfectly matched whereas the values in $(0,1)$
have a distance of less then 1% in each branch.

## References

  * (1) Brown, R.L., Lennox, W.S.: Ratemaking and Loss Reserving for Property and Casualty Insurance. ACTEX Publications (2015)
  * (2) Frees, E.W., Lee,G.Y., Yang,L.: Multivariate frequency-severity regression models in insurance. Risks 4(1), (2015)
  * (3) Frees, E.W., Lee,G.Y.: Rating endorsements using generalized linear models. Variance 10(1): 51–74 (2016)
  * (4) Lee,G.Y. General Insurance Deductible Ratemaking. North American Actuarial Journal 21(4), 620-638, (2017)
  * (5) Ospina, R., Ferrari, S.L.P.: Inflated beta distributions. Statistical Papers 51, 111 (2010)
  * (6) Rigby,R. A., Stasinopoulos, D. M. : Generalized additive models for location, scale and shape. Applied Statistics 54(3), 507-554 (2005)
  * (7) Stasinopoulos, D.M., Rigby, R.A., Akantziliotou, C.: Instructions on how to use the GAMLSS package in R. Accompanying documentation in the current GAMLSS help files, (see also http://www.gamlss.org/) (2006)
  * (8) Stasinopoulos, D.M., Rigby, R.A.: Generalized additive models for location scale and shape (GAMLSS) in R. Journal of Statistical Software 23(7), (2007)