Bayesian vs Limited Fluctuation Experience Analysis

Formulas

The Limited Fluctuation Method components: a credibility weight, \[Z\], an Observed Rate, and a Prior Rate.

\[ \text{Credibility-Weighted Rate} = Z × \text{Observed Rate} + (1 – Z) × \text{Prior Rate} \]

With probability equal to \(LC_p\) that the Observed Rate does not differ from the true rate by more than \(LC_r\).

LC_p = 0.9
LC_r = 0.05

0.05

\[ \text{Claims for full credibility} = \left(\frac{\text{Z-Score}}{\text{ratio}}\right)^{2} \]

LC_full_credibility = round(Int, (quantile(Normal(), 1 - (1 - LC_p) / 2) / LC_r)^2)

1082

Using the inputs above, we have a z-score corresponding to 0.9 of 1.6448536269514717, so: ( 1.6448536269514717 ÷ 0.05)² ≈ 1082. ““”

Atkinson goes on to nicely summarize the square root method which assigns full credibility, \[Z = 1\], when the number of actual claims equals or exceeds the full credibility threshold of $LC_full_credibility claims.

When the number of claims is less than full credibility:

\[ Z = \sqrt{\frac{\text{no. claims in group}}{\text{no. claims full credibility}}} \]

Issues with Limited Fluctuation

• Ignores available information:
  • Why does our \(Z\) not vary by exposures?
  • No mechanism for a related group to inform the credibility of another
  • Results in a single point estimate with only approximate measure of estimate uncertainty
• Relies on a number of assumptions/constraints:
  • That aggregated claims can be approximated by a normal distribution
  • The thresholds of $LC_p and $LC_r are arbitrary (e.g. all of the same arguments against p-value thresholds can apply to this approach) ““”

The formulation

Contrary to Limited Fluctuation, the Bayesian approach forces one to be explicit about the presumed structure of the probability model. The flexibility of the statistical model allows one to incorporate actuarial judgement in a quantitative way. For example, in this example we assume that the claims experience of each group informs a global (hyperparameter) prior distribution which we could use as a starting point for a new type of observation. More on this once the data is introduced.

Discussion of Bayesian model

Each group (make and model combination) has an expected claim rate \(μ_i\), which is informed by the global hyperparameter \(μ\) and variance \(\sigma^2\). Then, the observed claims by group are assumed to be distributed according to a Poisson distribution.

A complete overview of modern Bayesian models is beyond the scope of this notebook, but a few key points:

• We are forced with the Bayesian approach to be explicit about the assumptions (versus all of the implicit assumptions of alternative techniques like LF)
• We set priors which are the assumed distribution of model parameters before “learning” from the observations. With enough data, it can result in a posterior that the prior was very skeptical of beforehad.
• With the volume of data in the training set, the priors we select here are not that important. Some comments on why they were chosen though:
  • The rate of claims is linked with a logistic function, which looks like an integral sign of sorts. logistic(0.0) equals 0.5 while very negative inputs approach 0.0 and positive numbers approach 1.0. We do this so that the rate of claims is always constrained in the range \((0,1)\)
  • μ \sim Normal(-2,4) says that we expect the population of auto claims rate to be less than 0.5 (`logistic(-2) ≈ .12) but very wide range of possible values given the wide standard deviation.
  • σ \sim Exponential(0.25) says that we expect the standard devation of an individual group to be positive, but not super wide.
  • μ_i \sim Normal(μ,σ) is the actual prior for each group’s rate of claim and is informed by the above hyperparameters.
• @model is a Turing.jl macro which enables interpreting the nice ~ syntax, which makes the Julia code look very similar to the traditional mathematical notation.
• Note the use of broadcasting with the dot syntax (e.g. the . in .~). This tells Julia to vectorize and fuse the computation.
  • data.claims ~ product_distribution(Poisson.(data.n .* logistic.(μ_i))) means “each value in data.claims is a random outcome (~) distributed accroding to a corresponding Poisson distribution with \lambda =n \times \text{logistic}(\mu_i) where text{logistic}(\mu_i) is the average claims rate for each group.

@model function model_poisson(data, claims)
    # hyperparameter that informs the prior for each group
    μ ~ Normal(-2, 4)
    σ ~ Exponential(0.25)

    # the random variable representing the average claim rate for each group
    # filldist creates a set of random variable without needing to list out each one
    μ_i ~ filldist(Normal(μ, σ), length(unique(data.Blind_Model)))

    # use the poisson appproximation to the binomial claim outcome with a
    # logisitc link function to keep the probability between 0 and 1
    claims ~ product_distribution(Poisson.(data.n .* logistic.(μ_i)))
end

model_poisson (generic function with 2 methods)

Here’s what the prior distribution of claims looks like… and peeking ahead to what the posterior for a single group of claims looks like. See how even though our posterior was very wide (favoring small claims rates), it was dominated by the data to create a very narrow posterior average claims rate.

The model is combined with the data and Turing.jl is used to computationally arrive at the posterior distribution for the parameters in the statistical model, \(\mu\), \(\mu_i\), and \(\sigma\).

mp = let
    # combine the different years in the training set
    condensed = @chain train begin
        groupby(:Blind_Model)
        @combine begin
            :n = sum(:n)
            :claims = sum(:claims)
        end
    end
    model_poisson(condensed, condensed.claims)
end

DynamicPPL.Model{typeof(model_poisson), (:data, :claims), (), (), Tuple{DataFrame, Vector{Int64}}, Tuple{}, DynamicPPL.DefaultContext, false}(model_poisson, (data = 1238×3 DataFrame
  Row │ Blind_Model  n       claims 
      │ String7      Int64   Int64  
──────┼─────────────────────────────
    1 │ K.78          14058      94
    2 │ Q.22          71401     532
    3 │ AR.41         37587     370
    4 │ D.20          10019      50
    5 │ AJ.129       100824     693
    6 │ AQ.17         86702     541
    7 │ BW.3          75655     410
    8 │ BW.167        42737     294
  ⋮   │      ⋮         ⋮       ⋮
 1232 │ BG.8              1       0
 1233 │ CA.5              2       0
 1234 │ J.5               2       0
 1235 │ AQ.1              1       0
 1236 │ AJ.96             1       0
 1237 │ AJ.51             2       0
 1238 │ BQ.6              2       0
                   1223 rows omitted, claims = [94, 532, 370, 50, 693, 541, 410, 294, 755, 162  …  0, 0, 0, 0, 0, 0, 0, 0, 0, 0]), NamedTuple(), DynamicPPL.DefaultContext())

Bayesian Posterior Sampling

Here’s where recent advances in algorithms and computing power make Bayesian analyis possible. We can’t analytically compute the posterior distribution, but we can genererate samples from the posterior such that the frequency of the sampled result appears in proportion to the true posterior density. This uses a technique called Markov Chain Monte Carlo, abbreviated MCMC. There are different algorithms in this family, and we will use one called the No-U-Turn Sampler (NUTS for short).

The result of the sample function is a set of data containing data that is generated in proportion to the posterior distributions and we will use this data to make predictions and understand the distribution of our parameters.

Run the chain:

cp = MCMCChains.Chains(sample(mp, NUTS(), 500; check_model=false)) # this is the line that runs the MCMC sampler

# ╔═╡ 6b5d2186-78ac-48ae-8529-e0dcd51d0b00
let
    # sample from the priors before learning from any data
    ch_prior = MCMCChains.Chains(sample(mp, Prior(), 1500; check_model=false))

    f = Figure()
    ax = Axis(f[1, 1], title="Prior distribution of expected claims rate for a group")
    μ = vec(ch_prior["μ_i[1]"].data)
    hist!(ax, logistic.(μ); bins=50)
    ax2 = Axis(f[2, 1], title="Posterior distribution of expected claims rate for a group")
    μ = vec(cp["μ_i[1]"].data)
    hist!(ax2, logistic.(μ); bins=50)
    linkxaxes!(ax, ax2)
    f
end

Visualizing the Posterior Density

This is a plot of the posterior density for all of the many, many \(\mu\) parameters in our model. The black line shows the distribution which comes from our hyperparameters μ and σ. In the event of coming across a new group of interest (in our case, a new make of car), this is the prior distribution for the expected claims rate. The model has learned this from the data itself, and serves to regularize predictions.

"""
    dplot(chain,parameter_string)

plot the posterior density for model parameters matching the given string.
"""
function dplot(ch, param_string)
    f = Figure()
    ax = Axis(f[1, 1])

    # plot each group posterior
    for (ind, param) in enumerate(ch.name_map.parameters)
        if contains(string(param), param_string)
            v = vec(getindex(ch, param).data)
            density!(
                ax, logistic.(v), color=(:red, 0.0),
                strokecolor=(:red, 0.3), strokewidth=1, strokearound=false, label="Individual group posterior"
            )
        end
    end

    # Plot hyperparameters
    hyper_mean = mean(getindex(ch, :μ).data)
    hyper_sigma = mean(getindex(ch, :σ).data)
    d = density!(
        ax, logistic.(rand(Normal(hyper_mean, hyper_sigma), 1000)), color=(:black, 0.0),
        strokecolor=(:black, 0.3), strokewidth=3, strokearound=false, label="Hyper-parameter distribution"
    )


    xlims!(0.0, 0.03)
    hideydecorations!(ax)
    axislegend(ax, unique=true)
    return f
end

dplot(cp, "μ")

Predictions

Bayesian approach has lower error

Looking at the accumulated absolute error, the bayesian approach has about 16% lower error than the Limited Fluctuation approaches.

Total Actual to Expected

Here, offsetting errors happen to make the naive and LF year-by-year (the latter is a good proxy for the former) such that their total A/E is better than the Bayesian approach.

Given the lower error of the Bayesian approach, one would expect that over multiple years that it would produce more accurate predictions than the alternative methods.

Further thoughts

The Bayesian approach could be extended to improve its accuracy even further:

• Using vehicle make data to add more hierarchical structure to the statistical model. For example, one may observe that Porsches experience crashes at a higher rate than Volvos. LFC cannot embed that sort of overlapping hierarchy into its framework.
• The Bayesian hyperparameter provides a framework to think about “unseen” make-model combinations

Downsides to the Bayesian approach

• Computationally intensive. Complex @models can take very long time to run (many hours), compared to relatively quick frequentest methods like maximum likelihood estimation.

Further Reading

If this notebook has piqued your interest in Bayesian techniques, the following books are recommended learning resources (from easiest to most difficult):

• Statistical Rethinking
  • See also StatisticalRethinking.jl for Julia implementations of all of the book’s code and utility functions
• Bayes Rules!
• Bayeisan Data Analysis