Example Stochastic Reserving

Roger Hayne

9/6/2019

library(mvtnorm)
#library(MASS)
library(abind)
library(stochasticreserver)

Initialize Triangle

Input (B0) is a development array of cumulative averages with a the exposures (claims) used in the denominator appended as the last column. Assumption is for the same development increments as exposure increments and that all development lags with no development have # been removed. Data elements that are not available are indicated as such. This should work (but not tested for) just about any subset of an upper triangular data matrix.
Another requirement of this code is that the matrix contain no columns that are all zero.

Negative Loglikelihood Function to be Minimized

Note that the general form of the model has parameters in addition to those in the loss model, namely the power for the variance and the constant of proprtionality that varies by column. So if the original model has k parameters with size columns of data, the total objective function has k + size + 1 parameters

l.obj <- function(a, A) {
  npar <- length(a) - 2
  e <- g_obj(a[1:npar])
  v <- exp(-outer(logd[, 1], rep(a[npar + 1], size), "-")) * (e^2)^a[npar + 2]
  t1 <- log(2 * pi * v) / 2
  t2 <- (A - e) ^ 2 / (2 * v)
  sum(t1 + t2, na.rm = TRUE)
}
# Gradient of the objective function
l.grad <- function(a, A) {
  npar <- length(a) - 2
  p <- a[npar + 2]
  Av <- aperm(array(A, c(size, size, npar)), c(3, 1, 2))
  e <- g_obj(a[1:npar])
  ev <- aperm(array(e, c(size, size, npar)), c(3, 1, 2))
  v <- exp(-outer(logd[, 1], rep(a[npar + 1], size), "-")) * (e^2)^p
  vv <- aperm(array(v, c(size, size, npar)), c(3, 1, 2))
  dt <- rowSums(g_grad(a[1:npar]) * ((p / ev) + (ev - Av) / vv - p * 
                                      (Av - ev)^2 / (vv * ev)),
               na.rm = TRUE,
               dims = 1)
  yy <- 1 - (A - e) ^ 2 / v
  dk <- sum(yy / 2, na.rm = TRUE)
  dp <- sum(yy * log(e ^ 2) / 2, na.rm = TRUE)
  c(dt, dk, dp)
}

Hessian of the objective function

  • e is the expectated value matrix
  • v is the matrix of variances
  • A, e, v all have shape c(size, size)
  • The variables _v are copies of the originals to shape c(npar,size,size), paralleling the gradient of g.
  • The variables _m are copies of the originals to shape c(npar,npar,size,size), paralleling the hessian of g
l.hess <- function(a, A) {
  npar <- length(a) - 2
  p <- a[npar + 2]
  Av <- aperm(array(A, c(size, size, npar)), c(3, 1, 2))
  Am <- aperm(array(A, c(size, size, npar, npar)), c(3, 4, 1, 2))
  e <- g_obj(a[1:npar])
  ev <- aperm(array(e, c(size, size, npar)), c(3, 1, 2))
  em <- aperm(array(e, c(size, size, npar, npar)), c(3, 4, 1, 2))
  v <- exp(-outer(logd[, 1], rep(a[npar + 1], size), "-")) * (e ^ 2) ^ p
  vv <- aperm(array(v, c(size, size, npar)), c(3, 1, 2))
  vm <- aperm(array(v, c(size, size, npar, npar)), c(3, 4, 1, 2))
  g1 <- g_grad(a[1:npar])
  gg <- aperm(array(g1, c(npar, size, size, npar)), c(4, 1, 2, 3))
  gg <- gg * aperm(gg, c(2, 1, 3, 4))
  gh <- g_hess(a[1:npar])
  dtt <- rowSums(
    gh * (p / em + (em - Am) / vm - p * (Am - em) ^ 2 / (vm * em)) +
      gg * (
        1 / vm + 4 * p * (Am - em) / (vm * em) + p * (2 * p + 1) * (Am - em) ^ 2 /
          (vm * em ^ 2) - p / em ^ 2
      ),
    dims = 2,
    na.rm = TRUE
  )
  dkt <- rowSums((g1 * (Av - ev) + p * g1 * (Av - ev) ^ 2 / ev) / vv, na.rm = TRUE)
  dtp <- rowSums(g1 * (1 / ev + (
    log(ev ^ 2) * (Av - ev) + (p * log(ev ^ 2) - 1) * (Av - ev) ^ 2 / ev
  ) / vv),
  na.rm = TRUE)
  dkk <- sum((A - e) ^ 2 / (2 * v), na.rm = TRUE)
  dpk <- sum(log(e ^ 2) * (A - e) ^ 2 / (2 * v), na.rm = TRUE)
  dpp <- sum(log(e ^ 2) ^ 2 * (A - e) ^ 2 / (2 * v), na.rm = TRUE)
  m1 <- rbind(array(dkt), c(dtp))
  rbind(cbind(dtt, t(m1)), cbind(m1, rbind(cbind(dkk, c(
    dpk
  )), c(dpk, dpp))))
}

End of funciton specificaitons now on to the minimization

Minimization

Get starting values for kappa and p parameters, default 10 and 1

ttt <- c(10, 1)

For starting values use fitted objective function and assume variance for a cell is estimated by the square of the difference between actual and expected averages. Note since log(0) is -Inf we need to go through some machinations to prep the y values for the fit

E <- g_obj(a0)
yyy <- (A0 - E)^2
yyy <- logd + log(((yyy != 0) * yyy) - (yyy == 0))
sss <- na.omit(data.frame(x = c(log(E^2)), y = c(yyy)))
ttt <- array(coef(lm(sss$y ~ sss$x)))[1:2]
a0 <- c(a0, ttt)

set.seed(1) # to check reproducibility with original code
max <- list(iter.max = 10000, eval.max = 10000)

Actual minimization

mle <- nlminb(a0, l.obj, gradient = l.grad, hessian = l.hess, 
             scale = abs(1 / (2 * ((a0 * (a0 != 0)) + (1 * (a0 == 0))))),
             A = A0, control = max)

Model statistics

  • mean and var are model fitted values
  • stres is the standardized residuals
npar <- length(a0) - 2
p <- mle$par[npar + 2]
mean <- g_obj(mle$par[1:npar])
var <- exp(-outer(logd[, 1], rep(mle$par[npar + 1], size), "-")) * (mean ^
                                                                   2) ^ p
stres <- (A0 - mean) / sqrt(var)
g1 <- g_grad(mle$par[1:npar])
gg <- aperm(array(g1, c(npar, size, size, npar)), c(4, 1, 2, 3))
gg <- gg * aperm(gg, c(2, 1, 3, 4))
meanv <- aperm(array(mean, c(size, size, npar)), c(3, 1, 2))
meanm <- aperm(array(mean, c(size, size, npar, npar)), c(3, 4, 1, 2))
varm <- aperm(array(var, c(size, size, npar, npar)), c(3, 4, 1, 2))

Masks to screen out NA entries in original input matrix

s <- 0 * A0
sv <- aperm(array(s, c(size, size, npar)), c(3, 1, 2))
sm <- aperm(array(s, c(size, size, npar, npar)), c(3, 4, 1, 2))

Calculate the information matrix

  • Using second derivatives of the log likelihood function Second with respect to theta parameters
tt <- rowSums(sm + gg * (1 / varm + 2 * p ^ 2 / (meanm ^ 2)), dims = 2, na.rm = TRUE)

Second with respect to theta and kappa

kt <- p * rowSums(sv + g1 / meanv, na.rm = TRUE)

Second with respect to p and theta

tp <- p * rowSums(sv + g1 * log(meanv ^ 2) / meanv, na.rm = TRUE)

Second with respect to kappa

kk <- (1 / 2) * sum(1 + s, na.rm = TRUE)

Second with respect to p and kappa

pk <- (1 / 2) * sum(s + log(mean ^ 2), na.rm = TRUE)

Second with respect to p

pp <- (1 / 2) * sum(s + log(mean ^ 2) ^ 2, na.rm = TRUE)

Create information matrix in blocks

m1 <- rbind(array(kt), c(tp))
inf <- rbind(cbind(tt, t(m1)), cbind(m1, rbind(c(kk, pk), c(pk, pp))))

Variance-covariance matrix for parameters, inverse of information matrix

Simulation

Initialize simulation array to keep simulation results

sim <- matrix(0, 0, 11)
smn <- matrix(0, 0, 11)
spm <- matrix(0, 0, npar + 2)

Simulation for distribution of future amounts

Want 10,000 simulations, but exceeds R capacity, so do in batches of 5,000

Scatter plots of residuals & Distribution of Forecasts

## Summary From Simulation

Summary of mean, standard deviation, and 90% confidence interval from simulation, similar for one-period forecast

sumr <- matrix(0, 0, 4)
sumn <- matrix(0, 0, 4)

for (i in 1:11) {
  sumr <- rbind(sumr, c(mean(sim[, i]), sd(sim[, i]), quantile(sim[, i], c(.05, .95))))
  sumn <- rbind(sumn, c(mean(smn[, i]), sd(smn[, i]), quantile(smn[, i], c(.05, .95))))
}
sumr
##                                        5%       95%
##  [1,]         0.0        0.0         0.00         0
##  [2,]    642440.5   481028.8    -44194.36   1501201
##  [3,]   1252318.1   648086.8    234536.28   2350542
##  [4,]   3558346.6  1140237.1   1736139.47   5459264
##  [5,]   7343316.8  1685230.7   4618772.27  10135134
##  [6,]  17027739.1  2812746.4  12455686.23  21698143
##  [7,]  40274081.9  4840499.9  32377990.92  48358166
##  [8,]  74154813.9  7300843.4  62309918.99  86306505
##  [9,] 126342190.5 10639635.7 109032027.02 144116409
## [10,] 209675353.5 15905776.4 183850667.27 236026995
## [11,] 480270600.9 29296749.7 432905077.06 529557186
##                                        5%       95%
##  [1,]         0.0        0.0         0.00         0
##  [2,]    642440.5   481028.8    -44194.36   1501201
##  [3,]    527492.0   376363.7    -37645.58   1173939
##  [4,]   2230028.1   900647.6    790310.82   3758650
##  [5,]   3691095.9  1201650.8   1761276.09   5701639
##  [6,]   9561369.2  2183081.0   6055215.06  13228795
##  [7,]  20957144.4  3610684.4  15092490.22  26953460
##  [8,]  33428879.0  4931632.4  25433851.64  41583251
##  [9,]  46263909.9  6136563.8  36357479.01  56377102
## [10,]  59065904.7  7160890.5  47441520.52  70937231
## [11,] 176368263.7 12708050.9 155788022.10 197482374