2.1 Inverse distance weighting

A first possible approach to model combination is to intuitively build weights based on the distance between an expert’s projection about a variable of interest, or vector of variables, and the true value of this variable. This idea can be achieved through Inverse Distance Weighting (IDW). The advantage of this approach is its intuitiveness and ease of use.

Classically, IDW was used with Euclidean distance. In a geometric context, Shepard (Reference Shepard1968) used IDW to consider distance while taking angles into account. In a probabilistic setting, the challenge with this method is then to determine an appropriate distance measure. One such measure is the Wasserstein distance, which Kantorovitch & Rubinštein (Reference Kantorovitch and Rubinštein1958) first realised was applicable to probability distributions. This idea was expanded on by Givens & Shortt (Reference Givens and Shortt1984) and used recently by Pesenti et al. (Reference Pesenti, Bettini, Millossovich and Tsanakas2021) for sensitivity analysis. In the univariate case, the distance for expert e over time period $\vec{\tau}$ at location x is defined as

\begin{equation*} D^{(e)} = \int \Big|F_{Y_{\vec{\tau},x}^{(e)}}(y)-F_{Y_{\vec{\tau},x}^{\phantom{(e)}}}(y)\Big|dy, \end{equation*}

with $F_{Y_{\vec{\tau},x}}$ the real cumulative distribution function and $F_{Y_{\vec{\tau},x}^{(e)}}$ the expert’s CDF.

With this distance, the weight attributed to each expert’s projection is then

\begin{equation*} w_e = \frac{1/D^{(e)}}{\sum_{l=1}^{n}1/D^{(l)}}. \end{equation*}

2.2 Non-parametric calibration

From a literature-based approach, model combination can be approached from many angles. Cooke (Reference Cooke1991) offered a review of early expert combination methods. Clemen & Winkler (Reference Clemen and Winkler1999) further elaborated on this review, suggesting issues that need to be considered when combining expert opinions such as expert selection and the role of interaction between experts. Since then, Cooke & Goossens (Reference Cooke and Goossens2008) and Hammitt & Zhang (Reference Hammitt and Zhang2012) compared the performance of multiple combination methods, among which a classical approach which was first presented by Cooke (Reference Cooke1991).

This combination method uses desirable properties of scoring rules, namely that they should be coherent, strictly proper, and relevant (see Cooke, Reference Cooke1991 for details). A three-part method was established attributing weights to each expert distribution based on a relative information component, a calibration component, and an entropy component. This method has the advantage of being non-parametric, suggesting that an expert does not need to have a complete statistical model. Such a method can be appropriate for example in actuarial science, where an expert might reasonably provide an estimate for a small, medium, and large loss, but not a full loss distribution.

From the calibration variables $V_1$ to $V_K$ defined previously, we set $v_{k,0}$ and $v_{k,M+1}$ such that

\begin{equation*} v_{k,0} < v_{k,m}^{(e)} < v_{k,M+1} \text{ } \forall \text{ } m,e. \end{equation*}

We compare these selections and expert-provided values with the true observed values to find the proportion of calibration variables in each interquantile space. This forms an empirical distribution $\boldsymbol{{q}}=\{q_1, \ldots, q_{M+1}\}$ that we can compare to the theoretical proportion $\boldsymbol{{p}}=\{p_1, \ldots, p_{M+1}\}$ . As shown by Cooke (Reference Cooke1991), we can obtain the calibration and entropy components, C(e) and O(e) respectively, as

\begin{align*} C(e) = 1-\chi^2_{K-1} ((2K)I(q,p)),\end{align*}

where

\begin{align*} I(q,p) &= \sum_{k=1}^{M+1} q_k \ln\!\left(\frac{q_k}{p_k}\right)\end{align*}

is the relative information component, and

\begin{align*} O(e) = \frac{1}{K}\sum_{k=1}^{K} \left(\ln(v_{k,M+1}-v_{k,0})+\sum_{m=1}^{M+1} p_m \ln\!\left(\frac{p_m}{v_{k,m}^{(e)}-v_{k,m-1}^{(e)}}\right)\right). \end{align*}

It can readily be shown that the relative information component I(q, p) multiplied by 2K (i.e. twice the number of calibration variables) follows a chi-squared distribution. The calibration component uses this fact to measure the goodness of fit of each expert forecast, while the entropy component measures the distance of expert forecasts from a uniform distribution. The intuition for this component is that a uniform model provides very little useful information. From these, we finally obtain

\begin{equation*} w^{\prime}_e = C(e)O(e)I_{\{C(e)>\alpha\}} \end{equation*}

for a specified threshold $\alpha$ chosen by optimising the score of the combined distributions, where $0 < \alpha < 1$ . This $\alpha$ can be seen as a hyperparameter representing the minimal calibration level that each model needs to satisfy to receive weight. As such, a higher $\alpha$ means we give probability to less models. This also implies that the maximal value for $\alpha$ is the highest value of C(e). We can then recalibrate the weights to make their sum equal to 1 by dividing $w^{\prime}_e$ by the sum over all experts:

\begin{equation*} w_e= \frac{w^{\prime}_e}{\sum_{l=1}^n w^{\prime}_l}. \end{equation*}

These $w_e$ do not require the analyst to have a prior opinion of each expert’s projections. We will refer to this method as Cooke’s method for the sake of brevity. In the context of daily precipitation annual maxima, the corresponding calibration variable is then $\vec{Y}_{\vec{\tau},x}$ , where we consider K different sites x.

2.3 Bayesian model averaging

As an alternative to the previous approaches, one may seek to exploit their prior knowledge using Bayesian methods, updating a prior belief with observed data to obtain a posterior distribution more representative of recent data.

Bayesian Model Averaging (BMA) is a widely used tool for model combination. Recently, in the United States, BMA was used to study extreme rainfall density as well as daily mean rainfall by Zhu et al. (Reference Zhu, Forsee, Schumer and Gautam2013) and Massoud et al. (Reference Massoud, Lee, Gibson, Loikith and Waliser2020), respectively. First made popular by Raftery et al. (Reference Raftery, Madigan and Hoeting1997) in linear models, BMA uses observed data to update weights to different models based on their likeliness. This relies on the premise that any of the models could be right, but selecting only one model would fail to capture the uncertainty around this choice. This in turn leads to reducing overconfidence from ignoring a model’s uncertainty. BMA however implicitly relies on the assumption that one of the models must be right (Hoeting et al., Reference Hoeting, Madigan, Raftery and Volinsky1999). Note that the method presented in Cooke (Reference Cooke1991) relies on a similar assumption, given that the optimal $\alpha$ requires at least one model to be chosen.

Let $\mathcal{M}$ be a discrete variable representing this best model, with possible values $\{\mathcal{M}_1, \ldots, \mathcal{M}_n\}$ . An analyst has some prior belief about the probability that each expert’s model is right, $\Pr(\mathcal{M} = \mathcal{M}_e)$ , which we will denote $\Pr(\mathcal{M}_e)$ , normalised such that $\sum_{e=1}^n \Pr(\mathcal{M}_e) = 1$ . In the absence of prior information, then $\Pr(\mathcal{M}_e)=1/n, \,\forall \, e$ . Given data $\vec{y}_{\vec{\tau},x}$ , the analyst can update these probabilities through Bayesian updating, that is

\begin{equation*} \Pr\!\left(\mathcal{M}_e|\vec{y}_{\vec{\tau},x}\right) = \frac{\Pr\!\left(\vec{y}_{\vec{\tau},x}|\mathcal{M}_e\right)\Pr(\mathcal{M}_e)}{\sum_{l=1}^{n}\Pr(\vec{y}_{\vec{\tau},x}|\mathcal{M}_l)\Pr(\mathcal{M}_l)}, \end{equation*}

where $\Pr\!\left(\vec{y}_{\vec{\tau},x}|\mathcal{M}_e\right)$ is the probability of observing $\vec{y}_{\vec{\tau},x}$ under model $\mathcal{M}_e$ . Since we divide by $\sum_{l=1}^{n}\Pr\!\left(\vec{y}_{\vec{\tau},x}|\mathcal{M}_l\right)\Pr(\mathcal{M}_l)$ , it follows that $\sum_{e=1}^n \Pr\!\left(\mathcal{M}_e|\vec{y}_{\vec{\tau},x}\right) = 1$ , and posterior probabilities $\Pr\!\left(\mathcal{M}_e|\vec{y}_{\vec{\tau},x}\right)$ can therefore be considered as updated weights attributed to each expert. This supposes that all models are independent since we ignore possible interactions between models. This assumption is appropriate in this case since all experts rely on different approaches, but this will be discussed in section 4. There are different ways of calculating the expert-associated probabilities.

A first possibility is to use an expectation-maximisation (EM) algorithm, as shown by Darbandsari & Coulibaly (Reference Darbandsari and Coulibaly2019), where the residuals between the model projections $\vec{y}_{\vec{\tau},x}^{\,(e)}$ , representing an expert’s projection generated from model $\mathcal{M}_e$ about the variable $\vec{Y}_{\vec{\tau},x}$ , and actual data are assumed to follow a Gaussian distribution. This assumption allows for iterating through these residuals’ Gaussian likelihood while updating the weights attributed to each expert model until the difference between iterations is less than some threshold $\beta$ . The algorithm is outlined in Appendix A. The algorithm allows for projecting a posterior distribution for a period $\vec{\psi}=\{s^{\prime},s^{\prime}+1, \ldots, t^{\prime}\}$ , with $s^{\prime}\in \{t,t+1,\ldots,t^{\prime}\}$ , $t<t^{\prime} \le T$ . This approach must be used carefully as it can lead to overfitting. With a low threshold, EM will be optimised for training data, but will also learn the noise surrounding the signal. Because of this, the algorithm can then perform poorly on testing data. This limitation of the EM algorithm will be further explored in section 3.

The same hypothesis that residuals follow a normal distribution was used by Zhu et al. (Reference Zhu, Forsee, Schumer and Gautam2013), but with a different approach due to limited datasets, where the authors used bootstrapping under Generalised Likelihood Uncertainty Estimation (GLUE, see Beven & Freer, Reference Beven and Freer2001) to obtain the posterior likelihoods. The algorithm is presented in Algorithm 1, where $y_{\vec{\tau},x,m}$ is the $m{\text{th}}$ quantile of the vector $\vec{y}_{\vec{\tau},x}$ , $y_{\vec{\tau},x,m,b}$ is the $b{\text{th}}$ bootstrap resampling of this quantile with B resamplings, and $\Pr(Y_{\vec{\psi},x}=y|\mathcal{M}_e)$ is the probability distribution of extreme precipitation under model $\mathcal{M}_e$ for a future period $\vec{\psi}$ .

3.1 Non-equiprobable pooling

In the previous section, we saw different methods to attribute weights to expert opinions depending on the probability of each expert projection being accurate. We can incorporate these ideas into the pooling idea of Innocenti et al. (Reference Innocenti, Mailhot, Leduc, Cannon and Frigon2019). We use their pooling method as a baseline, where one may consider all expert-provided models as equally likely, which we will refer to as the equiprobable scenario. Instead of supposing that all model projections are equally likely ( $\Pr(\mathcal{M}_e)=1/n$ ), we can update our belief about the probability of each model with observed data. By defining

\begin{equation*} \Pr\!\left(Y_{\vec{\psi},x}=y\right) = \frac{\Pr\!\left(\mathcal{M}_e|\vec{y}_{\vec{\tau},x}\right)}{t^{\prime}-s^{\prime}+1}, \end{equation*}

with $y\in \vec{y}_{\vec{\psi},x}$ , we obtain a shifted distribution reflecting this updated belief, where $t^{\prime}-s^{\prime}+1$ is the number of years in the future projection period $\vec{\psi}$ , and $\vec{\tau}$ is the historical observed period.

3.2 Calculating areal reduction factors

We can now incorporate the model combination methods and pooling presented previously into ARFs to investigate their impact on extreme precipitation quantile and ARF projections.

Although there are slightly varying definitions of ARFs, we will focus on the one used by Le et al. (Reference Le, Davison, Engelke, Leonard and Westra2018), which can be thought of as a quantile of average areal precipitation over an average of point precipitation quantiles. This particular definition has the advantage of being applicable to any station within a region and not only one station. Starting from the notation introduced in section 2, let $Y_{\vec{\tau},A,R}$ and $Y_{\vec{\tau},x,R,}$ respectively, represent the areal and point rainfall for area A, point x, and frequency R over period $\vec{\tau}$ . The ARF based on daily precipitation is then

\begin{equation*} \text{ARF}_{(A,R,\vec{\tau})} = \frac{Y_{\vec{\tau},A,R}}{\frac{1}{\text{card}(X)}\sum_{x\in X} Y_{\vec{\tau},x,R}}, \end{equation*}

where there are a collection of sites $x\in X$ within area A.

An issue that arises when calculating ARFs with climate models is that expert projections are often not available at each point x, but rather at a grid scale. This issue can however be solved by assuming that scaling from point precipitation to grid average precipitation is time invariant. Li et al. (Reference Li, Sharma, Johnson and Evans2015) demonstrated the validity of this hypothesis, enabling the use of grid cells for ARF calculation, where we would have grid-to-area instead of point-to-area.

With this notion of time-constant scaling, we can thus consider the points x as grid cell coordinates instead of stations. This enables us to calculate ARFs using grid data, as made available by climate agencies such as Climate Data Canada and Copernicus Climate Change Service. Grid cells are available at a resolution of approximately 0.1 degrees of latitude and longitude and represent average precipitation over the grid cell. We consider zones of $6\times 4$ grid cells in the regions of Montreal and Quebec. We have access to 24 different climate models using historical data from 1951 to 2005 to project precipitation from 2006 to 2100. These models rely on three different Representative Concentration Pathways (RCP) emission scenarios: a low emissions scenario (RCP 2.6), a moderate emissions scenario (RCP 4.5) and a high emissions scenario (RCP 8.5). In keeping with Innocenti et al. (Reference Innocenti, Mailhot, Leduc, Cannon and Frigon2019), we will focus on the 8.5 scenario, corresponding to a 4.5 degree increase by 2100. We calibrate weights using data from 2001 to 2020, for which we have both real and projected precipitation. This allows us to compare quantiles for Bayesian Model Averaging, or interquantile space for Cooke’s method and inverse Distance Weighting, and so calibrate combination weights using each method. With the obtained weights, all future time periods are then forecasted. It is worth noting that this relies on the hypothesis that weights remain the same whether forecasting near or far future.

To use pooling, we need to have sufficient data for frequency analysis. Due to having 24 models instead of the 50 in Innocenti et al. (Reference Innocenti, Mailhot, Leduc, Cannon and Frigon2019), we consider 6-year periods, such as precipitation from 2016 to 2021, rather than 3-year periods to obtain a similar number of data points. Applying weights calculated using the different methods presented in section 2, we calculate shifted densities reflecting these adjusted weights, as can be observed in Figures 4 and 5. However, before using the BMA-EM algorithm, a threshold or number of iterations must be chosen to prevent overfitting. This is because too many iterations of the expectation-maximisation algorithm will lead to learning the signal as well as the noise in the training data. Figure 1 illustrates the average MSE resulting from splitting data from 2001 to 2020 into ten-year training and testing periods. Overfitting occurs passed 4 iterations of the expectation-maximisation algorithm, where we see that the testing sample MSE starts increasing significantly while the training sample MSE stabilises and even slightly increases. To prevent this overfitting, we choose to stop the algorithm after 4 iterations. This is a known issue of BMA (see e.g. Domingos, Reference Domingos2000), added to BMA tending to select only one model asymptotically, as BMA implicitly relies on the assumption that one of the models is true (Hoeting et al., Reference Hoeting, Madigan, Raftery and Volinsky1999). An $\alpha$ of 0.65 is also selected for Cooke’s method by optimising the error as shown in Figure 1.

Figure 1 Grid cell MSE of the expectation-maximisation algorithm (left) and Cooke’s method (right) in the Montreal region from 2001 to 2020.

We first note that different combination methods can yield very different weights attributed to each model. Figure 2 illustrates the difference in weights for the cities of Montreal and Quebec for a period from 2001 to 2020. Note that for the rest of the article, when we refer to Quebec, this will imply Quebec City and not the province. We see that for Montreal, the two BMA methods generally agree, whereas they do not for Quebec. On the other hand, both Cooke and IDW lead to relatively similar weights in both locations, but they differ from BMA results. These different weight attributions can lead to different projected quantiles.

Figure 2 Model weight by method for Montreal (left) and Quebec (right) for precipitation from 2001 to 2020.

One may seek to investigate the expert models with larger probability to ensure they agree with those models’ hypotheses. For example, in Montreal, the next to last model (MPI_MR) receives a large weight from the EM algorithm, but gets truncated by the calibration approach. This happens because the model has a jump in precipitation level around the $50{\text{th}}$ quantile, as illustrated in Figure 3. 7 observations out of 20 fall in the 45–50% interquantile space for model MPI_MR. This causes a poor fit in calibration in terms of Cooke’s method, but the quantile-to-quantile residuals are quite small, meaning that we still have a good fit in terms of low residuals when compared to real data, making the BMA methods give this model high weight. In similar fashion, one can gain additional insight by comparing the outputs of different combination approaches.

Figure 3 Cumulative distribution for model MPI_MR and real data in Montreal for a grid cell between 2001 and 2020.

Figures 4 and 5 illustrate the upper tail of the resulting empirical cumulative distribution functions under different possible combination methods for Montreal and Quebec, respectively. We see that the quantiles obtained from varying combination methods are substantially different depending on the weights attributed to each model. From a risk management perspective, such differences can alter conclusions reached by an analyst concerning risk level. As such, one would benefit from considering multiple combination methods, given that this would allow for better understanding of projection uncertainty.

Figure 4 Upper tail of empirical cumulative distribution functions of pooled annual maximum daily rainfall (mm) for Montreal from 2016 to 2021 with different weighting methods.

Figure 5 Upper tail of empirical cumulative distribution functions of pooled annual maximum daily rainfall (mm) for Quebec from 2016 to 2021 with different weighting methods.

Since different combination methods yield different results, one may be interested in the variability induced by attributing weights to each expert. Let $F_{\vec{\tau},x}^{(e)}$ be the cumulative distribution function corresponding to model $\mathcal{M}_e$ . We define the CDF of $Y_{\vec{\tau},A}$ as

\begin{equation*} F_{\vec{\tau},A}(y) = w_1 F_{\vec{\tau},A}^{(1)}(y) + \ldots + w_n F_{\vec{\tau},A}^{(n)}(y), \end{equation*}

where $w_1, \ldots, w_n$ are the weights attributed to each expert (which correspond to probabilities $\Pr(\mathcal{M}_e|\vec{y}_{\vec{\tau},A})$ ). It is easy to show that for a given return level, the boundaries for $Y_{\vec{\tau},A,R}$ will be the minimum and maximum of $\left\{Y_{\vec{\tau},A,R}^{(1)}, \ldots, Y_{\vec{\tau},A,R}^{(n)}\right\}$ . Indeed, we have

\begin{align*} Y_{\vec{\tau},A,R}&=\min\!\left(Y_{\vec{\tau},A}:\Pr(Y\le Y_{\vec{\tau},A}) \ge 1-1/R\right)\\ &=\min\!\left(Y_{\vec{\tau},A}:F_{\vec{\tau},A}\left(Y_{\vec{\tau},A}\right) \ge 1-1/R\right)\\ &=\min\!\left(Y_{\vec{\tau},A}:w_1 F_{\vec{\tau},A}^{(1)}\left(Y_{\vec{\tau},A}\right) + \ldots + w_n F_{\vec{\tau},A}^{(n)}\left(Y_{\vec{\tau},A}\right) \ge 1-1/R\right). \end{align*}

Now suppose $F_{\vec{\tau},A}^{(i)}\left(Y_{\vec{\tau},A}\right) \ge F_{\vec{\tau},A}^{(j)}\left(Y_{\vec{\tau},A}\right)$ for some $i \in \{1, \ldots, n\}$ and $\forall \, j \in \{1, \ldots, n\}$ . Then, it follows that $F_{\vec{\tau},A}^{(i)}\left(Y_{\vec{\tau},A}\right) \ge w_1 F_{\vec{\tau},A}^{(1)}\left(Y_{\vec{\tau},A}\right) + \ldots + w_n F_{\vec{\tau},A}^{(n)}\left(Y_{\vec{\tau},A}\right) \ge 1-1/R$ , provided that $w_1, \ldots, w_n \in [0,1]$ with $\sum w_i = 1$ , and so $F^{-1}_{\vec{\tau},A}\left(Y_{\vec{\tau},A}\right)$ must be the minimum for $Y_{\vec{\tau},A,R}$ for any combination of weights. Similarly, the reverse logic allows for stating that the lowest CDF must yield the maximum quantile.

From this reasoning, Figure 6 illustrates the CDF obtained with each combination method in Montreal between 2001 and 2020 compared to the minimum and maximum boundaries of quantiles, where the period is expanded to 20 years to allow for empirical quantiles from each expert in a short enough period that precipitation is not expected to change significantly. We notice that the combination methods are grouped within a much narrower range than the theoretical boundaries from the minimum and maximum projections.

Figure 6 Upper tail of empirical cumulative distribution functions of pooled annual maximum daily rainfall (mm) for Montreal from 2001 to 2020 with different weighting methods, and minimum and maximum boundaries.

We can suppose that the weights provided by the different combination methods will improve the variance around a quantile estimate compared to having no information about each expert. While we cannot obtain this variance mathematically, we can use bootstrap resampling to compare the quantile distribution under each scenario. Figure 7 illustrates the resulting density distributions for the $95{\text{th}}$ quantile in Montreal between 2001 and 2020. In keeping with intervals presented in Climate Data Canada, the 10% and 90% quantiles of the distribution supposing no information about experts are shown (corresponding to the equiprobable scenario), which can be thought of as the lower and upper bounds that a user with no evaluation of the expert models might consider as plausible. We notice that the two BMA methods differ largely from the other two methods, with modes lying outside the 10–90% boundaries, while the other methods are more similar to not evaluating experts, particularly for the $95{\text{th}}$ quantile.

Figure 7 Comparison of bootstrap densities under different combination methods for the 90th quantile (left) and 95th quantile (right) in Montreal between 2001 and 2020 for 10,000 iterations.

This difference is driven by the same phenomenon as the difference in weight attribution. BMA methods rely on the assumption that residuals between projections and real data follow a normal distribution, whereas Cooke’s method and IDW using Wasserstein distance use the distance between (cumulative) densities of the projections and real data. If expert distributions have jumps in their CDFs, this will cause aggregation for both Cooke and Wasserstein, leading to these models receiving little weight. Nonetheless, the residuals between these experts’ projections and real data might still be small, such that BMA methods will attribute larger weight to these models. These different weights cause the gap between quantile values of BMA methods compared to the other methods, as observed in Figure 6. Given the similarity in results between the non-parametric methods using densities, and the BMA methods using residuals as in Figure 7, it is natural to suppose that keeping only one method using densities and another using residuals provides sufficient information for analysis purposes.

Moreover, these combination methods allow for alternate confidence bounds based on an evaluation of expert models as opposed to supposing all expert projections are equally likely. Table 1 also highlights the reduction in variance for the $95{\text{th}}$ quantile in Montreal, while the much lower variance is similar for all methods in Quebec.

Table 1. Comparison of mean and variance of uniform weight attribution and model combination weights for Montreal and Quebec from 2001 to 2020 at the $95{\text{th}}$ quantile.

Applying the same exercise to multiple grid cells within the Montreal region, we can calculate the resulting ARF for each method. Given that we observe a 10% difference in $95{\text{th}}$ quantiles between methods, we can expect different weights to yield significantly different ARF curves.

It is worth noting that directly using quantiles found with model combination methods can yield non-sensical results when computing ARFs. This is because the spatiotemporal relation between the full region $Y_{\vec{\tau},A,R}$ and the underlying grid cells $Y_{\vec{\tau},x,R}$ for each expert’s projection is broken when comparing a weighted average of $y_{\vec{\tau},x,R}^{(e)}$ and $y_{\vec{\tau},A,R}^{(e)}$ , leading to ARFs possibly exceeding 1. From a point-to-area point of view, this would not make sense, seeing as a whole area cannot have more intense precipitation than its maximal component, limiting the applicability of such a method. This effect is lessened by using the same weights for all grid cells within an area.

From the significant variability in higher quantiles observed in the previous figures depending on the weights attributed to model projections, we choose to study percentage changes in ARF and quantiles because they yield more comparable information between the different combination methods than actual quantile and ARF values. Mathematically, the modelled quantile change for area A corresponds to $\Delta_{quant}=Y_{\vec{\psi},A,R}/Y_{\vec{\tau},A,R}$ , and the ARF change to $\Delta_{ARF} = ARF_{A,R,\vec{\psi}}/ARF_{A,R,\vec{\tau}}$ for future period $\vec{\psi}$ and current period $\vec{\tau}$ .

Using the quantile boundaries found previously, we can establish boundaries for possible quantile change by comparing the future maximum to the current minimum and vice versa for the minimum possible change. This exercise is not well-defined for ARFs, since the area value depends on the underlying grid cells, and so we cannot for example use the highest area quantile with the lowest grid quantiles, as this would not make sense from a rainfall perspective. Keeping the same 20-year period, we compare it to a near-future period of 2011–2030 and a far future of 2071–2090. The idea behind comparing two future periods is that the variability in near future should be lower than for a later period. Figures 8 and 9 show the change in quantiles and ARFs in Montreal for the near future and far future at a 1 in 20-year return level. While we observe the expected change in variability for quantiles, Figure 9 shows that change in ARF does not significantly vary between near and far projections. This could be explained by looking at the underlying composition of the ARF, where

\begin{align*} \Delta_{ARF} &= ARF_{A,R,\vec{\psi}}/ARF_{A,R,\vec{\tau}}= \frac{\left(\frac{Y_{\vec{\psi},A,R}}{\frac{1}{\text{card}(X)} \sum_{x\in X} Y_{\vec{\psi},x,R} }\right)}{\left(\frac{Y_{\vec{\tau},A,R}}{\frac{1}{\text{card}(X)} \sum_{x\in X} Y_{\vec{\tau},x,R} }\right)}\\ &= \frac{Y_{\vec{\psi},A,R}}{Y_{\vec{\tau},A,R}} \frac{\sum_{x\in X} Y_{\vec{\tau},x,R}}{\sum_{x\in X} Y_{\vec{\psi},x,R}} = \Delta_{quant}\frac{\sum_{x\in X} Y_{\vec{\tau},x,R}}{\sum_{x\in X} Y_{\vec{\psi},x,R}}, \end{align*}

such that the first ratio is the change in quantiles, but the second ratio has the current period and future period inverted, suggesting that it will be approximately inversely proportional to the quantile change. As such, the two ratios will cancel out, other than the random noise between different grid cell precipitation, which is what we observe in Figure 9. The fatter tails for Bayesian methods are induced by the distribution of quantile change, which is less centred around a mode, as seen in Figure 8.

Figure 8 Distribution of projected quantile change at a 1 in 20-year return level in Montreal between 2001–2020 and 2011–2030 (left) or 2071–2090 (right).

Figure 9 Distribution of projected ARF change at a 1 in 20-year return level in Montreal between 2001–2020 and 2011–2030 (left) or 2071–2090 (right).

The same idea is applied to Quebec in Appendix B, where all methods generally agree, and the Bayesian quantile change projections are more centred around their mode than for Montreal, such that the ARF change projection has smaller tails. The distributions resulting from different combination methods can provide valuable information about the uncertainty of projections, where for example in this case the confidence level is higher regarding Quebec projections than Montreal projections. Moreover, compared to the 10–90% confidence bounds usually presented, we see that the resulting distributions from combination methods provide alternate bounds based on an evaluation of expert projections. In an actuarial context, this could be very important as it can highlight whether a projection is too conservative or not conservative enough.

Figures 10 and 11 compare the mean percentage change in ARF and quantiles respectively for a 1 in 20-year return level for Montreal between Cooke’s method and BMA-EM, divided into approximately 24 km $\times$ 22 km areas. These two methods are chosen to illustrate the substantial variation between a density-based method and a residuals-based method. For example, both methods project increases in quantile, but one projects a 10% increase with little change to the ARF, while the other projects a 22% increase with a reduction to the ARF. From a risk management perspective, this would imply differing scenarios of a moderate increase with similar spatial distribution and a heavier increase with more localised precipitation, which can lead to different losses (see for example Cheng et al., Reference Cheng, Li, Li and Auld2012 and American Academy of Actuaries, 2020).

Figure 10 Percentage change in quantiles for a 1 in 20-year return level between 2001–2020 and 2071–2090 for the region of Montreal using Cooke’s method (left) and BMA-EM (right).

Figure 11 Percentage change in ARFs for a 1 in 20-year return level between 2001–2020 and 2071–2090 for the region of Montreal using Cooke’s method (left) and BMA-EM (right).

Flood losses provide a particular example of how Cooke’s method and Bayesian model averaging with expectation-maximisation would lead to different loss projections. While the link between extreme rainfall and flooding is complex, the difference in scenarios between Cooke’s method and BMA-EM allows for a theoretical discussion of its impact for an actuary. Through a combination of hydrological and hydraulic models such as Hydrotel (Fortin et al., Reference Fortin, Turcotte, Massicotte, Moussa, Fitzback and Villeneuve2001), HEC-RAS (Brunner, Reference Brunner2016) or the Hillslope Link Model (Demir & Krajewski, Reference Demir and Krajewski2013), one can produce discharge flood projections under different rainfall scenarios. Breinl et al. (Reference Breinl, Lun, Müller-Thomy and Blöschl2021) used elasticity to illustrate the relationship between extreme precipitation and flooding, where depending on ground dampness, an increase in precipitation will have an at least equivalent increase in river discharge, leading to increased flood severity. Supposing that the reduction in ARF will mitigate the impact of an increase in quantiles due to more localised rainfall, such that for example we have an approximately 7% and 19% increase under, respectively, the Cooke and BMA-EM scenarios, the relationship between discharge and rainfall would clearly imply a greater risk of increased flood losses in the latter case.

Using a hierarchical model such as the one used by Boudreault et al. (Reference Boudreault, Grenier, Pigeon, Potvin and Turcotte2020), flood intensities are associated with different levels of discharge, and their respective probabilities are established from frequency analysis. In their study, the second and third levels of flood intensities have discharges of 1,570 m $^3/$ s and of 1,740 m $^3/$ s, respectively, with occurrence probabilities of 0.01496 and 0.00842. This $10.8\%$ difference in discharge is lower than the projected increase in extreme precipitation using BMA-EM, which is not the case for Cooke’s method. All else being equal, the probability of observing more severe flooding in the BMA-EM scenario would increase relative to the Cooke scenario. This change in probability can then be used to calculate premiums and/or reserves for flooding, where BMA-EM would lead to a more conservative estimate than the other method in this case. In a changing climate perspective, the range of scenarios resulting from different combination methods becomes even more important to have a fuller understanding of the impact of climate change on insurable losses. An analyst using only one method would fail to obtain a complete picture of projection uncertainty and may find themselves being overconfident in the result of a single combination method.

Similar graphics are available in Appendix C for Quebec. Projections for this city are much more similar across methods, leading to smaller confidence intervals in this case.

In summary, we see that the different combination methods considered can yield varying sets of weights, or probabilities, assigned to each model, which impacts projected quantiles. From the similarities between methods using densities compared to methods using residuals, we see that one only needs to use one method from each approach to obtain a picture of the underlying projection uncertainty, and the difference between the approaches provides a measure of this uncertainty. In cases where methods agree, one could more confidently reach conclusions about the analysed data, but in cases where methods disagree, using only one method would fail to capture projection uncertainty. Moreover, combination methods can yield alternate confidence bounds based on an evaluation of expert models and offer an improved pooling projection over considering all expert projections as equally likely.