# Bootstrap Using Chain Ladder Method

This example shows how to apply a chain ladder bootstrap method to generate several developmentTriangle objects to estimate the ultimate claims.

Deterministic claim estimation methods produce point estimates of reserve values with no information about the uncertainty of these estimates. The goal of a stochastic claim estimation method is to assess the variability of estimated reserve values. The chain ladder bootstrapping approach is a simulation-based method to randomly modify the developmentTriangle data and produce a distribution of estimated reserves that represents the variability of the estimated reserve values. This example is based on the work of Wüthrich and Merz [1].

OriginYear    DevelopmentYear    ReportedClaims    PaidClaims
__________    _______________    ______________    __________

2010             12               3995.7          1893.9
2010             24                 4635          3371.2
2010             36               4866.8          4079.1
2010             48               4964.1            4487
2010             60               5013.7          4711.4
2010             72               5038.8          4805.6
2010             84                 5059          4853.7
2010             96               5074.1          4877.9

### Create developmentTriangle

Create a developmentTriangle object and use claimsPlot to visualize the developmentTriangle. For more information on unpaid claims estimation, see Overview of Claims Estimation Methods for Non-Life Insurance.

dTriangle = developmentTriangle(data);
dTriangleTable = view(dTriangle);
% visualize the development triangle
claimsPlot(dTriangle)

### Analyze the developmentTriangle

The developmentTriangle link ratios are estimated using the formula:

$\stackrel{ˆ}{{\mathit{f}}_{\mathit{j}}}=\frac{\sum _{\mathit{i}=0}^{\mathit{I}-\mathit{j}-1}{\mathit{C}}_{\mathit{i},\mathit{j}+1}}{\sum _{\mathit{i}=0}^{\mathit{I}-\mathit{j}-1}{\mathit{C}}_{\mathit{i},\mathit{j}}}$

Use linkRatios to calculate the age-to-age factors.

Use linkRatioAverages to calculate the averages of the age-to-age factors.

disp(averageFactorsTable);
12-24     24-36     36-48     48-60     60-72     72-84     84-96    96-108    108-120
______    ______    ______    ______    ______    ______    _____    ______    _______

Simple Average                        1.1767    1.0563    1.0249    1.0107    1.0054    1.0038    1.003    1.002      1.001
Simple Average - Latest 5              1.172     1.056    1.0268    1.0108    1.0054    1.0038    1.003    1.002      1.001
Simple Average - Latest 3               1.17    1.0533     1.027    1.0117    1.0057    1.0037    1.003    1.002      1.001
Medial Average - Latest 5x1           1.1733    1.0567    1.0267    1.0103     1.005     1.004    1.003    1.002      1.001
Volume-weighted Average               1.1766    1.0563     1.025    1.0107    1.0054    1.0038    1.003    1.002      1.001
Volume-weighted Average - Latest 5     1.172     1.056    1.0268    1.0108    1.0054    1.0038    1.003    1.002      1.001
Volume-weighted Average - Latest 3    1.1701    1.0534     1.027    1.0117    1.0057    1.0037    1.003    1.002      1.001
Geometric Average - Latest 4            1.17     1.055    1.0267     1.011    1.0055    1.0037    1.003    1.002      1.001

Display the selected age-to-age factors table and calculate the cumulative development factor (CDF) using cdfSummary.

dTriangle.TailFactor = 1;
selectedFactorsTable = cdfSummary(dTriangle);
disp(selectedFactorsTable);
12-24      24-36      36-48      48-60      60-72     72-84      84-96     96-108     108-120    Ultimate
_______    _______    _______    _______    _______    ______    _______    _______    _______    ________

Selected                    1.1766     1.0563      1.025     1.0107     1.0054    1.0038      1.003      1.002     1.001        1
CDF to Ultimate             1.3072      1.111     1.0518     1.0261     1.0153    1.0098      1.006      1.003     1.001        1
Percent of Total Claims    0.76501    0.90008    0.95075    0.97453    0.98496    0.9903    0.99402    0.99701     0.999        1

DIsplay the latest diagonal.

latestDiagonal = dTriangle.LatestDiagonal;

Compute the projected ultimate claims using ultimateClaims.

projectedUltimateClaims = ultimateClaims(dTriangle);

Display the full development triangle using fullTriangle.

fullTriangleTable = fullTriangle(dTriangle);
disp(fullTriangleTable);
12        24        36        48        60        72        84        96       108       120      Ultimate
______    ______    ______    ______    ______    ______    ______    ______    ______    ______    ________

2010    3995.7      4635    4866.8    4964.1    5013.7    5038.8      5059    5074.1    5084.3    5089.4     5089.4
2011      3968    4682.3    4963.2    5062.5    5113.1    5138.7    5154.1    5169.6    5179.9    5185.1     5185.1
2012      4217    5060.4      5364    5508.9    5558.4    5586.2    5608.6    5625.4    5636.7    5642.3     5642.3
2013    4374.2    5205.3    5517.7    5661.1    5740.4    5780.6    5803.7    5821.1    5832.7    5838.6     5838.6
2014    4499.7    5309.6    5628.2    5785.8    5849.4    5878.7    5900.8    5918.5    5930.3    5936.3     5936.3
2015    4530.2    5300.4    5565.4    5715.7    5772.8    5804.1    5825.9    5843.4    5855.1      5861       5861
2016    4572.6    5304.2    5569.5    5714.3    5775.4    5806.7    5828.6    5846.1    5857.7    5863.6     5863.6
2017    4680.6    5523.1    5854.4    6000.9    6065.1      6098    6120.9    6139.3    6151.6    6157.7     6157.7
2018    4696.7    5495.1    5804.4    5949.6    6013.3    6045.9    6068.6    6086.8      6099    6105.1     6105.1
2019    4945.9    5819.2    6146.7    6300.5    6367.9    6402.4    6426.5    6445.8    6458.7    6465.2     6465.2

Compute the total reserves using ultimateClaims.

IBNR = ultimateClaims(dTriangle) - dTriangle.LatestDiagonal;
IBNR = array2table(IBNR, 'RowNames', dTriangleTable.Properties.RowNames, 'VariableNames', {'IBNR'});
IBNR{'Total',1} = sum(IBNR{:,:});
disp(IBNR);
IBNR
______

2010          0
2011     5.1857
2012      16.89
2013     34.886
2014     57.583
2015     88.148
2016     149.34
2017     303.29
2018     609.99
2019     1519.3
Total    2784.6

To derive the resampling approaches, the Time Series Model of the distribution-free chain ladder (CL) model is defined as:

${\mathit{C}}_{\mathit{i},\mathit{j}+1}={\mathit{f}}_{\mathit{j}}{\mathit{C}}_{\mathit{i},\mathit{j}}+{\sigma }_{\mathit{j}}\sqrt{{\mathit{C}}_{\mathit{i},\mathit{j}}}{ϵ}_{\mathit{i},\mathit{j}+1}$

For the link ratio selected above, Wüthrich [1] and Mack [2] show that the standard deviation is estimated as:

${\stackrel{ˆ}{{\sigma }_{\mathit{j}}}}^{2}=\frac{1}{\mathit{I}-\mathit{j}-1}\text{\hspace{0.17em}}\sum _{\mathit{i}=0}^{\mathit{I}-\mathit{j}-1}{\mathit{C}}_{\mathit{i},\mathit{j}}{\left(\frac{{\mathit{C}}_{\mathit{i},\mathit{j}+1}}{{\mathit{C}}_{\mathit{i},\mathit{j}}}-\stackrel{ˆ}{{\mathit{f}}_{\mathit{j}}}\right)}^{2}$

${\stackrel{ˆ}{{\sigma }_{\mathit{J}-1}}}^{2}=\mathrm{min}\left\{\frac{{\stackrel{ˆ}{{\sigma }_{\mathit{J}-2}}}^{4}}{{\stackrel{ˆ}{{\sigma }_{\mathit{J}-3}}}^{3}}\text{\hspace{0.17em}};{\stackrel{ˆ}{{\sigma }_{\mathit{J}-3}}}^{2}\text{\hspace{0.17em}};{\stackrel{ˆ}{{\sigma }_{\mathit{J}-2}}}^{2}\right\}$

estimatedStandardDeviations = currentSelectedFactors;
for i=1:width(estimatedStandardDeviations)-1
estimatedStandardDeviations(1,i) = sqrt(sum(((factorsTable{1:end-i,i} - currentSelectedFactors(:,i)).^2).*dTriangleTable{1:end-i,i}) / (height(dTriangleTable)-i-1));
end
estimatedStandardDeviations(1,end) = sqrt(min([estimatedStandardDeviations(1,end-1)^4 / estimatedStandardDeviations(1,end-2)^2, estimatedStandardDeviations(1,end-2)^2, estimatedStandardDeviations(1,end-1)^2]));

disp(estimatedStandardDeviations);
0.8667    0.3699    0.2420    0.1310    0.0673    0.0361    0.0001    0.0001    0.0001

To apply the bootstrap method, you need to find the appropriate residuals that allow for the construction of the empirical distribution $\stackrel{ˆ}{{\mathit{F}}_{\mathit{n}}}$ to construct the bootstrap observations.

Consider the following residuals for $\mathit{i}+\mathit{j}\le \mathit{I},\mathit{j}\ge 1$.

$\stackrel{\sim }{{ϵ}_{\mathit{i},\mathit{j}}}=\frac{{\mathit{F}}_{\mathit{i},\mathit{j}}-\stackrel{ˆ}{{\mathit{f}}_{\mathit{j}-1}}}{{\sigma }_{\mathit{j}-1}{\mathit{C}}_{\mathit{i},\mathit{j}-1}^{-1/2}}$ where ${\mathit{F}}_{\mathit{i},\mathit{j}}=\frac{{\mathit{C}}_{\mathit{i},\mathit{j}}}{{\mathit{C}}_{\mathit{i},\mathit{j}-1}}$

Following Wüthrich [1], you can scale the residuals to adjust their variance upwards. Unscaled residuals tend to result in lighter tails in the simulated distribution.

Adjust the residuals such that the bootstrap distribution has an adjusted variance function.

${\mathit{Z}}_{\mathit{i},\mathit{j}}={\left(1-\frac{{\mathit{C}}_{\mathit{i},\mathit{j}-1}}{\sum _{\mathit{i}=0}^{\mathit{I}-\mathit{j}}{\mathit{C}}_{\mathit{i},\mathit{j}-1}}\right)}^{-\frac{1}{2}}\frac{{\mathit{F}}_{\mathit{i},\mathit{j}}-\stackrel{ˆ}{{\mathit{f}}_{\mathit{j}-1}}}{\stackrel{ˆ}{{\sigma }_{\mathit{j}-1}}{\mathit{C}}_{\mathit{i},\mathit{j}-1}^{-\frac{1}{2}}}$

You can apply the bootstrap algorithm using three different versions:

• Efron's nonparametric bootstrap for residuals $\stackrel{\sim }{{ϵ}_{\mathit{i},\mathit{j}}}$

• Efron's nonparametric bootstrap for scaled residuals ${\mathit{Z}}_{\mathit{i},\mathit{j}}$

• Parametric bootstrap under the assumption that the residuals have a standard Gaussian distribution, that is ${\mathit{Z}}_{\mathit{i},\mathit{j}}^{*}$ is resampled from $\mathit{N}\left(0,1\right)$

This example uses the second version (Efron's nonparametric bootstrap for scaled residuals) to calculate ${\mathit{Z}}_{\mathit{i},\mathit{j}}$.

% Create a copy of the factors table and modify it to create the
% residuals table
residuals = factorsTable.Variables;

colSums = sum(dTriangle.Claims,'omitnan');
for i=1:height(residuals)
for j=1:width(residuals)
residuals(i,j) = (1 - (dTriangleTable{i,j}/colSums(j)))^-0.5 * (factorsTable{i,j} - currentSelectedFactors(1,j)) / (estimatedStandardDeviations(1,j)*(dTriangleTable{i,j}^-0.5));
end
end

The residuals $\left\{{\mathit{Z}}_{\mathit{i},\mathit{j}},\mathit{i}+\mathit{j}\le \mathit{I}\right\}$ define a bootstrap distribution.

residualsVector = residuals(:);
residualsVector(isnan(residualsVector)) = [];
histogram(residualsVector,10)
title('Scaled Residuals')
xlabel('Residual Value')
ylabel('Frequency')

To simulate a new reserves scenario with the bootstrap method, follow these steps.

#### Step 1: Resample a triangle of residuals from the bootstrap distribution.

Resample the independent and identically distributed (i.i.d.) residuals$\left\{{{\mathit{Z}}^{*}}_{\mathit{i},\mathit{j}},\mathit{i}+\mathit{j}\le \mathit{I}\right\}$ from the bootstrap distribution.

resampledResiduals = residuals;

rng('default');
rng(1);

for i = 1:height(residuals)-1
for j = 1:width(residuals)-i+1
resampledResiduals(i,j) = datasample(residuals(~isnan(residuals)), 1);
end
end

disp(resampledResiduals);
-1.5522   -0.5120   -1.2668    0.7776   -1.3649    0.2799   -0.5495   -1.3146   -1.5364
-0.4041   -1.5522   -0.4784   -1.2189   -0.7591    0.2610   -0.4784   -1.5522       NaN
-0.4091   -1.3649   -0.5495   -1.6767   -0.8571   -1.3143   -0.4879       NaN       NaN
-0.7591    1.3226    1.0791    0.2610    0.2861   -0.7591       NaN       NaN       NaN
0.2799   -1.5522   -0.8571    0.3243   -0.4879       NaN       NaN       NaN       NaN
-1.3143   -0.4784    0.5556   -1.2668       NaN       NaN       NaN       NaN       NaN
1.9550         0    1.9550       NaN       NaN       NaN       NaN       NaN       NaN
0.7693    0.5169       NaN       NaN       NaN       NaN       NaN       NaN       NaN
0.2799       NaN       NaN       NaN       NaN       NaN       NaN       NaN       NaN
NaN       NaN       NaN       NaN       NaN       NaN       NaN       NaN       NaN

#### Step 2: Compute bootstrapped claims.

Define ${\mathit{C}}_{\mathit{i},0}^{*}={\mathit{C}}_{\mathit{i},0}$ and, for $\mathit{j}\ge 1$, assume that:

${\mathit{C}}_{\mathit{i},\mathit{j}}^{*}=\stackrel{ˆ}{{\mathit{f}}_{\mathit{j}-1}}{\mathit{C}}_{\mathit{i},\mathit{j}-1}^{*}+\stackrel{ˆ}{{\sigma }_{\mathit{j}-1}}\sqrt{{\mathit{C}}_{\mathit{i},\mathit{j}-1}^{*}}{\mathit{Z}}_{\mathit{i},\mathit{j}}^{*}$

This expression represents the new simulated claim values. Using the simulated claim values, you can create a new developmentTriangle to estimate new reserve values.

bootstrappedClaims = dTriangleTable.Variables;

for j = 2:width(bootstrappedClaims)
bootstrappedClaims(:,j) = currentSelectedFactors(1,j-1).*bootstrappedClaims(:,j-1) + estimatedStandardDeviations(1,j-1).*sqrt(bootstrappedClaims(:,j-1)).*resampledResiduals(:,j-1);
end

stackedClaims = reshape(bootstrappedClaims',100,1);
stackedClaims = stackedClaims(~isnan(stackedClaims));
newData = data;
newData.values = stackedClaims;
bootstrappedDevelopmentTriangle = developmentTriangle(newData,'Claims','values');

#### Step 3: Select a link ratio consistent with the model.

The volume-weighted average is the link ratio that is consistent with the model used in this bootstrap approach.

bootstrappedDevelopmentTriangle.TailFactor = 1;
bootstrappedSelectedFactorsTable = cdfSummary(bootstrappedDevelopmentTriangle);
disp(bootstrappedSelectedFactorsTable);
12-24      24-36     36-48      48-60      60-72     72-84      84-96     96-108     108-120    Ultimate
_______    _______    ______    _______    _______    ______    _______    _______    _______    ________

Selected                    1.1751      1.054    1.0253     1.0099     1.0048    1.0036      1.003      1.002     1.001        1
CDF to Ultimate              1.301     1.1072    1.0504     1.0245     1.0145    1.0096      1.006      1.003     1.001        1
Percent of Total Claims    0.76861    0.90321     0.952    0.97609    0.98572    0.9905    0.99403    0.99701     0.999        1

Use fullTriangle to display the full development triangle corresponding to the selected link ratio.

bootstrappedFullTriangle = fullTriangle(bootstrappedDevelopmentTriangle);
disp(bootstrappedFullTriangle);
12        24        36        48        60        72        84        96       108       120      Ultimate
______    ______    ______    ______    ______    ______    ______    ______    ______    ______    ________

2010    3995.7    4616.2    4863.2    4963.4    5023.7    5044.5    5064.1    5079.3    5089.5    5094.6     5094.6
2011      3968    4646.6      4869    4982.8    5024.8    5048.4    5068.1    5083.3    5093.4    5098.5     5098.5
2012      4217    4938.6    5181.1    5301.1    5341.9    5366.6    5383.3    5399.5    5410.2    5415.6     5415.6
2013    4374.2    5103.1    5425.3    5580.2    5642.5    5674.5    5693.8    5710.9    5722.3      5728       5728
2014    4499.7    5310.5    5567.5    5691.3    5755.4    5784.2    5804.8    5822.2    5833.8    5839.6     5839.6
2015    4530.2    5253.5    5536.3    5684.8    5733.2      5761    5781.5    5798.8    5810.4    5816.2     5816.2
2016    4572.6    5494.6    5803.9    5985.1    6044.2    6073.5    6095.1    6113.4    6125.6    6131.7     6131.7
2017    4680.6    5552.6    5879.4    6028.2    6087.7    6117.2      6139    6157.4    6169.7    6175.9     6175.9
2018    4696.7    5542.6      5842    5989.8    6048.9    6078.2    6099.9    6118.2    6130.4    6136.5     6136.5
2019    4945.9      5812      6126      6281      6343    6373.7    6396.4    6415.6    6428.4    6434.8     6434.8

#### Step 4: Compute the total reserves.

Compute the total reserves from the simulated developmentTriangle.

bootstrappedDevelopmentTriangleTable = view(bootstrappedDevelopmentTriangle);
bootstrappedIBNR = ultimateClaims(bootstrappedDevelopmentTriangle) - bootstrappedDevelopmentTriangle.LatestDiagonal;
bootstrappedIBNR = array2table(bootstrappedIBNR, 'RowNames', bootstrappedDevelopmentTriangleTable.Properties.RowNames, 'VariableNames', {'IBNR'});
bootstrappedIBNR{'Total',1} = sum(bootstrappedIBNR{:,:});
disp(bootstrappedIBNR);
IBNR
______

2010          0
2011     5.0881
2012     16.188
2013     34.197
2014     55.485
2015     83.048
2016     146.61
2017     296.45
2018     593.94
2019       1489
Total      2720

You can repeat the previous steps many times to genreate a full, simulated, distribution of reserves. The simulation produces reserves for each year and for the total reserves.

### Simulate Multiple Bootstrapped Scenarios

Create 1000 bootstrapped development triangles and calculate the incurred-but-not-reported (IBNR) for each developmentTriangle.

n = 1000;

simulatedIBNR = zeros(10,n);
for i = 1:n
simulatedResiduals = residuals;

for j = 1:height(residuals)-1
for k = 1:width(residuals)-j+1
simulatedResiduals(j,k) = datasample(residuals(~isnan(residuals)),1);
end
end

simulatedClaims = dTriangleTable.Variables;

for j = 2:width(simulatedClaims)
simulatedClaims(:,j) = currentSelectedFactors(1,j-1).*simulatedClaims(:,j-1) + estimatedStandardDeviations(1,j-1).*sqrt(simulatedClaims(:,j-1)).*simulatedResiduals(:,j-1);
end

simulatedClaims = reshape(simulatedClaims',100,1);
simulatedClaims = simulatedClaims(~isnan(simulatedClaims));
simulatedData = data;
simulatedData.ReportedClaims = simulatedClaims;
simulatedDevelopmentTriangle = developmentTriangle(simulatedData);

simulatedDevelopmentTriangle.TailFactor = 1;
simulatedLatestDiagonal = simulatedDevelopmentTriangle.LatestDiagonal;
simulatedProjectedUltimateClaims = ultimateClaims(simulatedDevelopmentTriangle);

simulatedIBNR(:,i) = simulatedProjectedUltimateClaims - simulatedLatestDiagonal;

end

simulatedIBNR(end+1,:) = sum(simulatedIBNR);

Select a year to plot the distribution of the IBNR, calculate the mean, and compare that mean to a calculated deterministic value.

originYear =5;
histogram(simulatedIBNR(originYear+1,:));
hold on;
plot(mean(simulatedIBNR(originYear+1,:)),0,'O','LineWidth',2)
plot(IBNR{originYear+1,1},0,'X','LineWidth',2);
legend('Simulated IBNR',['Simulated mean : ' num2str(round(mean(simulatedIBNR(originYear+1,:)),2))],['Deterministic IBNR : ' num2str(round(IBNR{originYear+1,1},2))]);
hold off;

Plot a histogram of the totals for IBNRs, simulated means, and deterministic values.

histogram(simulatedIBNR(11,:));
hold on;
plot(mean(simulatedIBNR(11,:)),0,'O','LineWidth',2)
plot(IBNR{11,1},0,'X','LineWidth',2);
legend('Simulated Total IBNR',['Simulated mean : ' num2str(round(mean(simulatedIBNR(11,:)),2))],['Deterministic Total IBNR : ' num2str(round(IBNR{11,1},2))]);
hold off;

### References

1. Wüthrich, Mario, and Michael Merz. Stochastic Claims Reserving Methods in Insurance. Hoboken, NJ: Wiley, 2008

2. Mack, Thomas. "Distribution-Free Calculation of the Standard Error of Chain Ladder Reserve Estimates." Astin Bulletin. Vol. 23, No. 2, 1993.