* •

Fact 1: The qu$d$it restriction of the bosonic loss-dephasing channel
$\pazocal{N}^{(d)}_{\lambda,\gamma}$ is anti-degradable if and only if
$\lambda\leq\lambda_{d}(\gamma)$

  * •

Fact 2: The quantity $\lambda_{d}(\gamma)$ is monotonically increasing in
$\gamma$

  * •

Fact 3: For $d=2$, it precisely holds that
$\lambda_{2}(\gamma)=\frac{1}{1+e^{-\gamma}}$

  * •

Fact 4: The quantity $\lambda_{d}(\gamma)$ is monotonically non-increasing in
$d$

  * •

Fact 5: It holds that
$\frac{1}{2}\leq\lambda_{d}(\gamma)\leq\frac{1}{1+e^{-\gamma}}$

  * •

Fact 6: When $d=3$ and $e^{-\gamma}\leq\sqrt{2}-1$ (or $\gamma\geq 0.881$), it
exactly holds that $\lambda_{3}(\gamma)=\frac{1}{1+e^{-\gamma}}$

###### Proof.

_Fact 1._ It suffices to show that for any $\gamma\geq 0$ and
$\lambda,\lambda^{\prime}\in[0,1]$ with $\lambda^{\prime}<\lambda$, if
$\pazocal{N}^{(d)}_{\lambda,\gamma}$ is anti-degradable, then so is
$\pazocal{N}^{(d)}_{\lambda^{\prime},\gamma}$. To show this, we exploit the
composition rule

$\displaystyle\pazocal{E}_{\lambda_{1}}\circ\pazocal{N}_{\lambda_{2},\gamma}=\pazocal{N}_{\lambda_{1}\lambda_{2},\gamma}\qquad\forall\,\lambda_{1},\lambda_{2}\in[0,1]\,,$
(S58)

as established by Lemma S25, implying that the channel
$\pazocal{N}^{(d)}_{\lambda^{\prime},\gamma}$ can be written as the
composition between $\pazocal{N}^{(d)}_{\lambda,\gamma}$ and another channel.
Consequently, Lemma S50 concludes the proof.

_Fact 2._ Analogously to Fact 1, it suffices to show that for any
$\lambda\in[0,1]$ and for any $\gamma^{\prime}\geq\gamma\geq 0$, if
$\pazocal{N}^{(d)}_{\lambda,\gamma}$ is anti-degradable, then so is
$\pazocal{N}^{(d)}_{\lambda,\gamma^{\prime}}$. This follows from the
composition rule

$\displaystyle\pazocal{D}_{\gamma_{1}}\circ\pazocal{N}_{\lambda,\gamma_{2}}=\pazocal{N}_{\lambda,\gamma_{1}+\gamma_{2}}\qquad\forall\,\gamma_{1},\gamma_{2}\geq
0\,,$ (S59)

as proved in Lemma S25. Furthermore, Lemma S50 concludes the proof.

_Fact 3._ This has already been proved in the Proof 3 of Theorem S35.

_Fact 4._ This follows from the observation that for all $d^{\prime}\geq d$,
if $\pazocal{N}_{\lambda,\gamma}^{(d^{\prime})}$ is anti-degradable, then so
is $\pazocal{N}_{\lambda,\gamma}^{(d)}$.

_Fact 5._ The upper bound $\lambda_{d}(\gamma)\leq\frac{1}{1+e^{-\gamma}}$
follows from Fact 3 and Fact 4. Moreover, since the pure-loss channel
$\pazocal{E}_{\lambda}$ is anti-degradable if and only if
$\lambda\leq\frac{1}{2}$, Fact 4 implies that $\lambda_{d}(0)\geq\frac{1}{2}$
(more specifically, one can also show that $\lambda_{d}(0)=\frac{1}{2}$).
Consequently, Fact 2 concludes the proof.

_Fact 6._ This proof relies on the equivalence between anti-degradability of a
channel and two-extendibility of its Choi state, as established in Lemma S49.
Let $\lambda$ and $\gamma$ be such that $\lambda=\frac{1}{1+e^{-\gamma}}$ and
$e^{-\gamma}\leq\sqrt{2}-1$, implying that $\lambda\geq\frac{1}{\sqrt{2}}$. By
using Lemma S22, we obtain the Choi state of
$\pazocal{N}^{(3)}_{\lambda,\gamma}$ as follows:

$\displaystyle\operatorname{id}_{A}\otimes\pazocal{N}_{\lambda,\gamma}(\outerproduct{\Phi_{3}}{\Phi_{3}})=\frac{1}{3}\sum_{m=0}^{2}\sum_{n=0}^{2}\sum_{\ell=0}^{\min(m,n)}e^{-\frac{\gamma}{2}(m-n)^{2}}\sqrt{\binom{m}{\ell}\binom{n}{\ell}}\lambda^{\frac{m+n}{2}-\ell}(1-\lambda)^{l}\outerproduct{m}{n}_{A}\otimes\outerproduct{m-\ell}{n-\ell}_{B},$

where $\Phi_{3}$ is the maximally entangled state of Schmidt rank $3$. We
define a two-extension $\tilde{\rho}_{AB_{1}B_{2}}$ of the Choi state by the
following conditions. First, $\tilde{\rho}_{AB_{1}B_{2}}$ satisfies the
following $B_{1}\leftrightarrow B_{2}$ symmetry for all
$i_{1},i_{2},i_{3},j_{1},j_{2},j_{3}\in\\{0,1,2\\}$:

$\displaystyle\bra{j_{1}}_{A}\bra{j_{2}}_{B_{1}}\bra{j_{3}}_{B_{2}}\tilde{\rho}_{AB_{1}B_{2}}\ket{i_{1}}_{A}\ket{i_{2}}_{B_{1}}\ket{i_{3}}_{B_{2}}$
$\displaystyle=\bra{j_{1}}_{A}\bra{j_{3}}_{B_{1}}\bra{j_{2}}_{B_{2}}\tilde{\rho}_{AB_{1}B_{2}}\ket{i_{1}}_{A}\ket{i_{2}}_{B_{1}}\ket{i_{3}}_{B_{2}}$
(S60)
$\displaystyle\bra{j_{1}}_{A}\bra{j_{2}}_{B_{1}}\bra{j_{3}}_{B_{2}}\tilde{\rho}_{AB_{1}B_{2}}\ket{i_{1}}_{A}\ket{i_{2}}_{B_{1}}\ket{i_{3}}_{B_{2}}$
$\displaystyle=\bra{j_{1}}_{A}\bra{j_{2}}_{B_{1}}\bra{j_{3}}_{B_{2}}\tilde{\rho}_{AB_{1}B_{2}}\ket{i_{1}}_{A}\ket{i_{3}}_{B_{1}}\ket{i_{2}}_{B_{2}}\,.$

Furthermore, if $i_{1}<\max(i_{2},i_{3})$ or $j_{1}<\max(j_{2},j_{3})$, or if
$i_{1}>i_{2}+i_{3}$ or $j_{1}>j_{2}+j_{3}$, then
$\bra{j_{1}}_{A}\bra{j_{2}}_{B_{1}}\bra{j_{3}}_{B_{2}}\tilde{\rho}_{AB_{1}B_{2}}\ket{i_{1}}_{A}\ket{i_{2}}_{B_{1}}\ket{i_{3}}_{B_{2}}=0\,.$
We can thus define $\tilde{\rho}_{AB_{1}B_{2}}$ by writing only the matrix
elements with respect the set
$\\{\ket{i_{1}}_{A}\ket{i_{2}}_{B_{1}}\ket{i_{3}}_{B_{2}}\\}$ with $2\geq
i_{1}\geq i_{2}\geq i_{3}\geq 0$ such that $i_{2}+i_{3}\geq i_{1}$. Hence, in
order to fully define $\tilde{\rho}_{AB_{1}B_{2}}$, it suffices to write the
matrix elements of $\tilde{\rho}_{AB_{1}B_{2}}$ with respect to the set $\\{$
$\ket{0}_{A}\ket{0}_{B_{1}}\ket{0}_{B_{2}}$,
$\ket{1}_{A}\ket{0}_{B_{1}}\ket{0}_{B_{2}}$,
$\ket{1}_{A}\ket{1}_{B_{1}}\ket{1}_{B_{2}}$,
$\ket{2}_{A}\ket{1}_{B_{1}}\ket{1}_{B_{2}}$ ,
$\ket{2}_{A}\ket{2}_{B_{1}}\ket{0}_{B_{2}}$,
$\ket{2}_{A}\ket{2}_{B_{1}}\ket{1}_{B_{2}}$ ,
$\ket{2}_{A}\ket{2}_{B_{1}}\ket{2}_{B_{2}}$ $\\}$. This gives rise to the
following $7\times 7$ matrix:

$\begin{array}[]{c|ccccccc}&0_{A}0_{B_{1}}0_{B_{2}}&1_{A}1_{B_{1}}0_{B_{2}}&1_{A}1_{B_{1}}1_{B_{2}}&2_{A}1_{B_{1}}1_{B_{2}}&2_{A}2_{B_{1}}0_{B_{2}}&2_{A}2_{B_{1}}1_{B_{2}}&2_{A}2_{B_{1}}2_{B_{2}}\\\
\hline\cr\\\
0_{A}0_{B_{1}}0_{B_{2}}&1&\sqrt{1-\lambda}&0&\sqrt{2}(1-\lambda)&\frac{(1-\lambda)^{2}}{\lambda}&0&0\\\
\\\
1_{A}1_{B_{1}}0_{B_{2}}&\sqrt{1-\lambda}&1-\lambda&0&\sqrt{2(1-\lambda)^{3}}&\sqrt{\frac{(1-\lambda)^{5}}{\lambda^{2}}}&0&0\\\
\\\
1_{A}1_{B_{1}}1_{B_{2}}&0&0&2\lambda-1&0&0&\frac{2\lambda-1}{\lambda}\sqrt{1-\lambda}&0\\\
\\\
2_{A}1_{B_{1}}1_{B_{2}}&\sqrt{2}(1-\lambda)&\sqrt{2(1-\lambda)^{3}}&0&2(1-\lambda)^{2}&\sqrt{2}\frac{(1-\lambda)^{3}}{\lambda}&0&0\\\
\\\
2_{A}2_{B_{1}}0_{B_{2}}&\frac{(1-\lambda)^{2}}{\lambda}&\sqrt{\frac{(1-\lambda)^{5}}{\lambda^{2}}}&0&\sqrt{2}\frac{(1-\lambda)^{3}}{\lambda}&(1-\lambda)^{2}&0&0\\\
\\\
2_{A}2_{B_{1}}1_{B_{2}}&0&0&\frac{2\lambda-1}{\lambda}\sqrt{1-\lambda}&0&0&2(1-\lambda)(2\lambda-1)&0\\\
\\\ 2_{A}2_{B_{1}}2_{B_{2}}&0&0&0&0&0&0&(2\lambda-1)^{2}\end{array}$

One can show by direct calculation that this matrix is positive semi-definite
if and only $\lambda\geq\frac{1}{\sqrt{2}}$. Note that the matrix is positive
semi-definite if and only if $\tilde{\rho}_{AB_{1}B_{2}}$ is positive semi-
definite. The latter follows from the following two simple facts: (i) A
$n\times n$ symmetric matrix with a duplicate column is positive semi-definite
if and only if the $(n-1)\times(n-1)$ matrix obtained by deleting one of the
two equal column and the corresponding row is positive semi-definite, and (ii)
A $n\times n$ symmetric matrix with a zero column is positive semi-definite if
and only if the $(n-1)\times(n-1)$ matrix obtained by deleting such a column
and the corresponding row is positive semi-definite. One can also verify
$\Tr_{B_{1}}\tilde{\rho}_{AB_{1}B_{2}}=\Tr_{B_{2}}\tilde{\rho}_{AB_{1}B_{2}}$,
and they are equal to the Choi state of the qutrit restriction with
$\lambda=\frac{1}{1+e^{-\gamma}}$. We therefore conclude that the curve
$\lambda=\frac{1}{1+e^{-\gamma}}$ provides a necessary and sufficient
condition for the anti-degradability of the qutrit channel
$\pazocal{N}_{\lambda,\gamma}^{(3)}$ when $\lambda\geq\frac{1}{\sqrt{2}}$, or
equivalently when $e^{-\gamma}\leq\sqrt{2}-1$. ∎

To numerically compute the quantity $\lambda_{d}(\gamma)$ given in (S57), we
utilise the equivalence between anti-degradability of a channel and two-
extendibility of its Choi state [42]. Specifically, for small values of $d$,
we can determine necessary and sufficient conditions for the anti-
degradability of $\pazocal{N}^{(d)}_{\lambda,\gamma}$ by numerically solving
the following _semi-definite program_ :

$\displaystyle\hskip 11.38092pt\min_{\rho_{AB_{1}B_{2}}}1$ (S61)
$\displaystyle\text{s.t. }\rho_{AB_{1}B_{2}}\geq 0\,,$
$\displaystyle\quad\Tr[\rho_{AB_{1}B_{2}}]=1,$
$\displaystyle\quad\Tr_{B_{2}}\left[\rho_{AB_{1}B_{2}}\right]=\operatorname{id}_{A}\otimes\pazocal{N}_{\lambda,\gamma}^{(d)}(\outerproduct{\Phi_{d}}{\Phi_{d}}_{AA^{\prime}})\,,$
$\displaystyle\quad\Tr_{B_{1}}\left[\rho_{AB_{1}B_{2}}\right]=\operatorname{id}_{A}\otimes\pazocal{N}_{\lambda,\gamma}^{(d)}(\outerproduct{\Phi_{d}}{\Phi_{d}}_{AA^{\prime}})\,.$

where $\ket{\Phi_{d}}$ is the maximally entangled state of schmidt rank $d$.
The channel $\pazocal{N}_{\lambda,\gamma}^{(d)}$ is anti-degradable if and
only if the semi-definite program admits a feasible solution. We compute the
quantity $\lambda_{d}(\gamma)$ defined in (S57) by numerically solving the
semi-definite program.

Figure 4: Each curve indicates necessary and sufficient conditions where the
qu$d$it restriction $\pazocal{N}_{\lambda,\gamma}^{(d)}$ of the bosonic loss-
dephasing channel $\pazocal{N}_{\lambda,\gamma}$ is anti-degradable. In the
region above the curve $\lambda_{5}(\gamma)$, the bosonic loss-dephasing
channel is never anti-degradable.

The results are plotted with respect to $e^{-\gamma}$ for various values of
$d$ in Fig. 4, showcasing the dependence of $\lambda_{d}(\gamma)$ on $\gamma$
for small values of $d$. Our numerical analysis reveals that when $\gamma$ is
sufficiently large, i.e. $e^{-\gamma}\lesssim 0.41$ or $\gamma\gtrsim 0.89$,
the value of $\lambda_{d}(\gamma)$ consistently equals
$\frac{1}{1+e^{-\gamma}}$ for all examined values of $d$. In particular, this
seems to suggest that within this range of the dephasing parameter,
$\pazocal{N}_{\lambda,\gamma}$ is anti-degradable if and only if
$\lambda\leq\frac{1}{1+e^{-\gamma}}$. Based on this numerical exploration, we
propose the following conjecture:

######  Conjecture S37.

If $\gamma$ is sufficiently large ($\gamma\gtrsim 0.89$), then the bosonic
loss-dephasing channel $\pazocal{N}_{\lambda,\gamma}$ is anti-degradable if
and only if $\lambda\leq\frac{1}{1+e^{-\gamma}}$.

Notably, from Fig. 4 we observe that if $\lambda$ and $\gamma$ satisfy
$\lambda=\frac{1}{1+e^{-\gamma}}$ with $e^{-\gamma}\gtrsim 0.41$ (or
$\gamma\lesssim 0.89$), then $\pazocal{N}_{\lambda,\gamma}$ is not anti-
degradable.

##  III Coherence preservation of the bosonic loss-dephasing channel

In this section we use the notation introduced in Section I.8.

######  Theorem S38.

Let $\gamma\geq 0$ and $\lambda\in(0,1]$. For any energy constraint $N_{s}>0$,
the energy-constrained two-way quantum and secret-key capacities of the
bosonic loss-dephasing channel are strictly positive,
$K(\pazocal{N}_{\lambda,\gamma},N_{s})\geq
Q_{2}(\pazocal{N}_{\lambda,\gamma},N_{s})>0$.

###### Proof.

We begin by assuming $N_{s}\in(0,1)$ and defining the two-mode state

$\displaystyle\ket{\Psi_{N_{s}}}_{AA^{\prime}}\coloneqq\sqrt{1-N_{s}}\ket{00}_{AA^{\prime}}+\sqrt{N_{s}}\ket{11}_{AA^{\prime}}\,,$
(S62)

where the mean photon number of $A^{\prime}$ system is equal to $N_{s}$. By
exploiting Lemma S22, one can observe that the state

$\displaystyle\rho_{AB}\coloneqq\operatorname{id}_{A}\otimes\pazocal{N}_{\lambda,\gamma}(\outerproduct{\Psi_{N_{s}}}{\Psi_{N_{s}}}_{AA^{\prime}})$
(S63)

is effectively a two-qubit state and its matrix with respect to the
computational basis $\\{\ket{00},\ket{01},\ket{10},\ket{11}\\}$ is given by:

$\displaystyle\rho_{AB}=$
$\displaystyle\left(\begin{matrix}1-N_{s}\quad&0&0&\sqrt{(1-N_{s})N_{s}e^{-\gamma}\lambda}\\\
0&0\quad&0\quad&0\\\ 0&0\quad&(1-\lambda)N_{s}\quad&0\\\
\sqrt{(1-N_{s})N_{s}e^{-\gamma}\lambda}&0&0&\lambda
N_{s}\quad\end{matrix}\right).$

If we perform partial transpose with respect to the system $B$, we find the
matrix

$\displaystyle(\rho_{AB})^{\intercal_{B}}=$
$\displaystyle\left(\begin{matrix}1-N_{s}\quad&0&0&0\\\
0&0\quad&\sqrt{(1-N_{s})N_{s}e^{-\gamma}\lambda}\quad&0\\\
0&\sqrt{(1-N_{s})N_{s}e^{-\gamma}\lambda}\quad&(1-\lambda)N_{s}\quad&0\\\
0&0&0&\lambda N_{s}\quad\end{matrix}\right)\,,$

whose eigenvalues are not all positive, i.e. the state $\rho_{AB}$ is not PPT
[52]. By exploiting the fact that any two-qubit state is distillable if and
only if it is not PPT [52], it follows that
$E_{d}\left(\operatorname{id}_{A}\otimes\pazocal{N}_{\lambda,\gamma}(\outerproduct{\Psi_{N_{s}}}{\Psi_{N_{s}}}_{AA^{\prime}})\right)>0\,$,
where $E_{d}$ is the distillable entanglement. On the other hand, from Lemma
S26 we have that

$\displaystyle K(\pazocal{N}_{\lambda,\gamma},N_{s})\geq
Q_{2}(\pazocal{N}_{\lambda,\gamma},N_{s})\geq
E_{d}\\!\left(\operatorname{id}_{A}\otimes\pazocal{N}_{\lambda,\gamma}(\outerproduct{\Psi_{N_{s}}}{\Psi_{N_{s}}}_{AA^{\prime}})\right)\,.$
(S64)

This concludes the proof for $N_{s}\in(0,1)$. Since the energy-constrained
capacities are monotonically non-decreasing in the energy constraint $N_{s}$,
the proof follows for any $N_{s}>0$. ∎

Since the state
$\operatorname{id}_{A}\otimes\pazocal{N}_{\lambda,\gamma}(\outerproduct{\Psi_{N_{s}}}{\Psi_{N_{s}}}_{AA^{\prime}})$
in (S63) is always entangled, it follows that the bosonic loss-dephasing
channel $\pazocal{N}_{\lambda,\gamma}$ is never entanglement breaking [33]. We
state this formally in the following theorem.

######  Theorem S39.

For all $\gamma\geq 0$ and all $\lambda\in(0,1]$, the bosonic loss-dephasing
channel $\pazocal{N}_{\lambda,\gamma}$ is not entanglement breaking.

##  IV Appendix

###  IV.1 Hadamard maps

Given an infinite matrix $A=(a_{mn})_{m,n\in\mathbb{N}},a_{mn}\in\mathbb{C}$,
we can introduce a superoperator $H$, recognised as the _Hadamard map_ , whose
action is defined as $H(\outerproduct{m}{n})=a_{n,m}\outerproduct{m}{n}$ for
all $n,m\in\mathbb{N}$. We are interested in establishing requirements for an
infinite matrix $A$ to ensure that the associated Hadamard map $H$ is a
quantum channel. We begin with some preliminaries. Let $\ell^{2}(\mathbb{N})$
be the space of square-summable complex-valued sequences defined as

$\displaystyle\ell^{2}(\mathbb{N})\coloneqq\left\\{x\coloneqq\\{x_{n}\\}_{n\in\mathbb{N}},x_{n}\in\mathbb{C}:\,\|x\|\coloneqq\sqrt{\sum_{n=0}^{\infty}|x_{n}|^{2}}<\infty\,\right\\}\,.$
(S65)

An infinite matrix $A\coloneqq(a_{mn})_{m,n\in\mathbb{N}},a_{mn}\in\mathbb{C}$
defines a linear operator on $\ell^{2}(\mathbb{N})$. The operator norm of $A$
is defined as follows:

$\displaystyle\|A\|_{\infty}\coloneqq\sup_{\begin{subarray}{c}x\in\ell^{2}(\mathbb{N})\\\
\|x\|=1\end{subarray}}\|Ax\|=\sup_{\begin{subarray}{c}\\{x_{n}\\}_{n\in\mathbb{N}},x_{n}\in\mathbb{C}\\\
\
\sum_{n=0}^{\infty}|x_{n}|^{2}=1\end{subarray}}\sqrt{\sum_{m=0}^{\infty}\left|\sum_{n=0}^{\infty}a_{mn}x_{n}\right|^{2}}.$

$A$ is said to be bounded if $\|A\|_{\infty}<\infty$. The following lemma,
referred to as _Schur test_ , gives a sufficent condition for an infinite
matrix to be bounded (e.g. [70, Page 24, Problem 45]).

######  Lemma S40.

Let $A\coloneqq(a_{mn})_{m,n\in\mathbb{N}},a_{mn}\in\mathbb{C}$, be an
infinite matrix. Suppose that there exist
$\\{p_{n}\\}_{n\in\mathbb{N}},p_{n}\in\mathbb{R}_{>0}$ and
$\\{q_{m}\\}_{m\in\mathbb{N}},q_{m}\in\mathbb{R}_{>0}$, and $\beta>0$, and
$\gamma>0$ such that

$\displaystyle\sum_{m=0}^{\infty}|a_{mn}|{p_{m}}\leq\beta
q_{n}\quad\text{and}\quad\sum_{n=0}^{\infty}|a_{mn}|{q_{n}}\leq\gamma
p_{m}\,,\quad\forall\,m,n\in\mathbb{N}\,.$

Then the matrix $A$ satisfies $\|A\|_{\infty}\leq\beta\gamma$. In particular,
$A$ is bounded.

By choosing $p_{n}=q_{n}=1$ and
$\gamma=\beta=\sup_{n\in\mathbb{N}}\sum_{m=0}^{\infty}|a_{mn}|$, we obtain the
following corollary:

######  Corollary S41.

Let $A=(a_{mn})_{m,n\in\mathbb{N}},a_{mn}\in\mathbb{C}$, be an infinite
Hermitian matrix. If $\sup_{n\in\mathbb{N}}\sum_{m=0}^{\infty}|a_{mn}|$ is
finite, then $A$ is bounded.

######  Lemma S42.

Let $A=(a_{mn})_{m,n\in\mathbb{N}},a_{mn}\in\mathbb{C}$, be a bounded
Hermitian infinite matrix. Then $A$ is positive semi-definite as an operator
on $\ell^{2}(\mathbb{N})$ if and only if
$A^{(d)}\coloneqq(a_{mn})_{m,n=0,1,\ldots,d-1}$ is positive semi-definite for
all $d\in\mathbb{N}$, where $A^{(d)}$ is the $d\times d$ top left corner of
$A$.

###### Proof.

Assume that $A^{(d)}$ is positive semi-definite for all $d\in\mathbb{N}$. Let
us pick an arbitrary $x\in\ell^{2}(\mathbb{N})$. It is known that for any
$\varepsilon>0$, there exists $d\in\mathbb{N}$ and
$y^{(d)}\coloneqq(y^{(d)}_{n})_{n\in\mathbb{N}},y^{(d)}_{n}\in\mathbb{C}$,
with $y^{(d)}_{n}=0$ for all $n>d$, such that $\|x-y^{(d)}\|<\varepsilon\,$.
Note that

$\displaystyle x^{\dagger}Ax$
$\displaystyle=\left(x-y^{(d)}\right)^{\dagger}A\,x+\left(y^{(d)}\right)^{\dagger}A\,(x-y^{(d)})+\left(y^{(d)}\right)^{\dagger}Ay^{(d)}$
$\displaystyle\stackrel{{\scriptstyle\mathclap{\scriptsize\mbox{(i)}}}}{{\geq}}-\|x-y^{(d)}\|\,\|A\|_{\infty}\,\left(\|x\|+\|y^{(d)}\|\right)+\left(y^{(d)}\right)^{\dagger}A^{(d)}\,y^{(d)}$
$\displaystyle\stackrel{{\scriptstyle\mathclap{\scriptsize\mbox{(ii)}}}}{{\geq}}-\varepsilon\,\|A\|_{\infty}\,\left(2\|x\|+\varepsilon\right)\,,$

where in (i) we applied Cauchy-Schwarz inequality twice and the definition of
infinity norm as follows

$\displaystyle|(x-y^{(d)})^{\dagger}A\,x|$
$\displaystyle\leq\|x-y^{(d)}\|\|A\,x\|\leq\|x-y^{(d)}\|\|A\|_{\infty}\|x\|\,,$
$\displaystyle|(y^{(d)})^{\dagger}A\,(x-y^{(d)})|$
$\displaystyle\leq\|y^{(d)}\|\|A\,(x-y^{(d)})\|\leq\|y^{(d)}\|\|A\|_{\infty}\|x-y^{(d)}\|\,,$

and, in (ii), we exploited triangular inequality to derive

$\displaystyle\|y^{(d)}\|\leq\|y^{(d)}-x\|+\|x\|\leq\varepsilon+\|x\|\,,$
(S66)

together with the fact that $A^{(d)}$ is positive semi-definite. Hence, since
$\varepsilon>0$ is arbitrary, we conclude that $x^{\dagger}Ax\geq 0$, meaning
that $A$ is positive semi-definite as an operator on $\ell^{2}(\mathbb{N})$. ∎

A square matrix is said to be diagonally dominant if

$\displaystyle\sum_{m\neq n}|a_{mn}|\leq|a_{nn}|,\quad\forall n.$

In words, a square matrix is said to be diagonally dominant if for every row
of the matrix, the absolute value of the diagonal entry in a row is larger
than or equal to the sum of the absolute values of all the other (non-
diagonal) entries in that row. Note that for Hermitian matrices one can
exchange row with column in this definition. We proceed to state a sufficient
condition for a matrix $A$ to be positive semi-definite as an operator on
$\ell^{2}(\mathbb{N})$. While the following sufficient condition is usually
stated for finite matrices, it can be generalised to infinite matrices as per
Lemma S42.

######  Lemma S43.

[44, Chapter 6] For any $d\in\mathbb{N}$, let
$A=(a_{mn})_{m,n=0,1,\ldots,d-1},a_{mn}\in\mathbb{C}$, be a $d\times d$
Hermitian matrix with $a_{nn}\in\mathbb{R}_{\geq 0}$ for all
$n\in\\{0,\ldots,d-1\\}$. $A$ is positive semi-definite if it is diagonally
dominant.

######  Lemma S44.

Let $A=(a_{mn})_{m,n\in\mathbb{N}},a_{mn}\in\mathbb{C}$ be an infinite
Hermitian matrix with $a_{nn}\in\mathbb{R}_{\geq 0},\forall n\in\mathbb{N}$.
Assume that $\sup_{n\in\mathbb{N}}a_{nn}$ is finite and that $A$ is diagonally
dominant. Then $A$ is bounded and positive semi-definite when seen as an
operator on $\ell^{2}(\mathbb{N})$.

###### Proof.

Since

$\displaystyle\sup_{n\in\mathbb{N}}\sum_{m=0}^{\infty}|a_{mn}|\leq
2\sup_{n\in\mathbb{N}}a_{n,n}<\infty\,,$

Corollary S41 implies that $A$ is bounded. The fact that $A$ is positive semi-
definite as an operator on $\ell^{2}(\mathbb{N})$ follows from Lemma S43
together with Lemma S42. ∎

We now present a lemma from the literature that establishes necessary and
sufficient conditions for an infinite matrix $A$ to ensure that the associated
Hadamard map $H$ is a quantum channel. We then use it to derive an explicit
sufficient condition for an infinite matrix $A$ to give rise to a cptp
Hadamard map.

######  Lemma S45.

[28, Lemma S4] Let $A\coloneqq(a_{mn})_{m,n\in\mathbb{N}},a_{mn}\in\mathbb{C}$
be a bounded infinite matrix. The following requirements establish the
necessary and sufficient conditions for the associated Hadamard map $H$ to
qualify as a quantum channel:

  1. (i)

$a_{nn}=1,\,\forall n\in\mathbb{N}$;

  2. (ii)

$A$ is positive semi-definite as on operator on $\ell^{2}(\mathbb{N})$.

As a consequence of Lemma S45 and Lemma S44, we obtain:

######  Lemma S46.

Let $A\coloneqq(a_{mn})_{m,n\in\mathbb{N}}$ be an infinite Hermitian matrix
that is diagonally dominant with $a_{n,n}=1$ for all $n\in\mathbb{N}$. In this
case, its associated Hadamard map $H$ is a quantum channel.

###  IV.2 Miscellaneous Lemmas

######  Lemma S47.

Let $\pazocal{H}_{A}$ and $\pazocal{H}_{B}$ be two Hilbert spaces and let
$\pazocal{N}:\pazocal{T}(\pazocal{H}_{A})\to\pazocal{T}(\pazocal{H}_{B})$ be a
quantum channel. For all normalised states
$\ket{\psi_{1}},\ket{\psi_{2}}\in\pazocal{H}_{A}$ and
$\ket{\phi_{1}},\ket{\phi_{2}}\in\pazocal{H}_{B}$ it holds that

$\displaystyle\left|\bra{\phi_{1}}\pazocal{N}(\outerproduct{\psi_{1}}{\psi_{2}})\ket{\phi_{2}}\right|\leq
1\,.$ (S67)

###### Proof.

It holds that

$\displaystyle\left|\bra{\phi_{1}}\pazocal{N}(\outerproduct{\psi_{1}}{\psi_{2}})\ket{\phi_{2}}\right|$
$\displaystyle\stackrel{{\scriptstyle\mathclap{\scriptsize\mbox{(i)}}}}{{\leq}}\|\pazocal{N}(\outerproduct{\psi_{1}}{\psi_{2}})\|_{\infty}$
(S68)
$\displaystyle\stackrel{{\scriptstyle\mathclap{\scriptsize\mbox{(ii)}}}}{{\leq}}\|\pazocal{N}(\outerproduct{\psi_{1}}{\psi_{2}})\|_{1}$
$\displaystyle\stackrel{{\scriptstyle\mathclap{\scriptsize\mbox{(iii)}}}}{{\leq}}\|\outerproduct{\psi_{1}}{\psi_{2}}\|_{1}$
$\displaystyle=1\,.$

Here, in (i), we exploited one of the definition of the operator norm in (S2).
In (ii), we exploited that the trace norm is always an upper bound on the
operator norm. Finally, in (iii), we leveraged the monotonicity of the trace
norm under quantum channels [60]. ∎

######  Lemma S48 ([62, 63]).

Let $\pazocal{H}_{A},\pazocal{H}_{A^{\prime}}$ be isomorphic Hilbert spaces,
possibly infinite dimensional. Let $\ket{\psi}_{A^{\prime}A}$ be a pure state
that satisfies
$\Tr_{A^{\prime}}\\!\left[\outerproduct{\psi}{\psi}_{AA^{\prime}}\right]>0$.
The generalised Choi–Jamiołkowski matrix defines an isomorphism between the
set of quantum channels from $\pazocal{H}_{A^{\prime}}$ to $\pazocal{H}_{B}$
and the set of bipartite states
$\sigma_{AB}\in\pazocal{P}\left(\pazocal{H}_{AB}\right)$ such that
$\Tr_{B}\sigma_{AB}=\Tr_{A^{\prime}}\\!\left[\outerproduct{\psi}{\psi}_{AA^{\prime}}\right]$.
Specifically, for any quantum channel $\pazocal{N}_{A^{\prime}\to
B}:\pazocal{T}(\pazocal{H}_{A^{\prime}})\to\pazocal{T}(\pazocal{H}_{B})$, it
holds that

$\displaystyle\pazocal{N}_{A^{\prime}\to
B}\\!\left(\outerproduct{e_{i}}{e_{j}}\right)=\frac{1}{\sqrt{\lambda_{i}\lambda_{j}}}\Tr_{A}\\!\left[\left(\outerproduct{e_{j}}{e_{i}}_{A}\otimes\mathbb{1}_{B}\right)\,\sigma_{AB}\right],\quad\forall\,i,j\in\mathbb{N}\,,$
(S69)

where $(\ket{e_{i}})_{i\in\mathbb{N}}$ and $(\lambda_{i})_{i\in\mathbb{N}}$
form a spectral decomposition of
$\Tr_{A^{\prime}}\\!\left[\outerproduct{\psi}{\psi}_{AA^{\prime}}\right]$,
i.e.
$\Tr_{A^{\prime}}\\!\left[\outerproduct{\psi}{\psi}_{AA^{\prime}}\right]=\sum_{i}\lambda_{i}\outerproduct{e_{i}}{e_{i}}_{A}\,$,
and where the state
$\sigma_{AB}\coloneqq\operatorname{id}_{A}\otimes\pazocal{N}_{A^{\prime}\to
B}(\outerproduct{\psi}{\psi}_{AA^{\prime}})$ is called the generalised Choi
state of $\pazocal{N}$. Eq. S69 is enough to specify the channel
$\pazocal{N}_{A^{\prime}\to B}$ completely, as the linear span of the
operators $\left(\outerproduct{e_{i}}{e_{j}}\right)_{i,j\in\mathbb{N}}$ (i.e.
the set of finite-rank operators) is dense in
$\pazocal{T}(\pazocal{H}_{A^{\prime}})$.

######  Lemma S49 ([42]).

Let $\pazocal{H}_{A},\pazocal{H}_{A^{\prime}},\pazocal{H}_{B}$ be isomorphic
Hilbert spaces, possibly infinite dimensional. Let $\pazocal{N}_{A^{\prime}\to
B}:\pazocal{T}(\pazocal{H}_{A^{\prime}})\to\pazocal{T}(\pazocal{H}_{B})$ be a
quantum channel. Let
$\ket{\psi}_{A^{\prime}A}\in\pazocal{H}_{A}\otimes\pazocal{H}_{A^{\prime}}$ be
a pure state such that the reduced state
$\Tr_{A^{\prime}}[\outerproduct{\psi}{\psi}_{AA^{\prime}}]$ is positive
definite. Then, $\pazocal{N}_{A^{\prime}\to B}$ is anti-degradable if and only
if the state $\operatorname{id}_{A}\otimes\pazocal{N}_{A^{\prime}\to
B}(\outerproduct{\psi}{\psi}_{AA^{\prime}})$ is two-extendible on $B$, meaning
that there exists a state
$\rho_{AB_{1}B_{2}}\in\pazocal{P}\left(\pazocal{H}_{A}\otimes\pazocal{H}_{B_{1}}\otimes\pazocal{H}_{B_{2}}\right)$
such that

$\displaystyle\Tr_{B_{2}}\left[\rho_{AB_{1}B_{2}}\right]$
$\displaystyle=\operatorname{id}_{A}\otimes\pazocal{N}_{A^{\prime}\to
B_{1}}(\outerproduct{\psi}{\psi}_{AA^{\prime}})\,,$ (S70)
$\displaystyle\Tr_{B_{1}}\left[\rho_{AB_{1}B_{2}}\right]$
$\displaystyle=\operatorname{id}_{A}\otimes\pazocal{N}_{A^{\prime}\to
B_{2}}(\outerproduct{\psi}{\psi}_{AA^{\prime}}),$

where $\pazocal{H}_{B_{1}}$ and $\pazocal{H}_{B_{2}}$ are Hilbert spaces that
are isomorphic to $\pazocal{H}_{B}$.

###### Proof.

Let $U_{A^{\prime}E\to BE}$ be a Stinespring dilation of the channel
$\pazocal{N}_{A^{\prime}\to B}$. Further assume that
$\pazocal{N}_{A^{\prime}\to B}$ is anti-degradable. By definition, there
exists a quantum channel $\pazocal{A}_{E\to B}$ such that $\pazocal{A}_{E\to
B}\circ\pazocal{N}^{\text{c}}_{A^{\prime}\to B}=\pazocal{N}_{A^{\prime}\to
B}$. Let us consider the tripartite state $\rho_{AB_{1}B_{2}}$, with
$B_{1},B_{2}$ being copies of $B$, defined as

$\displaystyle\rho_{AB_{1}B_{2}}=\operatorname{id}_{A}\otimes\operatorname{id}_{B_{1}}\otimes\pazocal{A}_{E\to
B_{2}}\left(U_{A^{\prime}E\to
B_{1}E}\,\big{(}\outerproduct{\psi}{\psi}_{AA^{\prime}}\otimes\outerproduct{0}{0}_{E}\big{)}\,U_{A^{\prime}E\to
B_{1}E}^{\dagger}\right)\,.$ (S71)

It holds that

$\displaystyle\Tr_{B_{2}}\left[\rho_{AB_{1}B_{2}}\right]$
$\displaystyle=\Tr_{B_{2}}\left[\operatorname{id}_{A}\otimes\operatorname{id}_{B_{1}}\otimes\pazocal{A}_{E\to
B_{2}}\left(U_{A^{\prime}E\to
B_{1}E}\,\big{(}\outerproduct{\psi}{\psi}_{AA^{\prime}}\otimes\outerproduct{0}{0}_{E}\big{)}\,U_{A^{\prime}E\to
B_{1}E}^{\dagger}\right)\right]$
$\displaystyle=\operatorname{id}_{A}\otimes\Tr_{E}\\!\left[U_{A^{\prime}E\to
B_{1}E}\big{(}\outerproduct{\psi}{\psi}_{AA^{\prime}}\otimes\outerproduct{0}{0}_{E}\big{)}U_{A^{\prime}E\to
B_{1}E}^{\dagger}\right]$
$\displaystyle=\operatorname{id}_{A}\otimes\pazocal{N}_{A^{\prime}\to
B_{1}}\left(\outerproduct{\psi}{\psi}_{AA^{\prime}}\right)\,,$

and

$\displaystyle\Tr_{B_{1}}\left[\rho_{AB_{1}B_{2}}\right]$
$\displaystyle=\Tr_{B_{1}}\left[\operatorname{id}_{A}\otimes\operatorname{id}_{B_{1}}\otimes\pazocal{A}_{E\to
B_{2}}\left(U_{A^{\prime}E\to
B_{1}E}\,\big{(}\outerproduct{\psi}{\psi}_{AA^{\prime}}\otimes\outerproduct{0}{0}_{E}\big{)}\,U_{A^{\prime}E\to
B_{1}E}^{\dagger}\right)\right]$
$\displaystyle=\operatorname{id}_{A}\otimes\pazocal{A}_{E\to
B_{2}}\left(\Tr_{B_{1}}\left[U_{A^{\prime}E\to
B_{1}E}\big{(}\outerproduct{\psi}{\psi}_{AA^{\prime}}\otimes\outerproduct{0}{0}_{E}\big{)}U_{A^{\prime}E\to
B_{1}E}^{\dagger}\right]\right)$
$\displaystyle=\operatorname{id}_{A}\otimes\pazocal{A}_{E\to
B_{2}}\circ\pazocal{N}^{\text{c}}_{A^{\prime}\to
B_{2}}\left(\outerproduct{\psi}{\psi}_{AA^{\prime}}\right)$
$\displaystyle=\operatorname{id}_{A}\otimes\pazocal{N}_{A^{\prime}\to
B_{2}}\left(\outerproduct{\psi}{\psi}_{AA^{\prime}}\right)\,.$

Now, let us establish the converse. Assume that there exists
$\rho_{AB_{1}B_{2}}$ which satisfies (S70). Let
$\ket{\Psi}_{AB_{1}B_{2}P}\in\pazocal{H}_{A}\otimes\pazocal{H}_{B_{1}}\otimes\pazocal{H}_{B_{2}}\otimes\pazocal{H}_{P}$
be a purification of $\rho_{AB_{1}B_{2}}$, with $\pazocal{H}_{P}$ being the
purifying Hilbert space. Note that both $\ket{\Psi}_{AB_{1}B_{2}P}$ and
$U_{A^{\prime}E\to
B_{1}E}\,\big{(}\ket{\psi}_{AA^{\prime}}\otimes\ket{0}_{E}\big{)}$ are
purifications of $\operatorname{id}_{A}\otimes\pazocal{N}_{A^{\prime}\to
B_{1}}\left(\outerproduct{\psi}{\psi}_{AA^{\prime}}\right)$, with
$\pazocal{H}_{B_{2}}\otimes\pazocal{H}_{P}$ and $\pazocal{H}_{E}$ being their
purifying Hilbert spaces, respectively. It follows that [60] there exists an
isometry $V_{E\to
B_{2}P}:\pazocal{H}_{E}\to\pazocal{H}_{B_{2}}\otimes\pazocal{H}_{P}$ such that
$V_{E\to B_{2}P}\,U_{A^{\prime}E\to
B_{1}E}\,\big{(}\ket{\psi}_{AA^{\prime}}\otimes\ket{0}_{E}\big{)}=\ket{\Psi}_{AB_{1}B_{2}P}\,.$
Hence, the quantum channel $\pazocal{A}_{E\to
B_{2}}:\pazocal{T}(\pazocal{H}_{E})\to\pazocal{T}(\pazocal{H}_{B_{2}})$,
defined by $\pazocal{A}_{E\to B_{2}}(\cdot)=\Tr_{P}\left[V_{E\to
B_{2}P}(\cdot)V_{E\to B_{2}P}^{\dagger}\right]$, satisfies that

$\displaystyle\operatorname{id}_{A}\otimes\pazocal{A}_{E\to
B_{2}}\circ\pazocal{N}^{\text{c}}_{A^{\prime}\to
B_{2}}(\outerproduct{\psi}{\psi}_{AA^{\prime}})$
$\displaystyle=\operatorname{id}_{A}\otimes\pazocal{A}_{E\to
B_{2}}\left(\Tr_{B_{1}}\left[U_{A^{\prime}E\to
B_{1}E}(\outerproduct{\psi}{\psi}_{AA^{\prime}}\otimes\outerproduct{0}{0}_{E})U_{A^{\prime}E\to
B_{1}E}^{\dagger}\right]\right)$
$\displaystyle=\Tr_{B_{1}P}\left[\outerproduct{\Psi}{\Psi}_{AB_{1}B_{2}P}\right]$
$\displaystyle=\Tr_{B_{1}}\left[\rho_{AB_{1}B_{2}}\right]$
$\displaystyle=\operatorname{id}_{A}\otimes\pazocal{N}_{A^{\prime}\to
B_{2}}\left(\outerproduct{\psi}{\psi}_{AA^{\prime}}\right)\,.$

Consequently, since the pure state $\ket{\psi}_{A^{\prime}A}$ satisfies
$\Tr_{A^{\prime}}\\!\left[\outerproduct{\psi}{\psi}_{AA^{\prime}}\right]>0$,
Lemma S48 implies that $\pazocal{A}_{E\to
B_{2}}\circ\pazocal{N}^{\text{c}}_{A^{\prime}\to
B_{2}}=\pazocal{N}_{A^{\prime}\to B_{2}}$, meaning that
$\pazocal{N}_{A^{\prime}\to B_{2}}$ is anti-degradable. ∎

######  Lemma S50.

Let
$\pazocal{N},\pazocal{M}:\pazocal{T}(\pazocal{H}_{S})\to\pazocal{T}(\pazocal{H}_{S})$
be quantum channels. If either $\pazocal{M}$ or $\pazocal{N}$ is anti-
degradable, then the composition $\pazocal{M}\circ\pazocal{N}$ is anti-
degradable. Specifically, let $E_{1}$ and $E_{2}$ be the Stinespring
environments of $\pazocal{N}$ and $\pazocal{M}$, respectively. If
$\pazocal{N}$ is anti-degradable with anti-degrading map
$\pazocal{A}_{E_{1}\to S}$, then $(\pazocal{M}\circ\pazocal{A}_{E_{1}\to
S})\otimes\Tr_{E_{2}}$ is an anti-degrading map of
$\pazocal{M}\circ\pazocal{N}$. Analogously, if $\pazocal{M}$ is anti-
degradable with anti-degrading map $\pazocal{A}_{E_{2}\to S}$, then
$\Tr_{E_{1}}\otimes\pazocal{A}_{E_{2}\to S}$ is an anti-degrading map of
$\pazocal{M}\circ\pazocal{N}$.

###### Proof.

Let $V^{S\to SE_{1}}$ and $W^{S\to SE_{2}}$ be Stinespring isometries
associated with $\pazocal{N}$ and $\pazocal{M}$, respectively. By considering
the following complementary channel of $\pazocal{M}\circ\pazocal{N}$,

$\displaystyle(\pazocal{M}\circ\pazocal{N})^{\text{c}}(\rho)=\Tr_{S}\\!\left[W^{S\to
SE_{2}}V^{S\to SE_{1}}\,\rho\,\left(V^{S\to
SE_{1}}\right)^{\dagger}\left(W^{S\to
SE_{2}}\right)^{\dagger}\right],\qquad\forall\,\rho\in\pazocal{T}(\pazocal{H}_{S})\,,$

one can easily check that if $\pazocal{N}$ is anti-degradable with anti-
degrading map $\pazocal{A}_{E_{1}\to S}$, then

$\displaystyle\left[(\pazocal{M}\circ\pazocal{A}_{E_{1}\to
S})\otimes\Tr_{E_{2}}\right]\circ(\pazocal{M}\circ\pazocal{N})^{\text{c}}=\pazocal{M}\circ\pazocal{N}\,.$
(S72)

Analogously, one can easily verify that if $\pazocal{M}$ is anti-degradable
with anti-degrading map $\pazocal{A}_{E_{2}\to S}$, then

$\displaystyle\left[\Tr_{E_{1}}\otimes\pazocal{A}_{E_{2}\to
S}\right]\circ(\pazocal{M}\circ\pazocal{N})^{\text{c}}=\pazocal{M}\circ\pazocal{N}\,.$
(S73)

∎


# Dynamic modeling of mortality via mixtures of skewed distribution functions

Emanuele Aliverti  _Department of Economics, University Ca’ Foscari Venezia_
Stefano Mazzuco  _Department of Statistical Sciences, University of Padova_
Bruno Scarpa  _Department of Statistical Sciences, University of Padova_
_Department of Mathematics “Tullio Levi-Civita”, University of Padova_

###### Abstract

There has been growing interest on forecasting mortality. In this article, we
propose a novel dynamic Bayesian approach for modeling and forecasting the
age-at-death distribution, focusing on a three-components mixture of a Dirac
mass, a Gaussian distribution and a Skew-Normal distribution. According to the
specified model, the age-at-death distribution is characterized via seven
parameters corresponding to the main aspects of infant, adult and old-age
mortality. The proposed approach focuses on coherent modeling of multiple
countries, and following a Bayesian approach to inference we allow to borrow
information across populations and to shrink parameters towards a common mean
level, implicitly penalizing diverging scenarios. Dynamic modeling across
years is induced through an hierarchical dynamic prior distribution that
allows to characterize the temporal evolution of each mortality component and
to forecast the age-at-death distribution. Empirical results on multiple
countries indicate that the proposed approach outperforms popular methods for
forecasting mortality, providing interpretable insights on its evolution.

Keywords: Bayesian inference; Dynamic modeling; Mixture model; Skew-Normal
distribution.

## 1 Introduction

Figure 1: Age-at-death distribution across Italy, France, Sweden and Norway
from $1970$ to $2015$.

In the last decades, changes in life expectancy and population growths have
stimulated growing interest in sophisticated models for mortality (e.g.,
Booth, , 2020; Zanotto et al., , 2020). These approaches characterize the
observed age-specific regularities of the mortality function and their trends
over time, providing a parsimonious representation of the components of
mortality and allowing to provide forecasts. Estimates of past mortality
trends provide insights on the global health of a population, describe
relevant aspects of its demographic structure and facilitate comparisons
across countries. Similarly, probabilistic forecasts are routinely employed by
researchers, practitioners and organizations to aid and support welfare
planning strategies such as pension funds, insurances and health systems,
among many others (e.g., Lutz and Samir, , 2010).

Human mortality has been generally characterized by several quantities, with
the hazard function, the survival function and the probability density
function being the most popular ones (e.g., Klein and Moeschberger, , 2006).
Although each representation can be directly obtained from the others, their
interpretations are complementary and highlight more clearly different aspects
of mortality (e.g., Basellini and Camarda, , 2019). In this article we will
focus on modeling the distribution of the ages at death, generally referred to
as “age-at-death distribution” and corresponding to the density function that
characterizes mortality (Keyfitz, , 2005). In the literature on modeling and
forecasting mortality, this quantity has received significantly less attention
than other quantities such as mortality rates (e.g., Lee and Carter, , 1992;
Li and Lee, , 2005). However, as outlined in Basellini and Camarda,  (2020)
and Pascariu et al.,  (2019), analysis of the age-at-death distribution can be
incredibly useful, since it provides a direct evaluation of longevity and
lifespan variability that cannot be directly inferred from other quantities
(see also Cheung et al., , 2005; Canudas-Romo, , 2010).

To further see that, Figure 1 illustrates the age-specific distribution of
deaths from four illustrative countries, over a time period that ranges from
$1970$ to $2015$. These data highlight crucial aspects of the temporal
evolution of mortality, such as the global shift of the curve to older
ages—due to the increase in life expectancy—and the drop of perinatal and
infant mortality. Figure 1 also suggests that such an evolution has
characterized countries differently, resulting in faster changes in some
countries (e.g., Italy) and slower in others (e.g., Sweden). Furthermore, the
age-specific distribution of deaths illustrates that several mortality trends
are converging, resulting in similar patterns in terms of modal age at death,
infant mortality and compression of old age mortality.

The first approaches for modeling the distribution of deaths date back to the
19-th century (e.g., Gompertz, , 1825; Makeham, , 1860; Pearson, , 1897), and
have been generalized in the following century (Siler, , 1983; Heligman and
Pollard, , 1980). The main interest of such methods has been on characterizing
mortality patterns across ages, focusing on a given country during a specific
year. According to these approaches, the mortality curve can be characterized
using different parametric specifications, accounting for distinct aspects of
the mortality curve; see also Dellaportas et al.,  (2001) and Mazzuco et al.,
(2018). One notable limitation of these approaches is the focus on single-
country modeling; when interest is on forecasting mortality for multiple
countries, independent country-specific analysis are required. This strategy
ignores the fact that in most developed nations similar trends are observed,
and this abundance of information can be appropriately included in the model
to improve the quality of the estimates. Joint modeling of multiple nations,
in fact, is expected to improve the estimation of country-specific parameters
and temporal trends, allowing to share information across different countries.

More recently, researchers have focused on explicitly modeling the temporal
evolution of the age-at-death distribution, forecasting the parameters using
time-series specifications (Pascariu et al., , 2019; Basellini and Camarda, ,
2019). These methods borrow ideas from death-rates forecasting, such as the
milestone approach introduced by Lee and Carter,  (1992) and sequent
developments (e.g., Hyndman and Ullah, , 2007). Potentially, these models can
be both used to describe the evolution of mortality in the past years and to
provide forecasts, although the focus of applications is more frequently on
forecasting only. However, separate forecasts often lead to unreasonable
results, with country-specific or sex-specific predictions that diverge
substantially across groups. Such divergent forecasts across sub-populations
are implausible, since industrialized countries have been showing global
convergence trends in terms of mortality levels, and this tendency is expected
to endure (Wilson, , 2001; Bergeron-Boucher et al., , 2017; Oeppen, , 2008).

Following the breakthrough approach of Li and Lee,  (2005), different authors
have focused on “coherent” mortality forecasts, defined as predictions that do
not diverge across groups. Such a strategy has constantly received attention
in the past years, being successfully implemented in several models for death
rates; for example, generalizing compositional-data techniques (Bergeron-
Boucher et al., , 2017) and functional-data models (Hyndman et al., , 2013)
for mortality. However, a critical issue with these methods concerns with the
choice of a population standard to be used as a reference. This choice is
substantially arbitrary, and a common consensus on which sub-population should
be used is still lacking, since results are often sensitive to different
references (Booth, , 2020; Kjærgaard et al., , 2016).

Motivated by the above discussion, in this article we propose a novel approach
for modeling and forecasting the age-at-death distributions of multiple
countries. We will focus on the age-distribution of deaths, combining the life
tables quantities with the total number of deaths. This procedure allows to
naturally take into account the population size of each country and its age-
structure, which are important determinants of the age-at-death distribution.
The proposed methods can also be applied, without modifications, directly to
the life table quantities, where we can assume that data are provided on
artificial populations with fixed population size. However, this would discard
an important source of uncertainty, which corresponds to the population size
of a specific nation. Indeed, imposing a fixed same sample size weights each
nation equally, implying that a relatively small country (e.g. Sweden) would
have the same influence as a much bigger country (e.g., Germany or United
Kingdom). Therefore, we prefer to use the actual sample size of each county to
have a proper quantification of uncertainty and genuinely account for
different population structures.

We propose a method based on a mixture distribution function of a Dirac mass,
a Gaussian distribution and a Skew-Normal distribution, in order to
characterize infant, adult and old-age mortality with an interpretable set of
parameters. Following a Bayesian approach to inference, we specify
hierarchical priors for the dynamic parameters and the country-specific
coefficients, relying on informative common prior specifications to improve
the borrowing of information across nations. Compared to the available
techniques, our method allows to model directly the temporal evolution and the
country-specific variations across multiple sub-populations, avoiding
sequential procedures and providing estimates shrunk towards common levels.
Under the proposed methods, coherent predictions are naturally obtained
without selecting a reference population or an external standard, but instead
estimating these quantities from the available data. We will refer to the
proposed dynamic model based on a skewed distribution function as dysm in the
sequel.

## 2 Methods

### 2.1 Model specification

Our interest is on providing a flexible representation of the age-at-death
distribution in multiple countries, dynamically across different years or time
periods. We will characterize this quantity through a continuous probability
density function (pdf) $f(x;\boldsymbol{\vartheta}_{jt})$ which parametrizes
the probability of dying at age $x$ – based on the individuals current living
at age $x$ – in country $j$ during year $t$ via a set of country- and time-
specific parameters $\boldsymbol{\vartheta}_{jt}$, with $x=0,\dots,110$,
$j=1,\dots,p$ countries and $t=1,\dots,T$ years. Although $f(x;\cdot)$ is
often specified as a continuous function, in practice we focus on
characterizing the age-at-death distribution on a discrete grid, corresponding
to the probability $p_{x\,jt}$ of dying at a specific age $x=0,\dots,110$.
Such an aim can be addressed discretizing $f(x;\boldsymbol{\vartheta}_{jt})$
over a set of disjoint intervals $(x-1/2,x+1/2]$ as

$p_{x\,jt}=\int_{x-1/2}^{x+1/2}f(z;\boldsymbol{\vartheta}_{jt})\,\mathop{}\\!\mathrm{d}z=F(x+1/2;\boldsymbol{\vartheta}_{jt})-F(x-1/2;\boldsymbol{\vartheta}_{jt}),\quad
x=0,\dots,110,$ (1)

where $F$ denotes the cumulative distribution function (cdf) associated with
$f$. This characterization of probabilities via a discretization of an
underlying pdf automatically guarantees that the vector
$\mathbf{p}_{jt}=\mathbf{p}_{jt}(\boldsymbol{\vartheta}_{jt})=[p_{0\,jt},\dots,p_{110\,jt}]$
with elements as in Equation 1 is a proper probability vector with
${p_{x\,jt}\in[0,1]}$ and ${\sum_{x}p_{x\,jt}=1}$, therefore avoiding further
restrictions (e.g., Oeppen, , 2008); see the left panel of Figure 2 for a
graphical representation of the discretization process. We denote as
${\mathbf{D}_{jt}=[D_{0\,jt},\dots,D_{110\,jt}]}$ the number of deaths and as
$\mathbf{p}_{jt}=[p_{0\,jt},\dots,p_{110\,jt}]$ the death probabilities, where
elements $D_{x\,jt}$ and $p_{x\,jt}$ denote, respectively, the number of
deaths and the probability of dying at age $x$ – based on the current
individuals living at age $x$ – in country $j$ during year $t$. Given the
total number of deaths ${n_{jt}=\sum_{x}D_{x\,jt}}$ and the probabilities
outlined in Equation 1, we model the vector $\mathbf{D}_{jt}$ as a multinomial
distribution with

$(\mathbf{D}_{jt}\mid
n_{jt},\mathbf{p}_{jt})\sim\mbox{{multinom}}(n_{jt},\mathbf{p}_{jt}),\quad
j=1,\dots,p,\quad t=1,\dots,T.$ (2)

A crucial aspect of this approach involves the specification $f(x;\cdot)$,
which induces a model on the vector of probabilities $\mathbf{p}_{jt}$. In
practice, different parametric models for $f(x;\cdot)$ are available in the
literature, with the Gompertz,  (1825) and Siler,  (1983) models being popular
options. The age-at-death distribution $f(x;\cdot)$ could be also modeled
semi-parametrically, using for example splines or other basis functions
(Hyndman et al., , 2013; Basellini and Camarda, , 2020). However, when
interest is explicitly on characterizing relevant aspects of mortality – such
as the impact of infant mortality or the shape of old-age mortality –
parametric models provide a much clearer interpretation and deeper insights on
the structure of the populations under investigation. In addition, mortality
can be naturally divided into different components, accounting for deaths at
different phases of life (e.g., Lexis, , 1879; Pearson, , 1897). Therefore, it
is desirable to rely on a specification which preserves such a decomposition,
for example using mixtures of distributions (Carriere, , 1992; Mazzuco et al.,
, 2018).

In this work, we draw inspiration from Mazzuco et al.,  (2018) and Zanotto et
al.,  (2020) and rely on a mixture of three components to characterize
mortality; see also Carriere,  (1992) for related arguments. Specifically, in
order to facilitate model’s interpretation, we decompose mortality as the
contribution of infant mortality, adult mortality and old-age mortality.
Infant mortality is accounted through a single Dirac mass at age $0$, while
adult mortality is modeled via a Gaussian distribution with mean
$\mu_{jt}\in\mathbb{R}$ and standard deviation $\sigma_{jt}\in\mathbb{R}^{+}$.
Old-age mortality, instead, is characterized using a Skew-Normal distribution
(e.g., Azzalini and Capitanio, , 2013) with location parameter
$\xi_{jt}\in\mathbb{R}$, scale parameter $\omega_{jt}\in\mathbb{R}^{+}$ and
shape parameter $\alpha_{jt}\in\mathbb{R}$. The impact of each component to
the overall mortality function is accounted assigning mixture weights
$\pi_{h\,jt}\in[0,1]$ such that ${\pi_{0\,jt}+\pi_{1\,jt}+\pi_{2\,jt}=1}$ in
order to guarantee that the resulting expression remains a proper pdf. This
choice leads to the following characterization of
$f(x;\boldsymbol{\vartheta}_{jt})$

$\displaystyle
f(x;\boldsymbol{\vartheta}_{jt})=\pi_{0\,jt}\mathds{1}[x=0]+\pi_{1\,jt}\frac{1}{\sigma_{jt}}\phi\left(\frac{x-\mu_{jt}}{\sigma_{jt}}\right)+\pi_{2\,2t}\frac{2}{\omega_{jt}}\phi\left(\frac{x-\xi_{jt}}{\omega_{jt}}\right)\Phi\left(\alpha_{jt}\frac{x-\xi_{jt}}{\omega_{jt}}\right),$
(3)

with
$\boldsymbol{\vartheta}_{jt}=\left(\pi_{1\,jt},\pi_{2\,jt},\mu_{jt},\sigma_{jt},\xi_{jt},\omega_{jt},\alpha_{jt}\right)$,
$\mathds{1}[\cdot]$ denoting the indicator function and $\phi(\cdot)$,
$\Phi(\cdot)$ denoting the pdf and the cdf of a standard Gaussian,
respectively; we omit the parameter $\pi_{0\,jt}$ from
$\boldsymbol{\vartheta}_{jt}$ since $\pi_{0\,jt}=1-\pi_{1\,jt}-\pi_{2\,jt}$
due to the constraint on the mixture weights; see the right panel of Figure 2
for a graphical representation of the proposed model for the age-at-death
distribution.

Figure 2: Proposed mixture density function. Panel (a) provides a sketch of
the discretization of the pdf into a set of probabilities $p_{j}$, with
$j=57,\dots,60$. Panel (b) highlights the three components of the mixture
density function.

Equation 3 allows to characterize the age-at-death distribution with
flexibility, while also retaining interpretability of the mortality structure.
The first component of Equation 3 accounts for infant and perinatal mortality,
characterizing deaths in the first year of life. As outlined, for example, in
Figure 1, perinatal mortality is notably decreased during recent years, mainly
due to technological advances (e.g., Wilson, , 2001; Zanotto et al., , 2020);
the mixture weight $\pi_{0\,jt}$ characterizes the impact of this component in
the overall mortality. The use of a Dirac mass is particularly adequate for
developed countries, where infant mortality occurs during the first year of
life and can be properly modeled with a single mass, avoiding more elaborate
specifications and additional parameters. However, dysm can be naturally
adapted to countries with high infant-mortality levels, replacing the Dirac
mass with an half-normal distribution (Zanotto et al., , 2020) and leveraging
on the same model architecture; in Section 3.4, we compare such a
specification with the proposed model. The Gaussian component, instead, is
assigned to a relative weight $\pi_{1\,jt}$ and characterizes adult mortality
via the location parameter $\mu_{jt}$ and the scale parameters $\sigma_{jt}$.
In addition, these quantities correspond also to the moments of the Gaussian
component, and therefore they can be interpreted directly as the mean and
standard deviation of the age at adult mortality.

The Skew-Normal component, lastly, accounts for the asymmetric shape of old-
age mortality, which has been frequently observed in developed countries
(e.g., Mazzuco et al., , 2018; Pascariu et al., , 2019). Since most deaths
occur at old ages, this component is expected to be associated to the largest
mixture weight $\pi_{2\,jt}$. Inference on the structure of old-age mortality
rely on such mixture weight and on the Skew-Normal parameters
$\xi_{jt},\omega_{jt}$ and $\alpha_{jt}$, providing insights on the evolution
of old-age mortality in terms of location, scale and shape of such a
component. In practice, it can be convenient to focus on the mean and standard
deviation of old-age mortality, which are more interpretable measures of the
shift and compression of the age-at-death distribution. According to the
proposed parametrization, these quantities can be conveniently expressed as

$\tilde{\mu}^{\mbox{ \tiny SN
}}_{jt}=\xi_{jt}+\omega_{jt}\delta_{jt}\sqrt{2/\pi},\quad\tilde{\sigma}^{\mbox{
\tiny SN }}_{jt}=\omega_{jt}\sqrt{1-2\delta^{2}_{jt}/\pi}\quad$ (4)

with $\delta_{jt}=\alpha_{jt}/\sqrt{1+\alpha_{jt}}$ (e.g., Azzalini and
Capitanio, , 2013, sec. 2.1.4). Alternatively, it is also possible to
parameterize directly the Skew-Normal distribution in terms of its cumulants
(e.g., Azzalini and Capitanio, , 2013, sec. 3.1.4). However, for simplicity in
implementation and exposition, we prefer to focus on the outlined parameters,
and eventually post-process the estimates to focus on the functionals of
interest.

### 2.2 Hierarchical dynamic model

The Multinomial model outlined in Expression 2, with probabilities given by
Equations 1 and 3, provides a flexible specification for the distribution of
the number of deaths occurred in a single country $j$ during a specific year
$t$. In order to account for the temporal variation of the age-at-death
distribution, we model explicitly the dynamic evolution of the parameters
$\boldsymbol{\vartheta}_{jt}$ through an hierarchical dynamic model. This
strategy will allow us to highlight what aspects of mortality have changed in
the past years, characterizing the rates of such changes and facilitating
comparison across different countries. In addition, the inclusion of a dynamic
component facilitates the developments of mortality forecasts extrapolating
parameters in time, obtaining predictions for the future age-at-death
distributions.

We specify the hierarchical model in an unconstrained domain, transforming the
parameters of Equation 3. This choice facilitates computation and
implementations, allowing to rely on Gaussian dynamic models in the
transformed space and avoiding restrictions on the support of the parameters.
For this purpose, we introduce the reparametrized vector
${\boldsymbol{\tilde{\vartheta}}}_{jt}$, with elements
$(\tilde{\vartheta}_{jt\,1},\dots,\tilde{\vartheta}_{jt\,7})$ corresponding to

$\widetilde{\boldsymbol{\vartheta}}_{jt}=(\tilde{\vartheta}_{jt\,1},\dots,\tilde{\vartheta}_{jt\,7})\coloneqq\left(\log\left(\frac{\pi_{1\,jt}}{1-\pi_{1\,jt}}\right),\log\left(\frac{\pi_{2\,jt}}{1-\pi_{1\,jt}-\pi_{2\,jt}}\right),\mu_{jt},\log(\sigma_{jt}),\xi_{jt},\log(\omega_{jt}),\alpha_{jt}\right).$
(5)

In Equation 5, the mixture weights $(\pi_{1\,jt},\pi_{2\,jt})$ are mapped into
$(\tilde{\vartheta}_{jt\,1},\tilde{\vartheta}_{jt\,2})$ leveraging a
cumulative logit transformation, while $\sigma_{jt}$ and $\omega_{jt}$ –
corresponding to adult and old-age mortality scale parameters – are
transformed into $\tilde{\vartheta}_{jt\,4}$ and $\tilde{\vartheta}_{jt\,6}$,
respectively, via a logarithmic transformations. Note that each
$\tilde{\vartheta}_{jt\,k}\in\mathbb{R}$ for $k=1,\dots,7$, and that the
vector $\widetilde{\boldsymbol{\vartheta}}_{jt}$ effectively corresponds to a
one-to-one reparametrization of ${\boldsymbol{\vartheta}}_{jt}$.

We model the dynamic evolution of the elements $\tilde{\vartheta}_{jt\,k}$
leveraging a random-walk model with drift. This choice leads to the following
specification.

$\begin{split}\tilde{\vartheta}_{j0\,k}&\sim\mathcal{N}(m_{k},s_{k}^{2})\\\
(\tilde{\vartheta}_{jt\,k}\mid\tilde{\vartheta}_{j\,t-1\,k},\beta_{j\,k},\eta_{j\,k}^{2})&\sim\mathcal{N}(\beta_{j\,k}+\tilde{\vartheta}_{j\,t-1\,k},\eta_{j\,k}^{2}),\quad
t=1,\dots,T,\end{split}$ (6)

In Equation 6, $\tilde{\vartheta}_{j0\,k}\in\mathbb{R}$ denotes the initial
condition of $\tilde{\vartheta}_{jt\,k}$, and it is assigned to a Gaussian
distribution with mean $m_{k}\in\mathbb{R}$ and standard deviation
$s_{k}\in\mathbb{R^{+}}$. These values correspond to the initial information
for the dynamic hierarchical model, and characterize a probabilistic
representation of our beliefs before the data are observed (e.g., West and
Harrison, , 1997). The drift parameter $\beta_{j\,k}\in\mathbb{R}$ measures
the expected difference across consecutive values of
$\tilde{\vartheta}_{jt\,k}$ and accounts for possible increasing or decreasing
linear trends in the evolution of the parameters. The inclusion of this
component is particularly relevant in our problem, where we expected to
observe clear and important trends in the evolution of mortality, accounted by
the parameters charactering the age-at-death distribution. Conditionally on
the value at time ${(t-1)}$, each $\tilde{\vartheta}_{jt\,k}$ is recursively
distributed as a Gaussian distribution with mean
$\beta_{j\,k}+\tilde{\vartheta}_{j\,t-1\,k}$, and standard deviation
$\eta_{j\,k}\in\mathbb{R^{+}}$. Overall, we can interpret dysm as a
Multinomial state-space model, where Equation 6 controls the evolution of the
latent states via random-walks with Gaussian innovations and, conditional on
the states, the observed vectors of deaths are modeled as independent
Multinomial variables with probabilities given by Equations 3 and 1; see
Figure 3 for a graphical representation of the proposed hierarchical dynamic
model.

Although this dynamic model might seem oversimplistic, we found that such a
specification leads to good results in a large variety of applications
involving developed countries; see Section 3 for further details. More
elaborate extensions are possible, for example including a moving-average
component, increasing the lag of the autoregressive terms or joint modeling
multiple components of $\tilde{\vartheta}_{jt\,k}$ via multivariate random-
walks. However, a random walk with drift on the latent states often
characterizes well many phenomena without over-parameterizing the dynamic
component, and is also general enough to be directly used, without further
modifications, in a large variety of settings; see Durbin and Koopman,  (2012)
and West and Harrison,  (1997) for related arguments and Section 3.4 for a
comparison with alternative specifications.

Additionally, we highlight that Equation 6 is valid for any equally-spaced
time grid, and therefore the proposed approach can be also applied to
different time scales. In fact, our description focuses on yearly data for
simplicity in exposition, but this requirement is not necessary for developing
the proposed model in practice. Several important settings involve analysis of
mortality across different time scales, such as monthly data or multi-year
periods; for example, when interest is on modeling short-term variations due
to extraordinary events (e.g., covid-19 pandemic) or large-scale trends,
respectively. The general structure of dysm facilitates its use in these
different settings, without modifying the model specification.

${\boldsymbol{\tilde{\vartheta}}}_{j1}$
${\boldsymbol{\tilde{\vartheta}}}_{j2}$ $\cdots$
${\boldsymbol{\tilde{\vartheta}}}_{jt}$ $\cdots$
${\boldsymbol{\tilde{\vartheta}}}_{j_{\,T-1}}$
${\boldsymbol{\tilde{\vartheta}}}_{j_{\,T}}$ $\mathbf{p}_{j1}$
$\mathbf{p}_{j2}$ $\cdots$ $\mathbf{p}_{jt}$ $\cdots$ $\mathbf{p}_{j_{\,T-1}}$
$\mathbf{p}_{j_{\,T}}$ $\mathbf{D}_{j1}$$\mathbf{D}_{j2}$$\cdots$
$\mathbf{D}_{jt}$$\cdots$ $\mathbf{D}_{j_{\,T-1}}$$\mathbf{D}_{j_{\,T}}$for
country $j=1,\ldots,p$ Figure 3: Graphical illustration of the dynamic
component of dysm. Solid circles, white squares and gray squares denote
vectors of dynamic parameters, probabilities characterizing the age-at-death
distribution and observed number of deaths, respectively.

### 2.3 Bayesian inference

As discussed in Section 1, a crucial aspect of mortality forecasting involves
the penalization of diverging scenarios via coherent modeling. This can be
obtained explicitly including a common trend across countries (e.g., Li and
Lee, , 2005; Li, , 2013), or choosing a reference population (e.g., Hyndman et
al., , 2013; Basellini and Camarda, , 2019); see also Booth,  (2020) for a
recent discussion. However, these choices are arbitrary and significantly
affect the final results (e.g., Kjærgaard et al., , 2016). Therefore, we
follow a different path and induce coherent predictions via a Bayesian
approach to inference, regularizing estimates toward a common mean trend. More
specifically, we introduce a common prior specification for the country-
specific dynamic hierarchical model defined in Equation 6, for each
coefficient $k=1,\dots,7$. Including the prior for
$\tilde{\vartheta}_{j0\,k}$, we let

$\begin{split}\tilde{\vartheta}_{j0\,k}&\sim\mathcal{N}(m_{k},s_{k}^{2})\\\
\beta_{j\,k}\overset{\mbox{\tiny
iid}}{\sim}\mathcal{N}(m_{\beta_{k}},s_{\beta_{k}}),&\quad\eta_{j\,k}^{2}\overset{\mbox{\tiny
iid}}{\sim}\mbox{Inv-
Gamma}(\mbox{a}_{k},\mbox{b}_{k})\quad\mbox{for}\,\,j=1,\dots,p,\end{split}$
(7)

where $m_{\beta_{k}},s_{\beta_{k}},\mbox{a}_{k}$ and $\mbox{b}_{k}$ denote
fixed hyper-parameters; see Figure 4 for an illustration of the hierarchical
dependence structure induced by the common prior specification.

The specification of a common prior distribution for the country-specific
parameters allows to share information across different nations, improving
estimates for the age-at-death distributions. In addition, the introduction of
a common prior shrinks estimates toward a common mean, providing an implicit
regularization on the parameters. Such a shrinkage is what makes dysm a
“coherent” model, using a data-driven reference mortality level without the
need of a user-defined benchmark.

$\beta_{1\,k},\eta_{1\,k},\left(\tilde{\vartheta}_{1t\,k}\right)_{t=0}^{T}$
$\cdots$
$\beta_{p\,k},\eta_{p\,k},\left(\tilde{\vartheta}_{pt\,k}\right)_{t=0}^{T}$
$m_{k},s_{k},m_{\beta_{k}},s_{\beta_{k}},a_{k},b_{k}$$k=1,\ldots,7$ Figure 4:
Graphical illustration of the hierarchical component of dysm, inducing
borrowing of information and coherent predictions across countries, with
$j=1,\dots,p$. Top level block denotes fixed hyper-parameters, while bottom
blocks denote model’s parameters.

### 2.4 Prior elicitation and posterior computation

The large abundance of research on the evolution of mortality, and the clear
interpretation of the parameters in Equation 3, facilitate informative
elicitation for the hyper-parameters of Equation 7. For example, European
adult mortality is expected to lie between $30$ and $70$ years with high
probability (e.g., Lewer et al., , 2020); therefore, we can assume
$\tilde{\vartheta}_{j0\,3}\sim\mathcal{N}(50,10)$ for the initial condition of
$\tilde{\vartheta}_{jt\,3}$, independently for $j=1,\dots,p$. Similarly,
informative specifications for the prior distribution on the drift parameter
$\beta_{j\,k}$ allow to include information on the direction and intensity of
the temporal trends. For instance, we could center the trend $\beta_{j\,3}$
around $m_{\beta_{3}}=2$, resulting in an expected increasing trend for the
age at adult mortality of $2$ years between consecutive time points. In
Section 3.4, we compare such informative prior specification with a non-
informative elicitation, confirming the importance of the proposed approach to
induce shrinkage across countries and obtain coherent predictions.

There is a large literature on simulation methods for non-Gaussian state-space
models; see, for example, Durbin and Koopman,  (2012, Chapter 13.3), Prado and
West,  (2010, Chapter 6) or Chopin and Papaspiliopoulos,  (2020) for recent
developments. A major complexity in our methods derive from the computation of
the age-at-death distribution, which links ${\boldsymbol{\vartheta}}_{jt}$ to
$\mathbf{p}_{jt}$ via Equation 1. Therefore, we conduct posterior inference
relying on the r package nimble (de Valpine et al., , 2020), that allows to
specify the main structure of dysm, compile it in c++, and implement an
adaptive Metropolis-within-Gibbs mcmc algorithm to simulate from the posterior
distribution of the model’s parameters. Specifically, the parameters of
Equation 7 are updated using Gibbs steps, while the remaining are updated with
multivariate Metropolis steps, due to the lack of conjugacy induced by
Equation 1. Pseudo-code illustrating the main steps to conduct posterior
inference is reported in the Supplementary materials. Forecasts for future
time points $t>T$ can be simulated recursively from the posterior predictive
distribution, continuing the conditional representation of Equation 6 for each
draw of the mcmc output.

In order to improve the computational performance, we included further
modifications to efficiently compute the cdf of the mixture density function
in Expression (3) using c++. This operation is non trivial since it requires
to compute the cdf of the Skew-Normal component, which relies on the Owen,
(1956) T-function; see also Azzalini and Capitanio,  (2013, Appendix B) for
the details and Patefield and Tandy,  (2000) for practical considerations.
Overall, posterior computations in a setting with $p=12$ countries and $T=20$
time points requires roughly $2$ minutes per $1000$ iterations on an $8$-core
intel I$7$-$7700$HQ, $2.8$ ghz processor running Linux, which is satisfactory
for our purposes.

## 3 European countries mortality

### 3.1 Data description

The focus of this Section is to evaluate the ability of dysm to characterize
the age-at-death distributions of different European countries. In Section
3.2, we are interested in comparing its performance with some state-of-the-art
methods for mortality forecast, focusing on rolling windows scenarios. In
Section 3.3 we use the proposed method to estimate mortality and to forecasts
for future years, illustrating in details how dysm can be used to draw
thoughtful inferential conclusions by making inference on different
quantities. Analysis will focus on the time range $1960-2016$, and on joint
modeling the following $p=12$ countries: Austria (aut), Switzerland (che),
England & Wales (gbr), Germany (deu), Denmark (dnk), Finland (fin), France
(fra), Italy (ita), Netherlands (nld), Norway (nor), Spain (esp) and Sweden
(swe), separately for the male and female populations. These countries were
selected according to the European countries analyzed in other approaches,
such as Li and Lee,  (2005). Mortality data are retrieved from the Human
Mortality Database,  (2020), and the year- and country-specific number of
deaths are obtained with the utilities of the r package demograpy (Hyndman et
al., , 2019).

### 3.2 Performance evaluation via rolling window

In order to quantitatively evaluate the performance of dysm in estimating and
forecasting mortality, we rely on rolling window scenarios. Specifically,
analysis focus on sub-periods of $20$ years of data to fit the models, and
subsequent periods of $10$ years of data to provide forecasts. Performance is
evaluated both in terms of in-sample accuracy—comparing the model fit with the
$20$ observed years—and in terms of forecasts ability—comparing predictions
with the out-of-sample block for the subsequent $10$ years. Multiple scenarios
are obtained rolling the window forward of steps of $1$ years, thereby
dividing the range $1960-2016$ in $30$ consecutive scenarios.

We compare dysm with different popular approaches for mortality forecasting:
the Lee-Carter model (Lee and Carter, , 1992), the coherent Li-Lee model (Li
and Lee, , 2005), the functional mortality model of Hyndman–Ullah (Hyndman and
Ullah, , 2007), the compositional Oeppen model (Oeppen, , 2008), and the mem5
model (Pascariu et al., , 2019). All methods are conveniently available in a
common interface within the r package MortalityForecast (Pascariu, , 2020).
Since these methods focus on modeling single population data, we perform the
analysis on each country separately. In order to compare the approaches on the
same basis, results are evaluated in terms of the quality of the fitted and
predicted age-at-death distributions and death rates. Under dysm, such
estimates are obtained post-processing the mcmc output, computing the
posterior distribution of Equation 1 via Monte Carlo integration. For the
remaining approaches, we follow standard practices and transform predictions
into death rates and age-at-death distributions (e.g., Pascariu, , 2020).

Posterior inference proceeds separately for the female and male populations,
relying on informative prior distributions. We elicited the initial conditions
$\tilde{\vartheta}_{j0\,k}$ for the hierarchical dynamic component considering
the expected level of mortality in the temporal window under investigation,
centering old-age mortality around $70$ years ($m_{5}=70$), adult mortality
around $50$ years $(m_{3}=50)$ and setting the remaining values of $m_{k}$ to
$0$. We also let ${s_{k}}=10$ for all $k$ to induce sufficient variability in
such initial prior guesses and cover the temporal window under investigation.
Lastly, we specified a shared standard Gaussian for the drift parameters
$\beta_{j\,k}$ (setting $m_{\beta_{k}}=0$ and $s_{\beta_{k}}=1$) and a weakly-
informative Inverse-Gamma on the variances of the dynamic innovations
$\eta_{j\,k}^{2}$, setting $\mbox{a}_{k}=0.01,\mbox{b}_{k}=0.01$ for all $k$.
In each scenario, posterior inference for dysm relies on $100000$ iterations
collected after a burn-in of $5000$, thinning one iteration every $5$.
Convergence was assessed in terms of mixing, autocorrelation of the chains and
Geweke diagnostics, which resulted satisfactory in all the windows considered,
with an Effective Sample Size larger than $18000$ out of $20000$ stored
iterations for all the parameters.

Performance is evaluated in terms of Mean-Squared-Error (mse) and Mean-
Absolute-Error (mae) between the fitted and predicted age-at-death
distributions or death rates and the ground truth. Since the elements of the
age-at-death distribution are generally very small numbers, we report results
in relative terms, dividing the performance of each approach by the
performance obtained by dysm in each scenario and country. Recalling that both
mse and mae are measure of discrepancy, results larger than $1$ indicate that
predictions under dysm achieve an error which is lower than the competitor,
and therefore we can conclude that dysm performs better; the opposite holds
for values lower than $1$.

| female | male
---|---|---
Age-at-death | mae | mse | mae | mse
dysm | 1.00 | 1.00 | 1.00 | 1.00
Hyndman–Ullah |  1.41 $[1.17,1.64]$ |  2.73 $[1.91,3.80]$ |  1.92 $[1.39,2.40]$ |  5.37 $[2.76,7.79]$
Lee-Carter |  1.40 $[1.17,1.62]$ |  2.74 $[1.90,3.79]$ |  1.93 $[1.39,2.36]$ |  5.34 $[2.77,7.77]$
Li-Lee |  1.23 $[0.94,1.45]$ |  2.19 $[1.22,2.97]$ |  1.75 $[1.30,2.28]$ |  4.71 $[2.53,6.99]$
mem5 |  1.22 $[1.07,1.31]$ |  2.01 $[1.55,2.55]$ |  1.86 $[1.56,2.19]$ |  4.20 $[3.01,6.61]$
Oeppen |  1.13 $[0.92,1.36]$ |  1.88 $[1.31,2.66]$ |  1.61 $[1.24,2.14]$ |  3.79 $[2.24,6.37]$

Death rates |  |  |  |
---|---|---|---|---
dysm | 1.00 | 1.00 | 1.00 | 1.00
Hyndman–Ullah |  0.86 $[0.76,0.95]$ |  1.02 $[0.86,1.11]$ |  0.84 $[0.79,0.92]$ |  1.02 $[0.85,1.19]$
Lee-Carter |  0.91 $[0.79,1.06]$ |  1.08 $[0.88,1.86]$ |  0.87 $[0.80,1.00]$ |  1.07 $[0.87,1.58]$
Li-Lee |  1.02 $[0.71,1.26]$ |  1.81 $[0.71,5.42]$ |  0.99 $[0.77,1.17]$ |  1.90 $[0.95,3.89]$
mem5 |  1.06 $[0.97,1.14]$ |  0.99 $[0.90,1.09]$ |  1.01 $[0.96,1.05]$ |  1.02 $[0.95,1.08]$
Oeppen |  1.27 $[1.03,4.71]$ |  3.87 $[1.10,6.51]$ |  1.09 $[0.89,2.98]$ |  1.32 $[0.81,7.01]$

Table 1: In-sample relative performance in estimating the age-at-death
distribution (top) and death rates (bottom), comparing predictions of dysm
against the competitors. Median across rolling windows and countries. First
and third quartiles are reported in squared brackets.

| females | males
---|---|---
Age-at-death | mae | mse | mae | mse
dysm | 1.00 | 1.00 | 1.00 | 1.00
Hyndman–Ullah |  1.38 $[1.14,1.69]$ |  1.99 $[1.39,2.88]$ |  1.91 $[1.12,2.48]$ |  3.64 $[1.42,6.36]$
Lee-Carter |  1.30 $[1.03,1.62]$ |  1.75 $[1.12,2.65]$ |  1.73 $[1.06,2.37]$ |  3.23 $[1.25,5.54]$
Li-Lee |  1.28 $[0.91,1.56]$ |  1.66 $[0.93,2.54]$ |  1.73 $[1.06,2.36]$ |  3.23 $[1.25,5.53]$
mem5 |  2.44 $[2.02,3.03]$ |  5.67 $[3.98,8.88]$ |  2.84 $[1.70,3.88]$ |  7.92 $[3.23,9.96]$
Oeppen |  1.29 $[0.94,1.82]$ |  2.30 $[1.29,9.34]$ |  1.66 $[1.06,2.36]$ |  3.86 $[1.45,8.11]$

Death rates |  |  |  |
---|---|---|---|---
dysm | 1.00 | 1.00 | 1.00 | 1.00
Hyndman–Ullah |  1.02 $[0.88,1.10]$ |  1.35 $[1.18,1.70]$ |  0.96 $[0.89,1.05]$ |  1.24 $[1.05,1.47]$
Lee-Carter |  1.12 $[0.95,2.66]$ |  1.73 $[1.22,8.96]$ |  1.06 $[0.94,2.08]$ |  1.51 $[1.14,9.24]$
Li-Lee |  1.12 $[0.95,2.66]$ |  1.73 $[1.22,7.95]$ |  1.06 $[0.94,2.08]$ |  1.51 $[1.14,9.24]$
mem5 |  1.06 $[0.89,1.19]$ |  1.03 $[0.76,1.48]$ |  1.32 $[1.03,1.58]$ |  1.05 $[0.83,1.41]$
Oeppen |  1.67 $[1.24,4.07]$ |  3.99 $[1.53,9.44]$ |  1.29 $[1.10,1.69]$ |  1.97 $[1.36,6.35]$

Table 2: Out-of-sample relative performance in estimating the age-at-death
distribution (top) and death rates (bottom), comparing predictions of dysm
against the competitors. Median across rolling windows and countries. First
and third quartiles are reported in squared brackets.

Results referring to in-sample evaluation are reported in Table 1, computing
the median, first and third quartiles of the relative performances across
countries and rolling windows scenarios. Current results suggest that dysm
outperforms the competitors both for the male and female populations, with a
performance that is better than the other approaches in all settings. These
results confirm the ability of dysm to characterize with flexibility age-at-
death distributions across multiple countries and time periods, providing
estimates that are consistent with the observed data and more accurate than
state-of-the art methods. When dysm is used for fitting death rates, a
relatively worse performance is observed, particularly in terms of mae. This
result can be explained considering that old-ages $(80-100)$ are assigned more
weight in the computation of death rates, compared to other quantities that
characterize mortality (e.g., Keyfitz, , 2005). Therefore, dysm performs
particularly well in characterizing premature and adult mortality, and this
leads to large improvements for modeling the age-at-death distribution, in
particular for male populations where the share of adult mortality is
relatively large.

Out-of sample performances are instead reported in Table 2. Results provide
further evidence in favor of the ability of dysm, which obtains out-of-sample
forecasts more accurate than the competitors. This important results suggest
that, overall, dysm is the most accurate methods for modeling and forecasting
the age-at-death distribution in the $12$ European countries under
investigation. Performances stratified by countries are reported in the
Supplementary materials.

### 3.3 Results for European mortality

In Section 3.2, we have demonstrated that dysm outperforms the competitors in
terms of the quality of the fitted and predicted age-at-death distributions.
In this Section we argue that the proposed method has concrete practical
benefits in terms of mortality modeling also from an interpretative point of
view. Indeed, posterior inference on the parameters characterizing dysm and
their functionals provides insights on the structure of mortality and its
future evolution, highlighting aspects that have larger impact in such changes
and explaining differences across countries.

Figure 5: Posterior medians (solid lines) and $90\%$ credible intervals
(dotted lines) for the parameters characterizing dysm and functionals of
interest. Black vertical dotted lines indicate the last observed period
(2017).

To illustrate this, we focus on the same European countries as in the previous
section, modeling the more recent period $1980-2017$ and providing forecasts
until $2037$. Posterior simulation for this application relies on the same
settings as in Section 3.2. Results in terms of mixing, autocorrelations,
Geweke diagnostic and effective sample size were satisfactory for all the
parameters considered, with an effective sample size larger than $19000$ out
of $20000$ iterations for all the parameter. We conduct inference on the set
of parameters
$(\pi_{0\,jt},\pi_{1\,jt},\pi_{2\,jt},\mu_{jt},\sigma_{jt},\alpha_{jt},\tilde{\mu}^{\mbox{
\tiny SN }}_{jt},\tilde{\sigma}^{\mbox{ \tiny SN }}_{jt})$ to further
facilitate interpretation. Results are reported in Figure 5, illustrating the
posterior medians and $90\%$ credible intervals for the estimated parameters
for male and female populations separately.

Current empirical findings are consistent with known aspects of the temporal
evolution of mortality of the past decades. For example, the impact of infant
mortality—measured via the parameters $\pi_{0j\,t}$—is notably decreased in
all countries in the past years, dropping from values around $0.01$ to
$0.005$. Although both male and female populations share such decreasing
trends, we also observe that infant mortality is larger for males than females
(e.g., Drevenstedt et al., , 2008). Another clear aspect involves the
evolution of the mean age of old-age mortality, captured via the Skew-Normal
mean $\tilde{\mu}^{\mbox{ \tiny SN }}_{jt}$, which has constantly increased in
all countries. The overall trend is similar between males and females, with
woman reporting older ages on average. Similarly, the decreasing trend of
$\tilde{\sigma}^{\mbox{ \tiny SN }}_{jt}$ is consistent with the empirical
evidence on compression of old-age mortality.

Under dysm, estimates and credible intervals for any complex functional of the
model’s parameters are easily obtained. More importantly, prediction intervals
are determined taking into account the uncertainty of the entire inferential
process, which is not the case for other coherent models, where all the
uncertainty related to choice of reference is generally not accounted for.
Some interesting functionals include death rates and life expectancies, which
can be obtained from the age-at-death distribution via deterministic
transformations; dysm allows to conduct coherent inference also on this
quantities, rigorously characterizing uncertainty post-processing the mcmc
output. As an illustrative example, Figure 6 and Figure 7 depict the life-
expectancy and the death-probability, respectively, computed at different ages
and considering male and female populations in Italy and Sweden. Results are
consistent with the conclusions outlines above, confirming regular patters of
mortality such as non-overlapping trends for males and females and an overall
improvement at all ages.

Figure 6: Life expectancy at ages $0$ (at birth), $20,40,60,70$ and $80$.
Italy and Sweden, female and male population. Estimates are obtaining post-
processing the mcmc for the age-at-death distribution. Solid and dashed lines
denote posterior medians and $90\%$ credible intervals, respectively. Figure
7: Logarithm of death probability at ages $50,60,70,80$ and $90$. Italy and
Sweden, female and male population. Estimates are obtaining post-processing
the mcmc for the age-at-death distribution. Solid and dashed lines denote
posterior medians and $90\%$ credible intervals, respectively.

### 3.4 Alternative parametric specifications

In this section we evaluate the robustness of the parametric assumptions
underlying dysm, focusing on the application outlined in Section 3.3. We
investigate different specifications for the mixture components characterizing
dysm and alternative prior specifications. Specifically, we modify dysm
focusing on

  * •

modeling old-age mortality with a rescaled Beta distribution, supported in the
interval $[75-110]$, instead of a Skew-Normal distribution;

  * •

modeling infant mortality with an Half-Normal distribution, as in Zanotto et
al.,  (2020), instead of a single Dirac mass at age $0$;

  * •

removing the Gaussian component accounting for adult mortality, setting
$\pi_{1\,jt}=0$;

  * •

characterizing the joint evolution of ${\boldsymbol{\tilde{\vartheta}}}_{jt}$
with a multivariate random walk, instead of independent innovations;

  * •

characterizing the evolution of ${\boldsymbol{\tilde{\vartheta}}}_{jt}$ via
random-walks with Student-$t$ innovations with $5$ degrees of freedom, instead
of Gaussian distributions;

  * •

replacing informative priors distributions in Equation 7 with flat uniform
priors.

All modifications are implemented and estimated with nimble, relying on the
same settings as in Section 3.3. Model evaluation is performed via marginal
likelihood – estimated via the harmonic mean estimator (e.g., Raftery et al.,
, 2007) – and comparing the mse and mae of predictions obtained via posterior
medians. The empirical findings are reported in Table 3, and suggest that the
current implementation for dysm is the best performing model in terms of
marginal likelihood and other measures. For example, modeling old-age
mortality with a scaled Beta or characterizing infant mortality with an Half-
normal distribution do not improve the fit of the current specification for
dysm; similar results are obtained with more complex specifications for the
dynamic prior distribution or non-informative elicitation, confirming the
importance of the proposed specifications. Interestingly, current results
indicate that the specification without adult mortality performs better for
female than for male mortality, confirming that the inclusion of a separate
component devoted to adult mortality is fundamental for modeling male
mortality, but also provides concrete benefits for female mortality.

| males | females
---|---|---
| $-\log p(\mathbf{D})$ | mae | mse | $-\log p(\mathbf{D})$ | mae | mise
dysm | ${55.65}$ | 0.10 | 10.19 | ${46.2}$ | ${0.08}$ | ${10.21}$
Beta (old-age mortality) | $453.55$ | 0.13 | 33.83 | $307.53$ | 0.92 | 98.72
Half-normal (infant mortality) | $176.03$ | 0.15 | 9.79 | $164.07$ | 0.95 | 98.72
Removing adult mortality | $264.20$ | 0.88 | 30.79 | $88.23$ | 0.74 | 83.74
Multivariate random-walk | $301.52$ | 8.78 | 116.76 | $205.52$ | 4.75 | 135.77
Student-$t$ innovations | $254.07$ | 0.90 | 10.26 | $178.23$ | 0.91 | 96.24
dysm, uniform prior | ${359.70}$ | 2.35 | 175.14 | ${466.52}$ | 2.01 | 121.61
Table 3: Evaluation of model fit: $-\log p(\mathbf{D})$ refer to the logarithm
of the marginal likelihood, changed of sign; mse and mae refer to Mean Squared
Error and Mean Absloute Error of the posterior medians, respectively. Smaller
values indicate better fit.

## 4 Conclusion

In this article we have introduced a dynamic model based on a Skewed
distribution function (dysm) for modeling and forecasting the age-at-death
distribution of multiple countries across time. The proposed method
automatically regularize diverging scenarios, and allows to obtain predictions
which are coherent across sub-populations. In addition, this strategy leads to
a borrowing of information which improved the quality of the fit and of the
forecasts, compared with state-of-the-art methods for mortality forecasting.
Another important advantage of the proposed method is that it provides
coherent forecasts without the need of choosing a reference population, which
is instead naturally derived from data. The choice of reference population, in
fact, is still an open issue of coherent forecast models and has received
notable attention from scholars in recent years (Booth, , 2020; Kjærgaard et
al., , 2016).

We have focused directly on the age-distribution for the number of deaths,
considering the total number of deaths as a given quantity and relying on a
Multinomial likelihood. This approach allows to focus directly on the
distributions of deaths across different ages, accounting for the
heterogeneity in the overall mortality trends. For future developments, it
might be useful to include a further dynamic level in the hierarchy and model
the process of death counts $n_{j\,t}$. Such an aim can be achieved, for
example, using an hierarchical Poisson or Negative-Binomial dynamic model.

## Acknowledgments

This work was supported by the grant miur–prin 2017 project 20177BR-JXS. The
authors would like to thank David Dunson, Sally Paganin and Tommaso Rigon for
their comments and the fruitful discussions.

## References

  * Azzalini and Capitanio, (2013) Azzalini, A. and Capitanio, A. (2013).  The skew-normal and related families, volume 3.  Cambridge University Press.
  * Basellini and Camarda, (2019) Basellini, U. and Camarda, C. G. (2019).  Modelling and forecasting adult age-at-death distributions.  Population Studies, 73(1):119–138.
  * Basellini and Camarda, (2020) Basellini, U. and Camarda, C. G. (2020).  A three-component approach to model and forecast age-at-death distributions.  In Mazzuco, S. and Keilman, N., editors, Developments in Demographic Forecasting, pages 105–129. Springer International Publishing, Cham.
  * Bergeron-Boucher et al., (2017) Bergeron-Boucher, M.-P., Canudas-Romo, V., Oeppen, J., and Vaupel, J. W. (2017).  Coherent forecasts of mortality with compositional data analysis.  Demographic Research, 37:527–566.
  * Booth, (2020) Booth, H. (2020).  Coherent mortality forecasting with standards: Low mortality serves as a guide.  In Mazzuco, S. and Keilman, N., editors, Developments in Demographic Forecasting, pages 153–178. Springer International Publishing, Cham.
  * Canudas-Romo, (2010) Canudas-Romo, V. (2010).  Three measures of longevity: Time trends and record values.  Demography, 47(2):299–312.
  * Carriere, (1992) Carriere, J. F. (1992).  Parametric models for life tables.  Transactions of the Society of Actuaries, 44:77–99.
  * Cheung et al., (2005) Cheung, S. L. K., Robine, J.-M., Tu, E. J.-C., and Caselli, G. (2005).  Three dimensions of the survival curve: Horizontalization, verticalization, and longevity extension.  Demography, 42(2):243–258.
  * Chopin and Papaspiliopoulos, (2020) Chopin, N. and Papaspiliopoulos, O. (2020).  An introduction to sequential Monte Carlo.  Springer International Publishing.
  * de Valpine et al., (2020) de Valpine, P., Paciorek, C., Turek, D., Michaud, N., Anderson-Bergman, C., Obermeyer, F., Wehrhahn Cortes, C., Rodrìguez, A., Temple Lang, D., and Paganin, S. (2020).  nimble: mcmc, particle filtering, and programmable hierarchical modeling.  R package version 0.9.1.
  * Dellaportas et al., (2001) Dellaportas, P., Smith, A. F., and Stavropoulos, P. (2001).  Bayesian analysis of mortality data.  Journal of the Royal Statistical Society: Series A (Statistics in Society), 164(2):275–291.
  * Drevenstedt et al., (2008) Drevenstedt, G. L., Crimmins, E. M., Vasunilashorn, S., and Finch, C. E. (2008).  The rise and fall of excess male infant mortality.  Proceedings of the National Academy of Sciences, 105(13):5016–5021.
  * Durbin and Koopman, (2012) Durbin, J. and Koopman, S. J. (2012).  Time series analysis by state space methods.  Oxford university press.
  * Gompertz, (1825) Gompertz, B. (1825).  On the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of life contingencies.  Philosophical Transactions of the Royal Society of London, (115):513–583.
  * Heligman and Pollard, (1980) Heligman, L. and Pollard, J. H. (1980).  The age pattern of mortality.  Journal of the Institute of Actuaries, 107(1):49–80.
  * Human Mortality Database, (2020) Human Mortality Database (2020).  University of California Berkeley, USA and Max Planck Institute for Demographic Research, Germany.  data downloaded on October 2020.
  * Hyndman et al., (2013) Hyndman, R. J., Booth, H., and Yasmeen, F. (2013).  Coherent mortality forecasting: the product-ratio method with functional time series models.  Demography, 50(1):261–283.
  * Hyndman and Ullah, (2007) Hyndman, R. J. and Ullah, M. S. (2007).  Robust forecasting of mortality and fertility rates: a functional data approach.  Computational Statistics & Data Analysis, 51(10):4942–4956.
  * Hyndman et al., (2019) Hyndman, R. J. w. c. f. H. B., Tickle, L., and Maindonald., J. (2019).  demography: Forecasting Mortality, Fertility, Migration and Population Data.  R package version 1.22.
  * Keyfitz, (2005) Keyfitz, N. (2005).  Applied mathematical demography.  Springer.
  * Kjærgaard et al., (2016) Kjærgaard, S., Canudas-Romo, V., and Vaupel, J. W. (2016).  The importance of the reference populations for coherent mortality forecasting models.  In European Population Conference.
  * Klein and Moeschberger, (2006) Klein, J. P. and Moeschberger, M. L. (2006).  Survival analysis: techniques for censored and truncated data.  Springer Science & Business Media.
  * Lee and Carter, (1992) Lee, R. D. and Carter, L. R. (1992).  Modeling and forecasting us mortality.  Journal of the American Statistical Association, 87(419):659–671.
  * Lewer et al., (2020) Lewer, D., Jayatunga, W., Aldridge, R. W., Edge, C., Marmot, M., Story, A., and Hayward, A. (2020).  Premature mortality attributable to socioeconomic inequality in england between 2003 and 2018: an observational study.  The Lancet Public Health, 5(1):e33–e41.
  * Lexis, (1879) Lexis, W. H. R. A. (1879).  Sur la durée normale de la vie humaine et sur la théorie de la stabilité des rapports statistiques.  Vve. F. Henry.
  * Li, (2013) Li, J. (2013).  A Poisson common factor model for projecting mortality and life expectancy jointly for females and males.  Population Studies, 67(1):111–126.
  * Li and Lee, (2005) Li, N. and Lee, R. (2005).  Coherent mortality forecasts for a group of populations: An extension of the Lee-Carter method.  Demography, 42(3):575–594.
  * Lutz and Samir, (2010) Lutz, W. and Samir, K. (2010).  Dimensions of global population projections: what do we know about future population trends and structures?  Philosophical Transactions of the Royal Society B: Biological Sciences, 365(1554):2779–2791.
  * Makeham, (1860) Makeham, W. M. (1860).  On the law of mortality and the construction of annuity tables.  Journal of the Institute of Actuaries, 8(6):301–310.
  * Mazzuco et al., (2018) Mazzuco, S., Scarpa, B., and Zanotto, L. (2018).  A mortality model based on a mixture distribution function.  Population Studies, 72(2):191–200.
  * Oeppen, (2008) Oeppen, J. (2008).  Coherent forecasting of multiple-decrement life tables: a test using Japanese cause of death data.  Compositional data analysis conference.
  * Owen, (1956) Owen, D. B. (1956).  Tables for computing bivariate normal probabilities.  The Annals of Mathematical Statistics, 27(4):1075–1090.
  * Pascariu, (2020) Pascariu, M. D. (2020).  MortalityForecast: Standard Tools to Compare and Evaluate Various Mortality Forecasting Methods.  R package version 0.8.0.
  * Pascariu et al., (2019) Pascariu, M. D., Lenart, A., and Canudas-Romo, V. (2019).  The maximum entropy mortality model: forecasting mortality using statistical moments.  Scandinavian Actuarial Journal, 2019(8):661–685.
  * Patefield and Tandy, (2000) Patefield, M. and Tandy, D. (2000).  Fast and accurate calculation of Owen T-function.  Journal of Statistical Software, 5(5):1–25.
  * Pearson, (1897) Pearson, K. (1897).  The Chances of Death, and Other Studies in Evolution, volume 2.  E. Arnold.
  * Prado and West, (2010) Prado, R. and West, M. (2010).  Time series: modeling, computation, and inference.  CRC Press.
  * Raftery et al., (2007) Raftery, A. E., Newton, M. A., Satagopan, J. M., and Krivitsky, P. N. (2007).  Estimating the integrated likelihood via posterior simulation using the harmonic mean identity.  Bayesian Statistics, 8:1–45.
  * Siler, (1983) Siler, W. (1983).  Parameters of mortality in human populations with widely varying life spans.  Statistics in Medicine, 2(3):373–380.
  * West and Harrison, (1997) West, M. and Harrison, J. (1997).  Bayesian forecasting and dynamic models.  Springer Science & Business Media.
  * Wilson, (2001) Wilson, C. (2001).  On the scale of global demographic convergence 1950–2000.  Population and Development Review, 27(1):155–171.
  * Zanotto et al., (2020) Zanotto, L., Canudas-Romo, V., and Mazzuco, S. (2020).  A mixture-function mortality model: illustration of the evolution of premature mortality.  European Journal of Population, pages 1–27.