This is a follow-up post on hierarchical compartmental reserving models using PK/PD models. It will show how differential equations can be used with Stan/ brms and how correlation for the same group level terms can be modelled.
PK/ PD is usually short for pharmacokinetic/ pharmacodynamic models, but as Eric Novik of Generable pointed out to me, it could also be short for Payment Kinetics/ Payment Dynamics Models in the insurance context.
In this post I’d like to discuss an extension to Jake Morris’ hierarchical compartmental reserving model, as described in his original paper (Morris (2016)) and my previous post, which allows for a time varying parameter \(k_{er}(t)\) describing the changing rate of earning and reporting and allows for correlation between \(RLR\), the reported loss ratio, and \(RRF\), the reserve robustness factor. A positive correlation would give evidence of a reserving cycle, i.e. in years with higher initial reported loss ratios a higher proportion of reported losses are paid.
I am very grateful for the support Jake Morris and Paul-Christian Bürkner have given me over the last few weeks, answering questions around the model and brms (Bürkner (2017)).
PK/ PD model
First of all let’s take a look at my new set of differential equations, describing the dynamics of exposure (\(EX\)), outstanding (\(OS\)) and paid claims (\(PD\)) over time. I replaced the constant parameter \(k_{er}\) with the linear function \(\beta_{er} \cdot t\).
\[
\begin{aligned}
dEX/dt &amp; = -\beta_{er} \cdot t \cdot EX \\
dOS/dt &amp; = \beta_{er} \cdot t \cdot RLR \cdot EX - k_p \cdot OS \\
dPD/dt &amp; = k_{p} \cdot RRF \cdot OS
\end{aligned}
\]
The dynamical system is no longer autonomous and initially I can’t be bothered to solve it analytically. Hence, I use an ODE solver instead, but I will get back to integrating the differential equations later.
Integrating ODEs with Stan
Fortunately, an ODE solver is part of the Stan language. The following code demonstrates how I can integrate the differential equations in Stan. You will notice that the model section is empty, as all I want to do for now is run the solver.
functions{
real[] claimsprocess(real t, real [] y, real [] theta, real [] x_r, int[] x_i){
real dydt[3];
dydt[1] = - theta[1] * t * y[1];
dydt[2] = theta[1] * t * theta[3] * y[1] - theta[2] * y[2];
dydt[3] = theta[2] * theta[4] * y[2];
return dydt;
}
}
data {
real theta[4];
int N_t;
real times[N_t]; real C0[3];
}
parameters{
}
transformed parameters{
real C[N_t, 3];
C = integrate_ode_rk45(claimsprocess, C0, 0, times, theta,
rep_array(0.0, 0), rep_array(1, 1));
}
model{
}
generated quantities{
real Exposure[N_t];
real OS[N_t];
real Paid[N_t];
Exposure = C[, 1];
OS = C[, 2];
Paid = C[, 3];
}
To run the Stan code I provide as an input a list with the parameters, the number of data points I’d like to generate, the time line and initial starting values. Note also that I have to set the algorithm = &quot;Fixed_param&quot; as I run a deterministic model and for that reason one iteration and one chain are sufficient.
library(rstan)
input_data