# AUTOCALIBRATION AND TWEEDIE-DOMINANCE FOR INSURANCE PRICING WITH MACHINE
LEARNING

Michel Denuit
Institute of Statistics, Biostatistics and Actuarial Science
UCLouvain
Louvain-la-Neuve, Belgium
Arthur Charpentier
Université du Québec à Montréal (UQAM)
Montreal, Quebec, Canada
Julien Trufin
Department of Mathematics
Université Libre de Bruxelles (ULB)
Brussels, Belgium

###### Abstract

Boosting techniques and neural networks are particularly effective machine
learning methods for insurance pricing. Often in practice, the sum of fitted
values can depart from the observed totals to a large extent. The possible
lack of balance when models are trained by minimizing deviance outside the
familiar GLM with canonical link setting has been documented in Wüthrich
(2019, 2020, 2021). The present paper aims to further study this phenomenon
when learning proceeds by minimizing Tweedie deviance. It is shown that
minimizing deviance involves a trade-off between the integral of weighted
differences of lower partial moments and the bias measured on a specific
scale. Hence, there is no guarantee that the sum of fitted values stays close
to observed totals if the latter bias term is dominated by the former one
entering deviance. Autocalibration is then proposed as a remedy. This new
method to correct for bias adds an extra local GLM step to the analysis with
the output of the first step as only predictor. Theoretically, it is shown
that it implements the autocalibration concept in pure premium calculation and
ensures that balance also holds on a local scale, not only at portfolio level
as with existing bias-correction techniques.
Keywords: Risk classification, Method of marginal totals, Tweedie distribution
family, Convex order.

##  1 Introduction and motivation

In the 1960s, North-American actuaries pioneered risk classification with the
help of minimum bias methods, after Bailey and Simon (1960) and Bailey (1963).
The central idea is that an acceptable set of premiums should reproduce the
experience within sub-portfolios corresponding to each level of meaningful
risk factors (like gender or age, for instance) and also the overall
experience, i.e. be balanced for each level and in total (leading to the
method of marginal totals, or MMT in short).

In the late 1980s, it turned out that for any GLM with canonical link function
and a score containing an intercept, there is an exact balance between fitted
and observed aggregated responses over the whole data set and for any level of
the categorical features. This formally related GLMs to minimum bias and MMT,
as documented in Mildenhall (1999). This connection greatly facilitated the
wide acceptance of GLMs in actuarial practice that has been particularly fast
after the development of powerful computer tools. The objective function used
for model training (often call loss function in machine learning) generally
corresponds to deviance or log-likelihood, in a model accounting for the
nature of insurance data under consideration: typically, Poisson for claim
counts, Gamma for average claim severities and compound Poisson sums with
Gamma-distributed terms for claim totals, all belonging to the Tweedie family
with power variance function.

Actuarial risk classification remained bridged to statistical regression
models and naturally followed their evolution from GLMs to GAMs, trees and
random forests, gradient and statistical boosting, projection pursuit and
neural networks, to name just a few. This evolution took place gradually,
following the availability of computational resources in statistical software.
The deviance established itself as the only objective function, even when the
interpretability of the likelihood equations in terms of balance was lost, and
the underlying MMT caution to GLMs has been progressively forgotten.

However, it turned out that tree-based boosting models and Neural Networks
trained to minimize deviance often violate total balance, even on the training
data set. This has been documented by Wüthrich (2019, 2020, 2021). Advanced
learning models are indeed able to produce scores that better correlate with
the response, as well as with the true premium compared to classical GLMs.
This comes from the additional freedom obtained by letting scores to depend in
a flexible way of available features, not only linearly. But breaking the
overall balance is the price to pay for this higher correlation. Because no
constraint on the replication of the observed total, or global balance is
imposed, machine learning tools are also able to substantially increase
overall bias. As pointed out by Wüthrich (2021), machine learning tools may
provide accurate fit at individual policy level but the global price level may
be completely wrong. It is therefore crucial to correct candidate premiums
produced by flexible machine learning tools before they can be used in
practice.

A natural remedy consists in restoring global balance at each step of the
iterative procedure used to optimize the loss function. This is easily
implemented by revising the intercept (under canonical link function) or by
constraining the optimization so that the observed total matches its fitted
counterpart. This simple solution restores global balance but does not ensure
that financial equilibrium also holds in meaningful sub-portfolios (remember
that GLMs with canonical link not only imposes global balance, at portfolio
level when the score comprises an intercept, but also at the level of risk
classes determined by binary features).

For this reason, we propose a new strategy based on the concept of
autocalibration (see, e.g., Kruger and Ziegel, 2020). This approach guarantees
global balance as well as local equilibrium in the spirit of the original MMT.
This simple and effective solution to the problem is implemented by adding an
extra step implementing MMT within a local GLM analysis. Specifically, after
the analysis has been performed with a method that does not necessarily
respect marginal totals, a local constant GLM fit is achieved in order to
restore the connection with MMT. Thus, as advocated by Wüthrich (2019, 2020,
2021), we also combine GLMs with advanced statistical learning tools but in a
different way. Wüthrich (2019, 2020, 2021) used Neural Networks to produce new
features in the last hidden layer, to be used in a GLM replacing the output
layer. In this approach, Neural Networks allow actuaries to perform feature-
engineering to feed the GLM score and marginal totals are respected by the use
of GLMs in the last step. Wüthrich (2019, 2020) then explains how to interpret
the working features generated by the last hidden layer of Neural Networks.
This approach can also be related to the polishing procedure proposed by Zumel
(2019) where random forests predictions are used to train a linear model in a
second step. As pointed out by this author, the extra step should not be
performed on the same data in order to avoid a potential source of overfitting
(called nested-model bias). The local GLM approach proposed in this paper
applies to any statistical learning model, not specifically to Neural Networks
or random forests. It consists in using the candidate premium produced in the
first step as unique feature in the second autocalibration step. This ensures
that the feature space reduces to the real line. By using a local constant, or
intercept-only GLM, the output of the first step defines optimal neighborhoods
to perform local averaging of observed losses.

The remainder of the paper is structured as follows. Section 2 provides the
reader with formal definitions and notation used throughout the text. Section
3 recalls Tweedie model and associated deviance. In Section 4, a mixture
representation of Tweedie deviance is derived. Precisely, Tweedie deviance is
decomposed into the sum of the integral of weighted differences of lower
partial moments and the bias measured on a specific scale. This decomposition
is used to understand the consequences of training the model by minimizing
Tweedie deviance. Special attention is devoted to the Poisson regression case.
The concept of Tweedie dominance is also introduced there, as a natural
counterpart to Bregman dominance, or forecast dominance discussed in Kruger
and Ziegel (2020). In Section 5, we propose a simple and powerful method to
restore balance at both global and local levels, based on the concept of
autocalibration. Tweedie dominance between autocalibrated predictors reduces
to the well-known convex order, or stop-loss order with equal means that has
been proposed to compare predictors by Denuit et al. (2019) and Kruger and
Ziegel (2020). A numerical study is provided in Section 6. The final Section 7
discusses the results and concludes.

##  2 Context, definition and notation

Let us now describe the notation used in this paper. We consider a response
$Y$ and a set of features $X_{1},\ldots,X_{p}$ gathered in the vector
$\boldsymbol{X}\in\mathcal{X}$ (classically,
$\mathcal{X}\subset\mathbb{R}^{p}$). In this paper, the response is typically
the number of claims reported to the insurance company by a given
policyholder, the average claim severity or the total claim amount in relation
with this contract. The dependence structure inside the random vector
$(Y,X_{1},\ldots,X_{p})$ is exploited to extract the information contained in
$\boldsymbol{X}$ about $Y$. In actuarial pricing, the aim is to evaluate the
pure premium as accurately as possible. This means that the target is the
conditional expectation $\mu(\boldsymbol{X})=\mathrm{E}[Y|\boldsymbol{X}]$ of
the response $Y$ (claim number or claim amount) given the available
information $\boldsymbol{X}$. Henceforth, $\mu(\boldsymbol{X})$ is referred to
as the true (pure) premium. Notice that in some applications,
$\mu(\boldsymbol{X})$ only refers to one component of the pure premium. For
instance, working in the frequency-severity decomposition of insurance losses,
$\mu(\boldsymbol{X})$ can be either the expected number of insured events or
the expected claim size, or severity.

The function
$\boldsymbol{x}\mapsto\mu(\boldsymbol{x})=\mathrm{E}[Y|\boldsymbol{X}=\boldsymbol{x}]$
is unknown to the actuary, and may exhibit a complex behavior in
$\boldsymbol{x}$. This is why this function is approximated by a (working, or
actual) premium $\boldsymbol{x}\mapsto\pi(\boldsymbol{x})$ with a simpler
structure. When the analyst is working in the frequency-severity decomposition
of insurance losses, $\pi(\boldsymbol{x})$ targets the expected number of
insured events or the mean claim severity, separately. Once fitted on the
training data set using an appropriate learning procedure, this produces
estimates $\widehat{\pi}(\boldsymbol{x})$ for $\mu(\boldsymbol{x})$, or fitted
values.

The developments in this paper apply in any setting where a global balance is
desirable, that is, where it is important that the sum of estimates does not
deviate too much from the sum of actual observations at both the entire
portfolio level and also more locally, in meaningful classes of policyholders.
The reason is obvious: the sum of the pure premiums must match the claim total
as accurately as possible so that the insurance company is able to indemnify
all third-parties and beneficiaries in execution of the contracts, without
excess nor deficit, by the very definition of pure premium (expense loadings
and cost-of-capital charges are added into the calculation at a later stage,
when moving to commercial premiums). This naturally translates into a global
balance constraint: considering that the total claim figures are
representative of next-year’s experience, it is important that the sum of
fitted premiums $\widehat{\pi}(\boldsymbol{x})$ match the sum of responses $Y$
taken as proxy for the total premium income (i.e., the sum of
$\mu(\boldsymbol{x})$), as closely as possible. But local equilibrium is also
essential to guarantee a competitive pricing.

The merits of a given pricing tool can be assessed using the pair
$\big{(}\mu(\boldsymbol{X}),\widehat{\pi}(\boldsymbol{X})\big{)}$ so that we
are back to the bivariate case even if there were thousands of features
comprised in $\boldsymbol{X}$. What really matters is the correlation between
$\widehat{\pi}(\boldsymbol{X})$ and $\mu(\boldsymbol{X})$ but as
$\mu(\boldsymbol{X})$ is unobserved the actuary can only use its noisy version
$Y$ to reveal the agreement of the true premium $\mu(\boldsymbol{X})$ with its
working counterpart $\widehat{\pi}(\boldsymbol{X})$. In insurance
applications, $\widehat{\pi}(\boldsymbol{X})$ is supposed to be used as a
premium so that correlation is important, but it is also essential that the
sum of predictions $\widehat{\pi}(\boldsymbol{X})$ matches the sum of actual
losses as closely as possible, as explained before. This is expressed by the
global balance condition and its local version.

To ease the exposition, we assume that predictor
$\widehat{\pi}(\boldsymbol{X})$ under consideration, as well as the
conditional expectation $\mu(\boldsymbol{X})$ are continuous random variables
admitting probability density functions. This is generally the case when there
is at least one continuous feature contained in the available information
$\boldsymbol{X}$ and the function $\widehat{\pi}$ is a continuously increasing
function of a real score built from $\boldsymbol{X}$. However, this rules out
predictions based on discrete features only, as well as piecewise constant
predictors, e.g., a single tree. Indeed, then $\widehat{\pi}(\boldsymbol{X})$
takes only a limited number of values. As actuarial pricing is nowadays based
on more sophisticated models (trees being combined into random forests, for
instance), this continuity assumption does not really restrict the generality
of the approach.

##  3 Tweedie model and deviance

In this paper, we assume that the response obeys a probability distribution
belonging to the Tweedie subclass of the Exponential Dispersion family.
Precisely, this means that we assume that the logarithm of the probability
mass function for a discrete response, or of the probability density function
for a continuous response, is of the form $\ln f=(y\theta-a(\theta))/\phi$ up
to a constant term, for some known dispersion parameter $\phi$ (that may
include a weight) and non-decreasing and convex cumulant function $a(\cdot)$
with $a^{\prime}(\theta)=\mathrm{E}[Y]=\mu$. The variance is then given by

$\mathrm{Var}[Y]=\phi a^{\prime\prime}(\theta)=\phi V(\mu)$

where the variance function $V(\cdot)$ corresponds to the second derivative of
the cumulant function whose argument $\theta$ has been replaced in terms of
$\mu$, that is, $\theta=(a^{\prime})^{-1}(\mu)$.

The Tweedie subclass corresponds to variance functions of the form
$V(\mu)=\mu^{\xi}$ for some power parameter $\xi$. Table 3.1 lists all Tweedie
distributions. Negative values of $\xi$ give continuous distributions on the
whole real axis. For $0<\xi<1$, there is no Exponential Dispersion
distribution with such variance function. Only the cases $\xi\geq 1$ are thus
interesting for applications in insurance. In the remainder of this paper, we
thus restrict our analysis to $\xi\geq 1$.

For $\xi\in(1,2)$, Tweedie distributions correspond to compound Poisson-Gamma
distributions, that is, compound Poisson sums with Gamma-distributed summands.
Starting from

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}\mu}\ln f$ $\displaystyle=$
$\displaystyle\frac{y-\mu}{\phi V(\mu)},$

we get with $V(\mu)=\mu^{\xi}$ that the probability density function over
$(0,\infty)$ is given by

$\displaystyle\ln f(y)$ $\displaystyle=$ $\displaystyle\int\frac{y-m}{\phi
m^{\xi}}\mathrm{d}m$ $\displaystyle=$
$\displaystyle\frac{1}{\phi}\left(y\frac{\mu^{1-\xi}}{1-\xi}-\frac{\mu^{2-\xi}}{2-\xi}\right)+\text{constant}$

with probability mass at zero
$\exp\left(-\frac{\mu^{2-\xi}}{\phi(2-\xi)}\right)$. Such compound Poisson-
Gamma distributions can be used for modeling annual claim amounts, having
positive probability at zero and a continuous distribution on the positive
real numbers. This offers an alternative to the decomposition of total losses
into claim numbers and claim severities, using Poisson distribution for
modeling claim counts and Gamma distribution for claim severities. We refer
the reader to Delong et al. (2021) for a through presentation of Tweedie
models.

| Type | Name
---|---|---
$\xi<0$ | Continuous | -
$\xi=0$ | Continuous | Normal
$0<\xi<1$ | Non existing | -
$\xi=1$ | Discrete | Poisson
$1<\xi<2$ | Mixed, non-negative | Compound Poisson sum
|  | with Gamma-distributed terms
$\xi=2$ | Continuous, positive | Gamma
$2<\xi<3$ | Continuous, positive | -
$\xi=3$ | Continuous, positive | Inverse Gaussian
$\xi>3$ | Continuous, positive | -
Table 3.1: Tweedie distributions and corresponding power parameters.

Let $\widehat{\pi}$ be the estimated mean response $Y$ built from some
training set (all formulas in this paper are meant given this training set).
The respective performances of competing models can then be assessed on the
basis of a validation set $\\{(Y_{i},\boldsymbol{X}_{i}),\hskip
5.69054pti=1,2,\ldots,n\\}$, that has not been used to obtain $\widehat{\pi}$.
Performances of $\widehat{\pi}$ are generally assessed with the help of out-
of-(training) sample deviance, also called predictive deviance and given by

$D_{n}(\xi,\widehat{\pi})=\frac{1}{n}\sum_{i=1}^{n}L\big{(}Y_{i},\widehat{\pi}(\boldsymbol{X}_{i})\big{)}$

where $L(\cdot,\cdot)$ is the loss function adopted to train the model. If $n$
is large enough then we can resort to the limiting value

$D_{n}(\xi,\widehat{\pi})\to
D(\xi,\widehat{\pi})=\mathrm{E}\big{[}L\big{(}Y^{\text{new}},\widehat{\pi}(\boldsymbol{X}^{\text{new}})\big{)}\big{]}\text{
as }n\to\infty,$

where $(Y^{\text{new}},\boldsymbol{X}^{\text{new}})$ is a new observation,
independent of, and distributed as those $(Y_{i},\boldsymbol{X}_{i})$
contained in the training set. In this paper, we compare models on the basis
of the large-sample version of the predictive deviance. This approach is
meaningful in insurance applications where the analyst is typically in a data-
rich situation.

Henceforth, we compare models on the basis of the predictive Tweedie deviance
that reduces to

$D(\xi,\widehat{\pi})=\left\\{\begin{array}[]{l}\mathrm{E}\left[\widehat{\pi}(\boldsymbol{X}^{\text{new}})-Y^{\text{new}}\ln\widehat{\pi}(\boldsymbol{X}^{\text{new}})\right]\text{
for $\xi=1$}\\\ \\\
\mathrm{E}\left[\ln\widehat{\pi}(\boldsymbol{X}^{\text{new}})+\frac{Y^{\text{new}}}{\widehat{\pi}(\boldsymbol{X}^{\text{new}})}\right]\text{
for $\xi=2$}\\\ \\\
\mathrm{E}\left[\frac{\widehat{\pi}(\boldsymbol{X}^{\text{new}})^{2-\xi}}{2-\xi}-\frac{Y^{\text{new}}\widehat{\pi}(\boldsymbol{X}^{\text{new}})^{1-\xi}}{1-\xi}\right]\text{
for $\xi>1$ and $\xi\neq 2$}\end{array}\right..$ (3.1)

In the remainder of the paper, we use (3.1) to assess the performances of a
given predictor $\widehat{\pi}$.

##  4 Tweedie dominance

Tweedie dominance is defined as dominance for every Tweedie deviances (3.1).
Precisely, $\widehat{\pi}_{2}$ outperforms $\widehat{\pi}_{1}$ in terms of
Tweedie dominance if the inequality $D(\xi,\widehat{\pi}_{2})\leq
D(\xi,\widehat{\pi}_{1})$ holds true for every power parameter $\xi\geq 1$.
Tweedie dominance appears to be a particular case of Bregman dominance, also
called forecast dominance defined as dominance for every Bregman loss
function. We refer the interested reader to Kruger and Ziegel (2020) and the
references therein for an extensive presentation of this concept. Tweedie
dominance is thus a particular stochastic order relation used to compare the
performances of two estimators $\widehat{\pi}_{1}$ and $\widehat{\pi}_{2}$ for
the conditional means. We refer the reader to Shaked and Shanthikumar (2007)
for a general presentation of stochastic order relations and to Denuit et al.
(2005) for applications to insurance.

The next result provides the actuary with a sufficient condition for a model
to outperform a competitor in terms of Tweedie dominance.

######  Proposition 4.1.

Define

$\psi_{\xi}(\pi)=\left\\{\begin{array}[]{l}\ln\pi\text{ for $\xi=2$}\\\ \\\
\frac{\pi^{2-\xi}}{2-\xi}\text{ else}.\end{array}\right.$ (4.1)

Then, $\widehat{\pi}_{2}$ outperforms $\widehat{\pi}_{1}$ in terms of Tweedie
dominance if

$\mathrm{E}\left[\psi_{\xi}(\widehat{\pi}_{1}(\boldsymbol{X}^{\mathrm{new}}))\right]\geq\mathrm{E}\left[\psi_{\xi}(\widehat{\pi}_{2}(\boldsymbol{X}^{\mathrm{new}}))\right]\text{
for all }\xi\geq 1$ (4.2)

and

$\mathrm{E}\big{[}Y^{\mathrm{new}}\mathrm{I}\left[\widehat{\pi}_{1}(\boldsymbol{X}^{\mathrm{new}})\leq
t\right]\big{]}\geq\mathrm{E}\big{[}Y^{\mathrm{new}}\mathrm{I}\left[\widehat{\pi}_{2}(\boldsymbol{X}^{\mathrm{new}})\leq
t\right]\big{]}\quad\text{for all }t\geq 0.$ (4.3)

###### Proof.

For $\xi=1$, $\widehat{\pi}_{2}$ is superior to $\widehat{\pi}_{1}$ if

$\displaystyle\mathrm{E}\left[\widehat{\pi}_{1}(\boldsymbol{X}^{\text{new}})\right]-\mathrm{E}\left[Y^{\text{new}}\ln\widehat{\pi}_{1}(\boldsymbol{X}^{\text{new}})\right]\geq\mathrm{E}\left[\widehat{\pi}_{2}(\boldsymbol{X}^{\text{new}})\right]-\mathrm{E}\left[Y^{\text{new}}\ln\widehat{\pi}_{2}(\boldsymbol{X}^{\text{new}})\right]$
$\displaystyle\Leftrightarrow$
$\displaystyle\mathrm{E}\left[\widehat{\pi}_{1}(\boldsymbol{X}^{\text{new}})\right]-\mathrm{E}\left[\widehat{\pi}_{2}(\boldsymbol{X}^{\text{new}})\right]+\mathrm{E}\left[Y^{\text{new}}\ln\widehat{\pi}_{2}(\boldsymbol{X}^{\text{new}})\right]-\mathrm{E}\left[Y^{\text{new}}\ln\widehat{\pi}_{1}(\boldsymbol{X}^{\text{new}})\right]\geq
0.$

Since the identity

$\displaystyle\mathrm{E}\left[Y^{\text{new}}\ln\widehat{\pi}(\boldsymbol{X}^{\text{new}})\right]$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\mathrm{E}\big{[}Y^{\text{new}}\mathrm{I}\left[\ln\widehat{\pi}(\boldsymbol{X}^{\text{new}})>t\right]\big{]}\mathrm{d}t-\int_{-\infty}^{0}\mathrm{E}\big{[}Y^{\text{new}}\mathrm{I}\left[\ln\widehat{\pi}(\boldsymbol{X}^{\text{new}})\leq
t\right]\big{]}\mathrm{d}t$ $\displaystyle=$
$\displaystyle\int_{1}^{\infty}\mathrm{E}\big{[}Y^{\text{new}}\mathrm{I}\left[\widehat{\pi}(\boldsymbol{X}^{\text{new}})>s\right]\big{]}\frac{1}{s}\mathrm{d}s-\int_{0}^{1}\mathrm{E}\big{[}Y^{\text{new}}\mathrm{I}\left[\widehat{\pi}(\boldsymbol{X}^{\text{new}})\leq
s\right]\big{]}\frac{1}{s}\mathrm{d}s$

holds true for any predictor $\widehat{\pi}$, we have

$\displaystyle\mathrm{E}\left[Y^{\text{new}}\ln\widehat{\pi}_{2}(\boldsymbol{X}^{\text{new}})\right]-\mathrm{E}\left[Y^{\text{new}}\ln\widehat{\pi}_{1}(\boldsymbol{X}^{\text{new}})\right]$
$\displaystyle=$
$\displaystyle\int_{1}^{\infty}\big{(}\mathrm{E}\big{[}Y^{\text{new}}\mathrm{I}\left[\widehat{\pi}_{2}(\boldsymbol{X}^{\text{new}})>s\right]\big{]}-\mathrm{E}\big{[}Y^{\text{new}}\mathrm{I}\left[\widehat{\pi}_{1}(\boldsymbol{X}^{\text{new}})>s\right]\big{]}\big{)}\frac{1}{s}\mathrm{d}s$
$\displaystyle-\int_{0}^{1}\big{(}\mathrm{E}\big{[}Y^{\text{new}}\mathrm{I}\left[\widehat{\pi}_{2}(\boldsymbol{X}^{\text{new}})\leq
s\right]\big{]}-\mathrm{E}\big{[}Y^{\text{new}}\mathrm{I}\left[\widehat{\pi}_{1}(\boldsymbol{X}^{\text{new}})\leq
s\right]\big{]}\big{)}\frac{1}{s}\mathrm{d}s$ $\displaystyle=$
$\displaystyle\int_{1}^{\infty}\big{(}\mathrm{E}\big{[}Y^{\text{new}}\left(1-\mathrm{I}\left[\widehat{\pi}_{2}(\boldsymbol{X}^{\text{new}})\leq
s\right]\right)\big{]}-\mathrm{E}\big{[}Y^{\text{new}}\left(1-\mathrm{I}\left[\widehat{\pi}_{1}(\boldsymbol{X}^{\text{new}})\leq
s\right]\right)\big{]}\big{)}\frac{1}{s}\mathrm{d}s$
$\displaystyle+\int_{0}^{1}\big{(}\mathrm{E}\big{[}Y^{\text{new}}\mathrm{I}\left[\widehat{\pi}_{1}(\boldsymbol{X}^{\text{new}})\leq
s\right]\big{]}-\mathrm{E}\big{[}Y^{\text{new}}\mathrm{I}\left[\widehat{\pi}_{2}(\boldsymbol{X}^{\text{new}})\leq
s\right]\big{]}\big{)}\frac{1}{s}\mathrm{d}s$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\big{(}\mathrm{E}\big{[}Y^{\text{new}}\mathrm{I}\left[\widehat{\pi}_{1}(\boldsymbol{X}^{\text{new}})\leq
s\right]\big{]}-\mathrm{E}\big{[}Y^{\text{new}}\mathrm{I}\left[\widehat{\pi}_{2}(\boldsymbol{X}^{\text{new}})\leq
s\right]\big{]}\big{)}\frac{1}{s}\mathrm{d}s$ $\displaystyle\geq$
$\displaystyle 0$

by (4.3). Hence, both conditions (4.2) and (4.3) ensure inequality (4).

Turning to the case $\xi=2$, $\widehat{\pi}_{2}$ is superior to
$\widehat{\pi}_{1}$ if

$\displaystyle\mathrm{E}\left[\ln\widehat{\pi}_{1}(\boldsymbol{X}^{\text{new}})\right]+\mathrm{E}\left[\frac{Y^{\text{new}}}{\widehat{\pi}_{1}(\boldsymbol{X}^{\text{new}})}\right]\geq\mathrm{E}\left[\ln\widehat{\pi}_{2}(\boldsymbol{X}^{\text{new}})\right]+\mathrm{E}\left[\frac{Y^{\text{new}}}{\widehat{\pi}_{2}(\boldsymbol{X}^{\text{new}})}\right]$
$\displaystyle\Leftrightarrow$
$\displaystyle\mathrm{E}\left[\ln\widehat{\pi}_{1}(\boldsymbol{X}^{\text{new}})\right]-\mathrm{E}\left[\ln\widehat{\pi}_{2}(\boldsymbol{X}^{\text{new}})\right]+\mathrm{E}\left[\frac{Y^{\text{new}}}{\widehat{\pi}_{1}(\boldsymbol{X}^{\text{new}})}\right]-\mathrm{E}\left[\frac{Y^{\text{new}}}{\widehat{\pi}_{2}(\boldsymbol{X}^{\text{new}})}\right]\geq
0.$

Hence, it suffices to notice that (4.3) implies
$\mathrm{E}\left[\frac{Y^{\text{new}}}{\widehat{\pi}_{1}(\boldsymbol{X}^{\text{new}})}\right]-\mathrm{E}\left[\frac{Y^{\text{new}}}{\widehat{\pi}_{2}(\boldsymbol{X}^{\text{new}})}\right]\geq
0$ since the identity

$\displaystyle\mathrm{E}\left[\frac{Y^{\text{new}}}{\widehat{\pi}(\boldsymbol{X}^{\text{new}})}\right]$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\mathrm{E}\left[Y^{\text{new}}\mathrm{I}\left[\frac{1}{\widehat{\pi}(\boldsymbol{X}^{\text{new}})}\geq
t\right]\right]\mathrm{d}t$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\mathrm{E}\left[Y^{\text{new}}\mathrm{I}\left[\widehat{\pi}(\boldsymbol{X}^{\text{new}})\leq\frac{1}{t}\right]\right]\mathrm{d}t$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\mathrm{E}\left[Y^{\text{new}}\mathrm{I}\left[\widehat{\pi}(\boldsymbol{X}^{\text{new}})\leq
s\right]\right]\frac{1}{s^{2}}\mathrm{d}s.$

is valid for every predictor $\widehat{\pi}$.

Finally, in the remaining cases, i.e. $\xi>1$ and $\xi\neq 2$,
$\widehat{\pi}_{2}$ is superior to $\widehat{\pi}_{1}$ if

$\displaystyle\mathrm{E}\left[\frac{\widehat{\pi}_{1}(\boldsymbol{X}^{\text{new}})^{2-\xi}}{2-\xi}\right]+\frac{1}{\xi-1}\mathrm{E}\left[\frac{Y^{\text{new}}}{\widehat{\pi}_{1}(\boldsymbol{X}^{\text{new}})^{\xi-1}}\right]\geq\mathrm{E}\left[\frac{\widehat{\pi}_{2}(\boldsymbol{X}^{\text{new}})^{2-\xi}}{2-\xi}\right]+\frac{1}{\xi-1}\mathrm{E}\left[\frac{Y^{\text{new}}}{\widehat{\pi}_{2}(\boldsymbol{X}^{\text{new}})^{\xi-1}}\right]$
$\displaystyle\Leftrightarrow$
$\displaystyle\mathrm{E}\left[\frac{\widehat{\pi}_{1}(\boldsymbol{X}^{\text{new}})^{2-\xi}}{2-\xi}\right]-\mathrm{E}\left[\frac{\widehat{\pi}_{2}(\boldsymbol{X}^{\text{new}})^{2-\xi}}{2-\xi}\right]+\frac{1}{\xi-1}\left(\mathrm{E}\left[\frac{Y^{\text{new}}}{\widehat{\pi}_{1}(\boldsymbol{X}^{\text{new}})^{\xi-1}}\right]-\mathrm{E}\left[\frac{Y^{\text{new}}}{\widehat{\pi}_{2}(\boldsymbol{X}^{\text{new}})^{\xi-1}}\right]\right)\geq
0.$

Whatever the predictor $\widehat{\pi}$, we can write

$\displaystyle\mathrm{E}\left[\frac{Y^{\text{new}}}{\widehat{\pi}(\boldsymbol{X}^{\text{new}})^{\xi-1}}\right]$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\mathrm{E}\left[Y^{\text{new}}\mathrm{I}\left[\widehat{\pi}(\boldsymbol{X}^{\text{new}})^{1-\xi}\geq
t\right]\right]\mathrm{d}t$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\mathrm{E}\left[Y^{\text{new}}\mathrm{I}\left[\widehat{\pi}(\boldsymbol{X}^{\text{new}})\leq\frac{1}{t^{\frac{1}{\xi-1}}}\right]\right]\mathrm{d}t$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\mathrm{E}\big{[}Y^{\text{new}}\mathrm{I}\left[\widehat{\pi}(\boldsymbol{X}^{\text{new}})\leq
s\right]\big{]}\frac{\xi-1}{s^{\xi}}\mathrm{d}s,$

The announced result then follows from

$\displaystyle\mathrm{E}\left[\frac{Y^{\text{new}}}{\widehat{\pi}_{1}(\boldsymbol{X}^{\text{new}})^{\xi-1}}\right]-\mathrm{E}\left[\frac{Y^{\text{new}}}{\widehat{\pi}_{2}(\boldsymbol{X}^{\text{new}})^{\xi-1}}\right]$
$\displaystyle=\int_{0}^{\infty}\Big{(}\mathrm{E}\big{[}Y^{\text{new}}\mathrm{I}\left[\widehat{\pi}_{1}(\boldsymbol{X}^{\text{new}})\leq
s\right]\big{]}-\mathrm{E}\big{[}Y^{\text{new}}\mathrm{I}\left[\widehat{\pi}_{2}(\boldsymbol{X}^{\text{new}})\leq
s\right]\big{]}\Big{)}\frac{\xi-1}{s^{\xi}}\mathrm{d}s.$

This ends the proof. ∎

Instead of Tweedie dominance, the actuary could select a specific parameter
$\xi$, only. This is typically the case with $\xi=1$ (Poisson regression for
counts), $\xi=2$ or $\xi=3$ (Gamma or Inverse Gaussian regression for claim
severities). In this case, condition (4.2) in Proposition 4.1 is imposed only
for that specific valued of $\xi$.

The main ingredient of the proof of Proposition 4.1 is the integral of the
difference of functions in (4.3) weighted by

$\frac{\xi-1+\mathrm{I}[\xi=1]}{s^{\xi}}=\frac{\xi-1+\mathrm{I}[\xi=1]}{V(s)}$

where $V(\cdot)$ is the variance function associated with the response
distribution. Weights thus appear to be inversely proportional to the
variability. Under the Tweedie variance function, $V(s)$ is a power of $s$ so
that smaller weights are assigned to the differences at larger values of the
response.

Proposition 4.1 shows that improving a predictor $\widehat{\pi}$ can be
achieved by

  * (i)

increasing the overall bias measured on a modified scale induced by the
auxiliary function $\psi_{\xi}$ defined in (4.1). This may act against the
conservation of observed totals.

  * (ii)

increasing the dependence between $\widehat{\pi}$ and the response, so to
decrease

$\displaystyle\mathrm{E}\Big{[}Y^{\text{new}}\mathrm{I}\big{[}\widehat{\pi}(\boldsymbol{X}^{\text{new}})\leq
t\big{]}\Big{]}$ $\displaystyle=$
$\displaystyle\mathrm{Cov}\Big{[}Y^{\text{new}},\mathrm{I}\big{[}\widehat{\pi}(\boldsymbol{X}^{\text{new}})\leq
t\big{]}\Big{]}+\mathrm{E}\Big{[}Y^{\text{new}}\Big{]}\mathrm{P}\Big{[}\widehat{\pi}(\boldsymbol{X}^{\text{new}})\leq
t\Big{]}$ $\displaystyle=$
$\displaystyle\mathrm{E}\Big{[}Y^{\text{new}}\Big{|}\widehat{\pi}(\boldsymbol{X}^{\text{new}})\leq
t\Big{]}\mathrm{P}\Big{[}\widehat{\pi}(\boldsymbol{X}^{\text{new}})\leq
t\Big{]}.$

The latter lower partial moment essentially depends on the correlation
structure of the pair
$\big{(}Y^{\text{new}},\widehat{\pi}(\boldsymbol{X}^{\text{new}})\big{)}$.
Notice that the quantities appearing in the decomposition above are closely
related to expectation dependence as defined by Wright (1987).

In an insurance ratemaking context, the lower partial moment can be
interpreted as best-profile premium income

$\displaystyle\mathrm{E}\Big{[}Y^{\text{new}}\mathrm{I}\big{[}\widehat{\pi}(\boldsymbol{X}^{\text{new}})\leq
t\big{]}\Big{]}$ $\displaystyle=$
$\displaystyle\mathrm{E}\Big{[}\mu(\boldsymbol{X}^{\text{new}})\mathrm{I}\big{[}\widehat{\pi}(\boldsymbol{X}^{\text{new}})\leq
t\big{]}\Big{]}$ $\displaystyle\approx$
$\displaystyle\frac{1}{n}\sum_{i|\widehat{\pi}(\boldsymbol{X}_{i})\leq
t}\mu(\boldsymbol{X}_{i}).$

Here, $\sum_{i|\widehat{\pi}(\boldsymbol{X}_{i})\leq
t}\mu(\boldsymbol{X}_{i})$ is the true premium income for the sub-portfolio
formed by gathering all policyholders with predicted premium at most equal to
$t$. These policyholders exhibit the best risk profiles according to the
candidate premium $\widehat{\pi}$.

Of course, there is a trade-off between these two goals when computing
$D(\xi,\widehat{\pi})$. If the model is very flexible then it can produce a
predictor that correlates a lot with the response and the lower partial moment
can be decreased to a large extent compared to models imposing a rigid form
for the predictor (like GLMs). There is in fact a fundamental trade-off
between (4.2) and (4.3). Since $\psi_{\xi}$ is increasing, (4.2) favors small
values of $\widehat{\pi}_{2}$. By contrast, (4.3) requires large values of
$\widehat{\pi}_{2}$ so that the indicator function on the right-hand side is
0. There is thus no guarantee that balance holds on a particular data set: if
one component dominates, we may end up with candidate premiums that are either
too small or too large compared to observed responses.

######  Example 4.2 (Poisson regression).

Poisson deviance is by far the most widely used one in insurance applications.
It applies for instance to claim counts in property and casualty insurance,
death counts in life insurance, and numbers of transitions in health
insurance. We assume that we deal with a response $Y$ obeying the Poisson
distribution. A vector $\boldsymbol{X}$ of features is available to predict
the mean response $\mu(\boldsymbol{X})$. To ease the presentation, we assume
unit exposures.

Considering (4.1)-(4.2) with $\xi=1$, the bias is thus measured on the
response scale. Condition (4.2) then favors $\widehat{\pi}_{2}$ if

$\mathrm{E}[\mu(\boldsymbol{X}^{\text{new}})]-\mathrm{E}[\widehat{\pi}_{2}(\boldsymbol{X}^{\text{new}})]\geq\mathrm{E}[\mu(\boldsymbol{X}^{\text{new}})]-\mathrm{E}[\widehat{\pi}_{1}(\boldsymbol{X}^{\text{new}})]$

or, equivalently,

$\mathrm{E}[Y^{\text{new}}]-\mathrm{E}[\widehat{\pi}_{2}(\boldsymbol{X}^{\text{new}})]\geq\mathrm{E}[Y^{\text{new}}]-\mathrm{E}[\widehat{\pi}_{1}(\boldsymbol{X}^{\text{new}})].$

A larger bias may thus be advantageous, depending on the fundamental trade-off
between (4.2) and (4.3) explained above. As a consequence, adopting Poisson
deviance as loss function outside GLMs may create a total premium income gap

$\displaystyle\mathrm{E}\Big{[}Y^{\text{new}}\Big{]}-\mathrm{E}\Big{[}\widehat{\pi}(\boldsymbol{X}^{\text{new}})\Big{]}$
$\displaystyle=$
$\displaystyle\mathrm{E}\Big{[}\mu(\boldsymbol{X}^{\text{new}})\Big{]}-\mathrm{E}\Big{[}\widehat{\pi}(\boldsymbol{X}^{\text{new}})\Big{]}$
$\displaystyle\approx$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\mu(\boldsymbol{X}_{i})-\frac{1}{n}\sum_{i=1}^{n}\widehat{\pi}(\boldsymbol{X}_{i})$

where $\sum_{i=1}^{n}\mu(\boldsymbol{X}_{i})$ is the true premium income for
the validation set whereas $\sum_{i=1}^{n}\widehat{\pi}(\boldsymbol{X}_{i})$
is the one obtained by adopting predictor $\widehat{\pi}$ for premium
calculation. This gap may become quite large with highly flexible models such
as Neural Networks or boosting.

##  5 Restoring balance at global and local scales

###  5.1 Autocalibration

Recall that a predictor $\widehat{\pi}$ is said to be autocalibrated if
$\widehat{\pi}(\boldsymbol{X})=\mathrm{E}[Y|\widehat{\pi}(\boldsymbol{X})]$.
We refer the reader to Kruger and Ziegel (2020) for a general presentation of
this concept. By Jensen inequality, autocalibration thus ensures that

$\mathrm{E}[g(\widehat{\pi}(\boldsymbol{X}))]\leq\mathrm{E}[g(Y)]\text{ for
every convex function }g,$

or equivalently, that

$\mathrm{E}[\widehat{\pi}(\boldsymbol{X})]=\mathrm{E}[Y]\text{ and
}\mathrm{E}[(\widehat{\pi}(\boldsymbol{X})-t)_{+}]\leq\mathrm{E}[(Y-t)_{+}]\text{
for all }t\geq 0.$

These inequalities correspond to the convex order between
$\widehat{\pi}(\boldsymbol{X})$ and $Y$. Thus, autocalibration implies that
the predictor is less variable than the response, in the sense of the convex
order.

Under mild technical requirement, a simple way to restore global balance
consists in switching from $\widehat{\pi}$ to its balance-corrected version
$\widehat{\pi}_{\text{BC}}$ defined as

$\widehat{\pi}_{\text{BC}}(\boldsymbol{X})=\mathrm{E}[Y|\widehat{\pi}(\boldsymbol{X})]$

that averages to $\mathrm{E}[Y]$, as shown in the next result.

######  Property 5.1.

If $s\mapsto\mathrm{E}[Y|\widehat{\pi}(\boldsymbol{X})=s]$ is continuously
increasing then the balance-corrected version $\widehat{\pi}_{\mathrm{BC}}$ of
the candidate premium $\widehat{\pi}$ satisfies the autocalibration property.

###### Proof.

If $s\mapsto\mathrm{E}[Y|\widehat{\pi}(\boldsymbol{X})=s]$ is continuously
increasing, that is, if $\widehat{\pi}_{\text{BC}}(\boldsymbol{X})$ is
continuously increasing in $\widehat{\pi}(\boldsymbol{X})$, then Lemma 2.2 in
Shaked et al. (2012) allows us to write

$\mathrm{E}[Y|\widehat{\pi}(\boldsymbol{X})]=\mathrm{E}[Y|\widehat{\pi}_{\text{BC}}(\boldsymbol{X})]=\widehat{\pi}_{\text{BC}}(\boldsymbol{X})$

so that the resulting $\widehat{\pi}_{\text{BC}}(\boldsymbol{X})$ is indeed
autocalibrated. This ends the proof. ∎

In insurance applications, autocalibration induces local balance and imposes
financial equilibrium not only at portfolio level but also in any sufficiently
large sub-portfolio. This concept thus appears to be particularly appealing in
a ratemaking context.

###  5.2 Autocalibrating a given predictor

We can restore global balance, or unbias the predictor by reconciling the
predicted and observed total on the training set. In this way, we recover the
global balance property imposed after the seminal work by Bailey and Simon
(1960). However, this does not ensure that balance holds locally. Indeed,
global balance is only one of the GLM likelihood equations under canonical
link, corresponding to the intercept. In order to extend the other likelihood
equations imposed in the GLM setting to general machine learning procedures,
we need to mimic the way local GLM proceeds for fitting, by defining
meaningful neighborhoods for statistical learning.

To this end, we work under canonical link function with data points
$(Y_{i},e_{i},\widehat{\pi}(\boldsymbol{x}_{i}))$ for some relevant exposure
$e_{i}$. An intuitively acceptable solution would consist in imposing marginal
constrains on local neighborhoods defined by mean of $\widehat{\pi}$. This
allows for some local transfers of claims and premiums from neighboring
policyholders and so implements local balance conditions in sub-portfolios
corresponding to these neighborhoods. This is in essence the local GLM
approach (see Loader, 1999, for a detailed account) that allows the actuary to
maintain the relationship with MMT.

In order to obtain an autocalibrated version $\widehat{\pi}_{\text{BC}}$ of
$\widehat{\pi}$, let us consider a specific risk profile $\boldsymbol{x}$. A
weight $\nu_{i}(\widehat{\pi}(\boldsymbol{x}))$ is assigned to each
$(Y_{i},e_{i},\widehat{\pi}(\boldsymbol{x}_{i}))$, $i=1,\ldots,n$, computed
from some weight function $\nu(\cdot)$ chosen to be continuous, symmetric,
peaked at 0 and defined on $[-1,1]$. These weights depend on the relative
distance of $\widehat{\pi}(\boldsymbol{x}_{i})$ with respect to
$\widehat{\pi}(\boldsymbol{x})$. A common choice for $\nu(\cdot)$ is the
tricube weight function but several alternatives are available, including
rectangular (or uniform), Gaussian or Epanechnikov, for instance. Here,
$\nu_{i}(\widehat{\pi}(\boldsymbol{x}))$ is larger for policyholders $i$ such
that $\widehat{\pi}(\boldsymbol{x}_{i})$ is close to
$\widehat{\pi}(\boldsymbol{x})$ and decreases when
$\widehat{\pi}(\boldsymbol{x}_{i})$ gets far away from
$\widehat{\pi}(\boldsymbol{x})$.

The local GLM likelihood equation

$\sum_{i=1}^{n}\nu_{i}(\widehat{\pi}(\boldsymbol{x}))y_{i}=\sum_{i=1}^{n}\nu_{i}(\widehat{\pi}(\boldsymbol{x}))e_{i}\widehat{\pi}_{\text{BC}}(\boldsymbol{x})$

thus matches MMT constraints: smoothing is ensured by transferring part of the
experience at neighboring $\widehat{\pi}$ values to obtain
$\widehat{\pi}_{\text{BC}}$. A local constant, or intercept-only GLM thus
provides the appropriate fitting procedure when balance must be respected.
Opting for a rectangular weight function complies with MMT: the weights
$\nu_{i}$ are constant within the smoothing window and zero otherwise so that
the sums reduce to observations comprised within this window and the uniform
weights factor out.

Our main message here is thus that a local, intercept-only GLM with
rectangular weights implements local balance, or MMT in a second step, within
sub-portfolios gathering policyholders with about the same predicted value
according to the first step. The rectangular weight function involved in the
statistical procedure optimally transfers part of claim experience between
neighboring policyholders. This approach implements smoothness from a
statistical point of view while remaining fully transparent and
understandable. Indeed, local averaging can just be seen as an application of
the mutuality principle at the heart of insurance.

##  6 Numerical illustrations

In this section, we first consider the freMTPL2freq dataset, from the
CASDataset R package of Charpentier (2014). The variable of interest is annual
claim frequency. Precisely, we run Poisson regression on response ClaimNb,
with exposure Exposure and various explanatory features (continuous and
categorical). Then we use the freMTPL2sev dataset, with the severity, and use
a Tweedie model to get an estimation of total annual losses.

###  6.1 Claim frequency

Four models are considered : a generalized linear model
$\widehat{\pi}^{\text{glm}}$, a generalized additive model
$\widehat{\pi}^{\text{gam}}$ where continuous features are transformed
nonlinearly using spline functions, a boosting algorithm
$\widehat{\pi}^{\text{bst}}$ and a Neural Network model
$\widehat{\pi}^{\text{nn}}$. For the boosting models, we use the h2o package,
with h2o.gbm. For the Neural Network approach, we use combined actuarial
neural net (CANN111see https://www.kaggle.com/floser/glm-neural-nets-and-
xgboost-for-insurance-pricing for a discussion and the codes used here.), from
Schelldorfer and Wüthrich (2019), using Keras. Notice that in the latter case,
even with the same seed, outputs of the Neural Nets procedure change.

In order to model our counting variable, we set distribution = "poisson" and
offset_column = "Exposure". For the boosting inference procedure, we use
ntrees = 30 (the impact of the number of iterations will be discussed later
on) and nfolds = 5222The codes used in this section can be found on the github
repository
https://github.com/freakonometrics/autocalibration. From the initial dataset,
with 678,013 rows, 70% are randomly chosen for our training dataset (474,609
rows) and used to construct $\widehat{\pi}$, and the remaining 30% (203,404
rows) are used as a validation dataset, to draw the following graphs.

On Figure 6.1, we can visualise the distribution of predictions
$\widehat{\pi}(\boldsymbol{x}_{i})$. The average value $\overline{\pi}$ for
the four models, on the training dataset, is given on the left of Table 6.1
(additional information is also given there, namely quantiles). As a
benchmark, if we run a simple Poisson regression of the intercept only, we
obtain $\widehat{\beta}_{0}=-2.2911$, corresponding to a baseline prediction
$\overline{\pi}=0.10115$. Observe that the boosting procedure globally
underestimates claims frequency, while neural nets globally overestimate.

Figure 6.1: Histogram of
$\\{\widehat{\pi}(\boldsymbol{x}_{1}),\cdots,\widehat{\pi}(\boldsymbol{x}_{n})\\}$
on the validation dataset. In this section,  corresponds to the standard
Poisson regression (GLM),  corresponds to the smooth additive regression
(GAM),  is the boosting model and  is the Neural Network.

| $\widehat{\pi}^{\text{glm}}$ | $\widehat{\pi}^{\text{gam}}$ | $\widehat{\pi}^{\text{bst}}$ | $\widehat{\pi}^{\text{nn}}$
---|---|---|---|---
average $\overline{\pi}$ | 0.1091 | 0.1091 | 0.0821 | 0.1230
10% quantile | 0.0602 | 0.0592 | 0.0494 | 0.0529
90% quantile | 0.1688 | 0.1729 | 0.1264 | 0.2051

| $\widehat{\pi}_{BC}^{\text{glm}}$ | $\widehat{\pi}_{BC}^{\text{gam}}$ | $\widehat{\pi}_{BC}^{\text{bst}}$ | $\widehat{\pi}_{BC}^{\text{nn}}$
---|---|---|---|---
average $\overline{\pi}$ | 0.1051 | 0.1059 | 0.1028 | 0.1051
10% quantile | 0.0573 | 0.0570 | 0.0518 | 0.0499
90% quantile | 0.1687 | 0.1806 | 0.1711 | 0.1776

Table 6.1: Summary statistics on
$\\{\widehat{\pi}(\boldsymbol{x}_{1}),\cdots,\widehat{\pi}(\boldsymbol{x}_{n})\\}$,
on the validation dataset (assuming an exposure of $1$ to provide annualized
predictions), for $\pi$ on the left, and the corrected version $\pi_{BC}$ on
the right.

Figure 6.2 displays $s\mapsto\mathrm{E}[Y|\widehat{\pi}(\boldsymbol{X})=s]$
when $s\in[0,0.2]$. For the Generalized Linear Model
$\mathrm{E}[Y|\widehat{\pi}(\boldsymbol{X})=s]\sim s$, while
$\mathrm{E}[Y|\widehat{\pi}(\boldsymbol{X})=s]>s$ for the boosting model. The
positive bias we observe on $\widehat{\pi}^{\text{bst}}$ means that this model
underestimates the true price of the risk almost everywhere; similarly, the
negative bias we observe on $\widehat{\pi}^{\text{nn}}$ almost eveywhere means
that this model overestimates the true price. As shown on Table 6.1, the bias
correction step solves this issue.

Figure 6.2: Evolution of
$s\mapsto\mathrm{E}[Y|\widehat{\pi}(\boldsymbol{X})=s]$ (the thin straight
line corresponds to $\mu=\widehat{\pi}$).

On Figure 6.3, we can compare $\widehat{\pi}$ and $\widehat{\pi}_{\text{BC}}$:
on top, we have the QQ plot of $\widehat{\pi}_{\text{BC}}$ against
$\widehat{\pi}$ and on the bottom, a scatterplot of
$\\{\widehat{\pi}(\boldsymbol{x}_{i}),\widehat{\pi}_{\text{BC}}(\boldsymbol{x}_{i})\\}$
(for a subset of the validation database). The distribution of
$\widehat{\pi}_{\text{BC}}$ can be visualized on Figure 6.4.

Figure 6.3: QQ plot of $\widehat{\pi}_{\text{BC}}$ against $\widehat{\pi}$
(plain line), on top, and scatterplot a subset of points
$\\{\widehat{\pi}(\boldsymbol{x}_{i}),\widehat{\pi}_{\text{BC}}(\boldsymbol{x}_{i})\\}$
from the validation dataset, below. Figure 6.4: Histogram of
$\\{\widehat{\pi}_{\text{BC}}(\boldsymbol{x}_{1}),\cdots,\widehat{\pi}_{\text{BC}}(\boldsymbol{x}_{n})\\}$
on the validation dataset.

The specification of bandwidth $h(s)$ is discussed in Loader (1999, Section
2.2.1). For a nearest neighbor bandwidth, compute all distances $|s-s_{i}|$
between the fitting point $s$ and the data points $s_{i}$, we choose $h(x)$ to
be the $k$th smallest distance, where $k=[n\alpha]$. In the R function locfit,
when alpha is given as a single number, it represents a nearest neighbor
fraction (the default smoothing parameter is $\alpha=70\%$). But a second
component can be added, $\alpha=(\alpha_{0},\alpha_{1})$. That second
component represents a constant bandwidth, and $h(s)$ will be computed as
follows: as previously, $k=[n\alpha_{0}]$, and if $d_{(i)}$ represents the
ordered statistics of $d_{i}=|s-s_{i}|$, $h(s)=\max\\{d_{(k)},\alpha_{1}\\}$.
The default value in R is alpha=c(0.7,0).

As we can see on Table 6.1, the boosting algorithm was globally
underestimating the average premium. Using a local regression can actually
lower the bias. On Figure 6.5, we used 60% of the original dataset to train
our three models (training dataset), then the local regression was performed
on 20% of the dataset (smoothing dataset) and finally various quantities (bias
and empirical loss) were computed on the remaining 20% (validation dataset).
On the left, we can visualize the evolution of the bias
$\displaystyle{\frac{1}{n}\sum_{i=1}^{n}(\widehat{\pi}^{\text{bst}}_{\text{BC}}(\boldsymbol{x}_{i})-y_{i})}$
as a function of $\alpha$. Here, smoothing has no impact on the overall bias
of GLM and GAM (here the two models have a (small) positive bias on the
validation dataset). Smoothing has an impact on the correction of
$\widehat{\pi}^{\text{bst}}$. We can observe that a small $\alpha$ can lead to
a much smaller bias (possibly null). Note that we used $\alpha=5\%$ for
previous graphs of this section. In the middle of Figure 6.5, we can visualize
the evolution of the empirical Poisson loss
$\displaystyle{\frac{1}{n}\sum_{i=1}^{n}(\widehat{\pi}^{\text{bst}}_{\text{BC}}(\boldsymbol{x}_{i})-y_{i}\ln\widehat{\pi}^{\text{bst}}_{\text{BC}}(\boldsymbol{x}_{i}))}$.
Note that again, a small $\alpha$ leads to a smaller loss. It is also
interesting to see that similar performances are obtained in case of
overfitting. Thus, autocalibration also corrects for overfitting by locally
averaging the initial predictors.

Figure 6.5: Evolution of the bias as a function a $\alpha$, from $2.5\%$ to
$97.5\%$, on the left, of the Poisson loss in the center, and joint
scatterplot of bias and empirical loss on the right. Horizontal light lines
correspond to original estimates $\widehat{\pi}$ while strong curves are
$\widehat{\pi}_{\text{BC}}$ for various smoothing parameter $\alpha$. As
previously  corresponds to the standard Poisson regression (GLM),  corresponds
to the smooth additive regression (GAM), while  is the boosting model with 30
trees (as previously) and  is the boosting model with 1000 trees
(overfitting).

As we have seen when moving from $\widehat{\pi}$ to
$\widehat{\pi}_{\text{BC}}$ the distribution of premiums changed substantially
(see Figures 6.1 and 6.4) and the transformation was not (perfectly) monotonic
(see Figure 6.3). Nevertheless, the rank correlation between $\widehat{\pi}$
and $\widehat{\pi}_{\text{BC}}$ is rather high (0.99 for all models), as
reported in Figure 6.6.

Figure 6.6: Spearman’s rank correlation between $\widehat{\pi}^{\text{glm}}$,
$\widehat{\pi}^{\text{gam}}$, $\widehat{\pi}^{\text{bst}}$,
$\widehat{\pi}^{\text{nn}}$ and their autocalibrated counterparts,
$\widehat{\pi}_{\text{BC}}^{\text{glm}}$,
$\widehat{\pi}_{\text{BC}}^{\text{gam}}$,
$\widehat{\pi}_{\text{BC}}^{\text{bst}}$ and
$\widehat{\pi}_{\text{BC}}^{\text{nn}}$, on the validation dataset

Finally, on Figure 6.7 we can visualize the concentration curves of the three
models, against $\widehat{\pi}^{\text{gam}}_{\text{BC}}$. We refer the reader
to Denuit et al. (2019) for a thorough presentation of this diagnostic tool.
Concentration curves clearly demonstrate here the improved performances of the
autocalibrated version of the predictor.

Figure 6.7: Concentration curves, with $\widehat{\pi}^{\text{glm}}$,
$\widehat{\pi}^{\text{gam}}$, $\widehat{\pi}^{\text{bst}}$, against
$\widehat{\pi}_{\text{BC}}^{\text{bst}}$, from the left to the right on top
and $\widehat{\pi}^{\text{glm}}-\widehat{\pi}_{\text{BC}}^{\text{bst}}$,
$\widehat{\pi}^{\text{gam}}-\widehat{\pi}_{\text{BC}}^{\text{bst}}$ and
$\widehat{\pi}^{\text{bst}}-\widehat{\pi}_{\text{BC}}^{\text{bst}}$ from the
left to the right, on the bottom.

###  6.2 Claims totals

Let us now consider aggregate losses and Tweedie deviance with $\xi\in(1,2)$.
Three models are considered: a generalized linear model
$\widehat{\pi}^{\text{glm}}$, a generalized additive model
$\widehat{\pi}^{\text{gam}}$ where continuous features are transformed
nonlinearly using spline functions, and a boosting algorithm
$\widehat{\pi}^{\text{bst}}$. For the boosting models, we used the TDboost
package. We did remove all observations with an exposure not close to 1, to
avoid major corrections when the exposure was too small (and additional
uncertainty).

On Figure 6.8, we can visualise the distribution of predictions
$\widehat{\pi}(\boldsymbol{x}_{i})$, with the three models. A grid search lead
to $\widehat{\xi}=1.6$. On Figure 6.9, we can visualize the evolution of
$s\mapsto\mathrm{E}[Y|\widehat{\pi}(\boldsymbol{X})=s]$ when $s$ varies
between 0 and 300. Observe that 300 is about 4 times the average value on the
training dataset ($\overline{y}=84.9$).

Figure 6.8: Histogram of
$\\{\widehat{\pi}(\boldsymbol{x}_{1}),\cdots,\widehat{\pi}(\boldsymbol{x}_{n})\\}$
on the validation dataset. In this section,  corresponds to the standard
Poisson regression (GLM),  corresponds to the smooth additive regression
(GAM), and  is the boosting model, with a Tweedie loss function. Figure 6.9:
Evolution of $s\mapsto\mathrm{E}[Y|\widehat{\pi}(\boldsymbol{X})=s]$ (the thin
straight line corresponds to $\mu=\widehat{\pi}$).

On Figure 6.10, we illustrate the impact of an underestimation of the power
index, with respectively 1.5 and 1.4 (instead of 1.6) for the GLM and GAM
procedures. We can see there that the results remain stable when the value of
$\xi$ changes. On Figure 6.11 we exclude policies with an annual loss
exceeding 10,000. Note that those represent 0.0539% of the training dataset
(10,000 is the 99% quantile of losses, for policies experiencing a loss).

Figure 6.10: Evolution of
$s\mapsto\mathrm{E}[Y|\widehat{\pi}(\boldsymbol{X})=s]$, when the Tweedie
power is 1.5 on top, and 1.4 below.  Figure 6.11: Evolution of
$s\mapsto\mathrm{E}[Y|\widehat{\pi}(\boldsymbol{X})=s]$, when online policies
with claims below 10,000 are kept in the training dataset.

##  7 Discussion

The main message of this paper is as follows. Advanced learning models are
able to produce scores that better correlate with the response, as well as
with the true premium compared to classical GLMs. This comes from the
additional freedom obtained by letting scores to depend in a flexible way of
available features, not only linearly. But breaking the overall balance is the
price to pay for this higher correlation. Because no constraint on the
replication of the observed total, or global balance is imposed, machine
learning tools are also able to substantially increase overall bias.

To prevent this to occur, the balance-corrected version of any predictor can
be obtained by local GLM, recognizing the nature of the response $Y$: local
Poisson GLM for claim counts, local Gamma GLM for average claim severities and
compound Poisson sums with Gamma-distributed terms for claim totals. The
canonical link function is adopted so that maximum-likelihood estimates
replicate observed totals. An intercept-only GLM is fitted locally, with
rectangular weight function, on a reduced set of observations consisting in
observed responses and candidate premiums, with proximity assessed with the
help of the candidate premium to be corrected.

In this approach, the supervised learning model is used to produce a real-
valued signal $\widehat{\pi}(\boldsymbol{X})$, reducing the high-dimensional
feature space to the real line. In fact, the score of the model is enough to
compute $\widehat{\pi}_{\text{BC}}$ so that the learning model is essentially
used to assess proximity among individuals, before performing local averaging.
In that respect, the proposed approach shares some similarities with $k$-NN,
or $k$ nearest-neighbors except that here, proximity is assessed with the help
of the real-valued $\widehat{\pi}$ produced by the learning model. And the
proposed extra autocalibration step also counteracts possible overfitting by
locally averaging the initial predictor.

The approach proposed in this paper reconciles minimum bias and method of
marginal totals, at the origin of insurance risk classification, with modern
learning tools. With minimum bias, the amount of premium is computed in order
to compensate insurance risks within meaningful sub-portfolios, maintaining an
overall balance. The proposed balance correction mechanism implements the very
same idea, using machine learning tools to define meaningful neighborhoods
where collected premiums must match claims to be compensated. This is in
accordance with the fundamental mutuality principle at the heart of insurance.
Autocalibration applies very generally, to any machine learning model and only
requires a statistical smoother like locfit.

## Acknowledgements

The authors thank two anonymous Referees and the Editor for their constructive
comments which greatly helped to improve this paper.

## References

  1. -

Bailey, R.A. (1963). Insurance rates with minimum bias. Proceedings of the
Casualty Actuarial Society 50, 4-11.

  2. -

Bailey, R.A., Simon, L.J. (1960). Two studies on automobile insurance
ratemaking. ASTIN Bulletin 1, 192-217.

  3. -

Charpentier, A. (2014). Computational Actuarial Science with R. Chapman and
Hall / CRC.

  4. -

Delong, L., Lindholm, M., Wüthrich, M. V. (2021). Making Tweedie’s compound
Poisson model more accessible. European Actuarial Journal, in press.

  5. -

Denuit, M., Dhaene, J., Goovaerts, M.J., Kaas, R. (2005). Actuarial Theory for
Dependent Risks: Measures, Orders and Models. Wiley, New York.

  6. -

Denuit, M., Sznajder, D., Trufin, J. (2019). Model selection based on Lorenz
and concentration curves, Gini indices and convex order. Insurance:
Mathematics and Economics 89, 128-139.

  7. -

Kruger, F., Ziegel, J.F. (2020). Generic conditions for forecast dominance.
Journal of Business & Economic Statistics, in press.

  8. -

Loader, C. (1999). Local Regression and Likelihood. Springer, New York.

  9. -

Mildenhall, S.J. (1999). A systematic relationship between minimum bias and
generalized linear models. Proceedings of the Casualty Actuarial Society 86,
393-487.

  10. -

Savage, L.J. (1971). Elicitation of personal probabilities and expectations.
Journal of the American Statistical Association 66, 783-810.

  11. -

Shaked, M., Shanthikumar, J.G. (2007). Stochastic Orders. Springer, New York.

  12. -

Shaked, M., Sordo, M. A., Suarez-Llorens, A. (2012). Global dependence
stochastic orders. Methodology and Computing in Applied Probability 14,
617-648.

  13. -

Schelldorfer, J., Wüthrich, M.V. (2019). Nesting classical actuarial models
into Neural Networks. Available at SSRN: https://ssrn.com/abstract=3320525

  14. -

Wright, R. (1987). Expectation dependence of random variables, with an
application in portfolio theory. Theory and Decision 22, 111-124.

  15. -

Wüthrich, M.V. (2019). From Generalized Linear Models to Neural Networks, and
back. Available at SSRN: https://ssrn.com/abstract=3491790

  16. -

Wüthrich, M.V. (2020). Bias regularization in neural network models for
general insurance pricing. European Actuarial Journal 10, 179-202.

  17. -

Wüthrich, M.V. (2021). The balance property in neural network modelling.
Statistical Theory and Related Fields, in press.

  18. -

Zumel, N. (2019). An ad-hoc method for calibrating uncalibrated models. Win-
Vector Blog, WordPress.


# Quantum algorithm for Wang-Landau sampling

Garrett T. Floyd, David P. Landau, Michael R. Geller  Department of Physics
and Astronomy, University of Georgia, Athens, Georgia 30602, USA and Center
for Simulational Physics, University of Georgia, Athens, Georgia 30602, USA

(August 28, 2024)

###### Abstract

It has been shown that the Metropolis algorithm can be implemented on quantum
computers in a way that avoids the sign problem. However, flat histogram
techniques are often preferred as they don’t suffer from the same limitations
that afflict Metropolis for problems of real-world interest and provide a host
of other benefits. In particular, the Wang-Landau method is known for its
efficiency and accuracy. In this work we design, implement, and validate a
quantum algorithm for Wang-Landau sampling, greatly expanding the range of
quantum many body problems solvable by Monte Carlo simulation.

##  I Introduction

Monte Carlo techniques are now ubiquitous in physics and related disciplines
because of their ability to calculate finite-temperature equilibrium
properties [1]. The simplest and most well-known of these techniques is the
“traditional”, Markov chain Metropolis algorithm. Nonetheless, traditional
Monte Carlo algorithms have found limited applicability to quantum systems due
to the infamous sign problem. On the other hand, quantum computers promise to
efficiently simulate a broad range of quantum models and provide vastly
increased speed and accuracy compared to classical algorithms.

It has been shown by Temme et al. [2] that the Metropolis algorithm can be
adapted for quantum computers while avoiding the sign problem. Since then
interest in the quantum simulation of thermal states has intensified. This has
given rise to a number of advances in quantum computer based Markov chain
methods [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14] as well as new techniques
such as those based on variational methods [15, 16, 17, 18, 19, 20, 21, 22,
23], and others [24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37]. As
an example, while no rigorous complexity bounds are known for the classical or
quantum Metropolis algorithms, many of the works based on Markov chain methods
have shown that an additional quadratic speedup of hitting times can be
achieved. This is often associated with an optimal choice of quantum walk
operator [13, 5, 8, 10, 29, 12, 14, 35, 36, 38, 39, 40].

Yet despite the popularity of the Metropolis and other simple Markov chain
based algorithms, these approaches suffer from considerable defects. It is
well known that Metropolis algorithms exhibit critical slowing down at second
order phase transitions and can become trapped at low temperatures when
sampling complex free energy landscapes. Metropolis algorithms also require
separate runs for each temperature investigated. Furthermore, not all
thermodynamic quantities of interest, e.g. entropy and free energy, can be
calculated directly from the raw data provided from Metropolis. These
limitations often necessitate the use of advanced post-Metropolis simulation
methods.

One such technique is the Wang-Landau Monte Carlo algorithm [41, 42]. The
Wang-Landau algorithm is known to sample complex free energy landscapes while
maintaining applicability and simplicity. The Wang-Landau algorithm works by
performing a random walk in energy space in a manner which depends on the
current estimate of the density of states (DOS), which is updated at each step
of the walk. As the walk transition probability discourages visiting energies
with higher DOS estimates, it is better at exploring the entire phase space
than the Metropolis method. Due to the fact that the algorithm samples a
temperature independent DOS, properties of the system can be estimated
simultaneously at all temperatures, avoiding trapping by free energy barriers.

Implementing a Wang-Landau algorithm on a quantum computer is challenging,
however. Unlike the Metropolis algorithm whose random walk depends solely on
the energies of the current and proposed wave functions, the Wang-Landau
algorithm is dependent on the current estimate of the DOS entries for the
current and proposed energies. It therefore requires storing, accessing, and
updating the DOS, $g(E)$. While the excellent performance of the Wang-Landau
algorithm has been convincingly demonstrated, at least on classical models,
rigorous bounds on the convergence and accuracy are not available.

To demonstrate the feasibility and advantages of Wang-Landau sampling with a
quantum computer, a proof of concept is critical. Given the considerably
increased number of qubits needed to incorporate the DOS into the Monte Carlo
algorithm, the extreme circuit depth required of Monte Carlo algorithms, and
the exponential scaling cost of simulating additional qubits on a classical
computer, providing a viable proof of concept requires some simplification.
Here we design and simulate a quantum Wang-Landau algorithm that overcomes the
sign problem and allows increased performance compared to Metropolis.

To begin any Wang-Landau algorithm, a Hamiltonian and state space is
specified. Here we denote the classical state of the system by a vector ${\bf
X}$ in ${\mathbb{R}}^{n}$ or a subset of ${\mathbb{R}}^{n}$, for a (typically
large) integer $n$ proportional to the number of degrees of freedom. For a
system of $N$ Ising spins, ${\bf X}\in\\{-1,+1\\}^{N}$. Next, a fixed set of
move operators $\\{C_{m}\\}_{m=1}^{M}$ is chosen, with $M$ independent of $n$.
The move operators are randomly chosen and applied to ${\bf X}$. The energy
range over which the system is simulated is then partitioned into bins which
are labeled by an index $i$. An initial guess for the (logarithm of the)
relative DOS $\ln(g(E))$ is then made, and a histogram is initialized with
$h(E_{i})=0\ \forall$ i. Next a modification factor $f>1$ and reduction factor
$0<\gamma<1$ are chosen. A criterion determining the acceptable flatness of
the histogram (flatness criterion) needed to move between iterations and the
total number of iterations (max number of iterations) to be performed in the
simulation is also defined.

For reference, the classical Wang-Landau algorithm is briefly summarized in
Algorithm 1. The output of any Wang-Landau algorithm is a $g(E)$ which
approximates the (absolute or relative) DOS of the model. Thermodynamic
quantities such as the internal energy and entropy are then obtained directly
from the $g(E)$.

Data: Energy function,

Set of move operators to be randomly selected $\\{C_{m}\\}_{m=1}^{M}$,

Number of energy bins $\ell$,

Initial guess for the logarithm of the relative DOS $\ln(g(E))$,

Initial histogram $h(E_{i})=0\ \forall$ i,

Modification factor $f>1$,

Modification reduction factor $0<\gamma<1$,

Histogram flatness criterion (flatness criterion),

Total number of rounds (max number of rounds).

Set round number = 0

while  _round number < max number of rounds_ do

for  _ a predetermined number of steps_ do

apply a randomly chosen move $C_{m}$

${\bf X}_{i}\xrightarrow{C_{m}}{\bf X}_{j}$

draw a random number $\in[0,1]$

if  _random number < min(1,$e^{\ln(g(E_{i}))-\ln(g(E_{j}))}$)_ then

Set $\ln(g(E_{j}))=\ln(g(E_{j}))+\ln(f)$

Set $h(E_{j})=h(E_{j})+1$

Set ${\bf X}_{i}={\bf X}_{j}$

else

Set $\ln(g(E_{i}))=\ln(g(E_{i}))+\ln(f)$

Set $h(E_{j})=h(E_{j})+1$

Set ${\bf X}_{i}={\bf X}_{i}$

end if

end for

Perform histogram check (Algorithm 2)

end while

Algorithm 1 Classical Wang-Landau Algorithm

Calculate average histogram entry ${\bar{h}}$

${\bar{h}}=\sum_{i=1}^{\ell}\frac{h(E_{i})}{\ell}$

if  _$\exists$ i : $h(E_{i})<$ (flatness criterion)${\bar{h}}$ or $h(E_{i})>$
(2-flatness criterion)${\bar{h}}$_ then

pass

else

$\forall$ i Set $h(E_{i})=0$

Set $f=f^{\gamma}$

Set round number = round number + 1

end if

Algorithm 2 Check Histogram Subalgorithm

With the classical algorithm as a starting point, we then take inspiration
from the method employed by Temme [2]. As a review, we note that in Ref. [2]
the authors perform Metropolis Monte Carlo in the following way: They begin by
allocating a quantum circuit with 4 registers. The first register is used to
contain the wavefunction of the system itself, the second and third are used
to contain the pre- and post-update energy estimations, respectively, and the
last register is a single qubit used to encode the acceptance/rejection
probability.

Starting from an initial eigenstate and energy, a random unitary move update
is chosen from a universal gate set $\\{C_{m}\\}_{m=1}^{M}$ that is
constructed to fulfill detailed balance. In general this will lead to the
state register being in a superposition of eigenstates. Phase estimation is
then used to get an estimate of the energy, which is used to compute and
encode the acceptance probability in the last qubit for each branch of the
superposition. The last qubit is then measured to see if the move is accepted
or not.

If accepted, the new energy is measured from register 3, collapsing the
wavefunction to a new eigenstate. The registers are then prepared for the next
step. If the move is rejected, a procedure is undertaken to return the system
wavefunction to an eigenstate which has the same energy as the pre-update
eigenstate. Following this prescription, physical quantities can be estimated
from the list of energies generated. However, from an efficiency standpoint
the rejection procedure is often lengthy and computationally expensive
compared to the rest of the algorithm. It has also faced criticism for the way
it uses inverse phase estimation [43].

To understand this, suppose that a step of the rejection procedure should use
inverse phase estimation to return the state of the third register to the all
zero state, for example. Note that the qubit used to store the acceptance-
rejection probabilities was entangled with the third register directly before
its measurement. This causes the probability amplitudes of each branch of the
superposition to be altered by the measurement outcome of the final qubit.
This, in turn, causes the probability of the inverse phase estimation
returning the third register to the all zero state to be less than one in the
case of finite-precision phase estimation. If this occurs it causes all
further phase estimations in the algorithm to fail, corrupting the algorithm.
This may possibly be alleviated in the limit of a large number of phase
estimation qubits, but it has also recently been argued that coherent
uncomputation of phase estimation after measurement cannot be generally
guaranteed to succeed [43]. Here we employ an alternative that both remedies
this issue and greatly improves the speed.

Consider the state prepared in step 2 of Yung and Aspuru-Guzik [13], which is
the infinite temperature thermal state of the Hamiltonian. Namely, we use the
result that

$\frac{1}{\sqrt{N}}\sum_{x=0}^{N-1}\ket{x}\ket{x}=\frac{1}{\sqrt{N}}\sum_{i=0}^{N-1}\ket{\psi_{i}}\ket{\bar{\psi_{i}}},$
(1)

where

$\ket{\bar{\psi_{i}}}=\sum_{x=0}^{N-1}\innerproduct{\psi_{i}}{x}\ket{x},$ (2)

with $\ket{x}$ a computational basis states and the $\ket{\psi_{i}}$ are
eigenstates of the simulated Hamiltonian. If phase estimation is performed on
the first register and measured, the effect is that there is an equal
probability of selecting any eigenstate of the Hamiltonian. Using preparation
of this state as a wavefunction update subroutine, we note that for any
initial eigenstate, a move to any other eigenstate is proposed with the same
probability. There is then no need to re-prepare the initial system
wavefunction in the case of rejection as the energy and all compatible
observables have already been recorded and the choice of move update has no
dependence on the initial wave function. Since this global move is able to
select any eigenstate and does so with equal probability, it also satisfies
detailed balance requirements. The need for uncomputation and wave function
re-initialization can then be removed at the cost of doubling the number of
qubits necessary for the system wavefunction.

It should also be noted that this approach allows the acceptance-rejection
step to be treated classically, and one of energy registers (and acceptance-
rejection qubit) to be dropped. It also has the the added benefit that
generalizing this approach to Wang-Landau then allows the DOS and histogram to
be treated classically as well. Finally, as the approach outlined above uses
phase estimation to estimate eigenstate energies in the same way as in Ref.
[2], it inherets its ability to avoid the sign problem. These considerations
lead to the quantum Wang-Landau sampling algorithm summarized in Algorithm 3.

Data: Hamiltonian,

Number of energy bins $\ell$,

Initial guess for the logarithm of the relative DOS $\ln(g(E))$,

Initial histogram $h(E_{i})=0\ \forall$ i,

Modification factor $f>1$,

Reduction factor $0<\gamma<1$,

Histogram flatness criterion,

Total number of rounds (max number of rounds)

Set round number = 0

Initialize the quantum computer to the state

$\ket{\psi}\ket{\bar{\psi_{i}}}\ket{E_{i}}$

Record $E_{i}$, wipe all registers.

while  _round number < max number of rounds_ do

for  _ A predetermined number of steps_ do

Prepare the state

$\frac{1}{\sqrt{N}}\sum_{x=0}^{N-1}\ket{x}\ket{x}=\frac{1}{\sqrt{N}}\sum_{i=0}^{N-1}\ket{\psi_{i}}\ket{\bar{\psi_{i}}}$

Perform phase estimation and measure $E_{j}$

$\ket{\psi_{i}}\ket{\bar{\psi_{i}}}\ket{E_{j}}$

Draw a random number

if  _random number < min(1,$e^{\ln(g(E_{i}))-\ln(g(E_{j}))}$)_ then

Set $\ln(g(E_{j}))=\ln(g(E_{j}))+\ln(f)$

Set $h(E_{j})=h(E_{j})+1$

Set $E_{i}=E_{j}$

Wipe all registers

else

Set $\ln(g(E_{i}))=\ln(g(E_{i}))+\ln(f)$

Set $h(E_{i})=h(E_{i})+1$

Set $E_{i}=E_{i}$

wipe all registers

end if

end for

Perform histogram check (Algorithm 2)

end while

Algorithm 3 Quantum Wang-Landau Algorithm Figure 1: Energy predicted by Wang-
Landau, Metropolis, and exact diagonalization vs. inverse temperature $\beta$.
Errors are shown in the inset. Figure 2: Specific heat predicted by Wang-
Landau, Metropolis, and exact diagonalization vs. inverse temperature $\beta$.
Figure 3: Free energy data predicted by Wang-Landau, Metropolis, and exact
diagonalization vs. inverse temperature $\beta$. Figure 4: Entropy predicted
by Wang-Landau, Metropolis, and direct Hamiltonian diagonalization vs. inverse
temperature $\beta$.

To assess the effectiveness of the technique, and its benefits with respect to
Metropolis sampling, we perform proof of concept simulations using a 1D
transverse-field Ising model with periodic boundary conditions. In this model,
a chain of quantum spins has the Hamiltonian

$H=J\sum_{i=1}^{N}\sigma^{z}_{i}\sigma^{z}_{i+1}+h\sum_{i=1}^{N}\sigma^{x}_{i}.$
(3)

Here the summation is over the spins in the chain and because of the periodic
boundary conditions the last term in the sum is
$\sigma^{z}_{N}\sigma^{z}_{1}$. Furthermore, $J$ and $h$ are constants which
can be changed to evaluate model behavior for different relative interaction
strengths. In this work we choose $J=2$ and $h=1$. The details for how the
Monte Carlo sampling is carried out are given in the Supplemental Material.

The internal energy, heat capacity, free energy, and Gibbs entropy for 9
system qubits (spins) and 11 phase estimation qubits obtained from simulated
quantum Wang-Landau and Metropolis sampling are given in Figures 1-4. The
differences between the values of each thermodynamic quantity predicted by
both Monte Carlo methods and the exact value determined by Hamiltonian
diagonalization are given in the insets. The error bars in each plot represent
one standard deviation. The time evolution interval for phase estimation was
chosen in accordance with the procedure given in the Supplemental Material.

The Wang-Landau data are obtained from the average of twenty runs, and
quantities are calculated directly from the relative DOS. The Metropolis
energy and heat capacity for each temperature is calculated from the average
of twenty runs with a constant number of Monte Carlo steps between runs and
temperatures. The total number of Monte Carlo steps per run per temperature
was chosen so that the combined number of Monte Carlo steps of all Metropolis
runs would be equal to or slightly greater than that of the Wang-Landau runs
after adjusting for an initial relaxation period. The first 5,000 steps of
each Metropolis run were cut from data analysis to allow the system to relax
to equilibrium. In the case of the data shown in Figures 1-4, 12,080,000
Metropolis steps were taken for each run, and all energies after 5000 steps
were used to calculate thermodynamic quantities.

The energy and heat capacity given by Metropolis are calculated from taking
the averages of the moments of energy. The entropy

$S=\int_{0}^{T_{f}}\frac{C_{v}}{T}\,dT=\int_{\beta_{f}}^{\infty}\frac{C_{v}}{\beta}\,d\beta$
(4)

is calculated by evaluating the second expression in Eq. (4). To evaluate that
integral it is necessary to truncate the upper limit; this is done by assuming
that the heat capacity drops to zero one inverse temperature interval after
the $\beta$ cutoff. The free energy $F=U-TS$ is calculated for each
temperature from the internal energy and entropy obtained above.

We find that there is a very small systematic error between the data given by
our simulations and the exact answer given by direct diagonalization of the
Hamiltonian across all observed quantities. However, excluding systematic
error due to integration of Metropolis data, the observed systematic error can
primarily be attributed to error in the phase estimation algorithm and can be
seen to decrease with an increasing number of phase estimation qubits.
Integration error in the data has the effect of systematically displacing the
Metropolis entropy curve compared to that given by Wang-Landau and results in
the erroneous sharp uptick in the Metropolis free energy curve at low inverse
temperature compared to Wang-Landau and that given by Hamiltonian
diagonalization. An analysis of the individual Wang-Landau sampling and phase
estimation errors is given in the Supplemental Material.

In conclusion, we have introduced and simulated an algorithm for performing
Wang-Landau sampling with a fault-tolerant quantum computer. This algorithm
represents an improvement over Metropolis sampling, avoids the sign problem,
and doesn’t become trapped in complex free energy landscapes or at low
temperature. For each of the thermodynamic quantities studied, Wang-Landau
sampling is found to significantly outperform Metropolis sampling in accuracy
when using the same total number of Monte Carlo steps. The improvement is
especially prominent in the entropy and free energy since these cannot be
directly obtained from Metropolis data and must be evaluated indirectly,
increasing sensitivity to noise. This algorithm brings the tools for
simulating quantum systems to the modern era of post-Metropolis Monte Carlo
sampling, enabling the simulation of complex and strongly correlated quantum
matter with fault-tolerant quantum computers.

## References

  * David P. Landau and Kurt Binder [2015] David P. Landau and Kurt Binder,  _A Guide to Monte Carlo Simulations in Statistical Physics_ , 4th ed. (Cambridge University Press, 2015).
  * Temme  _et al._ [2011] K. Temme, T. J. Osborne, K. G. Vollbrecht, D. Poulin, and F. Verstraete, Nature 471, 87 (2011).
  * Nelson  _et al._ [2022] J. Nelson, M. Vuffray, A. Y. Lokhov, T. Albash, and C. Coffrin, High-quality Thermal Gibbs Sampling with Quantum Annealing Hardware (2022), arXiv:2109.01690 [quant-ph] .
  * Somma  _et al._ [2008] R. D. Somma, S. Boixo, H. Barnum, and E. Knill, Phys. Rev. Lett. 101, 130504 (2008).
  * Lemieux  _et al._ [2020] J. Lemieux, B. Heim, D. Poulin, K. Svore, and M. Troyer, Quantum 4, 287 (2020), arXiv:1910.01659 [quant-ph] .
  * Moussa [2022] J. E. Moussa, Measurement-Based Quantum Metropolis Algorithm (2022), arXiv:1903.01451 [physics, physics:quant-ph] .
  * Terhal and DiVincenzo [2000] B. M. Terhal and D. P. DiVincenzo, Phys. Rev. A 61, 022301 (2000).
  * Wocjan  _et al._ [2009] P. Wocjan, C.-F. Chiang, D. Nagaj, and A. Abeyesinghe, Phys. Rev. A 80, 022340 (2009).
  * Chiang and Wocjan [2010] C.-F. Chiang and P. Wocjan, Quantum Algorithm for Preparing Thermal Gibbs States - Detailed Analysis (2010), arXiv:1001.1130 [physics, physics:quant-ph] .
  * Wild  _et al._ [2021a] D. S. Wild, D. Sels, H. Pichler, C. Zanoci, and M. D. Lukin, Phys. Rev. Lett. 127, 100504 (2021a).
  * Endo  _et al._ [2020] K. Endo, T. Nakamura, K. Fujii, and N. Yamamoto, Phys. Rev. Research 2, 043442 (2020).
  * Montanaro [2015] A. Montanaro, Proc. R. Soc. A. 471, 20150301 (2015).
  * Yung and Aspuru-Guzik [2012] M.-H. Yung and A. Aspuru-Guzik, Proc. Natl. Acad. Sci. U.S.A. 109, 754 (2012).
  * Arunachalam  _et al._ [2022] S. Arunachalam, V. Havlicek, G. Nannicini, K. Temme, and P. Wocjan, Simpler (classical) and faster (quantum) algorithms for Gibbs partition functions (2022), arXiv:2009.11270 [quant-ph] .
  * Motta  _et al._ [2020] M. Motta, C. Sun, A. T. K. Tan, M. J. O’Rourke, E. Ye, A. J. Minnich, F. G. S. L. Brandão, and G. K.-L. Chan, Nat. Phys. 16, 205 (2020).
  * Chen and Stoudenmire [2020] J. Chen and E. M. Stoudenmire, Phys. Rev. B 101, 195119 (2020).
  * Matsumoto  _et al._ [2022] K. Matsumoto, Y. Shingu, S. Endo, S. Kawabata, S. Watabe, T. Nikuni, H. Hakoshima, and Y. Matsuzaki, Jpn. J. Appl. Phys. 61, 042002 (2022), arXiv:2109.14880 [quant-ph] .
  * Verdon  _et al._ [2019] G. Verdon, J. Marks, S. Nanda, S. Leichenauer, and J. Hidary, Quantum Hamiltonian-Based Models and the Variational Quantum Thermalizer Algorithm (2019), arXiv:1910.02071 [quant-ph] .
  * Sagastizabal  _et al._ [2021] R. Sagastizabal, S. P. Premaratne, B. A. Klaver, M. A. Rol, V. Negîrneac, M. S. Moreira, X. Zou, S. Johri, N. Muthusubramanian, M. Beekman, C. Zachariadis, V. P. Ostroukh, N. Haider, A. Bruno, A. Y. Matsuura, and L. DiCarlo, npj Quantum Inf 7, 130 (2021).
  * Wu and Hsieh [2019] J. Wu and T. H. Hsieh, Phys. Rev. Lett. 123, 220502 (2019).
  * Chowdhury  _et al._ [2020] A. N. Chowdhury, G. H. Low, and N. Wiebe, A Variational Quantum Algorithm for Preparing Quantum Gibbs States (2020), arXiv:2002.00055 [quant-ph] .
  * Shingu  _et al._ [2021] Y. Shingu, Y. Seki, S. Watabe, S. Endo, Y. Matsuzaki, S. Kawabata, T. Nikuni, and H. Hakoshima, Phys. Rev. A 104, 032413 (2021).
  * Wang  _et al._ [2021] Y. Wang, G. Li, and X. Wang, Phys. Rev. Applied 16, 054035 (2021), arXiv:2005.08797 [cond-mat, physics:quant-ph] .
  * Lu  _et al._ [2021] S. Lu, M. C. Bañuls, and J. I. Cirac, PRX Quantum 2, 020321 (2021), arXiv:2006.03032 [quant-ph] .
  * Bassman  _et al._ [2021] L. Bassman, K. Klymko, D. Liu, N. M. Tubman, and W. A. de Jong, Computing Free Energies with Fluctuation Relations on Quantum Computers (2021), arXiv:2103.09846 [cond-mat, physics:quant-ph] .
  * Brandão and Kastoryano [2019] F. G. S. L. Brandão and M. J. Kastoryano, Commun. Math. Phys. 365, 1 (2019).
  * Cohn  _et al._ [2020] J. Cohn, F. Yang, K. Najafi, B. Jones, and J. K. Freericks, Phys. Rev. A 102, 022622 (2020).
  * Martyn and Swingle [2019] J. Martyn and B. Swingle, Phys. Rev. A 100, 032107 (2019).
  * Wild  _et al._ [2021b] D. S. Wild, D. Sels, H. Pichler, C. Zanoci, and M. D. Lukin, Phys. Rev. A 104, 032602 (2021b), arXiv:2109.03007 [cond-mat, physics:quant-ph] .
  * Ge  _et al._ [2016] Y. Ge, A. Molnár, and J. I. Cirac, Phys. Rev. Lett. 116, 080503 (2016).
  * Poulin and Wocjan [2009] D. Poulin and P. Wocjan, Phys. Rev. Lett. 103, 220502 (2009).
  * Wu and Wang [2021] Y. Wu and J. Wang, Estimating Gibbs partition function with quantum Clifford sampling (2021), arXiv:2109.10486 [quant-ph] .
  * Francis  _et al._ [2020] A. Francis, D. Zhu, C. H. Alderete, S. Johri, X. Xiao, J. K. Freericks, C. Monroe, N. M. Linke, and A. F. Kemper, Many Body Thermodynamics on Quantum Computers via Partition Function Zeros (2020), arXiv:2009.04648 [cond-mat, physics:quant-ph] .
  * Bilgin and Boixo [2010] E. Bilgin and S. Boixo, Phys. Rev. Lett. 105, 170405 (2010).
  * Chowdhury and Somma [2016] A. N. Chowdhury and R. D. Somma, Quantum algorithms for Gibbs sampling and hitting-time estimation (2016), arXiv:1603.02940 [quant-ph] .
  * van Apeldoorn  _et al._ [2020] J. van Apeldoorn, A. Gilyén, S. Gribling, and R. de Wolf, Quantum 4, 230 (2020), arXiv:1705.01843 [quant-ph] .
  * Alhambra and Cirac [2021] Á. M. Alhambra and J. I. Cirac, PRX Quantum 2, 040331 (2021), arXiv:2106.00710 [cond-mat, physics:quant-ph] .
  * An  _et al._ [2021] D. An, N. Linden, J.-P. Liu, A. Montanaro, C. Shao, and J. Wang, Quantum 5, 481 (2021), arXiv:2012.06283 [quant-ph, q-fin] .
  * Holbrook [2022] A. J. Holbrook, A quantum parallel Markov chain Monte Carlo (2022), arXiv:2112.00212 [quant-ph, stat] .
  * Herbert [2021] S. Herbert, Quantum Monte-Carlo Integration: The Full Advantage in Minimal Circuit Depth (2021), arXiv:2105.09100 [quant-ph] .
  * Wang and Landau [2001a] F. Wang and D. P. Landau, Phys. Rev. Lett. 86, 2050 (2001a).
  * Wang and Landau [2001b] F. Wang and D. P. Landau, Phys. Rev. E 64, 056101 (2001b).
  * Rall [2021] P. Rall, Quantum 5, 566 (2021), arXiv:2103.09717 [quant-ph] .


#  Lepton Flavor Violation in the Standard Model with dimension 6 operators

Saereh Najjari aaae-mail<EMAIL_ADDRESS>
We investigate Lepton Flavor Violation (LFV) in the general extension of the
Standard Model parametrized by the gauge invariant operators of dimension-6.
We discuss 3-body charged lepton decays and decay of $Z$ boson to lepton pair
and compare the obtained analytical expression with the current experimental
results. We derive the numerical bounds on the size of the Wilson coefficients
parameterizing the allowed size of New Physics effects.

##  1 Parametrization of the SM extensions in terms of effective operators

The renormalizable Standard Model (SM) is probably an effective theory valid
only up to some energy scale where a more fundamental theory could manifest
itself. At the electroweak scale the effects of New Physics (NP) can be
effectively parametrized by new interactions described by non-renormalizable
operators of higher mass dimensions. Coefficients of such operators (called
Wilson coefficients) are suppressed by heavy mass scale $(\Lambda)$ at which
New Physics should become effective. In general the dimension-4 renormalizable
SM Lagrangian can be extended as follows:

$\displaystyle{\cal L}_{\mathrm{SM}}={\cal
L}_{\mathrm{SM}}^{(4)}+\frac{1}{\Lambda}\sum_{k}C_{k}^{(5)}Q_{k}^{(5)}+\frac{1}{\Lambda^{2}}\sum_{k}C_{k}^{(6)}Q_{k}^{(6)}+\ldots\,.$
(1)

where $Q_{k}^{(n)}$ denote the higher-dimension operators and $C_{k}^{(n)}$
stand for the corresponding dimensionless Wilson coefficients.

Using such Lagrangian, the theoretical calculations of relevant physical
observables can be performed in a model-independent way, with final formulae
given in terms of Wilson coefficients. Having such expressions simplifies
significantly a comparison of various SM extension with the experimental
results, as now only the values Wilson coefficients of new operators need be
calculated within a given model of NP - this part of analysis is always model-
dependent.

In the extended SM with the neutrino mass term, the GIM mechanism makes the
branching ratio of the charged lepton flavor violating (CLFV) very small due
to smallness of the mass of the neutrino comparing to the mass of the heaviest
particle in the loop. Experimentally, the CLFV decays have never been observed
yet but there are many models beyond the SM predict sizable rates up to the
current experimental bounds. Here, we investigate the
$\ell_{i}\to\ell_{j}\ell_{j}\ell_{j}$ decay and $Z\to\ell_{f}^{+}\ell_{i}^{-}$
decays in the extension of the SM with dimension 6 operators.The full list of
operators of dimension 5 and 6 which can be constructed out of SM fields is
given in $\\!{}^{{\bf?}}$.

After simplification only 8 operators of dimension 6 give contribution to
these decays at tree level:

  * •

2 dipole-type operators, contributing to tree-level flavor non-diagonal $Z$
and photon couplings to leptons:

$Q_{eW}=(\bar{\ell}_{i}\sigma^{\mu\nu}e_{j})\tau^{I}\varphi W_{\mu\nu}^{I}$,
$Q_{eB}=(\bar{\ell}_{i}\sigma^{\mu\nu}e_{j})\varphi B_{\mu\nu}$.

  * •

3 operators modifying tree-level $Z$ and $W$ couplings:

$Q_{\varphi\ell}^{(1)}=(\mbox{${\varphi^{\dagger}i\,\raisebox{5.69054pt}{\boldmath${}^{\leftrightarrow}$}D_{\mu}\,\varphi}$})(\bar{\ell}_{i}\gamma^{\mu}\ell_{j})$,
$Q_{\varphi\ell}^{(3)}=(\mbox{${\varphi^{\dagger}i\,\raisebox{5.69054pt}{\boldmath${}^{\leftrightarrow}$}D_{\mu}^{\,I}\,\varphi}$})(\bar{\ell}_{i}\tau^{I}\gamma^{\mu}\ell_{j})$,
$Q_{\varphi
e}=(\mbox{${\varphi^{\dagger}i\,\raisebox{5.69054pt}{\boldmath${}^{\leftrightarrow}$}D_{\mu}\,\varphi}$})(\bar{e}_{i}\gamma^{\mu}e_{j})$

  * •

3 four-lepton contact operators:

$Q_{\ell
e}=(\bar{\ell}_{i}\gamma_{\mu}\ell_{j})(\bar{e}_{k}\gamma^{\mu}e_{l})$,
$Q_{\ell\ell}=(\bar{\ell_{i}}\gamma_{\mu}\ell_{j})(\bar{\ell}_{k}\gamma^{\mu}\ell_{l})$,
$Q_{ee}=(\bar{e}_{i}\gamma_{\mu}e_{j})(\bar{e}_{k}\gamma^{\mu}e_{l})$

In $\\!{}^{{\bf?}}$ we list the Feynman rules arising from dimension 6
operators.

##  2 Lepton Flavor violation in 3-body charged lepton and $Z$ to lepton pair
decays

In this section we compare the analytical results with experimental bounds on
the $\ell_{i}\to\ell_{j}\ell_{j}\ell_{j}$ decay and
$Z\to\ell_{f}^{+}\ell_{i}^{-}$ decays to constrain the Wilson
coefficients.Such decays can be generated already at tree level. Our results
for radiative lepton decays are given in $\\!{}^{{\bf?},{\bf?}}$. The
experimental bound on these decays are given in Tables 1 and 2.

Process | Experimental bound
---|---
${\cal B}\left[\tau^{-}\to\mu^{-}\mu^{+}\mu^{-}\right]$ |  $2.1\times 10^{-8}$ $\\!{}^{{\bf?}}$
${\cal B}\left[\tau^{-}\to e^{-}e^{+}e^{-}\right]$ |  $2.7\times 10^{-8}$ $\\!{}^{{\bf?}}$
${\cal B}\left[\mu^{-}\to e^{-}e^{+}e^{-}\right]$ |  $1.0\times 10^{-12}$ $\\!{}^{{\bf?}}$
Table 1: Experimental upper limits on the branching ratios of the three body charged lepton decays. Process | Experimental bound
---|---
$\mathrm{Br}\left[Z^{0}\to\mu^{\pm}e^{\mp}\right]$ |  $1.7\times 10^{-6}$ $\\!{}^{{\bf?}}$
$\mathrm{Br}\left[Z^{0}\to\tau^{\pm}e^{\mp}\right]$ |  $9.8\times 10^{-6}$ $\\!{}^{{\bf?}}$
$\mathrm{Br}\left[Z^{0}\to\tau^{\pm}\mu^{\mp}\right]$ |  $1.2\times 10^{-5}$ $\\!{}^{{\bf?}}$
Table 2: Experimental upper limits (95 % CL) on the lepton flavor violating
$Z$ decay rates.

In general, at the tree-level diagrams mediated by photon, $Z^{0}$ boson and
4-lepton contact interactions can contribute to the three body charged lepton
decays $\mu\to 3e$, $\tau\to 3e$ and $\tau\to 3\mu$. The general expression
for the Br($\ell_{i}\to 3\ell_{j}$) reads as:

$\displaystyle\mathrm{Br}(\ell_{i}\to\ell_{j}\ell_{j}\bar{\ell}_{j})$
$\displaystyle=$
$\displaystyle\frac{m_{\ell_{i}}^{5}}{12288\pi^{3}\Lambda^{4}\Gamma_{\ell_{i}}}\left(4\left(|C_{VLL}|^{2}+|C_{VRR}|^{2}+|C_{VLR}|^{2}+|C_{VRL}|^{2}\right)\right.$
(2) $\displaystyle+$ $\displaystyle|C_{SLR}|^{2}+|C_{SRL}|^{2}+48X_{\gamma}$

where $\Gamma^{\ell_{i}}$ is the total decay width of the initial lepton. The
$X_{\gamma}$ (photon contribution) and $C_{X}$ ($Z$ and contact interactions)
can be written in terms of Wilson coefficients of operators listed in previous
Section as:

$\displaystyle X_{\gamma}$ $\displaystyle=$
$\displaystyle-\frac{16ev}{m_{\ell_{i}}}\mathrm{Re}\left[\left(2C_{VLL}+C_{VLR}-\frac{1}{2}C_{SLR}\right)C_{\gamma
R}^{\star}+\left(2C_{VRR}+C_{VRL}-\frac{1}{2}C_{SRL}\right)C_{\gamma
L}^{\star}\right]$ (3) $\displaystyle+$
$\displaystyle\frac{64e^{2}v^{2}}{m_{\ell_{i}}^{2}}\left(\log\frac{m_{\ell_{i}}^{2}}{m_{\ell_{j}}^{2}}-\frac{11}{4}\right)(|C_{\gamma
L}|^{2}+|C_{\gamma R}|^{2})$

$\displaystyle C_{VLL}$ $\displaystyle=$ $\displaystyle
2\left((2s_{W}^{2}-1)\left(C_{\varphi\ell}^{(1)ji}+C_{\varphi\ell}^{(3)ji}\right)+C_{\ell\ell}^{jijj}\right)$
$\displaystyle C_{VRR}$ $\displaystyle=$ $\displaystyle
2\left(2s_{W}^{2}C_{\varphi e}^{ji}+C_{ee}^{jijj}\right)$ $\displaystyle
C_{VLR}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}C_{SRL}=\left(2s_{W}^{2}\left(C_{\varphi\ell}^{(1)ji}+C_{\varphi\ell}^{(3)ji}\right)+C_{le}^{jijj}\right)$
$\displaystyle C_{VRL}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}C_{SLR}=\left((2s_{W}^{2}-1)C_{\varphi
e}^{ji}+C_{\ell e}^{jjji}\right)$ $\displaystyle C_{\gamma L}^{ij}$
$\displaystyle=$ $\displaystyle C_{\gamma
R}^{ji\star}=2\sqrt{2}s_{W}\left(c_{W}C_{eB}^{ij\star}-s_{W}C_{eW}^{ij\star}\right)$
(4)

Knowing that $C_{\gamma}$ must be negligible, we neglect contribution from
$C_{\gamma}$ in order to constrain the Wilson coefficients by using the bounds
on $\ell_{i}\to\ell_{j}\ell_{j}\ell_{j}$ decays. Normalizing their branching
ratio to the limits in Table 1 we find following numerical equations
constraining the Wilson coefficients:

$\displaystyle{C_{\mu eee}}$ $\displaystyle\leq$ $\displaystyle
3.29\times{10^{-5}}{\left({\frac{\Lambda}{{1\;{\rm{TeV}}}}}\right)^{2}}\sqrt{\frac{{{\rm{Br}}\left[{\mu\to
eee}\right]}}{{1\times{{10}^{-12}}}}}\,,$ $\displaystyle{C_{\tau eee}}$
$\displaystyle\leq$ $\displaystyle
1.28\times{10^{-2}}{\left({\frac{\Lambda}{{1\;{\rm{TeV}}}}}\right)^{2}}\sqrt{\frac{{{\rm{Br}}\left[{\tau\to
eee}\right]}}{{2.7\times{{10}^{-8}}}}}\,,$ (5)
$\displaystyle{C_{\tau\mu\mu\mu}}$ $\displaystyle\leq$ $\displaystyle
1.13\times{10^{-2}}{\left({\frac{\Lambda}{{1\;{\rm{TeV}}}}}\right)^{2}}\sqrt{\frac{{{\rm{Br}}\left[{\tau\to\mu\mu\mu}\right]}}{{2.1\times{{10}^{-8}}}}}\,,$

with ${C_{{\ell_{i}}{\ell_{j}}{\ell_{j}}{\ell_{j}}}}$ given by

$\displaystyle\begin{array}[]{l}{C_{{\ell_{i}}{\ell_{j}}{\ell_{j}}{\ell_{j}}}}=\left({2{{\left|{{\rm{C}}_{\ell
e}^{jijj}-0.54\left({{\rm{C}}_{\varphi\ell}^{\left(1\right)ji}+{\rm{C}}_{\varphi\ell}^{\left(3\right)ji}}\right)}\right|}^{2}}}\par+2{{\left|{{\rm{C}}_{ee}^{jijj}+0.46{\rm{C}}_{\varphi
e}^{ji}}\right|}^{2}}\par\right.\\\
\left.{\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\par\par+{{\left|{{\rm{C}}_{\ell
e}^{jijj}}+0.46\left({{\rm{C}}_{\varphi\ell}^{\left(1\right)ji}+{\rm{C}}_{\varphi\ell}^{\left(3\right)fi}}\right)\right|}^{2}}\par+{{\left|{{\rm{C}}_{\ell
e}^{jjji}-0.54{\rm{C}}_{\varphi
e}^{ji}}\right|}^{2}}\par}\right)^{\frac{1}{2}}.\end{array}$ (8)

Similar analysis can be done for the LFV $Z\to\ell_{f}^{+}\ell_{i}^{-}$
decays. At the tree level the branching ratio is:

$\displaystyle\mathrm{Br}\left({Z^{0}\to\ell_{f}^{\pm}\ell_{i}^{\mp}}\right)=\frac{m_{Z}}{24\pi\Gamma_{Z}}\left[\frac{m_{Z}^{2}}{2}\left(\left|C_{fi}^{ZR}\right|^{2}+\left|C_{fi}^{ZL}\right|^{2}\right)+\left|\Gamma_{fi}^{ZL}\right|^{2}+\left|\Gamma_{fi}^{ZR}\right|^{2}\right]\,,$
(9)

with $\Gamma_{Z}\approx 2.495$ GeV being the total decay width of $Z$ boson.

where we included all tree-level contributions and

$\displaystyle\Gamma^{ZL}_{fi}$ $\displaystyle=$
$\displaystyle\dfrac{e}{2s_{W}c_{W}}\left(\dfrac{v^{2}}{\Lambda^{2}}\left(C_{\varphi
l}^{(1)fi}+C_{\varphi
l}^{(3)fi}\right)+\left(1-2s_{W}^{2}\right)\delta_{fi}\right)\,,$ (10)
$\displaystyle\Gamma^{ZR}_{fi}$ $\displaystyle=$
$\displaystyle\dfrac{e}{2s_{W}c_{W}}\left(\dfrac{v^{2}}{\Lambda^{2}}C_{\varphi
e}^{fi}-2s_{W}^{2}\delta_{fi}\right)\,,$ (11) $\displaystyle C^{ZR}_{fi}$
$\displaystyle=$ $\displaystyle C^{ZL\star}_{if}=-v\sqrt{2}C_{Z}^{fi}$ (12)

where $C_{Z}^{fi}$ is defined as

$\displaystyle C_{Z}^{fi}=\left(s_{W}C_{eB}^{fi}+c_{W}C_{eW}^{fi}\right)\,.$
(13)

In the experimental values for branching ratio are for the sum
$Z\to\ell_{f}^{\pm}\ell_{i}^{\mp}+\ell_{f}^{\mp}\ell_{i}^{\pm}$
$\\!{}^{{\bf?}}$ while in this equation branching ratio is for the decay
$Z\to\ell_{f}^{\pm}\ell_{i}^{\mp}$ or $\ell_{f}^{\mp}\ell_{i}^{\pm}$.
Therefore this equation must be multiply by a factor 2 in order to compare
into experimental value.

Again normalizing the formulae for the branching ratios to the current
experimental bound listed in Table 2 we derived the numerical equations
constraining the Wilson coefficients contributing to this decay:

$\displaystyle\sqrt{\left|C_{\varphi\ell}^{(1)12}+C_{\varphi\ell}^{(3)12}\right|^{2}+\left|C_{\varphi
e}^{12}\right|^{2}+\left|C_{Z}^{12}\right|^{2}+\left|C_{Z}^{21}\right|^{2}}$
$\displaystyle\leq
0.06\left(\frac{\Lambda}{1\;\rm{TeV}}\right)^{2}\sqrt{\dfrac{{\rm
Br}\left[Z^{0}\to\mu^{\pm}e^{\mp}\right]}{1.7\times 10^{-6}}}\,,$
$\displaystyle{}\sqrt{\left|C_{\varphi\ell}^{(1)13}+C_{\varphi
e}^{13}\right|^{2}+\left|C_{\varphi
e}^{13}\right|^{2}+\left|C_{Z}^{13}\right|^{2}+\left|C_{Z}^{31}\right|^{2}}$
$\displaystyle\leq
0.14\left(\frac{\Lambda}{1\;\rm{TeV}}\right)^{2}\sqrt{\dfrac{{\rm
Br}\left[Z^{0}\to\tau^{\pm}e^{\mp}\right]}{9.8\times 10^{-6}}}\,\,,$
$\displaystyle{}\sqrt{\left|C_{\varphi\ell}^{(1)23}+C_{\varphi\ell}^{(3)23}\right|^{2}+\left|C_{\varphi
e}^{23}\right|^{2}+\left|C_{Z}^{23}\right|^{2}+\left|C_{Z}^{32}\right|^{2}}$
$\displaystyle\leq
0.16\left(\frac{\Lambda}{1\;\rm{TeV}}\right)^{2}\sqrt{\dfrac{{\rm
Br}\left[Z^{0}\to\tau^{\pm}\mu^{\mp}\right]}{1.2\times 10^{-5}}}\,.$

These constraints are less stringent than the ones from
$\ell_{i}\to\ell_{j}\ell_{j}\ell_{j}$ decays.

##  3 Conclusions

We have calculated the rates of several lepton flavor violating decays in the
Standard Model extended with dimension 6 operators. We present numerical
expressions constraining the relevant Wilson coefficients, based on
experimental bounds for the $\ell_{i}\to\ell_{j}\ell_{j}\ell_{j}$ decay and
$Z\to\ell_{f}^{+}\ell_{i}^{-}$ decays. There are also some experiments
searching for charged lepton flavor violation and are going to upgrade the
sensitivity. Observation of charged lepton flavor violation at experiment
would be a clear hint for physics beyond the standard model. We show that the
bounds on lepton flavor violation couplings are already very strong if the
scale of New Physics is low.

## Acknowledgments

This work has been done in collaboration with A. Crivellin and J. Rosiek. I
acknowledges financial support from the organizing committee to attend the
conference.

## References

## References

  * [1] B. Grzadkowski, M. Iskrzynski, M. Misiak and J. Rosiek, JHEP 1010, 085 (2010) [arXiv:1008. 4884 [hep-ph]].
  * [2] A. Crivellin, S. Najjari and J. Rosiek, JHEP 1404, 167 (2014) [arXiv:1312.0634 [hep-ph]].
  * [3] S. Najjari, arXiv:1305.6819 [hep-ph].
  * [4] K. Hayasaka, K. Inami, Y. Miyazaki, K. Arinstein, V. Aulchenko, T. Aushev, A. M. Bakich and A. Bay et al., Phys. Lett. B 687, 139 (2010) [arXiv:1001.3221 [hep-ex]].
  * [5] U. Bellgardt et al., “Search for the Decay $\mu^{+}\to e^{+}e^{+}e^{-}$,” Nucl.Phys., vol. B299, p. 1, 1988.
  * [6] P. Abreu et al. [DELPHI Collaboration], Z. Phys. C 73, 243 (1997).
  * [7] J. I. Illana, M. Jack and T. Riemann, In *2nd ECFA/DESY Study 1998-2001* 490-524 [hep-ph/0001273].