Multi-population mortality modelling and forecasting: a hierarchical credibility regression approach

Abstract

This paper proposes a multi-level hierarchical credibility regression approach to model multi-population mortality data. Future mortality rates are derived using extrapolation techniques, while the forecasting performances between the proposed model, the original Lee–Carter model and two Lee–Carter extensions for multiple populations are compared for both genders of three northern European countries with small populations (Ireland, Norway, Finland). Empirical illustrations show that the proposed method produces more accurate forecasts than the Lee–Carter model and its multi-population extensions.

Appendices

Appendix

A Lee–Carter mortality modelling for multiple populations: a review of methods

We review the original Lee–Carter [35] model and its most popular and widely applied extensions for multiple populations, the joint-k model presented by Carter and Lee [11] and the augmented common factor model proposed by Li and Lee [38].

1.1 The Lee–Carter model

In its original form, the Lee–Carter model links the natural logarithm of the observed $\log m_{t,xgc}$ mortality rates with the following model predictor

$$\begin{aligned} \log m_{t,xgc} = \alpha _{xgc}^{(1)} + \alpha _{xgc}^{(2)} \kappa _{t,gc} + \epsilon _{t,xgc}, \end{aligned}$$

(62)

where $\alpha _{xgc}^{(1)}$ parameter reflects the average mortality at age x, gender g and country c, $\kappa _{t,gc}$ parameter indicates the general level of mortality for gender g and country c in year t and $\alpha _{xgc}^{(2)}$ parameter shows the deviation from the average mortality at age x, gender g and country c, as the general level of mortality changes. The errors $\epsilon _{t,xgc}$ are expected to be normally distributed, with zero mean and constant variance, reflecting effects not captured by the model. Thus, after assuming that errors are independent and homoscedastic with zero mean, Lee and Carter [35] suggested a close approximation to the SVD (Singular Value Decomposition) method, under the constraints $\textstyle \sum _{x = x_1}^{x_{k}} \alpha _{xgc}^{(2)} = 1$ and $\textstyle \sum _{t = t_1}^{t_{n}} \kappa _{t,gc} = 0$, to obtain the following parameter estimates

$$\begin{aligned} {{\widehat{\alpha }}}_{xgc}^{(1)}= & {} \frac{1}{t_{n} - t_{1} + 1} \sum _{t=t_1}^{t_{n}} \log m_{t,xgc}, \ \ {{\widehat{\kappa }}}_{t,gc} = \sum _{x=x_1}^{x_{k}}\left[ \log m_{t,xgc} - {{\widehat{\alpha }}}_{xgc}^{(1)}\right] ,\\ {{\widehat{\alpha }}}_{xgc}^{(2)}= & {} \frac{\sum _{t=t_1}^{t_{n}} \left[ \log m_{t,xgc} - {{\widehat{\alpha }}}_{xgc}^{(1)}\right] \widehat{\kappa }_{t,gc}}{\sum _{t=t_1}^{t_{n}} {{\widehat{\kappa }}}_{t,gc}^2} \ . \end{aligned}$$

The period estimates ${{\widehat{\kappa }}}_{t,gc}$ are extrapolated using time series methods. Lee and Carter [35] suggested a random walk with a drift parameter ${{\widehat{\theta }}}_{gc}$ to project period parameter for $h = 1, 2, \ldots , H$ years ahead, according to ${\widehat{\kappa }}_{t_{n} + h,gc} = {\widehat{\kappa }}_{t_{n},gc} + {{\widehat{\theta }}}_{gc} h$. Then, projected values of $\kappa _{t,gc}$ are utilized along with the estimates of age parameters $\alpha _x^{(1)}$ and $\alpha _x^{(2)}$ to obtain the following mortality forecasts

$$\begin{aligned} \log {\widehat{m}}_{t_{n} + h,xgc} = {{\widehat{\alpha }}}_{xgc}^{(1)} + {{\widehat{\alpha }}}_{xgc}^{(2)} {{\widehat{\kappa }}}_{t_{n} + h,gc} = \log {\widehat{m}}_{t_{n},xgc} + \left( {{\widehat{\alpha }}}_{xgc}^{(2)} \widehat{\theta }_{gc}\right) h. \end{aligned}$$

(63)

1.2 The joint-K model

Carter and Lee [11] proposed a model structure for multiple populations, in which mortality is jointly driven by a single period parameter $K_{t}$ as follows

$$\begin{aligned} \log m_{t,xgc} = \alpha _{xgc}^{(1)} + \alpha _{xgc}^{(2)} K_{t} + \epsilon _{t,xgc}, \end{aligned}$$

(64)

where $\alpha _{xgc}^{(1)}$ and $\alpha _{xgc}^{(2)}$ are defined as in the original Lee–Carter model, but now, $K_{t}$ is the common period parameter of the mortality level in year t for gender g and country c. Errors $\epsilon _{t,xgc}$ are assumed independent and identically distributed. Then, constraints $\textstyle \sum _{c = 1}^{C} \sum _{g = 1}^{G} \sum _{x = x_{1}}^{x_{k}} \alpha _{xgc}^{(2)} = 1$ and $\textstyle \sum _{t = t_{1}}^{t_{n}} K_{t} = 0$ lead to the following parameter estimates

$$\begin{aligned} {{\widehat{\alpha }}}_{xgc}^{(1)} = \frac{1}{t_{n} - t_{1} + 1} \sum _{t=t_{1}}^{t_{n}} \log m_{t,xgc}, \ {{\widehat{K}}}_{t} = \sum _{c = 1}^{C} \sum _{g = 1}^{G} \sum _{x=x_1}^{x_{k}} \left[ \log \widehat{m}_{t,xgc} - {{\widehat{\alpha }}}_{xgc}^{(1)}\right] \end{aligned}$$

and

$$\begin{aligned} {{\widehat{\alpha }}}_{xgc}^{(2)} = \frac{\sum _{t=t_{1}}^{t_{n}} \left[ \log m_{t,xgc} - {{\widehat{\alpha }}}_{xgc}^{(1)}\right] \widehat{K}_{t}}{\sum _{t=t_{1}}^{t_{n}} {{\widehat{K}}}_{t}^2} \ . \end{aligned}$$

The common period parameter ${{\widehat{K}}}_{t}$ is assumed to follow a random walk with a drift $\theta$, ${{\widehat{K}}}_{t} = {{\widehat{K}}}_{t-1} + \theta + \varepsilon _{t}$, where the time trend errors $\varepsilon _{t}$ are again assumed to be independent and identically distributed, and independent of the model errors $\epsilon _{t,xgc}$. The drift parameter is estimated by $\textstyle {{\widehat{\theta }}} = \frac{1}{t_{n}-t_{1}}\sum _{t=t_{1}+1}^{t_{n}}\left( {{\widehat{K}}}_{t} - {{\widehat{K}}}_{t-1}\right) = \frac{{{\widehat{K}}}_{t_{n}} - \widehat{K}_{t_{1}}}{n-1}$ and then it used to project period estimates $\textstyle {{\widehat{K}}}_{t_{n} + h} = {{\widehat{K}}}_{t_{n}} + \widehat{\theta }h$. Thus, projected mortality rates for $h = 1, 2, \ldots$ years ahead, for age x, gender g and country c are given by

$$\begin{aligned} \log {\widehat{m}}_{t_{n} + h,xgc} = {{\widehat{\alpha }}}_{xgc}^{(1)} + {{\widehat{\alpha }}}_{xgc}^{(2)} {{\widehat{K}}}_{t_{n} + h} = \log {\widehat{m}}_{t_{n},xgc} + \left( {{\widehat{\alpha }}}_{xgc}^{(2)} \widehat{\theta }\right) h \ . \end{aligned}$$

(65)

1.3 The augmented common factor model

To avoid long-run divergence in mean mortality forecasts for multiple countries, Li and Lee [38] modified the original Lee–Carter model by setting a common age parameter $\alpha _{x}^{(2)}$ and the same period parameter $K_{t}$ for all populations as follows

$$\begin{aligned} \log m_{t,xgc} = \alpha _{xgc}^{(1)} + \alpha _{x}^{(2)} K_{t} + \epsilon _{t,xgc}. \end{aligned}$$

(66)

Again, we use two constraints $\textstyle \sum _{c = 1}^{C} \sum _{g = 1}^{G} \sum _{x = x_{1}}^{x_{k}} w_{gc} \alpha _{x}^{(2)} = 1, \ \text {and} \ \sum _{t = t_{1}}^{t_{n}} K_{t} = 0$, where $w_{gc}$ is the weight for gender g in country c, set to be proportional to the total number of populations, i.e., $w_{gc} = 1/(G C)$. The model parameters are estimated by

$$\begin{aligned} {{\widehat{\alpha }}}_{xgc}^{(1)} = \frac{1}{t_{n} - t_{1} + 1} \sum _{t=t_{1}}^{t_{n}} \log m_{t,xgc}, \ {{\widehat{K}}}_{t} = \sum _{c = 1}^{C} \sum _{g = 1}^{G} \sum _{x=x_1}^{x_{k}} \ w_{gc} \left[ \log m_{t,xgc} - {{\widehat{\alpha }}}_{xgc}^{(1)} \right] \end{aligned}$$

and

$$\begin{aligned} {{\widehat{\alpha }}}_{x}^{(2)} = \frac{\sum _{c = 1}^{C} \sum _{g = 1}^{G} \sum _{t=t_{1}}^{t_{n}} \ w_{gc} \left[ \log m_{t,xgc} - \widehat{\alpha }_{xgc}^{(1)} \right] {{\widehat{K}}}_{t}}{\sum _{t=t_{1}}^{t_{n}} {{\widehat{K}}}_{t}^2} \ . \end{aligned}$$

To include the individual differences in the trends, Li and Lee [38] suggested an additional factor $\alpha _{xgc}^{(3)} \kappa _{t,gc}$ to form the augmented common factor model

$$\begin{aligned} \log m_{t,xgc} = \alpha _{xgc}^{(1)} + \alpha _{x}^{(2)} \ K_{t} + \alpha _{xgc}^{(3)} \ \kappa _{t,gc} + \epsilon _{t,xgc} \ . \end{aligned}$$

Assuming the extra constraint $\textstyle \sum _{x=x_{1}}^{x_{k}} \alpha _{xgc}^{(3)} = 1$, the additional parameters are estimated as

$$\begin{aligned} {{\widehat{\kappa }}}_{t,gc} = \sum _{x=x_1}^{x_{k}} \left[ \log m_{t,xgc} - {{\widehat{\alpha }}}_{xgc}^{(1)} - {{\widehat{\alpha }}}_{x}^{(2)} \widehat{K}_{t}\right] \end{aligned}$$

and

$$\begin{aligned} {{\widehat{\alpha }}}_{xgc}^{(3)} = \frac{\sum _{t=t_{1}}^{t_{n}} \left[ \log m_{t,xgc} - {{\widehat{\alpha }}}_{xgc}^{(1)} - \widehat{\alpha }_{x}^{(2)} {{\widehat{K}}}_{t}\right] \widehat{\kappa }_{t,gc}}{\sum _{t=t_{1}}^{t_{n}} {{\widehat{\kappa }}}_{t,gc}^{'2}}. \end{aligned}$$

Both period parameters ${{\widehat{K}}}_{t}$ and $\widehat{\kappa }_{t,gc}$ follow a random walk model with a drift, given by ${{\widehat{K}}}_{t} = {{\widehat{K}}}_{t-1} + \theta + \varepsilon _{t}$ and ${{\widehat{\kappa }}}_{t,gc} = {{\widehat{\kappa }}}_{t-1,gc} + \theta _{gc} + \varepsilon _{t,gc}$, respectively, where time trend errors $\varepsilon _{t}$ and $\varepsilon _{t,gc}$ are assumed to be independent and identically distributed, and independent of the model error $\epsilon _{t,xgc}$. The drift parameters can be estimated by

$$\begin{aligned} {{\widehat{\theta }}} =\ \frac{1}{t_{n}-t_{1}}\sum _{t=t_{1}+1}^{t_{n}}\left( {{\widehat{K}}}_{t} - {{\widehat{K}}}_{t-1}\right) = \frac{{{\widehat{K}}}_{t_{n}} - \widehat{K}_{t_{1}}}{n-1} \end{aligned}$$

and

$$\begin{aligned} {{\widehat{\theta }}}_{gc} =\ \frac{1}{t_{n}-t_{1}}\sum _{t=t_{1}+1}^{t_{n}}\left( {{\widehat{\kappa }}}_{t,gc} - {{\widehat{\kappa }}}_{t-1,gc}\right) = \frac{{{\widehat{\kappa }}}_{t_{n},gc} - {{\widehat{\kappa }}}_{t_{1},gc}}{n-1} \ . \end{aligned}$$

Then, projected mortality rates for $h = 1, 2, \ldots$ years ahead, for age x, gender g and country c are obtained by

$$\begin{aligned} \log {\widehat{m}}_{t_{n} + h,xgc}= & {} {{\widehat{\alpha }}}_{xgc} + {{\widehat{\alpha }}}_{x}^{(2)} {{\widehat{K}}}_{t_{n + h}} + \widehat{\alpha }_{xgc}^{(3)} {{\widehat{\kappa }}}_{t_{n + h},gc} \nonumber \\= & {} \log {\widehat{m}}_{t_{n},xgc} + \left( {{\widehat{\alpha }}}_{x}^{(2)} {{\widehat{\theta }}} + {{\widehat{\alpha }}}_{xgc}^{(3)} {{\widehat{\theta }}}_{gc}\right) h \ . \end{aligned}$$

(67)

We can easily observe that expressions (63), (65) and (67) are linear functions of the forecasting horizon h, where their intercept equal to the fitted rates of the last observed year and their slope is the product of the estimated age parameters with the drift terms.

B Proof of Theorem 1

Expressions (18)–(22) are notations. Expression (23) can be easily proved using assumptions (A6) and (4)

$$\begin{aligned} E\left[ {\widehat{{{\varvec{\beta }}}}}_{xgc}\bigg |\varTheta _{xgc}^{+123}\right]&= \left( {\varvec{Z}}_{xgc}^{'}{\varvec{W}}_{xgc}{\varvec{Z}}_{xgc}\right) ^{-1}{\varvec{Z}}_{xgc}^{'}{\varvec{W}}_{xgc} E\left( {\varvec{Y}}_{xgc}\bigg |\varTheta _{xgc}^{+123}\right) \\&=\left( {\varvec{Z}}_{xgc}^{'}{\varvec{W}}_{xgc}{\varvec{Z}}_{xgc}\right) ^{-1}{\varvec{Z}}_{xgc}^{'}{\varvec{W}}_{xgc}{\varvec{Z}}_{xgc} {{\varvec{\beta }}}\left( \varTheta _{xgc}^{+123}\right) \\&={{\varvec{\beta }}}\left( \varTheta _{xgc}^{+123}\right) . \end{aligned}$$

From (18), we can easily get (24)

$$\begin{aligned} E\left[ {\widehat{{{\varvec{\beta }}}}}_{xgc}\bigg |\varTheta _{gc}^{+23}\right] = E\left[ E\left[ {\widehat{{{\varvec{\beta }}}}}_{xgc}\bigg | \varTheta _{xgc}^{+123} \right] \bigg | \varTheta _{gc}^{+23} \right] = E\left[ {{\varvec{\beta }}}\left( \varTheta _{xgc}^{+123}\right) \bigg | \varTheta _{gc}^{+23}\right] = {{\varvec{\beta }}}\left( \varTheta _{gc}^{+23}\right) , \end{aligned}$$

while from (20) we can prove (25)

$$\begin{aligned} E\left[ {\widehat{{{\varvec{\beta }}}}}_{xgc}\bigg | \varTheta _{c}^{+3}\right] = E\left[ E\left[ {\widehat{{{\varvec{\beta }}}}}_{xgc}\bigg | \varTheta _{gc}^{+23} \right] \bigg | \varTheta _{c}^{+3} \right] = E\left[ {{\varvec{\beta }}}\left( \varTheta _{gc}^{+23}\right) \bigg | \varTheta _{c}^{+3}\right] \ = \ {{\varvec{\beta }}}\left( \varTheta _{c}^{+3}\right) . \end{aligned}$$

The first part of (26) can be proved by (22) and (25)

$$\begin{aligned} E\left( \widehat{\varvec{\beta }}_{xgc}\right) = E\left[ E\left( \widehat{\varvec{\beta }}_{xgc}\bigg | \varTheta _{c}^{+3}\right) \right] = E\left[ \varvec{\beta }\left( \varTheta _{c}^{+3}\right) \right] \ = \ \varvec{\beta }. \end{aligned}$$

Also, (6) gives

$$\begin{aligned} E\left( \widehat{\varvec{\beta }}_{gc}\right)&= \left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}\right) ^{-1} \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}E\left( \widehat{\varvec{\beta }}_{xgc}\right) = \varvec{\beta } \end{aligned}$$

and via (7), we can get

$$\begin{aligned} E\left( \widehat{\varvec{\beta }}_{c}\right)&= \left( \sum _{g=1}^{Gc}{\varvec{K}}_{gc}\right) ^{-1} \sum _{g=1}^{Gc}{\varvec{K}}_{gc}E\left( \widehat{\varvec{\beta }}_{gc}\right) \ = \ \varvec{\beta }. \end{aligned}$$

$\square$

C Proof of Theorem 2

Expression (27) is proved as follows

$$\begin{aligned}&{\text{ Cov }}\left[ \varvec{\beta }\left( \varTheta _{xgc}^{+123}\right) ,\widehat{\varvec{\beta }}_{x'g'c'}\right] = E\left\{ {\text{ Cov }}\left[ \varvec{\beta }\left( \varTheta _{xgc}^{+123}\right) , \widehat{\varvec{\beta }}_{x'g'c'}\bigg |\varTheta _{c}^{+3}\right] \right\} \\&\qquad + {\text{ Cov }}\left\{ E\left[ \varvec{\beta }\left( \varTheta _{xgc}^{+123}\right) \bigg |\varTheta _{c}^{+3}\right] ,E\left[ \widehat{\varvec{\beta }}_{x'g'c'}\bigg |\varTheta _{c}^{+3}\right] \right\} \\&\quad = \delta _{cc'} \delta _{gg'} E \left\{ {\text{ Cov }}\left[ \varvec{\beta }\left( \varTheta _{gc}^{+23}\right) \bigg |\varTheta _{c}^{+3}\right] \right\} \\&\qquad + \ E \left\{ E \left\{ {\text{ Cov }}\left\{ E\left[ \varvec{\beta }\left( \varTheta _{xgc}^{+123}\right) \bigg |\varTheta _{xgc}^{+123}\right] ,E\left[ \widehat{\varvec{\beta }}_{x'g'c'}\bigg |\varTheta _{xgc}^{+123}\right] \big |\varTheta _{gc}^{+23}\right\} \big |\varTheta _{c}^{+3}\right\} \right\} \\&\qquad + \ E \left\{ E \left\{ E \left\{ {\text{ Cov }}\left[ \varvec{\beta }\left( \varTheta _{xgc}^{+123}\right) , \widehat{\varvec{\beta }}_{x'g'c'}\bigg |\varTheta _{xgc}^{+123}\right] \bigg |\varTheta _{gc}^{+23}\right\} \bigg |\varTheta _{c}^{+3}\right\} \right\} + \ \delta _{cc'} {\text{ Cov }}\left[ \varvec{\beta }\left( \varTheta _{c}^{+3}\right) \right] \\&\quad = \delta _{cc'} \delta _{gg'} E \left\{ {\text{ Cov }}\left[ \varvec{\beta }\left( \varTheta _{gc}^{+23}\right) \bigg |\varTheta _{c}^{+3}\right] \right\} \\&\qquad + \delta _{cc'} \delta _{gg'} \delta _{xx'} E \left\{ E \left\{ {\text{ Cov }}\left[ \varvec{\beta }\left( \varTheta _{xgc}^{+123}\right) \bigg |\varTheta _{gc}^{+23}\right] \big |\varTheta _{c}^{+3} \right\} \right\} \\&\qquad + \ \delta _{cc'} {\text{ Cov }}\left[ \varvec{\beta }\left( \varTheta _{c}^{+3}\right) \right] \\&\quad = \delta _{cc'} \left[ \delta _{gg'} \left( \delta _{xx'} {\varvec{A}} + {\varvec{U}}\right) + \varvec{\varPsi }\right] . \end{aligned}$$

Similarly, we can obtain (28) and (29) as

$$\begin{aligned} {\text{ Cov }}\left[ \varvec{\beta }\left( \varTheta _{gc}^{+23}\right) ,\widehat{\varvec{\beta }}_{x'g'c'}\right]&= E \left\{ {\text{ Cov }}\left\{ E\left[ \varvec{\beta }\left( \varTheta _{gc}^{+23}\right) \bigg |\varTheta _{gc}^{+23}\right] , E\left[ \widehat{\varvec{\beta }}_{x'g'c'}\bigg |\varTheta _{gc}^{+23}\right] \bigg |\varTheta _{c}^{+3}\right\} \right\} \\&\qquad + \ E \left\{ E \left\{ {\text{ Cov }}\left[ \varvec{\beta }\left( \varTheta _{gc}^{+23}\right) ,\widehat{\varvec{\beta }}_{x'g'c'}\bigg |\varTheta _{gc}^{+23}\right] \bigg |\varTheta _{c}^{+3}\right\} \right\} \\&\qquad + \ \delta _{cc'} {\text{ Cov }}\left[ \varvec{\beta }\left( \varTheta _{c}^{+3}\right) \right] \\&= \delta _{cc'} \delta _{gg'} E \left\{ {\text{ Cov }}\left[ \varvec{\beta }\left( \varTheta _{gc}^{+23}\right) \bigg |\varTheta _{c}^{+3}\right] \right\} + \ \delta _{cc'} {\text{ Cov }}\left[ \varvec{\beta }\left( \varTheta _{c}^{+3}\right) \right] \\&= \delta _{cc'} \left( \delta _{gg'} {\varvec{U}} + \varvec{\varPsi }\right) \\ {\text{ Cov }}\left[ \varvec{\beta }\left( \varTheta _{c}^{+3}\right) ,\widehat{\varvec{\beta }}_{x'g'c'}\right] &= E\left\{ {\text{ Cov }}\left[ \varvec{\beta }\left( \varTheta _{c}^{+3}\right) , \widehat{\varvec{\beta }}_{x'g'c'}\bigg |\varTheta _{c}^{+3}\right] \right\} \\&\qquad + \ {\text{ Cov }}\left\{ E\left[ \varvec{\beta }\left( \varTheta _{c}^{+3}\right) \bigg |\varTheta _{c}^{+3}\right] ,E\left[ \widehat{\varvec{\beta }}_{x'g'c'} \bigg |\varTheta _{c}^{+3}\right] \right\} \\&= \delta _{cc'} {\text{ Cov }}\left[ \varvec{\beta }\left( \varTheta _{c}^{+3}\right) \right] \\&= \delta _{cc'} \varvec{\varPsi }. \end{aligned}$$

Expression (30) can be proved by

$$\begin{aligned}&{\text{ Cov }}\left( \widehat{\varvec{\beta }}_{xgc}, \widehat{\varvec{\beta }}_{x'g'c'}\right) \\&\quad = E\left[ {\text{ Cov }}\left( \widehat{\varvec{\beta }}_{xgc}, \widehat{\varvec{\beta }}_{x'g'c'}\bigg |\varTheta _{c}^{+3}\right) \right] + {\text{ Cov }}\left[ E\left( \widehat{\varvec{\beta }}_{xgc}\bigg |\varTheta _{c}^{+3}\right) ,E\left( \widehat{\varvec{\beta }}_{x'g'c'}\bigg |\varTheta _{c}^{+3}\right) \right] \\&\quad = \delta _{cc'} \delta _{gg'} E \left\{ {\text{ Cov }}\left[ \varvec{\beta }\left( \varTheta _{gc}^{+23}\right) \bigg |\varTheta _{c}^{+3}\right] \right\} \\&\qquad + \ E \left\{ E \left\{ {\text{ Cov }}\left[ E\left( \widehat{\varvec{\beta }}_{xgc}\bigg |\varTheta _{xgc}^{+123}\right) ,E\left( \widehat{\varvec{\beta }}_{x'g'c'}\bigg |\varTheta _{xgc}^{+123}\right) \bigg |\varTheta _{gc}^{+23} \right] \bigg |\varTheta _{c}^{+3}\right\} \right\} \\&\qquad + \ E \left\{ E \left\{ E \left[ {\text{ Cov }}\left( \widehat{\varvec{\beta }}_{xgc}, \widehat{\varvec{\beta }}_{x'g'c'}\bigg |\varTheta _{xgc}^{+123}\right) \bigg |\varTheta _{gc}^{+23} \right] \bigg |\varTheta _{c}^{+3}\right\} \right\} \\&\qquad + \delta _{cc'} \ {\text{ Cov }}\left[ \varvec{\beta }\left( \varTheta _{c}^{+3}\right) \right] \\&\quad = \delta _{cc'} \delta _{gg'} \ {\varvec{U}} + \ \delta _{cc'} \delta _{gg'} \delta _{xx'} \ {\varvec{A}} \\&\qquad + \ E \left\{ E \left\{ E \left[ \sigma ^2\left( \varTheta _{xgc}^{+123}\right) \left( {\varvec{Z}}_{xgc}^{'}{\varvec{W}}_{xgc}{\varvec{Z}}_{xgc}\right) ^{-1}\bigg |\varTheta _{gc}^{+23} \right] \bigg |\varTheta _{c}^{+3}\right\} \right\} \\&\qquad + \ \delta _{cc'} \ {\text{ Cov }}\left[ \varvec{\beta }\left( \varTheta _{c}^{+3}\right) \right] \\&\quad =\delta _{cc'} \left\{ \delta _{gg'}\left[ \delta _{xx'}\left( {\varvec{A}}+s^{2}\left( {\varvec{Z}}_{xgc}^{'}{\varvec{W}}_{xgc} {\varvec{Z}}_{xgc}\right) ^{-1}\right) +{\varvec{U}}\right] + \varvec{\varPsi }\right\} . \end{aligned}$$

For expression (31) we have

$$\begin{aligned}&{\text{ Cov }}\left( \widehat{\varvec{\beta }}_{gc}, \widehat{\varvec{\beta }}_{x'g'c'}\right) =\left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}\right) ^{-1}\left[ \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}{\text{ Cov }}\left( \widehat{\varvec{\beta }}_{xgc},\widehat{\varvec{\beta }}_{x'g'c'}\right) \right] \\&\quad =\left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}\right) ^{-1}\left\{ \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc} \ \delta _{cc'} \left\{ \delta _{gg'}\left[ \delta _{xx'}\right. \right. \right. \\&\quad \left. \left. \left. \left( {\varvec{A}}+s^{2}\left( {\varvec{Z}}_{xg}^{'}{\varvec{W}}_{xg} {\varvec{Z}}_{xg}\right) ^{-1}\right) +{\varvec{U}}\right] + \varvec{\varPsi }\right\} \right\} \\&\quad =\left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}\right) ^{-1} \left[ \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc} \ \delta _{cc'} \delta _{gg'} \delta _{xx'}({\varvec{A}}+s^{2}({\varvec{Z}}_{xg}^{'}{\varvec{W}}_{xg} {\varvec{Z}}_{xg})^{-1})\right] \\&\qquad + \ \left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}\right) ^{-1} \left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc} \ \delta _{cc'} \delta _{gg'} {\varvec{U}} \right) \\&\quad + \left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}\right) ^{-1} \left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc} \ \delta _{cc'} \varvec{\varPsi } \right) \\&\quad =\left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}\right) ^{-1} \left[ \delta _{cc'} \delta _{gg'} \ {\varvec{A}} \ \left[ {\varvec{A}}+ s^{2} \left( {\varvec{Z}}_{xgc}^{'}{\varvec{W}}_{xgc}{\varvec{Z}}_{xgc}\right) ^{-1} \right] ^{-1}\right. \\&\quad \left. \left( {\varvec{A}}+s^{2}\left( {\varvec{Z}}_{xg}^{'}{\varvec{W}}_{xg} {\varvec{Z}}_{xg}\right) ^{-1}\right) \right] \\&\qquad + \left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}\right) ^{-1} \left[ \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc} \left( \delta _{cc'} \delta _{gg'} {\varvec{U}} + \delta _{cc'} \varvec{\varPsi }\right) \right] \\&\quad =\delta _{cc'} \left\{ \delta _{gg'} \left[ \left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}\right) ^{-1} \ {\varvec{A}} + \ {\varvec{U}}\right] + \varvec{\varPsi }\right\} . \end{aligned}$$

Similarly, we can prove (32). Expression (33) is obtained as

$$\begin{aligned} {\text{ Cov }}\left( \widehat{\varvec{\beta }}_{gc}, \widehat{\varvec{\beta }}_{g'c'}\right)&=\left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}\right) ^{-1}\left[ \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}{\text{ Cov }}\left( \widehat{\varvec{\beta }}_{xgc},\widehat{\varvec{\beta }}_{g'c'}\right) \right] \\&=\delta _{cc'} \left\{ \delta _{gg'} \left[ {\varvec{A}} \ \left( \sum _{x'=x_1}^{x_{k_{gc}}}{\varvec{K}}'_{x'g'c'}\right) ^{-1} + {\varvec{U}}\right] + \varvec{\varPsi }\right\} \end{aligned}$$

and (34) can be proved in the same way. From expression (31), we can prove (35) as

$$\begin{aligned}&{\text{ Cov }}\left( \widehat{\varvec{\beta }}_{c},\widehat{\varvec{\beta }}_{x'g'c'}\right) = \left( \sum _{g=1}^{G_c}{\varvec{K}}_{gc}\right) ^{-1}\sum _{g=1}^{G_c}{\varvec{K}}_{gc} {\text{ Cov }}\left( \widehat{\varvec{\beta }}_{gc},\widehat{\varvec{\beta }}_{x'g'c'}\right) \\&\quad =\left( \sum _{g=1}^{G_c}{\varvec{K}}_{gc}\right) ^{-1} \left\{ {\varvec{U}}\left[ {\varvec{U}} + \left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}\right) ^{-1}{\varvec{A}}\right] ^{-1} \ \delta _{cc'} \left[ \left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}\right) ^{-1} \ {\varvec{A}} + {\varvec{U}} \right] \right\} \\&\qquad + \ \left( \sum _{g=1}^{G_c}{\varvec{K}}_{gc}\right) ^{-1} \delta _{cc'} \left( \sum _{g=1}^{G_c}{\varvec{K}}_{gc} \varvec{\varPsi } \right) \\&\quad = \delta _{cc'} \left[ \left( \sum _{g=1}^{G_{c}}{\varvec{K}}_{gc}\right) ^{-1}{\varvec{U}} + \varvec{\varPsi }\right] \end{aligned}$$

and similarly, we can get (36). Expression (37) is obtained by (33) as follows

$$\begin{aligned}&{\text{ Cov }}\left( \widehat{\varvec{\beta }}_{c}, \widehat{\varvec{\beta }}_{g'c'}\right) =\left( \sum _{g=1}^{G_c}{\varvec{K}}_{gc}\right) ^{-1}\sum _{g=1}^{G_c}{\varvec{K}}_{gc} {\text{ Cov }}\left( \widehat{\varvec{\beta }}_{gc},\widehat{\varvec{\beta }}_{g'c'}\right) \\&\quad =\left( \sum _{g=1}^{G_c}{\varvec{K}}_{gc}\right) ^{-1} {\varvec{U}}\left[ {\varvec{U}} + \left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}\right) ^{-1}{\varvec{A}}\right] ^{-1} \delta _{cc'} \left[ \left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}\right) ^{-1} \ {\varvec{A}} + \ {\varvec{U}}\right] \\&\qquad + \ \left( \sum _{g=1}^{G_c}{\varvec{K}}_{gc}\right) ^{-1} \sum _{g=1}^{G_c}{\varvec{K}}_{gc}\delta _{cc'} \varvec{\varPsi } \\&\quad =\delta _{cc'} \left[ \left( \sum _{g=1}^{G_c}{\varvec{K}}_{gc}\right) ^{-1} {\varvec{U}} + \varvec{\varPsi }\right] , \end{aligned}$$

while (38) yields similarly. Expression (39) can be proved via (38)

$$\begin{aligned} {\text{ Cov }}\left( \widehat{\varvec{\beta }}_{c}, \widehat{\varvec{\beta }}_{c'}\right)&=\left( \sum _{g=1}^{G_c}{\varvec{K}}_{gc}\right) ^{-1}\sum _{g=1}^{G_c}{\varvec{K}}_{gc} {\text{ Cov }}\left( \widehat{\varvec{\beta }}_{gc},\widehat{\varvec{\beta }}_{c'}\right) \\&=\delta _{cc'} \left[ {\varvec{U}} \left( \sum _{g'=1 }^{H_c}{\varvec{K}}_{g'c'}^{'}\right) ^{-1} + \varvec{\varPsi } \right] \end{aligned}$$

and similarly we can prove (40). Expression (41) can be obtained using (39) as follows

$$\begin{aligned} {\text{ Cov }}\left( \widehat{\varvec{\beta }}, \widehat{\varvec{\beta }}_{c'}\right)&= \left( \sum _{c=1}^{C}{\varvec{K}}_{c}\right) ^{-1} \left[ \sum _{c=1}^{C}{\varvec{K}}_{c} \ {\text{ Cov }}(\widehat{\varvec{\beta }}_{c},\widehat{\varvec{\beta }}_{c'})\right] \\&=\left( \sum _{c=1}^{C}{\varvec{K}}_{c}\right) ^{-1} \left\{ \sum _{c=1}^{C}{\varvec{K}}_{c} \ \delta _{cc'} \left[ \left( \sum _{g=1}^{G_c}{\varvec{K}}_{gc}\right) ^{-1} {\varvec{U}} + \varvec{\varPsi } \right] \right\} \\&=\left( \sum _{c=1}^{C}{\varvec{K}}_{c}\right) ^{-1} \left\{ \varvec{\varPsi }\left[ \left( \sum _{g=1}^{G_{c}}{\varvec{K}}_{gc}\right) ^{-1}{\varvec{U}} + \varvec{\varPsi }\right] ^{-1} \left[ \left( \sum _{g=1}^{G_c}{\varvec{K}}_{gc}\right) ^{-1} {\varvec{U}} + \varvec{\varPsi } \right] \right\} \\&=\left( \sum _{c=1}^{C}{\varvec{K}}_{c}\right) ^{-1} \varvec{\varPsi }, \end{aligned}$$

while (42) can be obtained as ${\text{ Cov }}\left( \widehat{\varvec{\beta }}_{c}, \widehat{\varvec{\beta }}\right) = {\text{ Cov }}\left( \widehat{\varvec{\beta }}, \widehat{\varvec{\beta }}_{c}\right) ^{'} = \varvec{\varPsi } \left( \sum _{c'=1}^{C}{\varvec{K}}_{c'}^{'}\right) ^{-1}$. Finally, expression (43) can be proved using (41)

$$\begin{aligned} {\text{ Cov }}\left( \widehat{\varvec{\beta }}, \widehat{\varvec{\beta }}\right)&= {\text{ Cov }}\left[ \widehat{\varvec{\beta }},\left( \sum _{c'=1}^{C}{\varvec{K}}_{c'}\right) ^{-1}\sum _{c'=1}^{C}{\varvec{K}}_{c'}\widehat{\varvec{\beta }}_{c'}\right] \\&=E\left[ \left( \widehat{\varvec{\beta }}-\varvec{\beta }\right) \sum _{c'=1}^{C}\left( \widehat{\varvec{\beta }}_{c'}-\varvec{\beta }\right) ^{'}{\varvec{K}}^{'}_{c'}\left( \sum _{c'=1}^{C}{\varvec{K}}_{c'}^{'}\right) ^{-1}\right] \\&=\left[ \sum _{c'=1}^{C} {\text{ Cov }}\left( \widehat{\varvec{\beta }},\widehat{\varvec{\beta }}_{c'}\right) {\varvec{K}}^{'}_{c'}\right] \left( \sum _{c'=1}^{C}{\varvec{K}}_{c'}^{'}\right) ^{-1} \\&=\left[ \sum _{c'=1}^{C} \left( \sum _{c=1}^{C}{\varvec{K}}_{c}\right) ^{-1} \varvec{\varPsi } {\varvec{K}}^{'}_{c'}\right] \left( \sum _{c'=1}^{C}{\varvec{K}}_{c'}^{'}\right) ^{-1} \\&=\left( \sum _{c=1}^{C}{\varvec{K}}_{c}\right) ^{-1} \varvec{\varPsi } \left( \sum _{c'=1}^{C} {\varvec{K}}^{'}_{c'}\right) \left( \sum _{c'=1}^{C}{\varvec{K}}_{c'}^{'}\right) ^{-1} \\&=\left( \sum _{c=1}^{C}{\varvec{K}}_{c}\right) ^{-1} \varvec{\varPsi }. \end{aligned}$$

$\square$

D Proof of Lemma 1

• (a) For country level it is sufficient to show that (53) satisfies (47) and (48) of Theorem 3. Expectation unbiasedness holds from (26) and (22) as follows

  $$\begin{aligned} E\left[ \varvec{\beta }_{c}^{Cred}\left( \varTheta _{c}^{+3}\right) \right]&= E\left[ {\varvec{K}}_{c} \ \widehat{\varvec{\beta }}_{c} + \left( {\varvec{I}}_{2}-{\varvec{K}}_{c}\right) \ \varvec{\beta }\right] \nonumber \\&= {\varvec{K}}_{c} \ E\left[ \widehat{\varvec{\beta }}_{c}\right] + \left( {\varvec{I}}_{2}-{\varvec{K}}_{c}\right) \ \varvec{\beta } \nonumber \\&= \varvec{\beta }, \end{aligned}$$

  (68)

  which is equal to $E\left[ \varvec{\beta }\left( \varTheta _{c}^{+3}\right) \right]$. The covariance condition is proved using expressions (35) and (29).

  $$\begin{aligned} {\text{ Cov }}\left[ \varvec{\beta }_{c}^{Cred}\left( \varTheta _{c}^{+3}\right) , \widehat{\varvec{\beta }}_{x'g'c'} \right)&= {\text{ Cov }}\left[ {\varvec{K}}_{c} \ \widehat{\varvec{\beta }}_{c}+\left( {\varvec{I}}_{2}-{\varvec{K}}_{c}\right) \ \varvec{\beta },\widehat{\varvec{\beta }}_{x'g'c'}\right] \nonumber \\&={\varvec{K}}_{c} \ {\text{ Cov }}\left( \widehat{\varvec{\beta }}_{c},\widehat{\varvec{\beta }}_{x'g'c'}\right) + ({\varvec{I}}_{2}-{\varvec{K}}_{c}) \ {\text{ Cov }}\left( \varvec{\beta } ,\widehat{\varvec{\beta }}_{x'g'c'}\right) \nonumber \\&= \varvec{\varPsi }\left[ \varvec{\varPsi } + \left( \sum _{g=1}^{G_{c}}{\varvec{K}}_{gc}\right) ^{-1}{\varvec{U}}\right] ^{-1} \delta _{cc'} \left[ \left( \sum _{g=1}^{G_{c}}{\varvec{K}}_{gc}\right) ^{-1}{\varvec{U}} + \varvec{\varPsi }\right] \nonumber \\&= \delta _{cc'} \varvec{\varPsi }, \end{aligned}$$

  (69)

  which is equal to ${\text{ Cov }}\left[ \varvec{\beta }\left( \varTheta _{c}^{+3}\right) ,\widehat{\varvec{\beta }}_{x'g'c'}\right]$.

• (b) Similarly, for gender level it is sufficient to show that (54) satisfies (49) and (50). The expectation condition can be proved by (68), (26) and (24) as

  $$\begin{aligned} E\left[ \varvec{\beta }_{gc}^{Cred}\left( \varTheta _{gc}^{+23}\right) \right]&= E\left[ {\varvec{K}}_{gc}\widehat{\varvec{\beta }}_{gc}+\left( {\varvec{I}}_{2}-{\varvec{K}}_{gc}\right) \ \varvec{\beta }^{Cred}\left( \varTheta _{c}^{+3}\right) \right] \nonumber \\&= {\varvec{K}}_{gc}E\left[ \widehat{\varvec{\beta }}_{gc}\right] +\left( {\varvec{I}}_{2}-{\varvec{K}}_{gc}\right) \ E\left[ \varvec{\beta }^{Cred}\left( \varTheta _{c}^{+3}\right) \right] \nonumber \\&= \varvec{\beta } = E\left( \widehat{\varvec{\beta }}_{xgc}\right) = E\left[ E\left( \widehat{\varvec{\beta }}_{xgc}\big |\varTheta _{gc}^{+23}\right) \right] = E\left[ {{\varvec{\beta }}}\left( \varTheta _{gc}^{+23}\right) \right] . \end{aligned}$$

  (70)

  Then, from (33), (69) and (28) we get

  $$\begin{aligned}&{\text{ Cov }}\left[ \varvec{\beta }_{gc}^{Cred}\left( \varTheta _{gc}^{+23}\right) , \widehat{\varvec{\beta }}_{x'g'c'}\right] \nonumber \\&\quad = {\text{ Cov }}\left[ {\varvec{K}}_{gc}\widehat{\varvec{\beta }}_{gc}+\left( {\varvec{I}}_{2}-{\varvec{K}}_{gc}\right) \ \varvec{\beta }^{Cred}\left( \varTheta _{c}^{+3}\right) ,\widehat{\varvec{\beta }}_{x'g'c'}\right] \nonumber \\&\quad = {\varvec{K}}_{gc}{\text{ Cov }}\left( \widehat{\varvec{\beta }}_{gc},\widehat{\varvec{\beta }}_{x'g'c'}\right) + \left( {\varvec{I}}_{2}-{\varvec{K}}_{gc}\right) \ {\text{ Cov }}\left[ \varvec{\beta }^{Cred}\left( \varTheta _{c}^{+3}\right) ,\widehat{\varvec{\beta }}_{x'g'c'}\right] \nonumber \\&\quad = {\varvec{U}}\left[ {\varvec{U}} + \left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}\right) ^{-1}{\varvec{A}}\right] ^{-1} \left\{ \delta _{cc'} \{\delta _{gg'}\left[ \left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}\right) ^{-1} \ {\varvec{A}} + \ {\varvec{U}}\right] + \varvec{\varPsi }\}\right\} \nonumber \\&\qquad + \left\{ {\varvec{I}}_{2}-{\varvec{U}}\left[ {\varvec{U}} +\left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}\right) ^{-1}{\varvec{A}}\right] ^{-1} \right\} \ \delta _{cc'} \varvec{\varPsi } \nonumber \\&\quad = \delta _{cc'} \left( \delta _{gg'} {\varvec{U}} + \varvec{\varPsi }\right) , \end{aligned}$$

  (71)

  which is equal to ${\text{ Cov }}[\varvec{\beta }(\varTheta _{gc}^{+23}),\widehat{\varvec{\beta }}_{x'g'c'}]$.

• (c) Finally, for the estimator (55) of age level, we have to prove (51) and (52). For the expectation condition, we use (26), (70) and (S1) as follows

  $$\begin{aligned} E\left[ \varvec{\beta }_{xgc}^{Cred}\left( \varTheta _{xgc}^{+123}\right) \right]&= E\left[ {\varvec{K}}_{xgc}\widehat{\varvec{\beta }}_{xgc}+\left( {\varvec{I}}_{2}-{\varvec{K}}_{xgc}\right) \ \varvec{\beta }^{Cred}\left( \varTheta _{gc}^{+23}\right) \right] \nonumber \\&= {\varvec{K}}_{xgc}E\left( \widehat{\varvec{\beta }}_{xgc}\right) +\left( {\varvec{I}}_{2}-{\varvec{K}}_{xgc}\right) \ E\left[ \varvec{\beta }^{Cred}\left( \varTheta _{gc}^{+23}\right) \right] \nonumber \\&= \varvec{\beta }, \end{aligned}$$

  (72)

  which is equal to $E\left[ \varvec{\beta }\left( \varTheta _{xgc}^{+123}\right) \right]$. For the covariance, the proof is given by using (30), (71) and (27).

  $$\begin{aligned}&{\text{ Cov }}\left[ \varvec{\beta }_{xgc}^{Cred}\left( \varTheta _{xgc}^{+123}\right) ,\widehat{\varvec{\beta }}_{x'g'c'}\right] \nonumber \\&\quad = {\text{ Cov }}\left[ {\varvec{K}}_{xgc} \ \widehat{\varvec{\beta }}_{xgc} + \left( {\varvec{I}}_{2}-{\varvec{K}}_{xgc}\right) \ \varvec{\beta }^{Cred}\left( \varTheta _{gc}^{+23}\right) , \widehat{\varvec{\beta }}_{x'g'c'}\right] \nonumber \\&\quad = {\varvec{K}}_{xgc}{\text{ Cov }}\left( \widehat{\varvec{\beta }}_{xgc},\widehat{\varvec{\beta }}_{x'g'c'}\right) \nonumber \\&\qquad + \left( {\varvec{I}}_{2}-{\varvec{K}}_{xgc}\right) \ {\text{ Cov }}\left[ \varvec{\beta }^{Cred}\left( \varTheta _{gc}^{+23}\right) ,\widehat{\varvec{\beta }}_{x'g'c'}\right] \nonumber \\&\quad = {\varvec{K}}_{xgc} \delta _{cc'} \left\{ \delta _{gg'}\left[ \delta _{xx'}\left( {\varvec{A}}+s^{2}\left( {\varvec{Z}}_{xgc}^{'}{\varvec{W}}_{xgc} {\varvec{Z}}_{xgc}\right) ^{-1}\right) +{\varvec{U}}\right] + \varvec{\varPsi }\right\} \nonumber \\&\qquad + \left( {\varvec{I}}_{2}-{\varvec{K}}_{xgc}\right) \ \delta _{cc'} \left( \delta _{gg'} {\varvec{U}} + \varvec{\varPsi }\right) \nonumber \\&\quad = \delta _{cc'} \left[ \delta _{gg'}\left( \delta _{xx'} {\varvec{A}} + {\varvec{U}}\right) + \varvec{\varPsi }\right] , \end{aligned}$$

  (73)

  which is equal to ${\text{ Cov }}\left[ \varvec{\beta }\left( \varTheta _{xgc}^{+123}\right) ,\widehat{\varvec{\beta }}_{x'g'c'}\right]$. $\square$

E Proof of Theorem 4

For the proof of expression (56) let

$$\begin{aligned} S^{\star }&=E\left[ \left( {\varvec{Y}}_{xgc}-{\varvec{Z}}_{xgc}\widehat{\varvec{\beta }}_{xgc}\right) ^{'}{\varvec{W}}_{xgc}\left( {\varvec{Y}}_{xgc}-{\varvec{Z}}_{xgc}\widehat{\varvec{\beta }}_{xgc}\right) \right] \\&=E\left[ {\varvec{Y}}_{xgc}^{'}\left( {\varvec{W}}_{xgc} - {\varvec{W}}_{xgc}{\varvec{Z}}_{xgc}\left( {\varvec{Z}}_{xgc}^{'}{\varvec{W}}_{xgc} {\varvec{Z}}_{xgc}\right) ^{-1} {\varvec{Z}}_{xgc}^{'}{\varvec{W}}_{xgc}\right) {\varvec{Y}}_{xgc}\right] \\&=E\left\{ tr\left[ {\varvec{Y}}_{xgc}^{'}\left( {\varvec{W}}_{xgc} - {\varvec{W}}_{xgc}{\varvec{Z}}_{xgc} \left( {\varvec{Z}}_{xgc}^{'}{\varvec{W}}_{xgc} {\varvec{Z}}_{xgc}\right) ^{-1} {\varvec{Z}}_{xgc}^{'}{\varvec{W}}_{xgc}\right) {\varvec{Y}}_{xgc}\right] \right\} \\&=tr\left[ \left( {\varvec{W}}_{xgc} - {\varvec{W}}_{xgc}{\varvec{Z}}_{xgc} \left( {\varvec{Z}}_{xgc}^{'}{\varvec{W}}_{xgc} {\varvec{Z}}_{xgc}\right) ^{-1} {\varvec{Z}}_{xgc}^{'}{\varvec{W}}_{xgc}\right) E\left( {\varvec{Y}}_{xgc}{\varvec{Y}}_{xgc}^{'}\right) \right] , \end{aligned}$$

where $E\left( {\varvec{Y}}_{xgc}{\varvec{Y}}_{xgc}^{'}\right) = s^2 \ {\varvec{W}}_{xgc}^{-1} + {\varvec{Z}}_{xgc} \varvec{\beta } \ \varvec{\beta }^{'} {\varvec{Z}}_{xgc}^{'} \ .$ By recalling the linearity of trace, we take

$$\begin{aligned} S^{\star }&= s^{2} \ tr\left( {\varvec{W}}_{xgc}^{'}{\varvec{W}}_{xgc}\right) + tr\left( {\varvec{W}}_{xgc} {\varvec{Z}}_{xgc} \ \varvec{\beta } \ \varvec{\beta }^{'} {\varvec{Z}}_{xgc}^{'} \right) \\&\qquad - \ s^2 \ tr\left[ \left( {\varvec{Z}}_{xgc}^{'}{\varvec{W}}_{xgc} {\varvec{Z}}_{xgc}\right) ^{-1} {\varvec{Z}}_{xgc}^{'}{\varvec{W}}_{xgc}{\varvec{Z}}_{xgc} \right] \\&\qquad \qquad - \ tr\left( {\varvec{W}}_{xgc} {\varvec{Z}}_{xgc} \ \varvec{\beta } \ \varvec{\beta }^{'} {\varvec{Z}}_{xgc}^{'} \right) \\&= s^{2} \ \left[ tr\left( {\varvec{I}}_{n}\right) - tr\left( {\varvec{I}}_{2} \right) \right] \\&= s^{2} \ (n - 2), \end{aligned}$$

that proves (56).

For the proof of (57), we use (30)–(33) as follows

$$\begin{aligned} E\left[ \widehat{{\varvec{A}}}^{\star }\right]&=E\left[ \sum _{c=1}^{C}\sum _{g=1}^{G_{c}}\sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc} \left( \widehat{\varvec{\beta }}_{xgc}-\widehat{\varvec{\beta }}_{gc} \right) \left( \widehat{\varvec{\beta }}_{xgc} -\widehat{\varvec{\beta }}_{gc}\right) ^{'}\right] \\&=\sum _{c=1}^{C}\sum _{g=1}^{G_{c}}\sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}\left[ {\text{ Cov }}\left( \widehat{\varvec{\beta }}_{xgc},\widehat{\varvec{\beta }}_{xgc}\right) \right. \\&\left. \quad -{\text{ Cov }}\left( \widehat{\varvec{\beta }}_{xgc},\widehat{\varvec{\beta }}_{gc}\right) -{\text{ Cov }}\left( \widehat{\varvec{\beta }}_{gc},\widehat{\varvec{\beta }}_{xgc}\right) +{\text{ Cov }}\left( \widehat{\varvec{\beta }}_{gc},\widehat{\varvec{\beta }}_{gc}\right) \right] \\&=\sum _{c=1}^{C}\sum _{g=1}^{G_{c}}\sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc} \left[ {\varvec{A}}+s^{2}\left( {\varvec{Z}}_{xgc}^{'}{\varvec{W}}_{xgc} {\varvec{Z}}_{xgc}\right) ^{-1} + {\varvec{U}} + \varvec{\varPsi } \right] \\&\qquad - \sum _{c=1}^{C}\sum _{g=1}^{G_{c}}\sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc} \left[ {\varvec{A}}\left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}^{'}\right) ^{-1} + {\varvec{U}} + \varvec{\varPsi } \right] \\&\qquad - \sum _{c=1}^{C}\sum _{g=1}^{G_{c}}\sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc} \left[ \left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}\right) ^{-1} {\varvec{A}} + {\varvec{U}} + \varvec{\varPsi } \right] \\&\qquad + \sum _{c=1}^{C}\sum _{g=1}^{G_{c}}\sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc} \left[ {\varvec{A}}\left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}^{'}\right) ^{-1} + {\varvec{U}} + \varvec{\varPsi } \right] \\&= \sum _{c=1}^{C}\sum _{g=1}^{G_{c}}\sum _{x=x_1}^{x_{k_{gc}}} {\varvec{K}}_{xgc} \left\{ \left[ {\varvec{A}}+s^{2}\left( {\varvec{Z}}_{xgc}^{'}{\varvec{W}}_{xgc} {\varvec{Z}}_{xgc}\right) ^{-1} + {\varvec{U}} + \varvec{\varPsi } \right] \right. \\&\left. \quad - \left[ \left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}\right) ^{-1} {\varvec{A}} + {\varvec{U}} + \varvec{\varPsi } \right] \right\} \\&= \sum _{c=1}^{C}\sum _{g=1}^{G_{c}}\sum _{x=x_1}^{x_{k_{gc}}} {\varvec{A}} - \sum _{c=1}^{C}\sum _{g=1}^{G_{c}} {\varvec{A}} = \sum _{c=1}^{C}\sum _{g=1}^{G_c}(x_{k_{gc}}-x_{1}) {\varvec{A}}. \end{aligned}$$

Similarly, (58) is proved via (34), (37), (38) and (39) as follows

$$\begin{aligned} E\left[ \widehat{{\varvec{U}}}^{\star }\right]&= E\left[ \sum _{c=1}^{C}\sum _{g=1}^{G_{c}}{\varvec{K}}_{gc} \left( \widehat{\varvec{\beta }}_{gc}-\widehat{\varvec{\beta }}_{c}\right) \left( \widehat{\varvec{\beta }}_{gc} -\widehat{\varvec{\beta }}_{c}\right) ^{'}\right] \\&= \sum _{c=1}^{C}\sum _{g=1}^{G_{c}} {\varvec{K}}_{gc}\left[ {\text{ Cov }}\left( \widehat{\varvec{\beta }}_{gc},\widehat{\varvec{\beta }}_{gc}\right) \right. \\&\left. \quad -{\text{ Cov }}\left( \widehat{\varvec{\beta }}_{gc},\widehat{\varvec{\beta }}_{c}\right) -{\text{ Cov }}\left( \widehat{\varvec{\beta }}_{c},\widehat{\varvec{\beta }}_{gc}\right) +{\text{ Cov }}\left( \widehat{\varvec{\beta }}_{c},\widehat{\varvec{\beta }}_{c}\right) \right] \\&= \sum _{c=1}^{C}\sum _{g=1}^{G_{c}} {\varvec{K}}_{gc} \left\{ \left[ \left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}\right) ^{-1}{\varvec{A}}+ {\varvec{U}} + \varvec{\varPsi }\right] \right. \\&\left. \quad - \left[ {\varvec{U}} \left( \sum _{g=1}^{G}{\varvec{K}}_{gc}^{'}\right) ^{-1} + \varvec{\varPsi } \right] \right. \\&\left. \qquad - \left[ \left( \sum _{g=1}^{G}{\varvec{K}}_{gc}\right) ^{-1} {\varvec{U}} + \varvec{\varPsi } \right] + \left[ {\varvec{U}} \left( \sum _{g=1}^{G}{\varvec{K}}_{gc}^{'}\right) ^{-1} + \varvec{\varPsi } \right] \right\} \\&= \sum _{c=1}^{C}\sum _{g=1}^{G_{c}} {\varvec{K}}_{gc} \left\{ \left[ \left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}\right) ^{-1}{\varvec{A}}+ {\varvec{U}} + \varvec{\varPsi }\right] - \left[ \left( \sum _{g=1}^{G}{\varvec{K}}_{gc}\right) ^{-1} {\varvec{U}} + \varvec{\varPsi } \right] \right\} \\&=\sum _{c=1}^{C}\sum _{g=1}^{G_{c}} {\varvec{U}}\left[ {\varvec{U}} + \left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}\right) ^{-1}{\varvec{A}}\right] ^{-1} \left[ \left( \sum _{x=x_1}^{x_{k_{gc}}}{\varvec{K}}_{xgc}\right) ^{-1}{\varvec{A}}+ {\varvec{U}} \right] - \sum _{c=1}^{C} {\varvec{U}} \\&=\sum _{c=1}^{C}(G_c -1) {\varvec{U}}. \end{aligned}$$

Finally, (59) is obtained by using (40)–(43) as follows

$$\begin{aligned} E\left[ \widehat{\varvec{\varPsi }}^{\star }\right]&=E\left[ \sum _{c=1}^{C}{\varvec{K}}_{c} \left( \widehat{\varvec{\beta }}_{c}-\widehat{\varvec{\beta }}\right) \left( \widehat{\varvec{\beta }}_{c} -\widehat{\varvec{\beta }}\right) ^{'}\right] \\&= \sum _{c=1}^{C} {\varvec{K}}_{c}\left[ {\text{ Cov }}\left( \widehat{\varvec{\beta }}_{c},\widehat{\varvec{\beta }}_{c}\right) -{\text{ Cov }}\left( \widehat{\varvec{\beta }}_{c},\widehat{\varvec{\beta }}\right) -{\text{ Cov }}\left( \widehat{\varvec{\beta }},\widehat{\varvec{\beta }}_{c}\right) +{\text{ Cov }}\left( \widehat{\varvec{\beta }},\widehat{\varvec{\beta }}\right) \right] \\&= \sum _{c=1}^{C} {\varvec{K}}_{c} \left\{ \left[ \left( \sum _{g=1}^{G}{\varvec{K}}_{gc}\right) ^{-1} {\varvec{U}} + \varvec{\varPsi } \right] - \left[ \varvec{\varPsi } \left( \sum _{c=1}^{C}{\varvec{K}}_{c}^{'}\right) ^{-1} \right] \right. \\&\qquad \left. - \ \left[ \left( \sum _{c=1}^{C}{\varvec{K}}_{c}\right) ^{-1} \varvec{\varPsi } \right] + \left[ \varvec{\varPsi } \left( \sum _{c=1}^{C}{\varvec{K}}_{c}^{'}\right) ^{-1} \right] \right\} \\&=\sum _{c=1}^{C} \varvec{\varPsi } \left[ \varvec{\varPsi } + \left( \sum _{g=1}^{G_{c}}{\varvec{K}}_{gc}\right) ^{-1}{\varvec{U}} \right] ^{-1} \left[ \left( \sum _{g=1}^{G}{\varvec{K}}_{gc}\right) ^{-1} {\varvec{U}} + \varvec{\varPsi } \right] \\&\qquad - \ \sum _{c=1}^{C}{\varvec{K}}_{c} \left[ \left( \sum _{c=1 }^{C}{\varvec{K}}_{c}\right) ^{-1}\varvec{\varPsi } \right] \\&=(C-1) \varvec{\varPsi }. \end{aligned}$$

$\square$