<figure><img src="image_1.png" /></figure>

HAL
open science

# Impacts of Climate Change on Mortality: An extrapolation of temperature effects based on time series data in France

Quentin Guibert, Gaëlle Pincemin, Frédéric Planchet

► To cite this version:

Quentin Guibert, Gaëlle Pincemin, Frédéric Planchet. Impacts of Climate Change on Mortality: An extrapolation of temperature effects based on time series data in France. 2024. hal-04596120

HAL Id: hal-04596120

https://hal.science/hal-04596120

Preprint submitted on 3 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Impacts of Climate Change on Mortality: An extrapolation of temperature effects based on time series data in France

Quentin Guibert, Gaëlle Pincemin, Frédéric Planchet

► To cite this version:

Quentin Guibert, Gaëlle Pincemin, Frédéric Planchet. Impacts of Climate Change on Mortality: An extrapolation of temperature effects based on time series data in France. International Journal of Forecasting, In press, 10.1016/j.ijforecast.2025.07.004. hal-04596120v4

HAL Id: hal-04596120

https://hal.science/hal-04596120v4

Submitted on 28 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Impacts of Climate Change on Mortality: An extrapolation of temperature effects based on time series data in France

Quentin Guibert*¹, Gaëlle Pincemin ², and Frédéric Planchet ³,⁴

¹CEREMADE, Université Paris-Dauphine, PSL University, CNRS, 75016 Paris, France

²Gallagher Re, 251 Boulevard Pereire, 75017 Paris, France

³Univ Lyon, Université Claude Bernard Lyon 1, Institut de Science Financière et d'Assurances (ISFA), Laboratoire SAF EA2429, F-69366, Lyon, France

⁴Prim'Act, 42 avenue de la Grande Armée, 75017 Paris, France

June 3, 2024

## Abstract

Most contemporary mortality models rely on extrapolating trends or past events. However, population dynamics will be significantly impacted by climate change, notably the influence of temperatures on mortality. In this paper, we introduce a novel approach to incorporate temperature effects on projected mortality using a multi-population mortality model. This method combines a stochastic mortality model with a climate epidemiology model, predicting mortality variations due to daily temperature fluctuations, be it excesses or insufficiencies. The significance of this approach lies in its ability to disrupt mortality projections by utilizing temperature forecasts from climate models and to assess the impact of this unaccounted risk factor in conventional mortality models. We illustrate this proposed mortality model using French data stratified by gender, focusing on past temperatures and mortality. Utilizing climate model predictions across various IPCC scenarios, we investigate gains and loss in life expectancy linked to temperatures and the additional mortality induced by extreme heatwaves, and quantify them by assessing this new risk factor in prediction intervals. Furthermore, we analyze the geographical differences across the Metropolitan France.

**Keywords:** Climate change, Temperature, Heatwave, Mortality forecasting, Multi-population mortality model

## 1 Introduction

The study of past trends since the 19th century has demonstrated unprecedented improvements in human mortality worldwide. In the fields of actuarial science and demography, the analysis of these trends has led to a multitude of models and approaches aimed at studying the dynamics of mortality, see Lee and Carter (1992), Brouhns et al. (2002), Barrieu et al. (2012), Plat (2009), and Dowd et al. (2020) among others. These dynamics are complex and arise from numerous factors that are more or less shared among different population groups, such as health factors, lifestyles, and the performance of health systems. These elements typically emerge through socioeconomic analyses that reveal mortality disparities within a region, see for example Bonnet et al. (2023) or Cairns et al. (2024) for recent results.

Climate change and the disruptions caused by humanity to the environment will have consequences on human health and mortality (Pörtner et al., 2022), especially in Europe (Weilnhammer et al., 2021). The WHO (2014) identify five main causes contributing to increase mortality in response to climate change: undernutrition, malaria, diarrhoeal disease, dengue, and heat. These causes will affect

<footnote>*Email: guibert@ceremade.dauphine.fr.</footnote>
# Impacts of Climate Change on Mortality: An extrapolation of temperature effects based on time series data in France

Quentin Guibert*¹, Gaëlle Pincemin † ², and Frédéric Planchet ‡ ³,⁴

¹CEREMADE, Université Paris-Dauphine, PSL University, CNRS, 75016 Paris, France

²Gallagher Re, 251 Boulevard Pereire, 75017 Paris, France

³Univ Lyon, Université Claude Bernard Lyon 1, Institut de Science Financière et d'Assurances (ISFA), Laboratoire SAF EA2429, F-69366, Lyon, France

⁴Prim'Act, 42 avenue de la Grande Armée, 75017 Paris, France

## Abstract

Most contemporary mortality models rely on extrapolating trends or past events. However, population dynamics will be significantly impacted by climate change, notably the influence of temperatures on mortality. In this paper, we introduce a novel approach to incorporate temperature effects on projected mortality using a multi-population mortality model. This method combines a stochastic mortality model with a climate epidemiology model, predicting mortality variations due to daily temperature fluctuations, be it excesses or insufficiencies. The significance of this approach lies in its ability to disrupt mortality projections by utilizing temperature forecasts from climate models and to assess the impact of this unaccounted risk factor in conventional mortality models. We illustrate this proposed mortality model using French data stratified by sex and age, focusing on past temperatures and mortality. Utilizing climate model predictions across various IPCC scenarios, we investigate gains and loss in life expectancy linked to temperature and the additional mortality induced by extreme heatwaves, and quantify them by assessing this new risk factor in prediction intervals. Furthermore, we analyze the geographical differences across the Metropolitan France.

**Keywords:** Climate change, Temperature, Heatwave, Mortality Forecasting, Multi-population mortality model

## 1 Introduction

The study of past trends since the 19th century has demonstrated unprecedented improvements in human mortality worldwide. In the fields of actuarial science and demography, the analysis of these trends has led to a multitude of models and approaches aimed at studying the dynamics of mortality, see Lee and Carter (1992), Brouhns et al. (2002), Barrieu et al. (2012), Plat (2009), and Dowd et al. (2020) among others. These dynamics are complex and arise from numerous factors that are more or less shared among different population groups, such as health factors, lifestyles, and the performance of health systems. These elements typically emerge through socioeconomic analyses that reveal mortality disparities within a region, see for example Bonnet et al. (2023) or Cairns et al. (2024) for recent results.

Climate change and the disruptions caused by humanity to the environment will have consequences on human health and mortality (Pörtner et al., 2022), including European countries (Weilnhammer et al., 2021). The WHO (2014) identify five main causes contributing to increase mortality in response

<footnote>*Email: guibert@ceremade.dauphine.fr.</footnote>

<footnote>†Email: gaelle.pincemin@gallagherre.com.</footnote>

<footnote>‡Email: frederic.planchet@univ-lyon1.fr.</footnote>
populations differently in the future depending on their level of education, access to healthcare, and age. For example, undernutrition primarily affects children under 15 years old, while heat mainly impacts mortality of the elderly. Malaria and diarrhoeal diseases are prevalent in children under 5 years, especially in South Asia and Africa. Heatwaves are expected to affect the entire world, with more than 68,000 heat-related deaths in Europe by 2030 without adaptation measures (Ballester et al., 2023). In these regions, the main mortality concerns among those listed by the WHO (2014) are actually heatwaves (Pottier et al., 2021) and vector-borne diseases (Semenza and Suk, 2018; Paz, 2020) with potential impacts on longevity and life expectancy. Recently, a significant number of studies have highlighted the impact of temperature fluctuations resulting from climate change on mortality (Gasparrini et al., 2015; Gronlund et al., 2018; Guo et al., 2018; Vicedo-Cabrera et al., 2018; Lee et al., 2020; Martínez-Solanas et al., 2021). With human activity, the rise in temperatures has accelerated since the pre-industrial era, reaching an average global increase of 1.1°C compared to the pre-industrial period (1850-1900) (Pörtner et al., 2022). These changes vary in magnitude according to latitude. For example, French population experienced exceptionally warm temperatures in 2022 (Haut Conseil pour le Climat, 2023) with consequences in terms of excess mortality (Ballester et al., 2023). According to IPCC projections, the direct impacts of temperature on human health are expected to intensify regardless of future climate trajectories. This could lead to changes in premature mortality due to temperature and adjustments in mortality trends related to the intra-annual death distribution. Under these conditions, it becomes crucial to study and incorporate the effects of climate change, particularly temperature, into the stochastic mortality models used for population projections. This is also a significant issue for insurers' and pension funds' models, as well as for supervisory authorities seeking to assess the vulnerability of financial institutions to climate change (ACPR, 2024). Therefore, the primary aim of this paper is to account for the effect of temperature variations on future projected mortality rates within the framework of multi-population stochastic mortality models.

The proposed model comprises two components, estimated separately. The first component is a conventional stochastic mortality model calibrated on historical trends. The second component is constructed using a climate epidemiology model. This last model is calibrated using daily temperature and mortality data allowing for the attribution of a fraction of deaths to the effects of temperature. By coupling these two models, our approach produces a stochastic mortality model accounting for the trend induced by temperatures. Our two-step approach provides significant flexibility to the mortality model, enabling projections that incorporate various climate scenarios concerning temperature evolution. Additionally, it allows for the measurement of gains or losses in life expectancy due to climate change. More precisely, we use daily temperature projections at different geographical levels and quantify the local and specific impact of temperature on mortality over specific periods of the year, e.g. isolating extreme temperatures. Furthermore, our approach also incorporates the uncertainty associated with various climate models.

The approach implemented in our article falls within the framework of modeling exogenous shocks to mortality rate dynamics. Starting from single-population mortality models, a relatively sparse literature has explored the introduction of stochastic shocks. Most of these models have been introduced to incorporate the effects of economic shocks (French and O'Hare, 2014) or pandemic shocks (Liu and Li, 2015; Robben et al., 2022; Goes et al., 2023) into the dynamics of mortality rates. In an insurance context, this may have substantial impacts on the amount of capital requirement, see e.g. (Planchet and Gautier de La Plaine, 2024) among others. Recently, Robben and Antonio (2024) introduce a model with mortality shocks based on a regime-switching model. Their approach considers a multi-population mortality model and mortality jumps that can occur in a regime of high volatility. So far, few stochastic mortality models have focused on integrating temperature-induced shocks. Seklecka et al. (2017) extend the framework of the Lee and Carter (1992) model by adding a component related to temperature-related mortality. This model is applied to data from the United Kingdom. Li and Tang (2022) implement an approach based on extreme value theory on US monthly data to model the dependence between extreme temperatures, both hot and cold, and excess mortality. For these different approaches, the generated mortality shocks are estimated from outliers observed generally in annual mortality rates. These outliers are generated by catastrophic events, such as pandemics, wars or other major events. In contrast, temperatures are responsible for lower-intensity mortality
to climate change: undernutrition, malaria, diarrhoeal disease, dengue, and heat. These causes will affect populations differently in the future depending on their level of education, access to healthcare, and age. For example, undernutrition primarily affects children under 15 years old, while heat mainly impacts mortality of the elderly. Malaria and diarrhoeal diseases are prevalent in children under 5 years, especially in South Asia and Africa. Heatwaves are expected to affect the entire world, with more than 68,000 heat-related deaths in Europe by 2030 without adaptation measures (Ballester et al., 2023). In these regions, the main mortality concerns among those listed by the WHO (2014) are actually heatwaves (Pottier et al., 2021) and vector-borne diseases (Semenza and Suk, 2018; Paz, 2020) with potential impacts on longevity and life expectancy. Recently, a significant number of studies have highlighted the impact of temperature fluctuations resulting from climate change on mortality (Gasparrini et al., 2015; Gronlund et al., 2018; Guo et al., 2018; Vicedo-Cabrera et al., 2018; Lee et al., 2020; Martínez-Solanas et al., 2021). With human activity, the rise in temperatures has accelerated since the pre-industrial era, reaching an average global increase of 1.1°C compared to the pre-industrial period (1850-1900) (Pörtner et al., 2022). These changes vary in magnitude according to latitude. For example, French population experienced exceptionally warm temperatures in 2022 (Haut Conseil pour le Climat, 2023) with consequences in terms of excess mortality (Ballester et al., 2023). According to IPCC projections, the direct impacts of temperature on human health are expected to intensify regardless of future climate trajectories. This could lead to changes in premature mortality due to temperature and adjustments in mortality trends related to the intra-annual death distribution. Under these conditions, it becomes crucial to study and incorporate the effects of climate change, particularly temperature, into the stochastic mortality models used for population projections. This is also a significant issue for insurers' and pension funds' models, as well as for supervisory authorities seeking to assess the vulnerability of financial institutions to climate change (ACPR, 2024). Therefore, the primary aim of this paper is to account for the effect of temperature variations on future projected mortality rates within the framework of multi-population stochastic mortality models.

The proposed model comprises two components, estimated separately. The first component is a conventional stochastic mortality model calibrated on historical trends. The second component is constructed using a climate epidemiology model. This last model is calibrated using daily temperature and mortality data allowing for the attribution of a fraction of deaths to the effects of temperature. By coupling these two models, our approach produces a stochastic mortality model accounting for the trend induced by temperatures. Our two-step approach provides significant flexibility to the mortality model, enabling projections that incorporate various climate scenarios concerning temperature evolution. Additionally, it allows for the measurement of gains or losses in life expectancy due to climate change. More precisely, we use daily temperature projections at different geographical levels and quantify the local and specific impact of temperature on mortality over specific periods of the year, e.g. isolating extreme temperatures. Furthermore, our approach also incorporates the uncertainty associated with various climate models.

The approach implemented in our article falls within the framework of modeling exogenous shocks to mortality rate dynamics. Starting from single-population mortality models such as the Lee and Carter (1992) model, see Basellini et al. (2022) for a review, a relatively sparse literature has explored the introduction of stochastic shocks. Most of these models have been introduced to incorporate the effects of economic shocks (French and O'Hare, 2014; Niu and Melenberg, 2014; Boonen and Li, 2017; Seklecka et al., 2019; Cupido et al., 2024; Ma and Boonen, 2024) or jumps in mortality rate dynamics, such as catastrophes, wars or pandemic shocks (Liu and Li, 2015; Robben et al., 2022; Goes et al., 2025) into the dynamics of mortality rates. In an insurance context, this may have substantial impacts on the amount of capital requirement, see e.g. Planchet and Gautier de La Plaine (2024) among others, or in pricing applications, see e.g. Chen and Cox (2009) and Li et al. (2023). Incorporating jumps whose effects gradually disappear over several years is part of model enhancement. Recently, Robben and Antonio (2024) develop a multi-population approach comprising a regime-switching model with a high-volatility state where Gaussian shocks on the mortality trend can occur. Their model involves maximum likelihood estimation on improvement rates in several steps and requires preprocessing of data labeled as outlines. Another recent example is provided by Goes et al. (2025) who use a Bayesian formulation and consider the impact of shocks that gradually diffuse over several years. So far, very
fluctuations that vary spatially, incorporate lagged effects, and are observed over specific periods of the year. Additionally, compensations between the effects of hot and cold temperatures can be observed over a calendar year. It is almost universally observed that the mortality response to hot and cold temperatures follows a U- or V-shaped curve depending on the population (Kunst et al., 1993; Guo et al., 2014). These effects are generally modeled using daily data through time-series regression models that capture exposure-response associations. The distributed lag non-linear generalized model (DLNM)(Gasparrini et al., 2010; Gasparrini, 2014) is a reference in this field and has been used in numerous studies, see among others Gasparrini et al. (2015), Gasparrini et al. (2017), and Vicedo-Cabrera et al. (2021). It estimates the number of deaths attributable to temperatures over specific periods, accounting for different lag effects depending on whether the temperatures are hot or cold. In a prospective approach, these temperature variations should be studied simultaneously, as the number of deaths induced by cold temperatures is expected to decrease while that induced by hot temperatures is expected to increase in Europe due to the climate change (Martínez-Solanas et al., 2021).

The numerical application developed in this paper focuses on the use of a mortality model for two populations using the augmented common factor approach introduced by Li and Lee (2005). Many multi-population mortality models have been developed in the literature in recent years. Most of these models are based on the framework proposed by Lee and Carter (1992) or the CDB models (Cairns et al., 2006), as reviewed by Villegas et al. (2017). These models typically assume coherence or semi-coherence (Li et al., 2017) among several sub-populations. In our paper, we impose a coherence assumption between women and men in France regarding mortality not affected by temperatures. Although the long-term convergence hypothesis is challenging to demonstrate for women and men, this modeling illustrates the possibility of coupling the DLNM climate epidemiology model with a standard multi-population mortality model. This issue is indeed essential given the importance of multi-population mortality models (Antonio et al., 2017) and the need to account for geographical factors when modeling temperature dynamics.

The contribution of our paper is fourfold. First, we develop a new specification of a multi-population mortality model that incorporates the effect of temperature variations on central mortality rates. These effects require a preliminary phase of estimating temperature-attributable deaths using a DLNM model. By subsequently modeling deaths unaffected by temperature through a multi-population model, our approach associates temperature-related and non-temperature-related deaths while preserving the Poisson assumption of the number of deaths. To the best of our knowledge, this is a novel approach in the literature enabling an easy coupling of most single-population or multi-population mortality models with temperature. We illustrate it using the Li and Lee (2005) model on French mortality data. Once calibrated, our modeling framework can project future mortality by incorporating various temperature-related effects. Second, the temperature effects are calibrated by gender and age groups, allowing us to integrate their specific sensitivities to temperatures into the projections. Third, our multi-population mortality model is associated to temperature projections from 12 climate models across three Representative Concentration Pathways (RCP) scenarios for the period 2020-2100. Then, we evaluate the loss or gain in life expectancy due to future temperature trends in these different climate scenarios in France, highlighting the specific impact of extreme heatwaves. Furthermore, this approach allows us to measure the impact of uncertainties caused by different climate models. Fourth, we examine the sensitivity of mortality projections to the local temperature dynamics of various French cities. This reveals divergences and a potentially new source of heterogeneity in mortality projections.

This paper is organized as follows. Section 2 presents a brief review of the literature on the relationship between mortality and temperature. Section 3 introduces the notations used and discusses the approach taken in this paper compared to existing models incorporating exogenous shocks. We then describe our multi-population mortality model excluding temperature effects. Finally, we explain how temperature-attributable deaths are modeled. Section 4 details the various estimation procedures adopted, starting with the DLNM model. We then explain how to calculate the temperature-attributable fractions annually, followed by the estimation of the mortality model. Section 5 describes
few mortality models incorporating a temperature-related effect have been developed, despite the abundance of epidemiological literature on the effects of heat, cold, or heatwaves on mortality, as discussed in Section 2. Seklecka et al. (2017) extend the framework of the Lee and Carter (1992) and Plat (2009) models by adding a component related to temperature-related mortality. Their model is applied to data from the United Kingdom. Li and Tang (2022) implement an approach based on extreme value theory on US monthly data to model the dependence between extreme temperatures, both hot and cold, and excess mortality. For these different approaches, the generated mortality shocks are estimated from outliers typically observed in annual mortality rates. These outliers are generated by catastrophic events, such as pandemics, wars or other major events. In contrast, temperatures are responsible for lower-intensity mortality fluctuations that vary spatially, incorporate lagged effects, and are observed over specific periods of the year. Additionally, compensations between the effects of hot and cold temperatures can be observed over a calendar year. For instance, heatwave-type shocks are acute and can be partly or totally compensate throughout the year due to harvest effects, i.e. the decrease in mortality observed after a peak of heat-related mortality. This effect has been observed notably in the French population (Pascal et al., 2018), but is not systematically found in all studies. For these reasons, approaches based on annual mortality data may encounter difficulties in being extended to model the effect of future temperatures on mortality. Another important aspect to consider is the increasing impact of climate change, particularly the consequences of rising temperatures and extreme heat events. This trend tends to limit the relevance of models in which shocks are estimated based on previous annual observations and suggests studying the temperature-mortality relationship.

It is almost universally observed that the mortality response to hot and cold temperatures follows a U- or V-shaped curve, depending on the population (Kunst et al., 1993; Guo et al., 2014). These effects are generally modeled using daily data through time-series regression models that capture exposure-response associations. From the aggregated temperature-mortality time series data, two traditional models are explored to capture the nonlinear correlation between temperature and mortality. Muggeo (2008) introduces the constrained segmented distributed lag (CSDL) model. This model assesses the effects of heat and cold using a fixed threshold within a Poisson regression framework. Recently, Pincemin (2021) implement it to measure the effect of temperatures on projected life expectancy of the French population. Alternatively, the distributed lag non-linear generalized model (DLNM)(Gasparrini et al., 2010; Gasparrini, 2014) is a reference in epidemiology and has been used in numerous studies, see among others Gasparrini et al. (2015), Gasparrini et al. (2017), and Vicedo-Cabrera et al. (2021). It estimates the number of deaths attributable to temperature effects over specific periods, accounting for different lag effects depending on whether the temperatures are hot or cold. In a prospective approach, these temperature variations should be studied simultaneously, as the number of deaths induced by cold temperatures is expected to decrease while that induced by hot temperatures is expected to increase in Europe due to the climate change (Martínez-Solanas et al., 2021).

The numerical application developed in this paper focuses on the use of a mortality model for two populations using the augmented common factor approach introduced by Li and Lee (2005). Many multi-population mortality models have been developed in the literature in recent years. Most of these models are based on the framework proposed by Lee and Carter (1992) or the CBD models (Cairns et al., 2006), as reviewed by Villegas et al. (2017). These models typically assume coherence or semi-coherence (Li et al., 2017) among several sub-populations. In our paper, we impose a coherence assumption between women and men in France regarding mortality not affected by temperatures. Although the long-term convergence hypothesis is challenging to demonstrate for women and men, this modeling illustrates the possibility of coupling the DLNM climate epidemiology model with a standard multi-population mortality model. This issue is indeed essential given the importance of multi-population mortality models (Antonio et al., 2017) and the need to account for geographical factors when modeling temperature dynamics.

The contribution of our paper is fourfold. First, we develop a new specification of a multi-population mortality model that incorporates the effect of temperature variations on central mortality rates. These effects require a preliminary phase of estimating temperature-attributable deaths using
the calibration of time series models used to project mortality trends. In addition, we present the projection of the attributable fractions based on climate scenarios and the life expectancy loss due to temperatures. In Section 6, we estimate and project the series of central death rates with and without the effect of temperatures. We then quantify the gains and losses in life expectancy induced by temperatures up to 2100 in Metropolitan France under different RCP scenarios and for different locations. Finally, Section 7 concludes the paper.

# 2 Background on temperature-mortality effects

This section presents a brief overview of the influence of temperatures on human mortality. These changes can have diverse consequences, directly impacting the body's function during extreme cold or heat situations, or indirectly affecting our environment-leading to issues such as disease outbreaks, food shortages, and deteriorating living conditions. Our specific focus lies in examining the effects of temperature fluctuations on mortality rates. For a comprehensive summary of other consequences of climate change on mortality, including aspects like vector-borne diseases, under-nutrition, and related factors, we refer to the sixth assessment report of the IPCC (Pörtner et al., 2022, TS.C.6).

## 2.1 The effects of heat and cold on human body

Human beings are homeotherms, equipped with various protective mechanisms to maintain stable internal temperatures. Prolonged exposure to high external temperatures can elevate internal heat, while cold temperatures can lower it (Beker et al., 2018). During periods of high heat, the body experiences several effects, including increased sweating, elevated heart rate, lower blood pressure, confusion, heat exhaustion, dehydration, heatstroke, heat rash, fainting, and heat edema, among other symptoms (WHO, 2018). To counter high temperatures, the body employs external mechanisms such as sweating and internal mechanisms that involve the production of specific proteins to safeguard internal cells (Bouchama and Knochel, 2002). However, these mechanisms have limitations and become inadequate with prolonged exposure. This can lead to more severe conditions like hyperthermia, identifiable through symptoms like rapid pulse, headaches, nausea, and loss of consciousness (Bouchama and Knochel, 2002). Severe dehydration is another consequence, carrying significant cardiovascular and pulmonary risks, potentially resulting in a rapid increase in body temperature and a heightened risk of fatality (WHO, 2004). The impact of heatstroke and dehydration during heatwaves can persist for up to three days post-exposure (Alsaqali et al., 2022). The risks escalate when high temperatures endure. Heatwaves pose greater threats when protective measures for vulnerable populations are lacking. Elderly individuals, those with existing health conditions, pregnant women, and outdoor workers face higher vulnerability. Consequently, targeted protection plans are implemented for these groups.

As with hot temperatures, cold temperatures elicit various physiological responses in individuals, see Ryti et al. (2016) for a review. Upon exposure to cold, initial symptoms include reduced blood circulation in extremities, arrhythmia, breathing challenges, muscle discomfort, and restricted mobility. Prolonged exposure to cold can lead to increasing levels of hypothermia. As the body temperature drops, the heart rate decreases, and severe bradycardia can occur when the internal temperature falls below 32°C, potentially resulting in fatality (Beker et al., 2018). Cold weather directly impacts heart and pulmonary health, notably affecting cardiovascular diseases (CVD). Particularly, mortality factors associated with cold weather include coronary heart disease, stroke and respiratory disorders. These health risks are notably prevalent in socioeconomically disadvantaged groups and regions with milder winters (Mercer, 2003).

## 2.2 Seasonal variations in mortality

The climate and seasonal variations significantly influence human mortality. Episodes of excess mortality can be observed during periods of heat, cold, or substantial temperature fluctuations, with
a DLNM model. By subsequently modeling deaths not attributable to temperatures through a multi-
population model, our approach associates temperature-related and non-temperature-related deaths
while preserving the Poisson assumption of the number of deaths. To the best of our knowledge, this
is a novel approach in the literature enabling an easy coupling of most single-population or multi-
population mortality models with temperature. We illustrate it using the Li and Lee (2005) model
on French mortality data. Once calibrated, our modeling framework can project future mortality by
incorporating various temperature-related effects. Second, the temperature effects are calibrated by
sex and age groups, allowing us to integrate their specific sensitivities to temperature into the projec-
tions. Third, our multi-population mortality model is associated to temperature projections from 12
climate models across three Representative Concentration Pathways (RCP) scenarios for the period
2020-2100. Then, we evaluate the loss or gain in life expectancy due to future temperature trends in
these different climate scenarios in France, highlighting the specific impact of extreme heatwaves. Fur-
thermore, this approach allows us to measure the impact of uncertainties caused by different climate
models. Fourth, we examine the sensitivity of mortality projections to the local temperature dynam-
ics of various French cities. This reveals divergences and a potentially new source of heterogeneity in
mortality projections.

This paper is organized as follows. Section 2 presents a brief review of the literature on the relationship between mortality and temperature. Section 3 introduces the notations and the modeling framework used. In particular, we describe our multi-population mortality model for all causes of death excluding temperature effects and explain how temperature-attributable mortality is incorporated through the DLNM model. Section 4 details the various estimation procedures adopted, starting with the DLNM model. We then explain how to calculate the total temperature-attributable fractions, followed by the estimation of the mortality model. Section 5 describes the calibration of time series models used to project mortality trends. In addition, we present the projection of the attributable fractions based on climate scenarios and the life expectancy loss due to temperature. In Section 6, we estimate and project the series of central death rates by adding the effect of future temperatures. We then quantify the gains and losses in life expectancy induced by temperatures up to 2100 in Metropolitan France under different RCP scenarios and for different locations. Finally, Section 7 concludes the paper.

2 Background on temperature-mortality effects

This section presents a brief overview of the influence of temperatures on human mortality. These
changes can have diverse consequences, directly impacting the body’s function during extreme cold
or heat situations, or indirectly affecting our environment-leading to issues such as disease outbreaks,
food shortages, and deteriorating living conditions. Our specific focus lies in examining the effects of
temperature fluctuations on mortality rates. For a comprehensive summary of other consequences of
climate change on mortality, including aspects like vector-borne diseases, under-nutrition, and related
factors, we refer to the sixth assessment report of the IPCC (Pörtner et al., 2022, TS.C.6).

## 2.1 The effects of heat and cold on human body

Human beings are homeotherms, equipped with various protective mechanisms to maintain stable internal temperatures. Prolonged exposure to high external temperatures can elevate internal heat, while cold temperatures can lower it (Beker et al., 2018). During periods of high heat, the body experiences several effects, including increased sweating, elevated heart rate, lower blood pressure, confusion, heat exhaustion, dehydration, heatstroke, heat rash, fainting, and heat edema, among other symptoms (WHO, 2018). To counter high temperatures, the body employs external mechanisms such as sweating and internal mechanisms that involve the production of specific proteins to safeguard internal cells (Bouchama and Knochel, 2002). However, these mechanisms have limitations and become inadequate with prolonged exposure. This can lead to more severe conditions like hyperthermia, identifiable through symptoms like rapid pulse, headaches, nausea, and loss of consciousness (Bouchama
characteristics varying by region. These differences in mortality patterns can be considerable across studied areas, yet the temperature-mortality relationship is generally assessed as having a U or V-shaped-relationship (Kunst et al., 1993; Guo et al., 2014). Both heat and cold increase mortality risk: cold effects can persist over several weeks, while heat impacts are immediate but shorter-lived.

Each region exhibits an optimal temperature, known as the Minimum Mortality Temperature (MMT), see Yin et al. (2019) for a current and projected overview across the planet. This temperature corresponds to the point at which the death rate is lowest, reflecting the population's acclimatization to its surroundings. The MMT varies significantly across countries, influenced by latitude, climate type, and average annual temperature, ranging from 14.2 to 31.1°C (Tobías et al., 2021). The observed MMT depends on the timeframe, serving as a crucial indicator of population adaptation. For instance, in France, the MMT is reported to have increased from 17.5°C (1968-1981) to 18.2°C (1996-2009) (Todd and Valleron, 2015). Globally, for every 1°C rise in average annual temperatures, the MMT increased by 0.8°C (Tobías et al., 2021). Of course, socioeconomic factors also influence these thresholds (Krummenauer et al., 2019).

Sensitivity to ambient temperatures varies by age, sex, health status, and social category (WHO, 2011). For example, the elderly (over 65 years) and young children (0-5 years) in tropical areas are particularly vulnerable to temperature fluctuations (Egondi et al., 2012). Among the elderly, women have a higher risk of heat-related mortality during heatwaves (Marí-Dell'Olmo et al., 2019; Yu et al., 2010). Geographic regions differ in their exposure to temperature variations and their capacity to respond to climate change, reflecting population acclimatization. Urban density within the same geographical area amplifies the risk of heat-related mortality, positively correlated with population density, air quality (e.g., fine particle concentration), or social inequality (Sera et al., 2019).

In countries with temperate and continental climates, mortality follows a distinct seasonal pattern (Madaniyazi et al., 2022). Mortality rates are higher in cold seasons compared to warm seasons, a trend observed globally but with varying seasonal amplitudes between regions. Temperature substantially explains this pattern, as stated in the literature (Pörtner et al., 2022), although the relation between mortality and temperatures can involve complex physiological aspects or indirect environmental and social aspects.

In tropical regions, besides temperature, rainfall plays also a crucial role. During monsoon periods, diseases like malaria, diarrheal diseases, or vector-borne illnesses spread more rapidly, significantly affecting women (Egondi et al., 2012). These diseases, identified by the World Health Organization (WHO) as major concerns related to climate change, contribute significantly to mortality (WHO, 2014).

Additionally, the human body's cooling mechanisms, primarily reliant on perspiration, can be compromised if ambient humidity levels are either too low or too high. From a physiological standpoint, scientific studies have highlighted the potential health hazards associated with extreme temperatures in highly humid environments (Davis et al., 2016). However, epidemiological evidence regarding the consequences of wet-hot events seems remain less conclusive (Armstrong et al., 2019), while Fang et al. (2023) indicating an increased risk of mortality in the presence of dry-hot events.

## 2.3 Temperature related deaths : past and expectation

Over the period 1998-2017, the WHO estimates the number of heatwave induced deaths at over 166,000 worldwide, including 70,000 deaths in Europe in the summer of 2003 alone (WHO, 2017). The number of people exposed to heatwaves is also expected to continue to rise. In Europe, it is 9.3 million and would reach almost 300 million in 2100 in a +3°C scenario (Naumann et al., 2020). The annual number of deaths linked to heatwaves is currently estimated at more than 61,000 in Europe during the summer in 2022 and would be more than 120,000 in 2050 (Ballester et al., 2023). However, it is important to note that its results strongly depend on the definition given for a heatwave, see Appendix A.1. In European countries and Japan, where the proportion of older people in the population is increasing, the damage from heat could be even more severe in the long term. Indeed, the average temperature
and Knochel, 2002). Severe dehydration is another consequence, carrying significant cardiovascular and pulmonary risks, potentially resulting in a rapid increase in body temperature and a heightened risk of fatality (WHO, 2004). The impact of heatstroke and dehydration during heatwaves can persist for up to three days post-exposure (Alsaiqali et al., 2022). The risks escalate when high temperatures endure. Heatwaves pose greater threats when protective measures for vulnerable populations are lacking. Elderly individuals, those with existing health conditions, pregnant women, and outdoor workers face higher vulnerability. Consequently, targeted protection plans are implemented for these groups.

As with hot temperatures, cold temperatures elicit various physiological responses in individuals, see Ryti et al. (2016) for a review. Upon exposure to cold, initial symptoms include reduced blood circulation in extremities, arrhythmia, breathing challenges, muscle discomfort, and restricted mobility. Prolonged exposure to cold can increase levels of hypothermia. As the body temperature drops, the heart rate decreases, and severe bradycardia can occur when the internal temperature falls below 32°C, potentially resulting in fatality (Beker et al., 2018). Cold weather directly impacts heart and pulmonary health, notably affecting cardiovascular diseases (CVD). Particularly, mortality factors associated with cold weather include coronary heart disease, stroke and respiratory disorders. These health risks are notably prevalent in socioeconomically disadvantaged groups, and increase when temperature falls below about 0°C. Paradoxically, mortality rates also rise sharply when ambient temperature falls in regions with milder winters (Mercer, 2003).

## 2.2 Seasonal variations in mortality

The climate and seasonal variations significantly influence human mortality. Episodes of excess mortality can be observed during periods of heat, cold, or substantial temperature fluctuations, with characteristics varying by region. These differences in mortality patterns can be considerable across studied areas, yet the temperature-mortality relationship is generally assessed as having a U or V-shaped-relationship (Kunst et al., 1993; Guo et al., 2014). Both heat and cold increase mortality risk: cold effects can persist over several weeks, while heat impacts are immediate but shorter-lived.

Each region exhibits an optimal temperature, known as the Minimum Mortality Temperature (MMT), see Yin et al. (2019) for a current and projected overview across the planet. This temperature corresponds to the point at which the death rate is the lowest, reflecting the population's acclimatization to its surroundings. The MMT varies significantly across countries, influenced by latitude, climate type, and average annual temperature, ranging from 14.2 to 31.1°C (Tobías et al., 2021). The observed MMT depends on the timeframe, serving as a crucial indicator of population adaptation. For instance, in France, the MMT is reported to have increased from 17.5°C (1968-1981) to 18.2°C (1996-2009) (Todd and Valleron, 2015). Globally, for every 1°C rise in average annual temperatures, the MMT increased by 0.8°C (Tobías et al., 2021). Of course, socioeconomic factors also influence these thresholds (Krummenauer et al., 2019).

Sensitivity to ambient temperatures varies by age, sex, health status, and social category (WHO, 2011). For example, the elderly (over 65 years) and young children (0-5 years) in tropical areas are particularly vulnerable to temperature fluctuations (Egondi et al., 2012). Among the elderly, women have a higher risk of heat-related mortality during heatwaves (Marí-Dell'Olmo et al., 2019; Yu et al., 2010). Geographic regions differ in their exposure to temperature variations and their capacity to respond to climate change, reflecting population acclimatization. Urban density within the same geographical area amplifies the risk of heat-related mortality, positively correlated with population density, air quality (e.g., fine particle concentration), or social inequality (Sera et al., 2019).

In countries with temperate and continental climates, mortality follows a distinct seasonal pattern (Madaniyazi et al., 2022). Mortality rates are higher in cold seasons compared to warm seasons, a trend observed globally but with varying seasonal amplitudes between regions. Temperature substantially explains this pattern, as stated in the literature (Pörtner et al., 2022), although the relation between mortality and temperatures can involve complex physiological aspects or indirect environmental and social aspects.
of the planet has already increased by +1.1°C (Pörtner et al., 2022, TS.C.6) since the pre-industrial
era and is expected to continue to rise under any scenario. Thus, extreme weather events such as
heatwaves are expected to become more frequent while at the same time the probability of cold waves
is expected to decrease (IPCC, 2023). The seasonality of deaths observed today in temperate climates
could strongly evolved. Climate change will reduce extreme cold spells and winter seasons period
(IPCC, 2023). This would conduct to reduce the cold exposure and cold fatalities (Madaniyazi et al.,
2021; Naumann et al., 2020) when most of the temperature related deaths was caused by cold period
(Gasparrini et al., 2015). In some parts of Europe, heat-associated health effects may even outweigh
the effects of cold (Gronlund et al., 2018; Martínez-Solanas et al., 2021).

So far, no real mitigation or change in the distribution of annual deaths has already been observed
in temperate countries (Cheng et al., 2019), which might explain a lesser focus in demographic or
actuarial literature.

3 A multi-population mortality model with temperature effects

In this section, we explain our approach to incorporating temperature change into our mortality model. After introducing the data and notations (Section 3.1), we describe existing mortality models with jumps in Section 3.2. Section 3.3 presents the proposed dynamics for central (virtual) mortality rates excluding temperature effects, followed by the Poisson formulation used to estimate the model parameters. Finally, Section 3.4 describes the DLNM model used to estimate temperature-attributable deaths.

**3.1 Notation**

In the following, some notations are introduced for modeling mortality. Let $\mu_{x,t}^{(g)}$, $E_{x,t}^{(g)}$, and $D_{x,t}^{(g)}$ represent, respectively, the force of mortality, observed exposure to risk, and observed number of deaths for gender $g \in \{\text{female, male}\} = \mathcal{G}$ within the age range $[x, x+1)$ and the calendar year $[t, t+1)$, where $x \in \mathcal{X} = \{x_{\min}, \dots, x_{\max}\}$ denotes the set of integer ages, and $t \in \mathcal{T}_y = \{y_{\min}, \dots, y_{\max}\}$ represents the set of calendar years considered as the calibrated period for mortality models. In this paper, two populations, women and men, are considered in Metropolitan France. However, our approach can be extended directly to other population groups, e.g. two or more countries. For clarity, bold notations are used for vectors and matrices.

Common stochastic mortality models such as the Lee-Carter model (Lee and Carter, 1992) rely on
the crude central death rate of mortality $\hat{m}_{x,t}^{(g)} = D_{x,t}^{(g)}/E_{x,t}^{(g)}$. It is equivalent to the estimated force of
mortality $\hat{\mu}_{x,t}^{(g)}$ assuming a piecewise constant force of mortality. Traditional mortality models may also
use the mortality rate $\hat{q}_{x,t}^{(g)} \approx 1 - \exp(-\hat{\mu}_{x,t}^{(g)})$ under the classical piecewise constant force of mortality
assumption.

This paper considers two distinct time scales. First, the mortality model component relating to
temperatures is fitted daily within the calibration period $\mathcal{T}_d = \{d_{\min}, \dots, d_{\max}\}$ across the range of
integer ages $\mathcal{X}$. Specifically, to evaluate the impact of temperature, we gather daily records of the
number of deaths over the calibration period. These specific details are sourced from the *Institut
national de la statistique et des études économiques* (INSEE, 2020), and we denote the number of
deaths in day $d$ of year $t$ for gender $g$ as $D_{x,t,d}^{(g)}$.

Second, the model forecasting annual mortality is calibrated using yearly data observed within the calibration period $\mathcal{T}_y$ as we as adjusted exposure to risk related to temperature effects. For that, we use annual mortality data derived from the series of observed numbers of deaths and observed exposure to risk, both available from the HMD (2024). It's noteworthy that the data series for observed exposure to risk are solely accessible on an annual basis. Further details on these data sources is provided in Section 6.1.
In tropical regions, besides temperature, rainfall plays also a crucial role. During monsoon periods, diseases like malaria, diarrheal diseases, or vector-borne illnesses spread more rapidly, significantly affecting women (Egondi et al., 2012). These diseases, identified by the World Health Organization (WHO) as major concerns related to climate change, contribute significantly to mortality (WHO, 2014).

Additionally, the human body's cooling mechanisms, primarily reliant on perspiration, can be compromised if ambient humidity levels are either too low or too high. From a physiological standpoint, scientific studies have highlighted the potential health hazards associated with extreme temperatures in highly humid environments (Davis et al., 2016). However, epidemiological evidence regarding the consequences of wet-hot events seems remain less conclusive (Armstrong et al., 2019), while Fang et al. (2023) indicating an increased risk of mortality in the presence of dry-hot events.

## 2.3 Temperature related deaths : past and expectation

Over the period 1998-2017, the WHO estimates the number of heatwave induced deaths at over 166,000 worldwide, including 70,000 deaths in Europe in the summer of 2003 alone (WHO, 2017). The number of people exposed to heatwaves is also expected to continue to rise. In Europe, it is 9.3 million and would reach almost 300 million in 2100 in a +3°C scenario (Naumann et al., 2020). The annual number of deaths linked to heatwaves is currently estimated at more than 61,000 in Europe during the summer in 2022 and would be more than 120,000 in 2050 (Ballester et al., 2023). However, it is important to note that its results strongly depend on the definition given for a heatwave, see Appendix A. In European countries and Japan, where the proportion of older people in the population is increasing, the damage from heat could be even more severe in the long term. Indeed, the average temperature of the planet has already increased by +1.1°C (Pörtner et al., 2022, TS.C.6) since the pre-industrial era and is expected to continue to rise under any scenario. Thus, extreme weather events such as heatwaves are expected to become more frequent while at the same time the probability of cold waves is expected to decrease (IPCC, 2023). The seasonality of deaths observed today in temperate climates could strongly evolved. Climate change will reduce extreme cold spells and winter seasons period (IPCC, 2023). This would conduct to reduce the cold exposure and cold fatalities (Madaniyazi et al., 2021; Naumann et al., 2020) when most of the temperature related deaths was caused by cold period (Gasparrini et al., 2015). In some parts of Europe, heat-associated health effects may even outweigh the effects of cold (Gronlund et al., 2018; Martínez-Solanas et al., 2021).

So far, no real mitigation or change in the distribution of annual deaths has already been observed in temperate countries (Cheng et al., 2019), which might explain a lesser focus in demographic or actuarial literature.

# 3 A multi-population mortality model with temperature effects

In this section, we explain our approach to incorporating temperature effects into our mortality model. After introducing the data and notations in Section 3.1, Section 3.2 presents the proposed dynamics for central mortality rates excluding temperature effects, followed by the Poisson formulation used to estimate the model parameters. Finally, Section 3.3 describes the DLNM model used to estimate temperature-attributable deaths.

## 3.1 Notation, death counts and temperature-attributable deaths

This paper considers two distinct time scales, daily and annual, as well as two human populations (women and men) in Metropolitan France. In the following, some notations are introduced for modeling mortality. For clarity, bold notations are used for vectors and matrices.

Let us introduce the notations for an annual mortality model that captures long-term trends by age and period in mortality rates. This model is estimated using yearly data observed during a calibration
In the following sections, our focus lies on the dynamics of the annual number of deaths, denoted as $D_{x,t}^{(g)}$. We assume that they can be decomposed into two components

$$D_{x,t}^{(g)} = \tilde{D}_{x,t}^{(g)} + \bar{D}_{x,t}^{(g)}, \qquad (1)$$

where $\tilde{D}_{x,t}^{(g)}$ represents the virtual number of deaths at age $x$ and for year $t$ without considering temperature effects, while $\bar{D}_{x,t}^{(g)}$ corresponds to the number of deaths attributable to temperature variations in year $t$ at age $x$ for gender $g$.

Additionally by dividing by $E_{x,t}^{(g)}$, we denote $\tilde{m}_{x,t}^{(g)}$ and $\bar{m}_{x,t}^{(g)}$ as the components of the crude central death rate excluding and including temperature effects, respectively, such as

$$\hat{m}_{x,t}^{(g)} = \tilde{m}_{x,t}^{(g)} + \bar{m}_{x,t}^{(g)}.$$

This relationship reveals two causes of mortality, with the mortality central death rate $\bar{m}_{x,t}^{(g)}$ associated with the temperature component of mortality. It captures a portion of the observed seasonal variations in daily deaths. We can also express it in a multiplicative form

$$\hat{m}_{x,t}^{(g)} = \frac{\tilde{m}_{x,t}^{(g)}}{1 - \theta_{x,t}^{(g)}}, \qquad (2)$$

where $\theta_{x,t}^{(g)}$ is the total attributable fraction related to temperatures in year $t$ at age $x$ for gender $g$

$$\theta_{x,t}^{(g)} = \frac{\bar{D}_{x,t}^{(g)}}{D_{x,t}^{(g)}}. \qquad (3)$$

This fraction will be calculated latter based on the estimation of deaths $\bar{D}_{x,t}^{(g)}$ attributable to temperatures.

## 3.2 Existing extrapolative mortality models with exogenous effects and shocks

The Lee and Carter (1992) model stands as the most famous model for extrapolating mortality rates, see Basellini et al. (2022) for a review. This model assumes that the logarithm of mortality rates follows a stochastic dynamic learned from past mortality trends. Using previously introduced notations and omitting gender for simplicity in this section, the central death rates are modeled as

$$\ln(\hat{m}_{x,t}) = \alpha_x + \beta_x \kappa_t + \epsilon_{x,t}, \qquad (4)$$

where $\alpha_x$ and $\beta_x$ represent age-specific effect parameters, $\kappa_t$ captures time-temporal factors reflecting the mortality trend in the logarithm of mortality rates over a past period, and $\epsilon_{x,t} \sim N(0, \sigma_{\epsilon}^2)$, $\sigma_{\epsilon} > 0$ denotes residual terms. To ensure model identifiability, the following constraints can, for instance, be imposed during model estimation

$$\sum_{t \in T_y} \kappa_t = 0 \text{ and } \sum_{x \in \mathcal{X}} \beta_x = 1. \qquad (5)$$

This model relies on a two-step estimation process: first estimating the parameters of the model (4), and then estimating a time series model to describe the dynamics of $\hat{\kappa}_t$. For the first step, Lee and Carter (1992) estimate the parameters $\alpha_x$ as the averages of log mortality central rates over the past period. Parameters $\beta_x$ and $\kappa_t$ can be calibrated using a singular value decomposition (SVD) method. Alternatively, Brouhns et al. (2002) assume that the number of deaths follows a Poisson distribution expressed as

$$D_{x,t} \sim \text{Pois}(E_{x,t} \exp(\alpha_x + \beta_x \kappa_t)). \qquad (6)$$
period. More precisely, we use annual mortality data derived from the series of observed numbers of deaths and observed exposure to risk, both available from the HMD (2024). Further details on these data sources are provided in Section 6.1.

Let $\mu_{x,t}^{(g)}$, $E_{x,t}^{(g)}$, and $D_{x,t}^{(g)}$ represent, respectively, the force of mortality, observed exposure to risk, and observed number of deaths for sex $g \in \mathcal{G} = \{\text{female, male}\}$ within the age range $[x, x+1)$ and the calendar year $[t, t+1)$, where $x \in \mathcal{X} = \{x_{\min}, \dots, x_{\max}\}$ denotes the set of integer ages, and $t \in \mathcal{T}_y = \{y_{\min}, \dots, y_{\max}\}$ represents the set of calendar years considered as the calibrated period for mortality models. Common stochastic mortality models such as the Lee-Carter model (Lee and Carter, 1992) rely on the crude central death rate of mortality $\hat{m}_{x,t}^{(g)} = D_{x,t}^{(g)}/E_{x,t}^{(g)}$. It is equivalent to the estimated force of mortality $\hat{\mu}_{x,t}^{(g)}$ assuming a piecewise constant force of mortality. Traditional mortality models may also use the mortality rate $\hat{q}_{x,t}^{(g)} \approx 1 - \exp(-\hat{\mu}_{x,t}^{(g)})$ under the classical piecewise constant force of mortality assumption.

In the following sections, our focus is on the dynamics of the annual number of deaths, denoted as $D_{x,t}^{(g)}$. We assume that the annual number of deaths can be decomposed into two components

$$D_{x,t}^{(g)} = \tilde{D}_{x,t}^{(g)} + \bar{D}_{x,t}^{(g)}, \qquad (1)$$

where $\tilde{D}_{x,t}^{(g)}$ represents the number of deaths at age $x$ for sex $g$ and for year $t$ not attributable to temperature effects, i.e. death counts for all causes of death excluding those attributed to temperature. The term $\bar{D}_{x,t}^{(g)}$ corresponds to the number of deaths attributable to temperature variations during the year. Additionally by dividing by $E_{x,t}^{(g)}$, we denote $\bar{m}_{x,t}^{(g)}$ and $\tilde{m}_{x,t}^{(g)}$ as the components of the crude central death rate attributable and not attributable to temperature, respectively, such as

$$\hat{m}_{x,t}^{(g)} = \tilde{m}_{x,t}^{(g)} + \bar{m}_{x,t}^{(g)}. \qquad (2)$$

This relationship reveals two causes of mortality. The central death rate $\bar{m}_{x,t}^{(g)}$ associated with the temperature component of mortality captures a portion of the observed seasonal variations in daily deaths. We can also express the all-cause central death rate in a multiplicative form

$$\hat{m}_{x,t}^{(g)} = \tilde{m}_{x,t}^{(g)} T_{x,t}^{(g)}, \qquad (3)$$

where $T_{x,t}^{(g)} = (1 - AF_{x,t}^{(g)})^{-1}$, and $AF_{x,t}^{(g)}$ is the total attributable fraction related to temperature in year $t$ at age $x$ for sex $g$

$$AF_{x,t}^{(g)} = \frac{\bar{D}_{x,t}^{(g)}}{D_{x,t}^{(g)}}. \qquad (4)$$

To use the decomposition (1) and then calculation of the total attributable fraction (4), it is necessary to establish a count of deaths attributable to temperature (heat, cold, or both). However, to the best of our knowledge, this categorization of deaths is not explicitly recorded. Two approaches are commonly employed in the literature to estimate temperature-attributable deaths. The first approach involves a detailed examination of causes of death, based notably on the International Classification of Diseases (ICD) established by the World Health Organization (WHO). This classification contains codes corresponding to diseases, symptoms, signs, or medical observations. The temperature-attributable fraction can then be estimated for different causes of death, following methodologies developed by the Global Burden of Disease Study, see e.g. Burkart et al. (2021) or Song et al. (2021). The second approach considers all-cause mortality and estimates the component attributable to temperature using a time-series regression model (Bhaskaran et al., 2013). In epidemiology, this method explains short-term variations in an outcome as a function of variations in other factors. It is particularly applied to analyze the relationship between temperature and mortality, based on the approach developed by Gasparrini et al. (2010). Their model serves as a reference methodological framework for numerous studies investigating the seasonality of deaths and the mortality response to external
The estimation can then be conducted via maximum likelihood. Another route to estimate the Lee-Carter model, as suggested by Haberman and Renshaw (2012), relies on the use of mortality improvement rates. This approach has recently garnered interest, as seen in for instance in Hunt and Villegas (2023), due to its better performance (Mitchell et al., 2013). The second step assumes a random walk with drift or an ARIMA model for the series of $\hat{\kappa}_t$ estimated in the first step.

Starting from this basic framework, several strategies have been devised in the literature to incorporate exogenous effects into mortality models. Indeed, it can be intriguing to showcase the specific effect that certain phenomena have had on past trends in order to disrupt projected trends.

First, one approach relies on econometric models using factor analysis. For instance, French and O'Hare (2014) describes the dynamics of the logarithm of death rates as

$$ \hat{m}_{x,t} = \lambda_x F_t + \epsilon_{x,t}, \qquad (7) $$

where $\lambda_x$ represents a vector of loading parameters associated with latent factor vectors $F_t$. This type of approach can integrate various latent factors that could be linked to economic or health indicators to explain mortality trends. Another line of literature focuses on the use of shocks or jumps in mortality rate dynamics. These models have been particularly explored for modeling catastrophes, wars, or epidemics effects on mortality, notably during the COVID-19 pandemic, see Robben et al. (2022) for a review. Among these approaches, Liu and Li (2015) considers adding jumps to the Lee-Carter model

$$ \ln(\hat{m}_{x,t}) = \alpha_x + \beta_x \kappa_t + N_t J_{x,t} + \epsilon_{x,t}, \qquad (8) $$

where $N_t$ is a Bernoulli variable identifying the occurrence of a temporary jump, i.e. over one year, at date $t$, and $J_{x,t}$ represents the intensity of a mortality jump at this date for age $x$ as a multivariate normal vector. This model is estimated by initially identifying outliers in mortality rate series to determine the occurrence and the severity of jumps.

Incorporating jumps whose effects gradually disappear over several years is part of model enhancement. Recently, Robben and Antonio (2024) develop a multi-population approach comprising a regime-switching model with a high-volatility state where Gaussian shocks on the mortality trend can occur. Their model involves maximum likelihood estimation on improvement rates in several steps and requires preprocessing of data labeled as outliers. Another recent example is provided by Goes et al. (2023) who use a Bayesian formulation and consider the impact of shocks that gradually diffuse over several years.

However, these approaches may encounter difficulties in being extended to model the effect of future temperatures on mortality. As described in section 2, temperature variations produce different responses to heat and cold, with effects lasting only a few days or weeks at most. The shocks caused are thus punctual and offset each other throughout the year, making it challenging to distinguish specific temperature outliers based solely on annual statistics. Additionally, heatwave-type shocks are acute and can be partly or totally compensate throughout the year due to harvest effects, i.e. the decrease in mortality observed after a peak of heat-related mortality. This effect has been observed notably in the French population (Pascal et al., 2018), but is not systematically found in all studies. Another important aspect to consider is the increasing impact of climate change, particularly the consequences of rising temperatures and extreme heat events. This trend tends to limit the relevance of models whose shocks are estimated based on previous annual observations.

Currently, very few mortality models incorporating a temperature-related effect have been developed despite the abundance of epidemiological literature on the effects of heat, cold, or heatwaves on mortality, as discussed in Section 2. However, it is worth mentioning the work of Seklecka et al. (2017), who expand the Lee-Carter model by integrating an additional temperature-related factor. Their approach describes the evolution of the logarithm of mortality rates following a logic similar to that proposed by Plat (2009)

$$ \ln(\hat{m}_{x,t}) = b_x^1 + k_t^1 + (\bar{x} - x)k_t^2 + (\bar{x} - x)^+ k_t^3 + ((a - x)^+ + \text{ct}_x(x - a)^+)^2 k_t^4 + \gamma_{t-x} + \epsilon_{x,t}, \qquad (9) $$

where $b_x^1$ represents the shape parameter at age $x$ of log-mortality rates, $k_t^1$, $k_t^2$, and $k_t^3$ are temporal effects applied respectively across all ages, a temporal effect varying by age, and a specific temporal
factors, as discussed in Section 2. It is important to note that in both cases, the number of deaths attributable to temperature remains an estimated quantity and not an observed one. The number of deaths $\tilde{D}_{x,t}^{(g)}$, which is also unobserved, can either be estimated with a model or calculated by the difference between observed all-cause death counts and the estimate of $\bar{D}_{x,t}^{(g)}$.

To precisely estimate the annual death counts $\bar{D}_{x,t}^{(g)}$, we focus on daily mortality and the impact of temperature on intra-annual seasonal variations in mortality. The daily number of temperature-attributable deaths is fitted within the daily calibration period $\mathcal{T}_d = \{d_{\min}, \dots, d_{\max}\}$, which includes all the days encompassed in the set $\mathcal{D}^* = \{1, 2, \dots, 365, (366)\}$ over the entire calibration period $\mathcal{T}_y$, as detailed in Section 4. Specifically, to evaluate the impact of temperature, we gather daily records of the death counts over the calibration period which are sourced from the *Institut national de la statistique et des études économiques* (INSEE, 2020). Hereto, we denote the number of deaths in day $d \in \mathcal{D}^*$ of year $t \in \mathcal{T}_y$ for sex $g \in \mathcal{G}$ as $D_{x,t,d}^{(g)}$. It can be decomposed between attributable and non-attributable death counts, respectively $\bar{D}_{x,t,d}^{(g)}$ and $\tilde{D}_{x,t,d}^{(g)}$, as in Equation (1). It is noteworthy that the data series for observed exposure to risk are typically not available on a daily basis.

## 3.2 Incorporating the effect of temperature into the Li-Lee model

Our approach involves introducing first a baseline stochastic mortality model for crude central death rates not attributable to temperature $\tilde{m}_{x,t}^{(g)}$. A wide range of mortality models is conceivable for a single population or multiple populations. However, these models primarily focus on all-cause mortality and do not differentiate temperature-related effects. Our approach assumes that long term mortality trends are consistently observed in central death rates not attributable to temperature effects. Among the different mortality models, we opt for the Li and Lee (2005) model, also known as the Augmented Common Factor (ACF) model. Here, the dynamics of mortality is defined for one country and two sexes.

Several reasons motivate this choice. Primarily, a model employing a biological coherence principle is favored here for projecting populations by sex. This choice allows us to observe whether the addition of a temperature component at a later stage diverges the mortality between men and women. In the context of multi-population mortality models, other authors emphasize a shared coherence assumption across multiple countries, see Dowd et al. (2011), Kleinow (2015), Li and Lu (2017), and Bergeron-Boucher et al. (2017) among others. The selection of the ACF model is justified by its status as a benchmark model in multi-population projections, building upon the Lee-Carter model. It's indeed one of the most commonly used mortality models across the literature, as highlighted by Antonio et al. (2017) for mortality projections in the Netherlands and Belgium. It is suitable for modeling mortality at all ages, whereas models like CBDX type models (Dowd et al., 2020) are generally better suited for retired populations. In addition, it should be noted that the approach subsequently developed to integrate the effect of temperatures is not dictated by the choice of the mortality model; other specifications based on forces of mortality or mortality rates are indeed entirely possible.

We introduce our multi-population mortality model based on the Li and Lee (2005) model for the crude central death rates not attributable to temperature effects as

$$ \ln(\tilde{m}_{x,t}^{(g)}) = A_x + B_x K_t + \alpha_x^{(g)} + \beta_x^{(g)} \kappa_t^{(g)}. \qquad (5) $$

This model imposes a common unisex trend $K_t$ in the logarithm of crude central death rates over time and considers specific dynamics by sex $\kappa_t^{(g)}$. The parameters $A_x$ and $\alpha_x^{(g)}$ represent level parameters at age $x$ for mortality, while the parameters $B_x$ and $\beta_x^{(g)}$ modulate the trend of the logarithm of death rates by age and by sex respectively.

The assumption of long-term coherence is ensured by the dynamics given to the series $K_t$ and $\kappa_t^{(g)}$

$$ K_t = \delta + K_{t-1} + e_t, \qquad (6) $$
effect for ages below the average age $\bar{x}$. The coefficient $\text{ctx}$ includes the effect of temperature fluctuations for elderly ages (ages over $a$). A specific temporal effect $k_t^4$ accounts for the temperature effect. Finally, $\gamma_{t-x}$ appears as a cohort effect. In this approach, the parameter $a$ is determined by experts, e.g. Seklecka et al. (2017) set $a = 50$ for the United Kingdom, and the coefficient $\text{ctx}_x$ is calculated as the Pearson's correlation coefficient between temperatures and mortality rates for age $x$.

## 3.3 The coherent Li-Lee model for the virtual central death rates

The aforementioned limitations lead us to consider a more flexible framework to isolate the effect of temperatures. Our approach involves introducing a baseline stochastic mortality model without temperature effects. As we see as in Section 4.2, its implementation requires processing daily data to estimate the effects of seasonal temperature fluctuations on annual mortality.

First, we focus on modeling the virtual death rates devoid of temperature effects. A wide range of mortality models is conceivable for a single population or multiple populations. However, these models primarily focus on overall mortality and do not differentiate temperature-related effects. Our approach assumes that long term mortality trends are consistently observed in annual central death rates without temperature effects. Among the different mortality models, we opt for the Li and Lee (2005) model, also known as the Augmented Common Factor (ACF) model. Here, the dynamics of mortality is defined for one country and two genders.

Several reasons motivate this choice. Primarily, a model employing a biological coherence principle is favored here for projecting populations by gender. This choice allows us to observe whether the addition of a temperature component at a later stage diverges the mortality between men and women. In the context of multi-population mortality models, other authors emphasize a shared coherence assumption across multiple countries, see Dowd et al. (2011), Kleinow (2015), Li and Lu (2017), and Bergeron-Boucher et al. (2017) among others. The selection of the ACF model is justified by its status as a benchmark model in multi-population projections, building upon the Lee-Carter model. It's indeed one of the most commonly used mortality models across the literature, as highlighted by Antonio et al. (2017) for mortality projections in the Netherlands and Belgium. It is suitable for modeling mortality at all ages, whereas models like CDBX type models (Dowd et al., 2020) are generally better suited for retired populations. In addition, it should be noted that the approach subsequently developed to integrate the effect of temperatures is not dictated by the choice of the mortality model; other specifications are indeed entirely possible.

Building upon the previously introduced notations, we introduce our multi-population mortality model based on the Li and Lee (2005) model for the crude central death rates $\tilde{m}_{x,t}^{(g)}$ excluding temperature effects

$$ \ln \tilde{m}_{x,t}^{(g)} = A_x + B_x K_t + \alpha_x^{(g)} + \beta_x^{(g)} \kappa_t^{(g)}. \quad (10) $$

This model imposes a common unisex trend $K_t$ in the logarithm of crude central death rates over time and considering specific dynamics by gender $\kappa_t^{(g)}$. The parameters $A_x$ and $\alpha_x^{(g)}$ represent a level parameter at age $x$ for mortality, while the parameters $B_x$ and $\beta_x^{(g)}$ modulate the trend of the logarithm of mortality rates by age and gender, respectively.

The assumption of long-term coherence is ensured by the dynamics given to the series $K_t$ and $\kappa_t^{(g)}$

$$ K_t = \delta + K_{t-1} + e_t, \quad (11) $$

$$ \kappa_t^{(g)} = c^{(g)} + \phi^{(g)} \kappa_{t-1}^{(g)} + r_t^{(g)}. \quad (12) $$

Here, we consider a first-order autoregressive (AR) model with drift for $\kappa_t^{(g)}$ with estimated coefficients $c^{(g)}$ and $\phi^{(g)}$, and a random walk with drift (RWD) for $K_t$ with a drift parameter $\delta$. As proposed by Antonio et al. (2017) and Antonio et al. (2020), the innovation errors $(e_t, r_t^{(f)}, r_t^{(m)})$ is assumed to be a vector of Gaussian white noises with a mean of zero and a variance-covariance matrix $\Sigma$. The
$$ \kappa_t^{(g)} = c^{(g)} + \phi^{(g)} \kappa_{t-1}^{(g)} + r_t^{(g)}. \qquad (7) $$

Here, we consider a first-order autoregressive (AR) model with drift for $\kappa_t^{(g)}$ with estimated coefficients $c^{(g)}$ and $\phi^{(g)}$, and a random walk with drift (RWD) for $K_t$ with a drift parameter $\delta$. As proposed by Antonio et al. (2017) and Antonio et al. (2020), the innovation errors $(e_t, r_t^{(f)}, r_t^{(m)})$ is assumed to be a vector of Gaussian white noises with a mean of zero and a variance-covariance matrix $\Sigma$. The parameters $c^{(g)}$ and $\phi^{(g)}$ are specific to sex $g$, and are estimated as described in Section 4.4. This specification ensures that the mortality of women and men does not diverge in the long run.

For calibrating this model, we apply the classical identifiability constraints for ensuring parameter uniqueness and efficient computations

$$ \sum_{t \in \mathcal{T}_y} K_t = 0 \quad \text{and} \quad \sum_{x \in \mathcal{X}} B_x^2 = 1. \qquad (8) $$

$$ \sum_{t \in \mathcal{T}_y} \kappa_t^{(g)} = 0 \quad \text{and} \quad \sum_{x \in \mathcal{X}} (\beta_x^{(g)})^2 = 1, \quad \text{for } g \in \mathcal{G}. \qquad (9) $$

Finally, building on the specification from Brouhns et al. (2002), we assume that the number of deaths not attributable to temperature at age $x$ and in year $t$ for sex $g$ follows a Poisson distribution

$$ \tilde{D}_{x,t}^{(g)} \sim \text{Pois} \left( E_{x,t}^{(g)} \tilde{m}_{x,t}^{(g)} \right). \qquad (10) $$

This model cannot be estimated directly because these death counts are not observed and depend on how the mortality-related temperature component is defined. Therefore, we reformulate the relation (10) using Equation (4) to reintroduce the (observable) all-cause death counts as $D_{x,t}^{(g)} = T_{x,t}^{(g)} \tilde{D}_{x,t}^{(g)}$. We deduce again a Poisson formulation with log-link function for $D_{x,t}^{(g)}$ as

$$ D_{x,t}^{(g)} \sim \text{Pois} \left( E_{x,t}^{(g)} T_{x,t}^{(g)} \tilde{m}_{x,t}^{(g)} \right). \qquad (11) $$

As we see in Section 4.2, the prior calculation of total attributable fractions requires processing daily data to estimate the effects of seasonal temperature fluctuations on mortality, and then summing these effects on an annual basis to obtain the terms $AF_{x,t}^{(g)}$ and $T_{x,t}^{(g)}$. From Equation (11), the model is estimated as a Poisson generalized linear model (GLM) through maximum likelihood estimation, as described in Section 4.3. We consider $\ln(E_{x,t}^{(g)} T_{x,t}^{(g)})$ as an offset term when we compute the predicted number of deaths based on the GLM structure

$$ \ln \left( \mathbb{E} \left[ D_{x,t}^{(g)} \right] \right) = \ln \left( E_{x,t}^{(g)} \right) + \ln \left( T_{x,t}^{(g)} \right) + \ln \left( \tilde{m}_{x,t}^{(g)} \right), $$

for each age $x$, sex $g$ and year $t$ based on the fitted model from Equation (5).

## 3.3 Modeling daily deaths with the distributed lag non-linear model

We now aim to model the influence of temperature on the daily number of deaths using the distributed lag non-linear model. We begin by describing the general framework of this model in Section 3.3.1, followed in Section 3.3.2 by the precise specification adopted in our study.

### 3.3.1 General framework

The so-called distributed lag models (DLMs) (Almon, 1965) are commonly employed by econometricians and epidemiologists to address the impact of delayed environmental factors on a response variable. Armstrong (2006) and Gasparrini et al. (2010) proposed a more versatile approach known as distributed lag non-linear generalized model (DLNM). This model, based on a time series regression
parameters $c^{(g)}$, $\phi^{(g)}$ and $\phi^{(g)}$ are specific to gender $g$, and are estimated as described in Section 4.3. This specification ensures that the mortality of women and men does not diverge in the long run.

For calibrating this model, we apply the classical identifiability constraints for ensuring parameter uniqueness and efficient computations

$$ \sum_{t \in \mathcal{T}} K_t = 0 \text{ and } \sum_{x \in \mathcal{X}} B_x^2 = 1. \qquad (13) $$

$$ \sum_{t \in \mathcal{T}_y} \kappa_t^{(g)} = 0 \text{ and } \sum_{x \in \mathcal{X}} (\beta_x^{(g)})^2 = 1, \text{ for } g \in \mathcal{G}. \qquad (14) $$

Finally, building on the specification from Brouhns et al. (2002), we assume that the number of virtual deaths at age $x$ and in year $t$ for gender $g$ follows a Poisson distribution

$$ \tilde{D}_{x,t}^{(g)} \sim \text{Pois} \left( E_{x,t}^{(g)} \exp \left( \tilde{m}_{x,t}^{(g)} \right) \right). \qquad (15) $$

This model cannot be estimated directly because virtual deaths are not observed and depend on how we define the mortality-related temperature component. Therefore, we reformulate the distribution of observed deaths knowing the attributable fraction $\theta_{x,t}^{(g)}$ defined in Equation (3). By noting $T_{x,t}^{(g)} = \exp(\tilde{m}_{x,t}^{(g)}) = 1 + \theta_{x,t}^{(g)}$, we deduce a Poisson formulation with log-link function for $D_{x,t}^{(g)}$ as

$$ D_{x,t}^{(g)} \sim \text{Pois} \left( E_{x,t}^{(g)} T_{x,t}^{(g)} \exp \left( \tilde{m}_{x,t}^{(g)} \right) \right). \qquad (16) $$

From Equation (16), the model is estimated as a Poisson generalized linear model through maximum likelihood estimation, as described in Section 4.3, using $\ln(E_{x,t}^{(g)} T_{x,t}^{(g)})$ as an offset term when we compute the predicted number of deaths

$$ \mathbb{E}[D_{x,t}^{(g)}] = \ln(E_{x,t}^{(g)}) + \ln(T_{x,t}^{(g)}) + \tilde{m}_{x,t}^{(g)}. $$

The term $T_{x,t}^{(g)}$ is then calculated by determining the mortality attributable to temperatures.

## 3.4 Modeling daily deaths with the distributed lag non-linear model

We now seek to model the influence of temperature on the daily number of deaths. The so-called distributed lag odels (DLMs) (Almon, 1965) are commonly employed by econometricians and epidemiologists to address the impact of delayed environmental factors on a response variable. From the aggregated temperature-mortality time series data, two traditional models is explored to capture the nonlinear correlation between temperature (U- or V-shaped) and mortality. Muggeo (2008) introduced the constrained segmented distributed lag (CSDL) model. This model assesses the effects of heat and cold using a fixed threshold within a Poisson regression framework. Recently, Pincemin (2021) implement it to measure the effect of temperatures on projected life expectancy of the French population.

Gasparrini et al. (2010) proposed a more versatile approach known as distributed lag non-linear generalized model (DLNM). This model, based on a quasi-Poisson regression, is adept at capturing nonlinear effects and has found extensive use in the epidemiological literature related to mortality and temperature. Notably, it has been employed in studies such as Gronlund et al. (2018), Guo et al. (2018), Vicedo-Cabrera et al. (2018), Lee et al. (2020), and Vicedo-Cabrera et al. (2021), among others. This model is able to reproduce lagged impacts induced by temperatures persist even after several days have elapsed. When a heatwave occurs, the most vulnerable, particularly the elderly, are the first to be affected. This group experiences an immediate increase in mortality, followed by a period of reduced mortality. Conversely, the effects of cold temperatures usually manifest after a delay of several days. Thus, it's crucial to utilize a statistical model that can effectively account for both immediate and delayed effects on a daily basis.
framework, has the ability to capture the short-term and nonlinear effects of explanatory variables on a response variable, as well as the seasonal and long-term trend effects. It has found extensive use in the epidemiological literature related to mortality and temperature, see Gronlund et al. (2018), Guo et al. (2018), Vicedo-Cabrera et al. (2018), Lee et al. (2020), and Vicedo-Cabrera et al. (2021) among others.

The model employs an exposure–lag–response function that captures the lagged impacts of exposure variables, which may persist several days after exposure. Indeed, immediate and delayed effects on a daily basis are critical to consider in the temperature-mortality relationship. Specifically, during a heatwave, the most vulnerable populations, particularly the elderly, are the first to be affected. This group experiences an immediate increase in mortality, possibly followed by a period of reduced mortality called the harvesting effect. Conversely, the effects of cold temperatures usually manifest after a delay of several days.

The DLNM is designed to examine the impact of exposure variables, such as temperature, humidity or air pollution, on univariate time series data of the outcome, e.g. the number of deaths, considering delayed effects and nonlinear relationships. This model describes the relationship between a time series ($Y_d, d \in T$) for a set of indices $T$ as follows

$$ \rho(\mathbb{E}[Y_d]) = \eta + \sum_{j=1}^{J} s_j(x_{d,j}, L; \theta_j) + \sum_{m=1}^{M} r_m(u_{d,m}; \gamma_m), + \sum_{p=1}^{P} h_p(z_{d,p}; \zeta_p), \quad (12) $$

where $\rho(\cdot)$ is a monotonic link function and the outcome variable $Y$ is assumed to belong to the exponential family of distributions (McCullagh and Nelder, 1989). The mean $\mathbb{E}[Y_d]$ is modeled based on different explanatory variables ($x_{d,1}, \dots, x_{d,J}$), ($u_{d,1}, \dots, u_{d,M}$), and ($z_{d,1}, \dots, z_{d,P}$) for $J, M, P \in \mathbb{N}^*$. The terms $\eta, \theta_j, \gamma_m$ and $\zeta_p$ represent the unknown parameters to be estimated.

The functions $s_j(\cdot, \cdot)$, $j \in \{1, \dots, J\}$, represent smooth bi-dimensional functions of the parameters $\theta_j$ that capture delayed effects of past exposures $q_{x_{d,j},d} = (x_{d,j}, x_{d-1,j}, \dots, x_{d-L,j})^\top$ with a lag $L \in \mathbb{N}$. Gasparrini et al. (2010) and Gasparrini (2014) specify these smooth functions using the so-called basis transformation (Wood, 2006), which refers to a space of known functions, e.g. polynomial functions, splines or B-splines. The transformation of a variable $x_{d,j}$ using a basis function is called a basis variable. Following this methodology and omitting the index $j$ to simplify the notation, each function $s_j(\cdot, \cdot)$ can be expressed as a generic function $s(\cdot, \cdot)$

$$ s(x_d, L; \boldsymbol{\theta}) = \int_0^L f \cdot w(x_{d-l}, l; \boldsymbol{\theta}) dl \approx \sum_{l=0}^L f \cdot w(x_{d-l}, l; \boldsymbol{\theta}), \quad (13) $$

where $f \cdot w(\cdot, \cdot)$ is a bi-dimensional integrable function, called the exposure–lag–response function. This function is constructed using the concept of cross-basis function, which corresponds to the bi-dimensional function space obtained by combining two independent sets of basis functions, $f(\cdot)$ and $w(\cdot)$, with respective dimensions $v_x \in \mathbb{N}^*$ and $v_l \in \mathbb{N}^*$. Formally, the function $s(\cdot, \cdot)$ can also be expressed as a tensor product

$$ s(x_d, L; \boldsymbol{\theta}) \approx \boldsymbol{w}_{x,d}^{\top} \boldsymbol{\theta} = (\boldsymbol{1}_{v_x \cdot v_l}^{\top} \mathbf{A}_{x,d}) \boldsymbol{\theta}, \quad (14) $$

where $\mathbf{1}_v$ is a $v$-dimensional vector of 1's, and $v_x \cdot v_l$ represents the dimension of the cross-basis function $f \cdot w$. It is important to note that an intercept term should not be included in $f(\cdot)$ to ensure identifiability. The matrix $\mathbf{A}_{x,d}$ is defined as

$$ \mathbf{A}_{x,d} = (\mathbf{1}_{v_l}^{\top} \otimes \mathbf{R}_{x,d}) \odot (\mathbf{C} \otimes \mathbf{1}_{v_x}^{\top}), $$

where $\mathbf{C}$ is a $(L+1) \times v_l$ matrix of basis variables from $w(\cdot)$, i.e. the transformation of the lag vector $(0, \dots, L)^\top$ using $w(\cdot)$, and $\mathbf{R}_{x,d}$ is a $(L+1) \times v_x$ matrix obtained by applying $f \cdot w(\cdot, \cdot)$ to the vector of exposures $q_{x,d}$. The symbols $\otimes$ and $\odot$ respectively denote the Kronecker product and the Hadamard product. We introduce the matrix of cross-basis functions $\mathbf{W}$ of dimension $|T| \times (v_x \times v_l)$ composed of the terms $w_{x,d}$ from Equation (14).
The DLNM is designed to examine the impact of exposure variables, e.g. temperature or air pollution, on univariate time series data of outcome, e.g. number of deaths, considering delayed effects and nonlinear relationships. However, it isn't explicitly designed to incorporate both age and time dimensions simultaneously. Hence, to address this limitation, we adopt an approach proposed by Vicedo-Cabrera et al. (2019) and introduce age and gender effects by fitting a model for age buckets and gender. To achieve this, we partition the age range $\mathcal{X}$ into $K \in \mathbb{N}^*$ distinct strata $X_k = [x_{k-1}, x_k), k \in \{1, \dots, K\}$. We aggregate the observed number of deaths for gender $g$ within each stratum $X_k$ and denote it as $D_{k,t,d}^{(g)}$ for $k \in \{1, \dots, K\}$ on day $d \in D^* = \{1, 2, \dots, 365, (366)\}$ of year $t$ over the period $\mathcal{T}_d$.

Let $\lambda_{k,t,d}^{(g)} = \mathbb{E}[D_{k,t,d}^{(g)}]$ represent the expected number of deaths for day $d$, $\vartheta_{d,t}$ denote the daily mean temperature of day $d$ for year $t$, $\mathbf{z}_d$ be a vector comprising $P$ explanatory variables such as year, day of the week, or month, and $\mathbf{u}_d$ be a vector of additional confounding variables like pollutants. We assume that $D_{k,t,d}^{(g)}$ follows a Poisson distribution with overdispersion and a canonical log-link. Utilizing these notations, the DLMN is implemented within each sub-group $k \in \{1, \dots, K\}$ as follows

$$ \ln(\lambda_{k,t,d}) = \eta_k + s(\vartheta_{d,t}, L; \theta_k) + \mathbf{u}_d \gamma_k^\top + \sum_{p=1}^{P} h_p(z_{d,p}; \zeta_k), \quad (17) $$

where $s(\vartheta_{d,t}; l, \beta_k)$ denotes a cross-basis (non-linear) function that simultaneously specifies the association between $\vartheta_{d,t}$ and its lag structure over a maximum of $L$ days. The functions $h_p(z_{d,p}; \zeta_k)$ are natural cubic spline capturing the effect of explanatory variables to account for seasonal variations. Confounding variables $\mathbf{u}_d$ are specified using linear effects. The terms $\eta_k, \theta_k, \gamma_k$ and $\zeta_k$ represent the parameters to be estimated. It's important to note that this specification does not incorporate the exposure to risk as an offset factor within the usual Poisson regression model as in Brouhns et al. (2002).

# 4 Estimation procedure

In this section, we employ a four-step methodology to integrate temperature variations into our estimation procedure. The initial stage of our approach draws from the established methodological framework commonly employed in environmental epidemiology, as seen in studies such as Gasparrini et al. (2017) or Vicedo-Cabrera et al. (2019). This step relies on a statistical approach that enables the association between temperature fluctuations and mortality from daily time series data. The subsequent step involves quantifying the surplus mortality and deaths directly attributable to temperature variations for past periods. These deaths are then considering for calculating an attributable fraction that we use to recalibrate the baseline mortality model in the third step. Finally, a time series model is introduced.

## 4.1 Estimating the DLNM model

We estimate the model given by (17) using the methodology described by Gasparrini (2014). For that, we consider a bi-dimensional spline function $s$

$$ s(\vartheta_{d,t}; L, \theta_k) = \int_0^L f \cdot w(\vartheta_{d-l,t}, l; \theta_k) dl \approx \sum_{l=0}^L f \cdot w(\vartheta_{d-l,t}, l; \theta_k), \quad (18) $$

where $f \cdot w$ is a bi-dimensional integrable function, known as the exposure-lag-response function. For instance, this function can be a natural cubic B-splines or B-splines of various degrees. Equation (18) can be expressed as a tensor product using the concept of cross-basis function, as detailed in Gasparrini et al. (2010). Hence, function $s$ represents a surface of relative risk (RR) indicating the increase or decrease in mortality in response to a rise or fall in temperature relative to a certain threshold.
The functions $r_m(\cdot)$, $m \in \{1, \dots, M\}$, in Equation (12) are smooth univariate functions that capture the effects of confounding daily variables $u_{d,m}$. Linear effects can also be incorporated. For instance, Gasparrini et al. (2010) note that air pollution is usually modeled using a linear relationship in studies focusing on the mortality-temperature association.

The outcome variable is likely influenced by exposure variables and daily confounding variables, but it may also experience short-term residual seasonal fluctuations and long-term effects. In the context of time series regression, it is therefore crucial to finely control for the effects of seasonal phenomena, such as work and holiday periods, weekends, and long-term patterns, in order to properly isolate the short-term exposure-outcome association of interest, as emphasized by Bhaskaran et al. (2013). Various approaches can be adopted to control for these effects. The authors notably highlight the advantages of using spline functions of time instead of time stratified model or Fourier terms since they often fail to adequately capture non-seasonal components. Thus, we introduce functions $h_p(\cdot)$, $p \in \{1, \dots, P\}$, as smooth univariate functions of categorical time variables $z_{d,p}$, e.g. day, day of the week, month or year. They capture residual seasonal effects, demographic shifts, or other long-term trends not accounted for by the other covariates.

### 3.3.2 The DLNM model for daily mortality

We now examine the modeling of daily deaths using the DLNM framework presented above. In our study, we use daily average temperature as the environmental factor, along with its delayed effects.

The DLNM is not explicitly designed to incorporate both age and time dimensions simultaneously. Hence, to address this limitation, we adopt an approach proposed by Vicedo-Cabrera et al. (2019) and introduce age and sex effects by fitting a model for age buckets and sex. To achieve this, we partition the age range $\mathcal{X}$ into $K \in \mathbb{N}^*$ distinct strata $\mathcal{X}_k = [x_{k-1}, x_k)$, $k \in \{1, \dots, K\}$. We aggregate the observed number of deaths for sex $g$ within each stratum $\mathcal{X}_k$, and denote it as $D_{k,t,d}^{(g)}$ for $k \in \{1, \dots, K\}$ on day $d \in \mathcal{D}^*$ of year $t \in \mathcal{T}_y$

$$D_{k,t,d}^{(g)} = \sum_{x \in \mathcal{X}_k} D_{x,t,d}^{(g)}$$

Let $\lambda_{k,t,d}^{(g)} = \mathbb{E}[D_{k,t,d}^{(g)}]$ represent the expected number of deaths for day $d$, year $t$, sex $g$ and sub-group $k$. $\vartheta_{d,t}$ denotes the daily average temperature of day $d$ for year $t$. We assume that $D_{k,t,d}^{(g)}$ follows a Poisson distribution with overdispersion and a canonical log-link. Using the general formulation of the DLNM model in Equation (12), we specify the following quasi-Poisson regression model within each subgroup $k \in \{1, \dots, K\}$ and sex $g \in \mathcal{G}$ as follows

$$\ln(\lambda_{k,t,d}^{(g)}) = \eta_k^{(g)} + s(\vartheta_{d,t}, L; \theta_k^{(g)}) + \sum_{p=1}^{P} h_p(z_{d,p}; \zeta_{k,p}^{(g)}), \quad (15)$$

where $\eta_k^{(g)}$, $\theta_k^{(g)}$ and $\zeta_{k,p}^{(g)}$ represent the parameters to be estimated.

The function $s(\vartheta_{d,t}, L; \theta_k^{(g)})$ denotes a cross-basis non-linear function that simultaneously specifies the association between $\vartheta_{d,t}$ and its lag structure over a maximum of $L$ days. This function captures the short-term cumulative association between mortality and temperature over the period $L$. Based on the work of Gasparrini et al. (2017), Lee et al. (2020), and Martínez-Solanas et al. (2021), we consider for the exposure-response curve $f(\cdot)$ a natural cubic B-spline with three internal knots that are placed at the 10th, 75th, and 90th percentiles of the daily average temperature distribution within the observational period. This choice of a cubic spline and knots allows for easier extrapolation of results to extreme temperatures. We also use a natural cubic B-spline for the lag-response curve $w(\cdot)$ with an intercept and three internal knots placed at evenly spaced intervals on the log scale. A maximum of $L = 21$ days is considered to account for the short-term delayed effects of both cold and hot temperatures.
The uncertainty linked to the estimation of parameters $\theta_k$ in the spline function $s$ can be measured through its variance-covariance matrix $V[\theta_k]$. Subsequently, this estimated uncertainty is incorporated into our projections using a parametric bootstrap technique based on the assumption of a multivariate normal distribution for the estimated coefficients $\theta_k$ (Vicedo-Cabrera et al., 2019).

In this paper, we refer to recent literature in environmental epidemiology to measure the temperature-mortality relationship and calibrate this model. For the sake of comparability, our choice for the lag $L$ and the function $s$ is done consistently with this literature, rather than searching for hyperparameters best suited to our data. Based on works such as Gasparrini et al. (2017), Lee et al. (2020), and Martínez-Solanas et al. (2021), the association between mortality and temperature is modeled using a natural cubic spline. Three internal knots are placed at the 10th, 75th, and 90th percentiles of the daily temperature distribution within the observational period. This choice of a cubic spline allows for easier extrapolation of results to extreme temperatures. We also use a natural cubic spline with an intercept and three internal knots placed at evenly spaced intervals on the log scale. A maximum of 21 days is considered to account for delayed effects of both cold and hot temperatures. Finally, we retain two explanatory variables in the model (17): $z_{d,1}$ representing time, fitted as a natural cubic spline $h_1$ with 8 degrees of freedom per year to accommodate seasonality and long-term trends, and $z_{d,2}$ indicating the day of the week as a categorical variable with a linear relationship for $h_2$. The study does not consider any confounding variables $u_d$.

## 4.2 Excess of mortality and attributable risk

Based on the previously estimated model using temperature and daily death series, we calculate the historical count of deaths attributed to temperatures for each day of the year $d \in T_d$, and age $x \in X_k$, $k \in \{1, \dots, K\}$ using the formula

$$ \bar{D}_{x,t,d}^{(g)} = (1 - \exp(-s(\vartheta_d; L, \hat{\theta}_k))) \times D_{x,t,d}^{(g)}, \quad (19) $$

where $\hat{\theta}_k$ is estimated as described in Section 4.1. The model's setup allows for a range of predictions that precisely gauge the influence of heat, cold, or specific events occurring on each days of subset $D \subseteq D^*$. It's important to note that the temperature effect is shared across all ages within the stratum $X_k$.

To determine the total attributed deaths due to temperatures for a calendar year $t \in T_y$, we sum of the contributions from each day of a subperiod $D \subseteq D^*$ for a year $t \in T_y$

$$ \bar{D}_{x,t}^{(g)} = \sum_{d \in D} \bar{D}_{x,t,d}^{(g)} \mathbb{1}_{\{d \in D\}}. $$

In practice, the period $D$ considered could correspond to the entire year, certain targeted months, heatwaves, or days above or below a certain temperature threshold. We can also assess heatwave effects by considering only the four hottest months of the year or on specific heatwave dates. In this paper, we consider $D = D^*$, i.e. deaths attributable to both cold and hot temperatures. The total attributable fraction is then estimated for each age $x \in X_k$ and year $t \in T_y$ as

$$ \theta_{x,t}^{(g)} = \frac{\bar{D}_{x,t}^{(g)}}{D_{x,t}^{(g)}}. \quad (20) $$

## 4.3 Estimating the mortality model for virtual deaths

Now, let's delve into our estimation approach for the stochastic mortality model (10) under constraints (13) and (14) as introduced in Section 3.3. After estimating the attributable fraction (20) over $D$, we calibrate our model using sets of exposure to risk adjusted by temperature effects and death count data

$$ \mathcal{E} = \{ E_{x,t}^{(g)} T_{x,t}^{(g)}, x \in \mathcal{X}, t \in T_y, g \in \mathcal{G} \} $$
In studies focusing on daily mortality, seasonal variations in mortality are observed and well-documented in the literature (Madaniyazi et al., 2022). This seasonality is only partially explained by temperature effects. Moreover, it is also crucial to consider trends in improving life expectancy and changes in exposure to risk when fitting the model using data observed on many years. To control these effects, we consider the day of the week indicator and the day number within the year as the time categorical variables $z_{d,p}$. In our specification, we define the functions $h_p(\cdot)$ following the approaches of Gasparrini et al. (2015), Gasparrini et al. (2017), and Vicedo-Cabrera et al. (2021). Specifically, we assume a linear relationship for the day of the week variable and employ a natural cubic B-splines for the day of the year variable with 8 degrees of freedom per year in the calibration period. Thus, we can reformulate the component related to the residual seasonal and long-term effects in Equation (15) as

$$ \sum_{p=1}^{P} h_p(z_{d,p}; \zeta_{k,p}^{(g)}) = z_{d,1}\zeta_{k,1}^{(g)} + \sum_{t \in T_y} h_t(z_{d,t}; \zeta_{k,t}^{(g)}), $$

where $z_{1,d}$ is the day of the week indicator related to day $d$, and $z_{d,t}$ is the day number within year $t$. This specification captures the effects of certain days of the week, such as weekends, as well as the residual seasonal mortality patterns observed in each year. Since our regression model (15) does not include confounding variables, their effects are implicitly captured within the component related to seasonal and long-term effects. Moreover, it does not explicitly incorporate the exposure to risk as an offset factor within the usual Poisson regression model as in Brouhns et al. (2002). Hence, the effects of variations in exposure to risk are also captured through $\eta_k^{(g)}$ and $h_p(\cdot)$ functions.

# 4 Estimation procedure

In this section, we employ a four-step methodology to integrate temperature variations into our estimation procedure. The initial stage of our approach draws from the established methodological framework commonly employed in environmental epidemiology, as seen in studies such as Gasparrini et al. (2017) or Vicedo-Cabrera et al. (2019). This step relies on a statistical approach that enables the association between temperature fluctuations and mortality from daily time series data. The subsequent step involves quantifying the excess mortality and death counts attributable to temperature variations for past periods. These deaths are then considered for calculating a total attributable fraction that we use to calibrate the mortality model in the third step. Finally, a time series model is introduced.

## 4.1 Estimating the DLNM model

We estimate the model given in Equation (15) using the methodology described by Gasparrini (2014). We first note that the matrix $W$ can be integrated into the design matrix of the regression model specified in Equation (15). Consequently, the estimation of parameters $\theta_k^{(g)}$ for $k \in \{1, \dots, K\}$ and $g \in G$, along with parameters $\eta_k^{(g)}, \zeta_{k,1}^{(g)},$ and $\zeta_{k,t}^{(g)}$ for $t \in T_y$, can be performed using a maximum likelihood method (Wood, 2006). It is important to note that this procedure does not impose constraints on the function $s(\cdot, \cdot)$. As a result, this smooth function may exhibit negative values over limited time intervals, which could lead to the prediction of negative death counts attributable to temperature. A potential solution could be to constrain the model estimation by imposing non-negative values. However, this solution is complex to implement when the model is non-linear and is outside the scope of this paper (Gasparrini, 2014). It is therefore preferable to adopt a specification for the function $s(\cdot, \cdot)$ that minimizes the occurrence of negative values and to verify that the predicted number of negative values remains negligible once the model is estimated.

The uncertainty linked to the estimation of parameters $\theta_k^{(g)}, k \in \{1, \dots, K\}, g \in G$, in the smooth function $s(\cdot, \cdot)$ can be measured through its variance-covariance matrix $V[\theta_k^{(g)}]$. Subsequently, this estimated uncertainty is incorporated into our projections using a parametric bootstrap technique
$$ \mathcal{D} = \{ D_{x,t}^{(g)}, x \in \mathcal{X}, t \in \mathcal{T}_y, g \in \mathcal{G} \}. $$

The estimation of the chosen model involves maximizing the following expression derived from the Poisson log-likelihood

$$ \max_{A_x, B_x, K_t, \alpha_x^{(g)}, \beta_x^{(g)}, \kappa_t^{(g)}} \sum_{x \in \mathcal{X}} \sum_{t \in \mathcal{T}_y} \sum_{g \in \{f, m\}} \left( D_{x,t}^{(g)} \ln(\tilde{m}_{x,t}^{(g)}) - E_{x,t}^{(g)} T_{x,t}^{(g)} \tilde{m}_{x,t}^{(g)} \right), $$

where $\tilde{m}_{x,t}^{(g)} = \exp(A_x + B_x K_t + \alpha_x^{(g)} + \beta_x^{(g)} \kappa_t^{(g)})$. The parameter estimation in the model through maximum likelihood can be calculated following an iterative procedure described, for instance, in Li (2013) or Chen and Millossovich (2018). In practice, we employ a simpler alternative suggested by Brouhns et al. (2002) based on a the conditional maximum-likelihood approach, see e.g. Li (2013) or Robben et al. (2022). Hence, the parameters are determined in two steps. Firstly, we estimate under constraints (13)

$$ \max_{A_x, B_x, K_t} \sum_{x \in \mathcal{X}} \sum_{t \in \mathcal{T}_y} \left( D_{x,t}^{\text{agg}} \ln(\tilde{m}_{x,t}^{\text{agg}}) - E_{x,t}^{\text{agg}} \tilde{m}_{x,t}^{\text{agg}} \right), $$

where $D_{x,t}^{\text{agg}} = D_{x,t}^{(f)} + D_{x,t}^{(m)}$, $E_{x,t}^{\text{agg}} = E_{x,t}^{(f)} T_{x,t}^{(f)} + E_{x,t}^{(m)} T_{x,t}^{(m)}$ and $\tilde{m}_{x,t}^{\text{agg}} = \exp(A_x + B_x K_t)$.

Once these parameters are estimated and using the predicted values of $\tilde{m}_{x,t}^{\text{agg}}$, we secondly estimate the gender-specific parameters again by maximizing a Poisson log-likelihood under the constraints (14)

$$ \forall g \in \mathcal{G}, \max_{\alpha_x^{(g)}, \beta_x^{(g)}, \kappa_t^{(g)}} \sum_{x \in \mathcal{X}} \sum_{t \in \mathcal{T}_y} \left( D_{x,t}^{(g)} \ln(\tilde{m}_{x,t}^{(g)}) - E_{x,t}^{(g)} T_{x,t}^{(g)} \tilde{m}_{x,t}^{(g)} \right). $$

As presented in Section 6.3, the parameters estimated for the mortality model described above on temperature-reprocessed data are quite similar to those estimated on non-adjusted exposure to risk.

## 4.4 Estimating the time series model

The dynamics of our model comprise two components: temperature dynamics and a time series model for forecasting the vector $(K_t, \kappa_t^{(f)}, \kappa_t^{(m)})$. In this article, we assume that temperature dynamics are exogenous information, with projected temperature series provided by climate models presented in Section 6.1. Here, we focus on calibrating the time series model described in Equations (11) and (12), rewritten in matrix form for all $t \in \mathcal{T}_y$ as

$$ Y_t = D + \Phi Y_{t-1} + E_t, \qquad (21) $$

where

$$ Y_t = \begin{pmatrix} K_t \\ \kappa_t^{(f)} \\ \kappa_t^{(m)} \end{pmatrix}, D = \begin{pmatrix} \delta \\ c^{(f)} \\ c^{(m)} \end{pmatrix}, \Phi = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \phi^{(f)} & 0 \\ 0 & 0 & \phi^{(m)} \end{pmatrix} \text{ and } E_t = \begin{pmatrix} e_t \\ r_t^{(f)} \\ r_t^{(m)} \end{pmatrix}. $$

Drawing from Antonio et al. (2017), the parameters $D$, $\Phi$ and $\Sigma$ are estimated through maximum likelihood since $E_t \sim N_3(0, \Sigma)$. This estimation can be performed using the R-package Multi-MoMo (Antonio et al., 2022) by maximizing the three-dimensionnal Gaussian log-likelihood related to $(Y_t)_{t \in \mathcal{T}_y}$.

# 5 Mortality forecast with temperature effect

In this section, we now aim to incorporate the effects of projected mortality. Indeed, daily projections of temperatures from climate models encompass changes in the intensity in both cold and hot periods, as well as in the severity and duration of heatwaves.

We consider the projection period on annual basis $\mathcal{T}_y^{\text{for}} = \{y_{\text{max}} + 1, \dots, y_{\text{max}} + H\}$ where $y_{\text{max}}$ is the last year of the calibration set $\mathcal{T}_y$ and $H \in \mathbb{N}^*$ is the projection horizon.
based on the assumption of a multivariate normal distribution for the estimated coefficients $\theta_k^{(g)}$ (Vicedo-Cabrera et al., 2019).

Appendix C.2 includes an analysis of the model's goodness of fit. As described in Section 3.3.2, we refer to recent literature in environmental epidemiology to measure the temperature-mortality relationship, set some hyperparameters and calibrate this model. More precisely, our choice of lag $L$, the function $s(\cdot, \cdot)$, and $h_p(\cdot)$ is made consistently with the existing literature to ensure comparability, rather than optimizing hyperparameters specifically for our dataset. In particular, the lag parameter $L$ modulates the short-term seasonality attributable to temperatures and is challenging to determine. If the lag is too short, the temperature effects may fail to capture the harvesting effect resulting from certain intense heatwaves, as well as the effects of cold waves. Conversely, if the lag is too long, the function $s(\cdot, \cdot)$ may capture other seasonal effects unrelated to the short-term impact of temperatures. Therefore, this hyperparameter is conventionally fixed based on references from the literature to ensure a comparable definition of temperature-attributable mortality. Additional analyses to verify the relevance of hyperparameter choices and check the robustness of our results are provided in Appendix C.3.

## 4.2 Death counts, excess of mortality and attributable risk

Based on the previously estimated model given by Equation (15) using daily average temperature and daily death series on period $T_d$, the estimated daily death counts for each day $d \in D^*$ of year $t \in T_y$, and for each age $x \in \mathcal{X}_k$, $k \in \{1, \dots, K\}$ and $g \in G$ is given by

$$ \hat{D}_{k,t,d}^{(g)} = \exp \left( s(\vartheta_{d,t}, L; \hat{\theta}_{k}^{(g)}) \right) \exp \left( \hat{\eta}_{k}^{(g)} + z_{d,1} \hat{\zeta}_{k,1}^{(g)} + \sum_{t \in T_y} h_t(z_{t,p}; \hat{\zeta}_{k,t}^{(g)}) \right) \quad (16) $$

$$ = \exp \left( s(\vartheta_{d,t}, L; \hat{\theta}_{k}^{(g)}) \right) \hat{D}_{k,t,d}^{(g)}, \quad (17) $$

where $\hat{D}_{k,t,d}^{(g)}$ corresponds to the estimated number of deaths not attributable to temperature, and $\hat{\theta}_k^{(g)}$, $\hat{\zeta}_{k,1}^{(g)}$, $\hat{\zeta}_{k,t}^{(g)}$ and $\hat{\eta}_k^{(g)}$ are estimated as described in Section 4.1. The term $\exp(s(\vartheta_{d,t}, L; \hat{\theta}_k^{(g)})$ corresponds to the classical concept of Relative Risk (RR) and can be calculated using Equation (13). We deduce the estimated attributable fraction for day $d$ as

$$ 1 - \exp \left( - \sum_{l=0}^{L} f \cdot w(\vartheta_{d-l,t}, l; \hat{\theta}_{k}^{(g)}) \right), $$

since the attributable fraction satisfies the classical relationship $AF = (RR - 1)/RR$. It's important to note that the temperature effect is shared across all ages within the stratum $\mathcal{X}_k$.

This formulation corresponds to an estimate of attributable risk from a backward perspective, i.e. assessing the current risk based on exposure to temperatures on past dates. Following Gasparini and Leone (2014), we adopt an alternative estimator of this attributable fraction, based on a forward perspective. It is defined for day $d$ as

$$ \widehat{AF}_{k,t,d}^{(g)} = 1 - \exp \left( - \sum_{l=0}^{L} f \cdot w(\vartheta_{d,t}, l; \hat{\theta}_{k}^{(g)}) \right), \quad (18) $$

and can be interpreted on a given date $d$ as the cumulative future effects resulting from a temperature variation on that date. With this alternative, the death counts attributable to temperature effects is estimated based on the mean of daily observed death counts $D_{x,t,d}^{(g)}$ for future dates in period $t, \dots, t+L$

$$ \hat{D}_{x,t,d}^{(g)} = \widehat{AF}_{k,t,d}^{(g)} \times \sum_{l=0}^{L} \frac{D_{x,t,d+l}^{(g)}}{L+1}. \quad (19) $$
## 5.1 Forecasting the virtual mortality without temperature effect

Firstly, we describe the approach adopted to simulate a sample of central death rate trajectories $\tilde{m}_{x,t}^{(g)}$ for $x \in \mathcal{X}$, $g \in \mathcal{G}$ et $t \in \mathcal{T}_y^{\text{for}}$. For this purpose, we use Equation (10) with the estimated parameters $\hat{A}_x$, $\hat{B}_x$, $\hat{K}_t$, $\hat{\alpha}_x^{(g)}$, $\hat{\beta}_x^{(g)}$ and $\hat{\kappa}_t^{(g)}$ obtained as described in Section 4.3.

The mortality model for virtual deaths also involves the projection of the vector $\mathbf{Y}_t$, derived from Equation (21) using the estimated parameters $\hat{D}$, $\hat{\Phi}$ et $\hat{\Sigma}$, see Section 4.4. The dynamics of $\mathbf{Y}_t$ are iteratively deduced by simulating the innovation errors $\mathbf{E}_t$ from a multivariate Gaussian vector with zero mean and covariance matrix $\hat{\Sigma}$. Thus, we can obtain a Monte Carlo sample of trajectories $\tilde{m}_{x,t}^{(g)}$, allowing us to derive prediction intervals for any date $t \in \mathcal{T}_y^{\text{for}}$.

## 5.2 Forecasting temperature-Related Central Death Rates

Secondly, we focus on adding shocks to mortality resulting from daily temperature variations. Based on the usual methodology in ecological epidemiology (Vicedo-Cabrera et al., 2019), the excess of deaths can be deduced by using scenarios of temperatures available on a daily basis, see Section 6.1. Thus, considering the daily temperature trajectory ($\vartheta_d$) given by an external climate scenario, we deduce the daily number of deaths attributed to temperatures for age $x \in X_k$, $k \in \{1, \dots, K\}$, gender $g \in \mathcal{G}$ and for any day $d \in D^*$ of the year $t \in \mathcal{T}_y^{\text{for}}$ using Equation (19). We rewrite it simply as

$$ \bar{D}_{x,t,d}^{(g)} = \tilde{D}_{x,t,d}^{(g)} \left( \exp \left( s(\vartheta_{d,t}; L, \hat{\theta}_k) \right) - 1 \right). $$

Thus, the cumulative number of deaths over a subperiod $D \subseteq D^*$ for a year $t \in \mathcal{T}_y^{\text{for}}$ is expressed as

$$ \bar{D}_{x,t}^{(g)} = \sum_{d \in D} \bar{D}_{x,t,d}^{(g)} \mathbb{1}_{\{d \in D\}}. $$

It is worth noting that a period $D$ different to $D^*$ can be considered when forecasting central mortality rates. This allows us to specifically measure the projected effects of hot and cold temperatures on mortality, distinguishing the compensations that occur throughout the year.

By using Equation (1), we obtain the following expression for the projected central mortality rates with temperature effects measured over period $D$ of year $t \in \mathcal{T}_y^{\text{for}}$

$$ \begin{aligned} \hat{m}_{x,t}^{(g)} &= \tilde{m}_{x,t}^{(g)} + \bar{m}_{x,t}^{(g)} \\ &= \tilde{m}_{x,t}^{(g)} + \frac{\tilde{D}_{x,t}^{(g)}}{E_{x,t}} \sum_{d \in D} \frac{\bar{D}_{x,t,d}^{(g)}}{\tilde{D}_{x,t}^{(g)}} \mathbb{1}_{\{d \in D\}} \\ &= \tilde{m}_{x,t}^{(g)} \left[ 1 + \sum_{d \in D} \omega_{x,t,d}^{(g)} \left( \exp \left( s(\vartheta_{d,t}; L, \hat{\theta}_k) \right) - 1 \right) \mathbb{1}_{\{d \in D\}} \right]. \end{aligned} \quad (22) $$

In Equation (22), the term $\tilde{m}_{x,t}^{(g)}$ corresponds to the forecasted central mortality rate without temperature effects. These effects are then added multiplicatively. The term $w_{x,t,d}^{(g)} = \tilde{D}_{x,t,d}^{(g)}/\tilde{D}_{x,t}^{(g)}$ corresponds to the weight associated with the distribution of virtual deaths over the period $D$ of year $t$. These weights are unknown for future periods, but can be fixed by introducing a distribution assumption of deaths throughout the year. As the seasonality in death counts is mainly captured by temperature effects, it is reasonable to consider a uniform distribution for $\tilde{D}_{x,t,d}^{(g)}$. Under this usual assumption in actuarial literature, the hazard rate function is constant over the year and the weight can be fixed as $w_{x,t,d} = 1/D^*$. Other assumptions about the distribution of deaths could be considered, e.g. the historical distribution of virtual deaths observed empirically over a past period.

Furthermore, we recognize in expression (22) the projected total attributable fraction to temperature related to period $D$, which we denote

$$ \theta_{x,t}^{(g)}(D) = \sum_{d \in D} \omega_{x,t,d}^{(g)} \left( \exp \left( s(\vartheta_{d,t}; L, \hat{\theta}_k) \right) - 1 \right) \mathbb{1}_{\{d \in D\}}. \quad (23) $$
This forward definition offers several advantages. It provides a simpler interpretation of the results and significantly reduces the occurrence of daily negative deaths attributable to temperature <fnref n="1" />. Moreover, it enables an additive representation of the risk attributable to different ranges of temperatures. However, this specification introduces a slight negative bias in the estimation of attributable fractions.

Our model's setup allows for a range of predictions that precisely gauge the influence of heat, cold, or specific events occurring on each day of subset $D_t \subseteq D^*$ of year $t$. To determine the total attributable number of deaths to temperatures for a calendar year $t \in T_y$, we sum of the contributions from each day of a subperiod $D_t \subseteq D^*$

$$ \hat{D}_{x,t}^{(g)} = \sum_{d \in D_t} \hat{D}_{x,t,d}^{(g)} 1_{\{d \in D_t\}} $$

Finally, we estimate the total attributable fraction as defined in Equation (4) for each age $x \in \mathcal{X}_k$ and year $t \in T_y$ as

$$ \widehat{\mathrm{AF}}_{x,t}^{(g)}(\mathcal{D}_t) = \frac{\hat{D}_{x,t}^{(g)}}{\sum_{d \in \mathcal{D}_t} D_{x,t,d}^{(g)}}, \qquad (20) $$

as well as the temperature adjustment

$$ \hat{T}_{x,t}^{(g)}(\mathcal{D}_t) = \left(1 - \widehat{\mathrm{AF}}_{x,t}^{(g)}(\mathcal{D}_t)\right)^{-1}. $$

This adjustment can be interpreted as a cumulative RR over period $D_t$.

In practice, the period $D_t$ considered could correspond to the entire year, certain targeted months, heatwaves, or days above or below a certain temperature threshold. We can also assess heatwave effects by considering only the four hottest months of the year or on specific heatwave dates. To estimate the stochastic mortality model in Section 4.3, we consider $D_t = D^*$ so that $\tilde{m}_{x,t}^{(g)}$ corresponds to the central death rate not attributable to temperature effects (both hot and cold). Conversely, if $D_t$ is strictly included in $D^*$, the central death rate would partially capture some of the temperature effects.

In Equation (20), it is worth noting that all quantities are derived from daily data. While annual death counts can be obtained by summing daily deaths, daily exposure to risk is generally unavailable. Thus, an additional annual dataset containing exposures to risk and death counts is required to calculate annual crude central death rates. Before estimating the stochastic mortality model, it is essential to verify that the annual and daily data sources align and provide consistent death counts by age, sex, and year.

## 4.3 Estimating the stochastic mortality model

We now detail our estimation approach for the stochastic mortality model under constraints (8) and (9) as presented in Section 3.2. The dynamics of the Li-Lee model (5) rely on unobserved central death rates, which are fitted, according to Equation (11), using the total attributable fraction (20). After estimating these attributable fractions, we calibrate our model using sets of exposure to risk adjusted for temperature effects and death counts from annual data

$$ \begin{aligned} \mathcal{E} &= \{ E_{x,t}^{(g)} \hat{T}_{x,t}^{(g)}(D^*), x \in \mathcal{X}, t \in T_y, g \in \mathcal{G} \}, \\ \mathcal{C} &= \{ D_{x,t}^{(g)}, x \in \mathcal{X}, t \in T_y, g \in \mathcal{G} \}. \end{aligned} $$

<footnote n="1">¹It is worth noting that negative deaths attributable to temperature are not necessarily anomalous. As explained by Gasparrini and Leone (2014), the prediction of negative deaths may occur in response to certain events, such as heatwaves, and can be interpreted as a "harvesting paradox."</footnote>
It is worth noting that this quantity is not necessarily positive. Indeed, the reduction in cold-related mortality in a context of overall temperature increase can lead to a decrease in mortality during certain periods.

Finally, using Equation (22), we easily obtain an expression for the total mortality rates with and without temperature effects

$$q_{x,t}^{(g)} = 1 - \exp(-\hat{m}_{x,t}^{(g)}), \quad \tilde{q}_{x,t}^{(g)} = 1 - \exp(-\tilde{m}_{x,t}^{(g)}),$$

for all $x \in X_k, g \in G, t \in T_y^{\text{for}}$. Additionally, the periodic life expectancy at age $x$ truncated at age $t_{\text{max}}$, is given by

$$e_{x,t}^{(g)} = \int_{x}^{t_{\text{max}}} e^{-\int_{x}^{t} \mu_{x,u}^{(g)} du} dt \approx \sum_{k=1}^{t_{\text{max}}} \prod_{j=0}^{k-1} (1 - q_{x,t}^{(g)}).$$

In the presence of two causes of mortality, we examine the number of years of life expectancy lost (or gained) due to temperatures for a person of age $x$ at date $t$ due to the temperature effect

$$\Delta e_{x,t}^{(g)} = \int_{x}^{t_{\text{max}}} e^{-\int_{x}^{t} \tilde{\mu}_{x,u}^{(g)} du} dt - \int_{x}^{t_{\text{max}}} e^{-\int_{x}^{t} \mu_{x,u}^{(g)} du} dt \approx \sum_{k=1}^{t_{\text{max}}} \left[ \prod_{j=0}^{k-1} (1 - \tilde{q}_{x,t}^{(g)}) - \prod_{j=0}^{k-1} (1 - q_{x,t}^{(g)}) \right], \quad (24)$$

with the force of mortality $\mu_{x,t}^{(g)} = \tilde{\mu}_{x,t}^{(g)} + \bar{\mu}_{x,t}^{(g)}$.

# 6 Assessing the impact of future temperatures on French mortality projections

Numerical illustrations presented in this case study have been carried out in the R statistical software (R Core Team, 2023). We use the DLNM model from package dlnm (Gasparrini, 2011). This package proposes a set of functions allowing the construction of the crossbasis transformation presented in Section 4.1, as well as many basis plot commands. The mortality projection is performed by adapting the R code from the MultiMoMo package (Antonio et al., 2022).

## 6.1 Data sources

In this section, we provide an illustration of our model using French mortality data for both women and men, both on a daily and annual basis. Daily mortality data are compared to past data on average national temperatures. Additionally, to project the future effect of temperatures, we utilize temperature scenarios from various climate models.

The calibration period chosen for this study corresponds to the range $T_y = \{1980, \dots, 2019\}$. This selection allows us to focus on relatively recent temperature event. This calibration period ensures adaptability conditions comparable to those experienced by populations at the beginning of the forecasting period. Additionally, we choose to use temperature and mortality data up to the end of the year 2019 to avoid any undesirable effects of the Covid-19 pandemic on the mortality results. We concentrate our analysis on the age group $\mathcal{X} = \{0, \dots, 105\}$ available in the HMD. Older age groups are excluded due to the lower volumes of available data, which could compromise the stability of the results.

### 6.1.1 Temperature dataset for past periods

The temperature records come from the GHCN database (Climate Data Record (CDR) program, 2018). This database gathers more than 20 different sources and contains the main climatological information that can be observed at the surface of the globe, station by station. It undergoes several testing phases and validation processes to ensure the accuracy of its data.
The estimation of the chosen model involves maximizing the following expression derived from the Poisson log-likelihood

$$ \max_{A_x, B_x, K_t, \alpha_x^{(g)}, \beta_x^{(g)}, \kappa_t^{(g)}} \sum_{x \in \mathcal{X}} \sum_{t \in \mathcal{T}_y} \sum_{g \in \{f, m\}} \left( D_{x,t}^{(g)} \ln(\tilde{m}_{x,t}^{(g)}) - E_{x,t}^{(g)} \hat{T}_{x,t}^{(g)}(\mathcal{D}^*) \tilde{m}_{x,t}^{(g)} \right), $$

where $\tilde{m}_{x,t}^{(g)} = \exp(A_x + B_x K_t + \alpha_x^{(g)} + \beta_x^{(g)} \kappa_t^{(g)})$. The parameter estimation in the model through maximum likelihood can be calculated following an iterative procedure described, for instance, in Li (2013) or Chen and Millossovich (2018). In practice, we employ a simpler alternative suggested by Brouhns et al. (2002) based on a the conditional maximum-likelihood approach, see e.g. Li (2013) or Robben et al. (2022). Hence, the parameters are determined in two steps. Firstly, we estimate under constraints (8)

$$ \max_{A_x, B_x, K_t} \sum_{x \in \mathcal{X}} \sum_{t \in \mathcal{T}_y} \left( D_{x,t}^{\text{agg}} \ln(\tilde{m}_{x,t}^{\text{agg}}) - E_{x,t}^{\text{agg}} \tilde{m}_{x,t}^{\text{agg}} \right), $$

where $D_{x,t}^{\text{agg}} = D_{x,t}^{(f)} + D_{x,t}^{(m)}$, $E_{x,t}^{\text{agg}} = E_{x,t}^{(f)} \hat{T}_{x,t}^{(f)}(\mathcal{D}^*) + E_{x,t}^{(m)} \hat{T}_{x,t}^{(m)}(\mathcal{D}^*)$ and $\tilde{m}_{x,t}^{\text{agg}} = \exp(A_x + B_x K_t)$.

Once these parameters are estimated and using the predicted values of $\tilde{m}_{x,t}^{\text{agg}}$, we secondly estimate the sex-specific parameters again by maximizing a Poisson log-likelihood under the constraints (9)

$$ \forall g \in \mathcal{G}, \max_{\alpha_x^{(g)}, \beta_x^{(g)}, \kappa_t^{(g)}} \sum_{x \in \mathcal{X}} \sum_{t \in \mathcal{T}_y} \left( D_{x,t}^{(g)} \ln(\tilde{m}_{x,t}^{(g)}) - E_{x,t}^{(g)} \hat{T}_{x,t}^{(g)}(\mathcal{D}^*) \tilde{m}_{x,t}^{(g)} \right). $$

As presented in Section 6.3, the parameters estimated for the mortality model described above on temperature-reprocessed data are quite similar to those estimated on non-adjusted exposure to risk.

## 4.4 Estimating the time series model

The dynamics of our model comprise two components: temperature dynamics and a time series model for forecasting the vector $(K_t, \kappa_t^{(f)}, \kappa_t^{(m)})$. In this article, we assume that temperature dynamics are exogenous information, with projected temperature series provided by climate models presented in Section 6.1. Here, we focus on calibrating the time series model described in Equations (6) and (7), rewritten in matrix form for all $t \in \mathcal{T}_y$ as

$$ \mathbf{Y}_t = \mathbf{\Upsilon} + \mathbf{\Phi} \mathbf{Y}_{t-1} + \mathbf{E}_t, \qquad (21) $$

where

$$ \mathbf{Y}_t = \begin{pmatrix} K_t \\ \kappa_t^{(f)} \\ \kappa_t^{(m)} \end{pmatrix}, \mathbf{\Upsilon} = \begin{pmatrix} \delta \\ c^{(f)} \\ c^{(m)} \end{pmatrix}, \mathbf{\Phi} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \phi^{(f)} & 0 \\ 0 & 0 & \phi^{(m)} \end{pmatrix} \text{ and } \mathbf{E}_t = \begin{pmatrix} e_t \\ r_t^{(f)} \\ r_t^{(m)} \end{pmatrix}. $$

Drawing from Antonio et al. (2017), the parameters $\mathbf{\Upsilon}$, $\mathbf{\Phi}$ and $\mathbf{\Sigma}$ are estimated through maximum likelihood since $\mathbf{E}_t \sim \mathcal{N}_3(0, \mathbf{\Sigma})$. This estimation can be performed using the R-package **MultiMoMo** (Antonio et al., 2022) by maximizing the three-dimensional Gaussian log-likelihood related to $(\mathbf{Y}_t)_{t \in \mathcal{T}_y}$.

# 5 Mortality forecast with temperature effects

In this section, we now aim to incorporate the effects of projected temperatures on mortality. Our framework relies on coupling the estimates of the DLNM model with the Li-Lee model. Although the DLNM model can capture both long-term and seasonal effects, it encounters challenges when used to project mortality. Indeed, the functions $h_p(\cdot)$ are B-spline functions with numerous parameters, making it challenging to extract clear trends for projection purposes. Additionally, unlike the Li-Lee model, this approach is purely data-driven. Thus, it is complex to enforce consistency across the population groups, whether for the K age groups or by sex.
The oldest data in the GHCN database can be up to 175 years old, and in the case of temperatures this poses various problems. The weather stations, over the years, have evolved technically, migrated from one location to another or may have been recalibrated. This is why internal recalibration algorithms and processes have been applied to this database. Finally, from this database, it is possible to select different meteorological stations of the globe.

In this paper, we select 14 stations related to major French cities: Bordeaux, Brest, Caen, Clermont-Ferrand, Dijon, Lille, Lyon, Marseille, Nantes, Paris, Perpignan, Strasbourg, Toulouse and Tours, to cover the French territory and the main part of the population over the calibration period. For each of them, daily records give access to the average, minimum and maximum temperature of each city. These tree indicators are then aggregated to obtain their average at the scale of Metropolitan France, allowing for a comparison with the daily death count series collected at an aggregated level for metropolitan France. Figure 1 presents the distribution of average daily temperatures over the observation period. It illustrates the temperature variations throughout the year as well as the range of values taken for each month of the year. We describe in Appendix A a historical overview of recent heatwaves in France and the definition used for their counting.

<figure><img src="image_2.png" /><figcaption>Figure 1: Distribution of average daily temperatures for each month of the year. The average daily temperatures are calculated city by city, and then an average of 14 cities is considered to derive average daily temperatures for Metropolitan France. Temperatures are expressed in degrees Celsius. Negative temperatures are represented in a blue color range, while positive temperatures are represented in a yellow color range. Calibration period 1980-2019, GHCN data.</figcaption></figure>

## 6.1.2 Mortality datasets

**Daily data.** Our daily mortality data come from a specific inquiry made to the Quetelet-Prodego Diffusion network (INSEE, 2020), a French network of data centers for social sciences. This portal serves as a valuable resource for researchers and students. The data is acquired through a connection with the INSEE and is provided by the *Archives des Données Issues de la Statistique Publique* team.

Our request involves aggregating deaths from the civil registry. Various variables were accessible, including gender, age at death, socio-professional category, urban area, department, and region of death. The selected deaths are limited to metropolitan France, covering the years 1980 to 2019 and are provided on a daily basis. We select the age-period observation $\mathcal{X} = \{0, \dots, 105\}$ $\mathcal{T}_y = \{1980, \dots, 2019\}$ and the death counts at national level.

Figure 2 depicts the evolution of daily death counts over the period and compares them with the
These challenges explain why we do not rely solely on a model based on daily data for mortality projections. Instead, we select a Li-Lee-type model applied to annual data for mortality projections. By adjusting for exposure to risk, we remove the total annual effect attributable to temperature when we train the Li-Lee model. As a result, the terms $K_t$ and $\kappa_t^{(g)}$ estimated in Section 4.3 are no longer influenced by the short-term seasonality of temperature. Finally, the future effects of temperature are projected separately using the DLNM model and subsequently added to the projected central mortality rates. For that, we use daily temperature projections from climate models, which reflect changes in the intensity of both cold and hot periods, as well as variations in the severity and duration of heatwaves.

We consider the projection period on annual basis $T_y^{\text{for}} = \{y_{\text{max}} + 1, \dots, y_{\text{max}} + H\}$ where $y_{\text{max}}$ is the last year of the calibration set $T_y$ and $H \in \mathbb{N}^*$ is the projection horizon.

## 5.1 Forecasting central death rates not attributable to temperatures

Firstly, we describe the approach adopted to simulate a sample of central death rate trajectories $\tilde{m}_{x,t}^{(g)}$ for $x \in \mathcal{X}$, $g \in \mathcal{G}$ and $t \in T_y^{\text{for}}$. For this purpose, we use Equation (5) with the estimated parameters $\hat{A}_x$, $\hat{B}_x$, $\hat{K}_t$, $\hat{\alpha}_x^{(g)}$, $\hat{\beta}_x^{(g)}$ and $\hat{\kappa}_t^{(g)}$ obtained as described in Section 4.3.

The mortality model for central death rates not attributable to temperature effects also involves the projection of the vector $\mathbf{Y}_t$, derived from Equation (21) using the estimated parameters $\hat{\mathbf{Y}}$, $\hat{\mathbf{\Phi}}$ and $\hat{\mathbf{\Sigma}}$, see Section 4.4. The dynamics of $\mathbf{Y}_t$ are iteratively deduced by simulating the innovation errors $\mathbf{E}_t$ from a multivariate Gaussian vector with zero mean and covariance matrix $\hat{\mathbf{\Sigma}}$. Thus, we can obtain a Monte Carlo sample of trajectories $\tilde{m}_{x,t}^{(g)}$, allowing us to derive prediction intervals for any date $t \in T_y^{\text{for}}$.

## 5.2 Forecasting central death rates with temperature effects

Secondly, we focus on adding shocks to mortality resulting from daily temperature variations, based on scenarios of temperatures available on a daily basis, see Section 6.1. Using the usual methodology in ecological epidemiology (Vicedo-Cabrera et al., 2019), our approach involves projecting the annual central mortality rates not attributable to temperatures, and then adjusting them for the impact of future temperatures. This effect can be assessed on a period $\mathcal{D}_t$ different from $\mathcal{D}^*$ to specifically measure the projected impacts of hot and cold temperatures on mortality, distinguishing the compensations that occur throughout the year.

To do this, we first establish a relationship between the annual central mortality rates and the daily death counts. Indeed, the relationship in Equation (4) and the estimator in Equation (20) cannot be used for projection purposes, as future deaths are not observed. We denote $\text{AF}_{x,d,t}^{(g)} = \bar{D}_{x,t,d}^{(g)}/D_{x,t,d}^{(g)}$ for age $x \in \mathcal{X}_k$, $k \in \{1, \dots, K\}$, sex $g \in \mathcal{G}$ and for any day $d \in \mathcal{D}_t$ of the year $t \in T_y^{\text{for}}$. Then, we obtain

$$ 
\begin{aligned}
\tilde{m}_{x,t}^{(g)} &= \tilde{m}_{x,t}^{(g)} + \bar{m}_{x,t}^{(g)} \\
&= \tilde{m}_{x,t}^{(g)} + \frac{\tilde{D}_{x,t}^{(g)}}{E_{x,t}^{(g)}} \sum_{d \in \mathcal{D}_t} \frac{\bar{D}_{x,t,d}^{(g)}}{\tilde{D}_{x,t}^{(g)}} \mathbb{1}_{\{d \in \mathcal{D}_t\}} \\
&= \tilde{m}_{x,t}^{(g)} \left[ 1 + \sum_{d \in \mathcal{D}_t} \frac{\tilde{D}_{x,t,d}^{(g)} \text{AF}_{x,d,t}^{(g)} (1 - \text{AF}_{x,d,t}^{(g)})^{-1}}{\tilde{D}_{x,t}^{(g)}} \mathbb{1}_{\{d \in \mathcal{D}_t\}} \right] \\
&= \tilde{m}_{x,t}^{(g)} \left[ 1 + \sum_{d \in \mathcal{D}_t} \omega_{x,t,d}^{(g)} \text{AF}_{x,d,t}^{(g)} (1 - \text{AF}_{x,d,t}^{(g)})^{-1} \mathbb{1}_{\{d \in \mathcal{D}_t\}} \right].
\end{aligned}
\quad (22) $$

In Equation (22), the term $\omega_{x,t,d}^{(g)} = \tilde{D}_{x,t,d}^{(g)}/\tilde{D}_{x,t}^{(g)}$ corresponds to the weight associated with the distribution of death counts not attributable to top temperature effects over the period $\mathcal{D}_t$ of year $t$. These weights are unknown for future periods, but can be estimated by introducing a distribution
daily average temperatures, for both women and men. This figure highlights the relationship between mortality and temperature by revealing a classical "U" shape: the daily death count tends to decrease with temperature until reaching a minimal level near 20°C. Beyond this average threshold over a day, mortality tends to increase with rising temperatures. The figure highlights a more pronounced association between heat and mortality for women.

<figure><img src="image_3.png" /><figcaption>Figure 2: Representation of the daily death count according to the average daily temperature in France for women and men. The blue trend curve is obtained by smoothing using a generalized additive model with a confidence level of a 95%. The average daily temperatures are calculated city by city, and then an average of 14 cities is considered to derive average daily temperatures for Metropolitan France. Temperatures are expressed in degrees Celsius. Observation period: 1980-2019, based on GHCN data.</figcaption></figure>

This relationship between mortality and temperature leads to observing seasonal cycles in the number of deaths. Figure 3 illustrates the distribution of deaths per month over the observation period for women et men. It highlights both the mortality associated with cold during the winter months, which can be marked by some extreme events, e.g., flu epidemics, and the peaks in mortality associated with heatwaves, notably an extreme event related to the summer heatwave of 2003. This figure underscores the greater variability in the number of deaths among women during the winter and in the height of summer.

**Annual data.** The annual mortality data considered are extracted from the Human Mortality Database (HMD, 2024). It enables the representation of the crude central death rate $\hat{m}_{x,t}^{(g)}$ for each age $x \in \mathcal{X}$, year $t \in \mathcal{T}_y$, and gender $g \in \mathcal{G}$.

### 6.1.3 Climate scenario dataset

In this paper, our daily future temperature projections are derived from twelve climate models for which we examine three climate change scenarios (Representative Concentration Pathway [RCP]2.6, RCP4.5, and RCP8.5) presented in the framework of the 5th IPCC Assessment Report (IPCC, 2014). These pathways depict the evolution of radiative forcing, which represents the difference between incoming solar radiation and outgoing infrared radiation. Radiative forcing is measured in $W.m^{-2}$.

The international CORDEX project aim to coordinate global research in terms of climate projections, its European branch EUROCORDEX made various simulations available in 2014. The outputs of these simulations have been improved and corrected in 2020 and are freely available through the DRIAS (2020) portal. The data available for Europe have a resolution of 7 to 10 km depending on the geographical areas, with different daily variables such as mean, minimum and maximum temper-
assumption of deaths throughout the year. If the seasonality in death counts is primarily driven by temperature effects, it is reasonable to assume that $\hat{D}_{x,t,d}^{(g)}$ follows a uniform distribution. Under this common assumption in actuarial literature, the hazard rate function remains constant throughout the year, and the weight can be fixed as $\omega_{x,t,d}^{(g)} = 1/|\mathcal{D}^*|$. However, if residual seasonality in death counts, i.e. excluding temperature effects, is more significant, alternative assumptions about the distribution of deaths could be considered. For instance, the daily distribution of deaths not attributable to temperature effects estimated over a previous period can be used. However, this implies that the residual seasonality is assumed to remain constant over time (Madaniyazi et al., 2024).

The effect of temperatures along a given trajectory of projected temperatures ($\hat{\vartheta}_{d,t}$) from an external climate scenario is determined by applying Equation (18). Then, using the central death rates not attributable to temperatures, $\hat{\tilde{m}}_{x,t}^{(g)}$ projected as described in Section 5.1, we derive the projected central death rates as

$$ \hat{\tilde{m}}_{x,t}^{(g)} = \hat{\tilde{m}}_{x,t}^{(g)} \left[ 1 + \sum_{d \in \mathcal{D}_t} \omega_{x,t,d}^{(g)} \hat{\mathrm{AF}}_{x,d,t}^{(g)} (1 - \hat{\mathrm{AF}}_{x,d,t}^{(g)})^{-1} \mathbb{1}_{\{d \in \mathcal{D}_t\}} \right]. \quad (23) $$

Furthermore, we recognize in Equation (23) a multiplicative form similar to that of Equation (3), which implies that the projected total attributable fraction to temperature for the period $\mathcal{D}_t$ can be calculated as

$$ \hat{\mathrm{AF}}_{x,t}^{(g)}(\mathcal{D}_t) = \frac{\sum_{d \in \mathcal{D}_t} \omega_{x,t,d}^{(g)} \hat{\mathrm{AF}}_{x,d,t}^{(g)} (1 - \hat{\mathrm{AF}}_{x,d,t}^{(g)})^{-1} \mathbb{1}_{\{d \in \mathcal{D}_t\}}}{1 + \sum_{d \in \mathcal{D}_t} \omega_{x,t,d}^{(g)} \hat{\mathrm{AF}}_{x,d,t}^{(g)} (1 - \hat{\mathrm{AF}}_{x,d,t}^{(g)})^{-1} \mathbb{1}_{\{d \in \mathcal{D}_t\}}}. \quad (24) $$

It is worth noting that this quantity is not necessarily positive. Indeed, the reduction in cold-related mortality in a context of overall temperature increase can lead to a decrease in mortality during certain periods.

Finally, using Equation (23), we easily obtain an expression for the projected all-cause mortality rates and the projected mortality rates not attributable to temperatures

$$ \hat{q}_{x,t}^{(g)} = 1 - \exp(-\hat{\tilde{m}}_{x,t}^{(g)}), \quad \hat{\tilde{q}}_{x,t}^{(g)} = 1 - \exp(\hat{\tilde{m}}_{x,t}^{(g)}), $$

for all $x \in \mathcal{X}_k$, $k \in \{1, \dots, K\}$, sex $g \in \mathcal{G}$ and $t \in \mathcal{T}_y^{\text{for}}$. Additionally, the periodic life expectancy at age $x$ truncated at age $t_{\text{max}}$, is estimated based on

$$ \hat{e}_{x,t}^{(g)} = \sum_{k=1}^{t_{\text{max}}} \prod_{j=0}^{k-1} (1 - \hat{q}_{x,j}^{(g)}). $$

In the presence of two causes of mortality, we examine the number of years of life expectancy lost (or gained) due to temperatures for a person of age $x$ at date $t$ due to the temperature effect

$$ \Delta \hat{e}_{x,t}^{(g)} = \sum_{k=1}^{t_{\text{max}}} \left[ \prod_{j=0}^{k-1} (1 - \hat{q}_{x,j}^{(g)}) - \prod_{j=0}^{k-1} (1 - \hat{\tilde{q}}_{x,j}^{(g)}) \right]. \quad (25) $$

# 6 Assessing the impact of future temperatures on French mortality projections

Numerical illustrations presented in this case study have been carried out in the R statistical software (R Core Team, 2023). We use the DLNM model from package dlnm (Gasparrini, 2011). This package proposes a set of functions allowing the construction of the crossbasis transformation presented in Section 3.3, as well as many basis plot commands. The mortality projection is performed by adapting the R code from the MultiMoMo package (Antonio et al., 2022).
<figure><img src="image_4.png" /><figcaption>Figure 3: Probability density of the number of deaths by month of the year. The mean plus and minus the standard deviation is displayed as a pointrange. Daily average temperatures are calculated for each city, and then an average of 14 cities is used to derive the daily average temperatures for Metropolitan France. Temperatures are measured in degrees Celsius. Observation period: 1980-2019, based on GHCN data.</figcaption></figure>

ature, total precipitation or wind speed. In our study, we extract temperature scenarios generated using global circulation models (GCMs) and regional circulation models (RCMs) described in Table 1. Temperature data is available in France on a daily time step over the period 2006-2100. It undergoes bias correction and statistical downscaling using the ADAMONT method (Verfaillie et al., 2017) to enable projections on an 8 km resolution grid (SAFRAN).

<table><caption>Table 1: List of climate models provided by the DRIAS</caption><thead><tr><th>Short Name</th><th>GCM</th><th>RCM</th><th>RCPs available</th><th>Period</th></tr></thead><tbody><tr><td>CNRM-CM5 / ALADIN63</td><td>CNRM-CM5</td><td>ALADIN63</td><td>RCP8.5, RCP4.5, RCP2.6</td><td>2006-2100</td></tr><tr><td>MPI-ESM / CCLM4-8-17</td><td>MPI-ESM</td><td>CCLM4-8-17</td><td>RCP8.5, RCP4.5, RCP2.6</td><td>2006-2100</td></tr><tr><td>HadGEM2 / RegCM4-6</td><td>HadGEM2</td><td>RegCM4-6</td><td>RCP8.5, RCP2.6</td><td>2006-2099</td></tr><tr><td>EC-EARTH / RCA4</td><td>EC-EARTH</td><td>RCA4</td><td>RCP8.5, RCP4.5, RCP2.6</td><td>2006-2100</td></tr><tr><td>IPSL-CM5A / WRF381P</td><td>IPSL-CM5A</td><td>WRF381P</td><td>RCP8.5, RCP4.5</td><td>2006-2100</td></tr><tr><td>NorESM1 / REMO2015</td><td>NorESM1</td><td>REMO2015</td><td>RCP8.5, RCP2.6</td><td>2006-2100</td></tr><tr><td>MPI-ESM / REMO2009</td><td>MPI-ESM</td><td>REMO2009</td><td>RCP8.5, RCP4.5, RCP2.6</td><td>2006-2100</td></tr><tr><td>HadGEM2 / CCLM4-8-17</td><td>HadGEM2</td><td>CCLM4-8-17</td><td>RCP8.5, RCP4.5</td><td>2006-2099</td></tr><tr><td>EC-EARTH / RACMO22E</td><td>EC-EARTH</td><td>RACMO22E</td><td>RCP8.5, RCP4.5, RCP2.6</td><td>2006-2100</td></tr><tr><td>IPSL-CM5A / RCA4</td><td>IPSL-CM5A</td><td>RCA4</td><td>RCP8.5, RCP4.5</td><td>2006-2100</td></tr><tr><td>CNRM-CM5 / RACMO22E</td><td>CNRM-CM5</td><td>RACMO22E</td><td>RCP8.5, RCP4.5, RCP2.6</td><td>2006-2100</td></tr><tr><td>NorESM1 / HIRHAM5 v3</td><td>NorESM1</td><td>HIRHAM5 v3</td><td>RCP8.5, RCP4.5</td><td>2006-2100</td></tr></tbody></table>

To remain consistent with the observed mortality trends in 14 major cities described in Section 6.1.2, we extract daily averages, minimums, and maximums of temperatures for the points of the SAFRAN grid corresponding to these major cities. Subsequently, we calculate daily averages of these temperature indicators to reconstruct national averages. This simplification is adopted to have aggregated temperature trends enabling a connection with the selected stochastic mortality model at the scale of Metropolitan France. It is important to note that this simplification results in smoothing out certain regional temperature variations. It is conducted because we do not have access to historical daily death records by city or other areas.

Figure 4 depicts temperature trajectories from the three RCP2.6, RCP4.5, and RCP8.5 scenarios, as well as the distribution of indicators measuring the duration, cumulative severity, and intensity of
## 6.1 Data sources

In this section, we provide an illustration of our model using French mortality data for both women and men, both on a daily and annual basis. Daily mortality data are compared to past data on average national temperatures. Additionally, to project the future effect of temperatures, we utilize temperature scenarios from various climate models.

The calibration period chosen for this study corresponds to the range $T_y = \{1980, \dots, 2019\}$. The period $T_d$ includes all the days corresponding to the years in $T_y$. This selection allows us to focus on relatively recent temperature event. This calibration period ensures adaptability conditions comparable to those experienced by populations at the beginning of the forecasting period. Additionally, we choose to use temperature and mortality data up to the end of the year 2019 to avoid any undesirable effects of the Covid-19 pandemic on the mortality results. We concentrate our analysis on the age group $\mathcal{X} = \{0, \dots, 105\}$ available in the HMD. Older age groups are excluded due to the lower volumes of available data, which could compromise the stability of the results.

### 6.1.1 Temperature dataset for past periods

The temperature records come from the GHCN database (Menne et al., 2012a; Menne et al., 2012b). This database gathers more than 20 different sources and contains the main climatic information that can be observed at the surface of the globe, station by station. It undergoes several testing phases and validation processes to ensure the accuracy of its data.

The oldest data in the GHCN database can be up to 175 years old, and in the case of temperatures this poses various problems. The weather stations, over the years, have evolved technically, migrated from one location to another or may have been recalibrated. This is why internal recalibration algorithms and processes have been applied to this database. Finally, from this database, it is possible to select different meteorological stations of the globe.

In this paper, we select 14 stations related to major French cities: Bordeaux, Brest, Caen, Clermont-Ferrand, Dijon, Lille, Lyon, Marseille, Nantes, Paris, Perpignan, Strasbourg, Toulouse and Tours, to cover the French territory and the main part of the population over the calibration period. For each of them, daily records give access to the average, minimum and maximum temperature of each city. These three indicators are then aggregated to obtain their average at the scale of Metropolitan France, allowing for a comparison with the daily death count series collected at an aggregated level for Metropolitan France.<fnref n="2" /> Figure 1 presents the distribution of average daily temperatures over the observation period. It illustrates the temperature variations throughout the year as well as the range of values taken for each month of the year. We describe in Appendix A a historical overview of recent heatwaves in France.

### 6.1.2 Mortality datasets

**Daily data.** Our daily mortality data come from a specific inquiry made to the Quetelet-Prodego Diffusion network (INSEE, 2020), a French network of data centers for social sciences. This portal serves as a valuable resource for researchers and students. The data is acquired through a connection with the INSEE and is provided by the Archives des Données Issues de la Statistique Publique team.

Our request involves aggregating deaths from the civil registry. Various variables were accessible, including sex, age at death, socio-professional category, urban area, department, and region of death. The selected deaths in our inquiry are limited to metropolitan France, covering the years 1980 to

<footnote n="2">²To our knowledge, there is no existing national index of average temperatures weighted by observed mortality data. As an alternative, we have considered a series of daily mean temperatures weighted by the population size in the departments corresponding to the 14 cities studied. Since this series is nearly identical to the one obtained by calculating a simple average of the temperatures across the 14 cities, both for past and projected periods (assuming a constant population), we opted for the simple average in this study for the sake of simplicity.</footnote>
a heatwave for each decade, as defined in Appendix A.1. Figure 4.a illustrates a gradual increase in average annual temperatures, exacerbated for the RCP8.5 scenario from the 2040s-2050s. This overall temperature rise impacts not only the number but also the characteristics of heatwaves, which gain both in duration, severity, and intensity in the RCP4.5 and RCP8.5 scenarios. Specifically, we observe a similar trajectory of these characteristics for the 2020s and 2030s for all three scenarios, followed by an increase in the dangerosity of heatwaves starting from the 2040s in the RCP8.5 scenario. For these indicators (Figure 4.b to Figure 4.d), the increase in median values in the RCP4.5 and RCP8.5 scenarios is also accompanied by greater model volatility over time beyond the 2040s. The RCP2.6 scenario presents a relatively stable trajectory between 2020 and 2100. Additionally, we can observe a slight decrease in median values in the RCP2.6 scenario over the decades 2020 to 2040, with no clear trend thereafter. These observations are consistent with the findings on hot days presented by the IPCC (2022). Indeed, the regions identified as the most affected by this warming are central and eastern North America, central and southern Europe, and the Mediterranean region.

Regarding the RCP2.6 scenario, the HadGEM2/RegCM4-6 model predicts the most heatwaves, with 459 heatwaves recorded over the period 2020-2100. According to this model, the most severe heatwave occurs in 2047, lasting for 30 consecutive days with a cumulative severity of 209.4. In 2097, it predicts a heatwave with an intensity of 16.3. In scenario RCP4.5, the most pessimistic model in terms of heatwaves is the HadGEM2/CCLM4-8-17 model with a total of 1243 heatwaves. According to this model, the most severe heatwave occurs in 2062, lasting for 49 consecutive days with a cumulative severity of 572.0. In 2097, it predicts a heatwave with an intensity of 22.1. Finally, for scenario 8.5, this model is again the most pessimistic and predicts the most severe heatwave in 2088 (lasting for 89 consecutive days with a cumulative severity of 1099.5), as well as the most intense heatwave in 2094 with an intensity of 27.2. In comparison with these thresholds, the 2003 heatwave lasted 12 consecutive days, with a cumulative severity of 92 and an intensity of 9.2.

## 6.2 Estimation of temperature-mortality association with the DLNM

Using the approach outlined in Section 4.1, we calibrate model (18) over the period $T_y$ using daily temperature and death data with a lag of $L = 21$ days. In our study, we consider $K = 4$ age buckets to account for the differing sensitivities of women's and men's ages based on their specific MMTs: 0-64, 65-74, 75-84, and 85+. Sensitivity analysis regarding the parameter $L$ is conducted in the appendix B.2, and does not challenge the obtained results for the body of the temperature distribution. However, we decide to observe responses over 21 days rather than 14 days in order to better capture long-term responses to heat and cold stress.

As explained in Section 4.1, we choose to adopt reference hyperparameters for the calibration of our DLNM models, as they provide a satisfactory fit to the data and facilitate comparisons with the literature. Alternatively, the search for the most accurate hyperparameters given to the data could be easily conducted by examining the Log-likelihood, the Bayes information criterion (BIC), and the Akaike information criterion (AIC). The goodness of fit of our DLNM models is assessed by comparing the distributions of observations and predictions, as well as the residuals of the Deviance, see Appendix B.1.

Figure 5 depicts the cumulative temperature-mortality association estimated with the DLNM model for both women and men, and for each age group, along with the corresponding 95% confidence intervals (CIs). These curves reveal a non-linear relationship between temperature and the Relative Risk (RR) of mortality over a 21-day period. The shapes of the curves are typical, showing sensitivity to both extreme cold and heat variations, above and below an optimal temperature range, i.e. the MMT, typically between 17 and 21°C, with little variation by gender. We observe that this sensitivity increases notably with age and exhibits differences between genders. For individuals over 85 years old, the responses to temperature variations are generally similar between women and men for both extreme cold and heat temperatures. For ages between 65 and 84 years, the RR for cold temperatures at the 1% percentile of the coldest temperatures exceeds the response to heat at the 1% percentile of the hottest temperatures. However, this risk significantly increases during the most extreme heatwaves.
<figure><img src="image_2.png" /><figcaption>Figure 1: Distribution of average daily temperatures for each month of the year. The average daily temperatures are calculated city by city, and then an average of 14 cities is considered to derive average daily temperatures for Metropolitan France. Temperatures are expressed in degrees Celsius. Negative temperatures are represented in a blue color range, while positive temperatures are represented in a yellow color range. Calibration period 1980-2019, GHCN data.</figcaption></figure>

2019 and are provided on a daily basis. We select the age-period observation $\mathcal{X} = \{0, \dots, 105\}$
$\mathcal{T}_y = \{1980, \dots, 2019\}$ and the death counts at national level.

Figure 2 depicts the evolution of daily death counts over the period and compares them with the daily average temperatures, for both women and men. This figure highlights the relationship between mortality and temperature by revealing a classical "U" shape: the daily death count tends to decrease with temperature until reaching a minimal level near 20°C. Beyond this average threshold over a day, mortality tends to increase with rising temperatures. The figure highlights a more pronounced association between heat and mortality for women.

This relationship between mortality and temperature leads to observing seasonal cycles in the number of deaths. Figure 3 illustrates the distribution of deaths per month over the observation period for women and men. It highlights both the mortality associated with cold during the winter months, which can be marked by some extreme events, e.g., flu epidemics, and the peaks in mortality associated with heatwaves, notably an extreme event related to the summer heatwave of 2003. This figure underscores the greater variability in the number of deaths among women during the winter and in the height of summer.

**Annual data.** The annual mortality data considered are extracted from the Human Mortality Database (HMD, 2024). It enables the representation of the crude central death rate $\hat{m}_{x,t}^{(g)}$ for each age $x \in \mathcal{X}$, year $t \in \mathcal{T}_y$, and sex $g \in \mathcal{G}$. We verify the consistency of annual and daily death counts in Appendix B.

6.1.3 Climate scenario dataset

In this paper, our daily future temperature projections are derived from twelve climate models for which we examine three climate change scenarios (Representative Concentration Pathway [RCP]2.6, RCP4.5, and RCP8.5) presented in the framework of the 5th IPCC Assessment Report (IPCC, 2014). These pathways depict the evolution of radiative forcing, which represents the difference between incoming solar radiation and outgoing infrared radiation. Radiative forcing is measured in W.m⁻².
<figure><img src="image_5.png" /><figcaption>Figure 4: Projection of temperatures and heatwaves by RCP scenario in Metropolitan France over the period 2020-2100. (a) Annual average of daily temperatures by RCP scenario. Each curve is calculated as the average across available climate models. The shaded area corresponds to the 95% confidence intervals across all available climate models. (b) Distribution of heatwave duration per decade and per RCP scenario. (c) Distribution of cumulative heatwave severity per decade and per RCP scenario. (d) Distribution of heatwave intensity per decade and per RCP scenario. For each RCP scenario and each decade, the boxplots represent the distribution of duration, cumulative severity, and intensity calculated over the sample consisting of all years of the decade and across all available climate models. Daily average temperatures are calculated for each city, and then an average of 14 cities is used to derive the daily average temperatures for Metropolitan France. Temperatures are measured in degrees Celsius. Based on DRIAS data.</figcaption></figure>

We also note a greater sensitivity of women to high temperatures and, conversely, a higher sensitivity of men to the coldest temperatures in these age groups. Regarding the age group 0-64 years, the curve for response to extreme heat is steeper for men.

Over the period from 1980 to 2019, we finally represent in Figure 6 the temperature attributable fraction for both women and men, defined in Equation (3), and aggregated over all ages to facilitate visualization. This fraction globally oscillates between 6% and 9% over the considered period, except for women where a peak appears for the 2003 extreme heatwave (Charpentier, 2011). This fraction is decomposed into different temperature effects, depending on the observed daily average temperatures.
<figure><img src="image_3.png" /><figcaption>Figure 2: Representation of the daily death counts according to the average daily temperature in France for women and men. The blue trend curve is obtained by smoothing using a generalized additive model with a confidence level of a 95%. The average daily temperatures are calculated city by city, and then an average of 14 cities is considered to derive average daily temperatures for Metropolitan France. Temperatures are expressed in degrees Celsius. Observation period: 1980-2019, based on GHCN and INSEE data.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Probability density of the number of deaths by month of the year. The mean plus and minus the standard deviation is displayed as a pointrange. Daily average temperatures are calculated for each city, and then an average of 14 cities is used to derive the daily average temperatures for Metropolitan France. Temperatures are measured in degrees Celsius. Observation period: 1980-2019, based on GHCN and INSEE data.</figcaption></figure>

The international CORDEX project aim to coordinate global research in terms of climate projections, its European branch EUROCORDEX made various simulations available in 2014. The outputs of these simulations have been improved and corrected in 2020 and are freely available through the DRIAS (2020) portal. The data available for Europe have a resolution of 7 to 10 km depending on the geographical areas, with different daily variables such as mean, minimum and maximum temperature, total precipitation or wind speed. In our study, we extract temperature scenarios generated using global circulation models (GCMs) and regional circulation models (RCMs) described in Table 1, see
<figure><img src="image_6.png" /><figcaption>Figure 5: Cumulative relative risk of mortality over a 21-day period in Metropolitan France calculated for the years 1980-2019 for women (red) and men (blue) across age groups 0-64, 65-74, 75-84, and 85+, with their 95% confidence intervals (shaded areas). The dashed vertical lines represent the quantiles of daily mean temperatures at 1% (blue), 5% (cyan), 95% (orange), and 99% (red) calculated over the period 1980-2019. Daily average temperatures are calculated for each city, and then an average of 14 cities is used to derive the daily average temperatures for Metropolitan France. Temperatures are measured in degrees Celsius. Based on GHCN data.</figcaption></figure>

Hot or cold temperature days corresponds to days *D* with average temperature above or below the MMT. Thus, similarly to Martínez-Solanas et al. (2021), we define moderate cold days (respectively moderate hot days) as cold days where the average temperature is below the 2.5% quantile (respectively above the 97.5% quantile)<fnref n="1" /> of observed average temperatures over the period. Extreme cold (respectively extreme hot) days are those where the average temperature is below the 2.5% quantile (respectively above the 97.5% quantile). The most significant contribution to the attributable fraction for both women and men comes from moderate cold, while the effects of moderate hot temperatures are generally constant and very low. Extreme cold and hot events produce peaks in excess mortality in certain years, such as in 1985 and 2010 for cold effect, or in 2003 for the heatwave effect. We observe from around 2015 what appears to be the beginning of an upward trend in the extreme hot effect for

<footnote n="1">¹Over the observation period, the 2.5% quantile of the average temperature serie is 0.86°C and the 97.5% quantile is 23.44°C.</footnote>
Flato et al. (2013) for a comprehensive presentation of climate models and their evaluations. Temperature data is available in France on a daily time step over the period 2006-2100. It undergoes bias correction and statistical downscaling using the ADAMONT method (Verfaillie et al., 2017) to enable projections on an 8 km resolution grid (SAFRAN).

<table><caption>Table 1: List of climate models provided by the DRIAS</caption><thead><tr><th>Short Name</th><th>GCM</th><th>RCM</th><th>RCPs available</th><th>Period</th></tr></thead><tbody><tr><td>CNRM-CM5 / ALADIN63</td><td>CNRM-CM5</td><td>ALADIN63</td><td>RCP8.5, RCP4.5, RCP2.6</td><td>2006-2100</td></tr><tr><td>MPI-ESM / CCLM4-8-17</td><td>MPI-ESM</td><td>CCLM4-8-17</td><td>RCP8.5, RCP4.5, RCP2.6</td><td>2006-2100</td></tr><tr><td>HadGEM2 / RegCM4-6</td><td>HadGEM2</td><td>RegCM4-6</td><td>RCP8.5, RCP2.6</td><td>2006-2099</td></tr><tr><td>EC-EARTH / RCA4</td><td>EC-EARTH</td><td>RCA4</td><td>RCP8.5, RCP4.5, RCP2.6</td><td>2006-2100</td></tr><tr><td>IPSL-CM5A / WRF381P</td><td>IPSL-CM5A</td><td>WRF381P</td><td>RCP8.5, RCP4.5</td><td>2006-2100</td></tr><tr><td>NorESM1 / REMO2015</td><td>NorESM1</td><td>REMO2015</td><td>RCP8.5, RCP2.6</td><td>2006-2100</td></tr><tr><td>MPI-ESM / REMO2009</td><td>MPI-ESM</td><td>REMO2009</td><td>RCP8.5, RCP4.5, RCP2.6</td><td>2006-2100</td></tr><tr><td>HadGEM2 / CCLM4-8-17</td><td>HadGEM2</td><td>CCLM4-8-17</td><td>RCP8.5, RCP4.5</td><td>2006-2099</td></tr><tr><td>EC-EARTH / RACMO22E</td><td>EC-EARTH</td><td>RACMO22E</td><td>RCP8.5, RCP4.5, RCP2.6</td><td>2006-2100</td></tr><tr><td>IPSL-CM5A / RCA4</td><td>IPSL-CM5A</td><td>RCA4</td><td>RCP8.5, RCP4.5</td><td>2006-2100</td></tr><tr><td>CNRM-CM5 / RACMO22E</td><td>CNRM-CM5</td><td>RACMO22E</td><td>RCP8.5, RCP4.5, RCP2.6</td><td>2006-2100</td></tr><tr><td>NorESM1 / HIRHAM5 v3</td><td>NorESM1</td><td>HIRHAM5 v3</td><td>RCP8.5, RCP4.5</td><td>2006-2100</td></tr></tbody></table>

To remain consistent with the observed mortality trends in 14 major cities described in Section 6.1.2, we extract daily averages, minimums, and maximums of temperatures for the points of the SAFRAN grid corresponding to these major cities. Subsequently, we calculate daily averages of these temperature indicators to reconstruct national averages. This simplification is adopted to have aggregated temperature trends enabling a connection with the selected stochastic mortality model at the scale of Metropolitan France. It is important to note that this simplification results in smoothing out certain regional temperature variations. It is conducted because we do not have access to historical daily death records by city or other areas.

To illustrate the impact of global warming on heatwaves, we introduce several descriptive indicators to characterize each heatwave period. As explained in Appendix A, we define a heatwave in France for a given day $d \in D^*$ as follows

$$ \bar{\vartheta}_{d}^{\max }=\frac{1}{3} \sum_{j=1}^{3} \vartheta_{d-j}^{\max }, \quad \bar{\vartheta}_{d}^{\min }=\frac{1}{3} \sum_{j=1}^{3} \vartheta_{d-j}^{\min } $$

where $\vartheta_{d}^{\max}$ and $\vartheta_{d}^{\min}$ are the aggregated maximum and minimum temperatures in Metropolitan France for day $d$. This 3-day moving average is sometimes referred to as a biometeorological indicator. Thus, a day is classified as a heatwave if the 3-day moving averages of the minimum temperature $\bar{\vartheta}_{d}^{\min}$ and maximum temperature $\bar{\vartheta}_{d}^{\max}$ exceed the respective minimum $r^{\min}$ and maximum $r^{\max}$ thresholds. These thresholds vary by geographical area. At a national scale, we consider $r^{\min} = 17^{\circ}C$ and $r^{\max} = 30^{\circ}C$. For an heatwave period $D_t \subset D^*$ in year $t$, we then introduce:

* the duration corresponds to the number of consecutive days above the alert thresholds, i.e., $|D_t|$,

* the cumulative severity corresponds to the sum of the daily severities $sev_d$ over the heatwave period $D_t$ where

$$ sev_{d,t} = |\vartheta_{d,t}^{\min} - r^{\min}| + |\vartheta_{d,t}^{\max} - r^{\max}|, $$

* the intensity of a heatwave corresponds to

$$ \text{int} = \max_{d \in D_t} \{|\vartheta_d^{\min} - r^{\min}|\} + \max_{d \in D_t} \{|\vartheta_d^{\max} - r^{\max}|\}. $$

These indicators allow for triggering the appropriate level of alert in line with the danger associated with a heatwave, and for evaluating and ranking heatwaves among themselves.
both women and men.

<figure><img src="image_7.png" /><figcaption>Figure 6: Temperature attributable fraction in Metropolitan France, calculated for the years 1980-2019 for both women and men across all age groups. These values, along with their 95% confidence intervals (shaded areas) are calculated through 1,000 Monte Carlo simulations of the DLNM coefficients. Each line corresponds to the average of these simulations. Attributable fraction is expressed in % for all effect, decomposed into moderate cold, moderate hot, extreme cold and extreme hot.</figcaption></figure>

## 6.3 Calibrating and projecting the Li-Lee model

In this section, we present the estimates obtained for the Li-Lee model described in Section 4.3, in order to project the predicted values of the virtual central death rate $\tilde{m}_{x,t}^{(g)}$. To assess the impact of adjusting exposure to risk, we present both the results of the models (15) and (16). In both cases, our estimates of the Li-Lee model are obtained based on annual mortality data over the period $T_y$.

Figure 7 illustrates the dynamics of the parameters of model (10). It first appears the close proximity between the estimated parameters with and without reprocessing exposures to temperature-attributable deaths. In particular, the temporal trend $\hat{K}_t$ is not disrupted by this reprocessing. Although the impact is moderate, the most affected trend is visually that related to $\hat{\kappa}^{(m)}$.

Figure 8 examines the goodness of fit of the estimated Li-Lee model for model (16) by presenting the Pearson residuals of the fitted Poisson model. Although we have chosen a relatively simple mortality model, it adequately captures age and period effects for both female and male populations. We observe that some cohort effects are less well taken into account, especially for the male population. This limitation could potentially be addressed by using a model with a cohort component. However, including this effect generally introduces numerous other issues, such as less robustness than an Lee-Carter model and difficulties in projecting model parameters. Additionally, cohort effect modeling should be justified by arguments supporting the existence of real cohort effects in the data (Hunt and Blake, 2021). This situation is not straightforward concerning French data. For example, Boumezoued et al. (2020) shows a significant reduction in isolated cohort effects on the French population by incorporating information on monthly birth counts. Another reason is the diagonal trends in the female population appear relatively diffuse and limited. Since the issue of cohort effects is beyond the scope of this paper, we choose to retain a model capturing only age and period effects. It is also worth noting that our chosen specification for handling temperature-attributable deaths could easily be adapted to another stochastic mortality model.

The joint dynamics of the parameters $(\hat{K}_t, \hat{\kappa}_t^{(f)}, \hat{\kappa}_t^{(m)})$ are determined using the time series model described in Section 4.4, employing respectively a RWD process and two AR(1) processes. The
Figure 4 depicts temperature trajectories from the three RCP2.6, RCP4.5, and RCP8.5 scenarios, as well as the distribution of indicators measuring the duration, cumulative severity, and intensity of a heatwave for each decade. Figure 4.a illustrates a gradual increase in average annual temperatures, exacerbated for the RCP8.5 scenario from the 2040s-2050s. This overall temperature rise impacts not only the number but also the characteristics of heatwaves, which gain both in duration, severity, and intensity in the RCP4.5 and RCP8.5 scenarios. Specifically, we observe a similar trajectory of these characteristics for the 2020s and 2030s for all three scenarios, followed by an increase in the danger of heatwaves starting from the 2040s in the RCP8.5 scenario. For these indicators (Figure 4.b to Figure 4.d), the increase in median values in the RCP4.5 and RCP8.5 scenarios is also accompanied by greater model volatility over time beyond the 2040s. The RCP2.6 scenario presents a relatively stable trajectory between 2020 and 2100. Additionally, we can observe a slight decrease in median values in the RCP2.6 scenario over the decades 2020 to 2040, with no clear trend thereafter. These observations are consistent with the findings on hot days presented by the IPCC (2022). Indeed, the regions identified as the most affected by this warming are central and eastern North America, central and southern Europe, and the Mediterranean region.

Regarding the RCP2.6 scenario, the HadGEM2/RegCM4-6 model predicts the most heatwaves, with 459 heatwaves recorded over the period 2020-2100. According to this model, the most severe heatwave occurs in 2047, lasting for 30 consecutive days with a cumulative severity of 209.4. In 2097, it predicts a heatwave with an intensity of 16.3. In scenario RCP4.5, the most pessimistic model in terms of heatwaves is the HadGEM2/CCLM4-8-17 model with a total of 1,243 heatwaves. According to this model, the most severe heatwave occurs in 2062, lasting for 49 consecutive days with a cumulative severity of 572.0. In 2097, it predicts a heatwave with an intensity of 22.1. Finally, for scenario RCP8.5, this model is again the most pessimistic and predicts the most severe heatwave in 2088 (lasting for 89 consecutive days with a cumulative severity of 1,099.5), as well as the most intense heatwave in 2094 with an intensity of 27.2. In comparison with these thresholds, the 2003 heatwave lasted 12 consecutive days, with a cumulative severity of 92 and an intensity of 9.2.

## 6.2 Estimation of temperature-mortality association with the DLNM model

Using the approach outlined in Section 4.1, we calibrate model (15) over the period $T_y$ using daily temperature and death data with a lag of $L = 21$ days. In our study, we consider $K = 4$ age buckets to account for the differing sensitivities of women's and men's ages based on their specific MMTs: 0-64, 65-74, 75-84, and 85+.

Figure 5 depicts the cumulative temperature-mortality association estimated with the DLNM model for both women and men, and for each age group, along with the corresponding 95% confidence intervals (CIs).

After centering the values of $\vartheta_{d,t}$ on the MMT threshold, it represents a cumulative RR surface that quantifies the increase or decrease in mortality in response to deviations in temperature above or below this threshold. These curves reveal a non-linear relationship between temperature and the cumulative RR of mortality over a 21-day period. The shapes of the curves are typical, showing sensitivity to both extreme cold and heat variations, above and below an optimal temperature range, i.e. the MMT, typically between 17 and 21°C, with little variation by sex. We observe that this sensitivity increases notably with age and exhibits differences between sexes. For individuals over 85 years old, the responses to temperature variations are generally similar between women and men for both extreme cold and heat temperatures. For ages between 65 and 84 years, the RR for cold temperatures at the 1% percentile of the coldest temperatures exceeds the response to heat at the 1% percentile of the hottest temperatures. However, this risk significantly increases during the most extreme heatwaves. We also note a greater sensitivity of women to high temperatures and, conversely, a higher sensitivity of men to the coldest temperatures in these age groups. Regarding the age group 0-64 years, the curve for response to extreme heat is steeper for men.

Over the period from 1980 to 2019, we finally represent in Figure 6 the temperature attributable fraction for both women and men, defined in Equation (20), and aggregated over all ages to facilitate
<figure><img src="image_8.png" /><figcaption>Figure 7: Estimated parameters ($\hat{A}_x$, $\hat{B}_x$, $\hat{K}_t$, $\hat{\alpha}_x^{(f)}$, $\hat{\beta}_x^{(f)}$, $\hat{\kappa}_t^{(f)}$, $\hat{\alpha}_x^{(m)}$, $\hat{\beta}_x^{(m)}$, $\hat{\kappa}_t^{(m)}$) of the Li-Lee model for the calibration period 1980-2019 and ages between 0-105 for the entire population of Metropolitan France (Common), females (Female), and males (Male). We distinguish between estimated parameters based on all observed deaths $D_{x,t}^{(g)}$ represented by black solid lines and virtual deaths $\tilde{D}_{x,t}^{(g)}$ reprocessed to account for temperature-attributable deaths in blue dashed lines.</figcaption></figure>

estimated parameters from this model are presented in Table 2, both for the calibration performed on the observed data and for the data virtually reprocessed to exclude temperature-attributable deaths. Table 3 describes the correlation matrix of the residuals for these two sets of parameters.

To ensure compliance with the consistency constraint specified by the Lee-Li model, the dynamics of the processes $\hat{\kappa}_t^{(f)}$ and $\hat{\kappa}_t^{(m)}$ must remain stationary, i.e. the absolute value of $\phi^{(g)}$, $g \in \{f, m\}$, must be strictly lower than 1. This constraint is met during the maximum likelihood estimation of the model (21) regarding $\kappa_t^{(m)}$ and $\kappa_t^{(f)}$ for observed deaths. However, it is observed that the values of the coefficients $\phi^{(m)}$ and $\phi^{(f)}$ are very close to 1. Regarding $\phi^{(f)}$ when adjusting exposure to risk to the temperature effect $T_{x,t}^{(g)}$, the parameter values are slightly above 1. This undesirable situation would thus lead to projecting mortality rates for women diverging from the common trend $K_t$ in the long term. To overcome this issue, Antonio et al. (2020) proposes reducing the calibration period of the model or increasing the lag of the autoregressive processes. Other strategies could be proposed, such as using a mortality model based on mortality increments (Hunt and Villegas, 2023) or employing a
<figure><img src="image_5.png" /><figcaption>Figure 4: Projection of temperatures and heatwaves by RCP scenario in Metropolitan France over the period 2020-2100. (a) Annual average of daily temperatures by RCP scenario. Each curve is calculated as the average across available climate models. The shaded area corresponds to the 95% confidence intervals across all available climate models. (b) Distribution of heatwave duration per decade and per RCP scenario. (c) Distribution of cumulative heatwave severity per decade and per RCP scenario. (d) Distribution of heatwave intensity per decade and per RCP scenario. For each RCP scenario and each decade, the boxplots represent the distribution of duration, cumulative severity, and intensity calculated over the sample consisting of all years of the decade and across all available climate models. Daily average temperatures are calculated for each city, and then an average of 14 cities is used to derive the daily average temperatures for Metropolitan France. Temperatures are measured in degrees Celsius. Based on DRIAS data.</figcaption></figure>

visualization. This fraction globally oscillates between 5% and 9% over the considered period, except for women where a peak appears for the 2003 extreme heatwave (Charpentier, 2011). This fraction is decomposed into different temperature effects, depending on the observed daily average temperatures. Hot or cold temperature days corresponds to days $D_t$ with average temperature above or below the MMT. Thus, similarly to Martínez-Solanas et al. (2021), we define moderate cold days (respectively moderate hot days) as cold days where the average temperature is below the 2.5% quantile (respectively above the 97.5% quantile)<fnref n="3" /> of observed average temperatures over the period. Extreme cold

<footnote n="3">³Over the observation period, the 2.5% quantile of the average temperature series is 0.86°C and the 97.5% quantile</footnote>
<figure><img src="image_9.png" /><figcaption>Figure 8: Pearson residuals of the Li-Lee model for the calibration period 1980-2019 and ages between 0-105 for the female and male populations of Metropolitan France. The model is fitted on observed deaths and temperature-adjusted risk exposures.</figcaption></figure>

data-driven model, as seen in Guibert et al. (2019). However, this topic is beyond the scope of this paper, and we choose to address it more simply by constraining the coefficient $\phi^{(g)}$ to be strictly less than 1 during maximum likelihood maximization.

<table><caption>Table 2: Calibrated parameters in the Li-Lee model with observed deaths with and without temperature effects on the period 1980-2019 and age range 0-105.</caption><thead><tr><th>Calibration data</th><th>δ</th><th>c<sup>(m)</sup></th><th>φ<sup>(m)</sup></th><th>c<sup>(f)</sup></th><th>φ<sup>(f)</sup></th></tr></thead><tbody><tr><td>Without temp. effects</td><td>-0.234</td><td>-0.0066</td><td>0.9692</td><td>-0.0266</td><td>0.9963</td></tr><tr><td>With temp. effects</td><td>-0.235</td><td>-0.0034</td><td>0.9705</td><td>-0.0258</td><td>0.9999</td></tr></tbody></table>

<table><caption>Table 3: Correlation matrix of residuals in the Li-Lee model observed deaths with and without temperature effects on the period 1980-2019 and age range 0-105.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Without temp. effects</th><th colspan="3">With temp. effects</th></tr><tr><th>e<sub>t</sub></th><th>r<sub>t</sub><sup>(m)</sup></th><th>r<sub>t</sub><sup>(f)</sup></th><th>e<sub>t</sub></th><th>r<sub>t</sub><sup>(m)</sup></th><th>r<sub>t</sub><sup>(f)</sup></th></tr></thead><tbody><tr><td>e<sub>t</sub></td><td>1.0000</td><td>-0.7393</td><td>0.9449</td><td>1.0000</td><td>-0.5745</td><td>0.9342</td></tr><tr><td>r<sub>t</sub><sup>(m)</sup></td><td>-0.7393</td><td>1.0000</td><td>-0.7152</td><td>-0.5745</td><td>1.0000</td><td>-0.5537</td></tr><tr><td>r<sub>t</sub><sup>(f)</sup></td><td>0.9449</td><td>-0.7152</td><td>1.0000</td><td>0.9342</td><td>-0.5537</td><td>1.0000</td></tr></tbody></table>

Figure 9 finally presents the projection of parameters $\hat{K}_t$, $\hat{\kappa}_t^{(f)}$, and $\hat{\kappa}_t^{(m)}$ over the period 2020-2100, along with the 95% prediction intervals obtained by 1,000 Monte-Carlo simulation of innovation errors. These projections are shown for both models with and without temperature effects. More precisely, the trends without temperature effects correspond to the dynamic of virtual deaths. We first observe that the common trend $K_t$ is virtually identical in terms of median for both parameter estimations. However, as expected, the projected uncertainty for the model built on observed mortality data without temperature corrections is wider than the trend for virtual central death rates. This is explained by adjustments in exposure to risk related to mortality peaks induced by temperature variations. Regarding the parameters $\hat{\kappa}_t^{(f)}$ and $\hat{\kappa}_t^{(m)}$, we notice the similarity between the forecasted series, and the same phenomenon in terms of uncertainty.

## 6.4 Simulating temperatures effects

We now focus on the simulation of attributable fractions to temperatures $\theta_{x,t}^{(g)}(D)$, as defined in equation (23). These attributable fractions depend on a period $D$ that allows for distinguishing the effects of heat and cold over the year. These quantities are also calculated for each age group, gender and calendar year. To simplify notation, we do not include the index corresponding to the Monte-Carlo
<figure><img src="image_6.png" /><figcaption>Figure 5: Cumulative relative risk of mortality over a 21-day period in Metropolitan France calculated for the years 1980-2019 for women (red) and men (blue) across age groups 0-64, 65-74, 75-84, and 85+, with their 95% confidence intervals (shaded areas). Daily average temperatures are calculated for each city, and then an average of 14 cities is used to derive the daily average temperatures for Metropolitan France. Temperatures are measured in degrees Celsius. Based on GHCN and INSEE data.</figcaption></figure>

(respectively extreme hot) days are those where the average temperature is below the 2.5% quantile (respectively above the 97.5% quantile). The most significant contribution to the attributable fraction for both women and men comes from moderate cold, while the effects of moderate hot temperatures are generally constant and very low. Extreme cold and hot events produce peaks in excess mortality in certain years, such as in 1985 and 2010 for cold effect, or in 2003 for the heatwave effect. We observe from around 2015 what appears to be the beginning of an upward trend in the extreme hot effect for both women and men.

Appendix C.1 includes an empirical analysis presenting the predicted all-cause death counts derived from the DLNM model. We visually assess its ability to reproduce the seasonality of daily deaths for women and men and for each age group. This appendix also details the decomposition of seasonal mortality and illustrates that the temperature-attributable component aligns with the evolution of daily average temperatures. As explained in Section 4.1, we choose to adopt reference hyperparameters

<footnote>is 23.44°C.</footnote>
<figure><img src="image_10.png" /><figcaption>Figure 9: Projection of trend parameters $\hat{K}_t$, $\hat{\kappa}_t^{(f)}$, and $\hat{\kappa}_t^{(m)}$ over the period 2020-2100 for the Li-Lee model for observed central death rates (with temperature effects) and central virtual death rates (without temperature effects) for the calibration period 1980-2019 and ages between 0-105 for the entire population of Metropolitan France. The solid lines represent from 2020 the median of the projections. The 95% prediction intervals are constructed based on 1,000 Monte-Carlo simulations.</figcaption></figure>

simulation numbers.

To facilitate visual analysis, we calculate an aggregated attributable fraction to temperatures using the distribution of deaths known at the end of 2019, as follows

$$ \theta_t(D) = \sum_{g \in \mathcal{G}} \sum_{x \in \mathcal{X}} \theta_{x,t}^{(g)}(D) \frac{D_{x,2019}^{(g)}}{D_{2019}}, $$

where $D_{2019} = \sum_{g \in \mathcal{G}} \sum_{x \in \mathcal{X}} D_{x,2019}^{(g)}$. The details of the attributable fractions to temperatures $\theta_{x,t}^{(g)}(D)$ per age bucket and gender are provided in Appendix C.1.

Figure 10 breaks down the evolution of attributable fractions related to hot and cold temperatures according to the RCP scenarios for the period 2020-2100 and all population. The presented 95% confidence intervals combine, for each RCP scenario, the estimation error associated with the coefficients $\theta_k$, measured by 1,000 bootstrap replications, and the uncertainty related to different climate models listed in 1. Similar to Figure 6, we examine the evolution of components associated with moderate cold and hot days and extreme cold and hot days. Firstly, we observe a relatively stable overall effect of mean temperatures until the end of the century in the RCP2.6 scenario and a slightly decreasing effect in the RCP4.5 scenario for women and men. The RCP4.5 scenario is also characterized by a decrease in the effect of moderate cold and a progressive increase in the effect of extreme heat. It is noteworthy that there is an increase in the uncertainty of the effect of extreme heat over time.

The RCP8.5 scenario exhibits a different behavior. Specifically, it shows a more pronounced decrease over time than the RCP4.5 scenario in the effect of moderate cold temperatures, which is the main component of temperature-attributable mortality for both women and men. Concurrently, a marked increase in the effect of extreme heat is observed starting from the 2050s. However, as highlighted in Appendix C.1, this increase is not sufficient to reverse the overall declining trend in the effect of temperatures over the century for the males. For the female population, we observe a resurgence in the effect of temperatures towards the end of the century, as well as greater uncertainty in temperature-attributable mortality.

It is also noteworthy that the effects of temperatures vary according to geographical locations (Martínez-Solanas et al., 2021). To illustrate this pronounced effect between the north and south of France, we simulate the trajectory of attributable fractions for different French cities across the country, namely Brest, Marseille, Paris, Perpignan, and Strasbourg, using data from DRIAS. Detailed results are presented in Appendix C.2. For this purpose, we employ the DLNM model estimated on the daily mortality data of Metropolitan France as presented in Section 6.2. This is an approximation that could be improved by having access to time-series of death counts per city. Indeed, each of the
<figure><img src="image_7.png" /><figcaption>Figure 6: Temperature attributable fraction in Metropolitan France, calculated for the years 1980-2019 for both women and men across all age groups. These values, along with their 95% confidence intervals (shaded areas) are calculated through 1,000 Monte Carlo simulations of the DLNM coefficients. Each line corresponds to the average of these simulations. Attributable fraction is expressed in % for all effect, decomposed into moderate cold, moderate hot, extreme cold and extreme hot.</figcaption></figure>

for the calibration of our DLNM model, as they provide a satisfactory fit to the data and facilitate comparisons with the literature. The goodness of fit of our DLNM models is assessed by comparing the distributions of observations and predictions, as well as the residuals of the deviance, see Appendix C.2. Sensitivity analysis to model hyperparameters is conducted in Appendix C.3 to assess its robustness. Additionally, we check that our model predicts an insignificant number of negative temperature-attributable deaths.

## 6.3 Calibrating and projecting the Li-Lee model

In this section, we present the estimates obtained for the Li-Lee model described in Section 4.3, in order to project the predicted values of the central death rate not attributable to temperatures $\tilde{m}_{x,t}^{(g)}$. To assess the impact of adjusting exposure to risk, we compare the results of the model (5), which incorporates observed death counts and adjusted exposure to risk according to Equation (11), with those of the original Li-Lee model for $\hat{m}_{x,t}^{(g)}$, which is based on observed death counts and unadjusted exposure to risk, i.e.

$$D_{x,t}^{(g)} \sim \text{Pois} \left( E_{x,t}^{(g)} \exp \left( A_x + B_x K_t + \alpha_x^{(g)} + \beta_x^{(g)} \kappa_t^{(g)} \right) \right). \quad (26)$$

In both cases, our estimates of the Li-Lee model are obtained based on annual mortality data over the period $T_y$.

Figure 7 illustrates the dynamics of the parameters of both Li-Lee models. It first appears the close proximity between their estimated parameters. In particular, the temporal trend $\hat{K}_t$ is not disrupted by adjusting exposure to risk. Although the impact is moderate, the most affected trend is visually that related to $\hat{K}^{(m)}$.

Figure 8 examines the goodness of fit of the estimated Li-Lee model (5) by presenting the Pearson residuals of the fitted Poisson model. Although we have chosen a relatively simple mortality model, it adequately captures age and period effects for both female and male populations. We observe that some cohort effects are less well taken into account, especially for the male population. This limitation could potentially be addressed by using a model with a cohort component. However, including this effect generally introduces numerous other issues, such as less robustness than an Lee-Carter model
<figure><img src="image_11.png" /></figure>

Figure 10: Temperature attributable fraction in Metropolitan France, simulated for the years 2020-2100 for both women and men across. These values, along with their 95% confidence intervals (shaded areas) are calculated through 1,000 Monte Carlo simulations of the DLNM coefficients and the ensemble of climate models. Each line corresponds to the average of these simulations. Attributable fraction is expressed in % for all effect, decomposed into moderate cold, moderate hot, extreme cold and extreme hot.

sub-populations residing there may have specific adaptation capacities to hot and cold temperatures, e.g. air conditioning or thermal insulation of buildings, which we do not account for.

These simulations reveal different sensitivities of each of these cities to extreme hot temperatures and moderate cold temperatures. For the city of Brest for example, the contribution of extreme hot temperatures to mortality remains very limited in the latter part of the century, even under the RCP8.5 scenario. On the other hand, the contribution of cold temperatures decreases gradually, which could tend to decrease the contribution of temperatures to mortality in this region. In contrast, the city of Perpignan appears significantly more exposed to extreme hot temperatures in the RCP4.5 scenario and in the RCP8.5 scenario.

## 6.5 Impacts on life-years lost due to temperature

Finally, we examine the effect of temperatures on the loss of life expectancy at birth. Figure 11 shows the loss of life expectancy at birth, as defined in Equation (24), for women and men between 2020 and 2100 due to temperatures in Metropolitan France. We specifically distinguish the overall temperature effects, which account for both hot and cold temperature-related mortality, and the sole effects of extreme hot temperatures. The different trajectories presented correspond to RCP2.6, RCP4.5 and RCP8.5 scenarios. These assessments are conducted under the assumption that populations and healthcare systems do not adapt to temperature changes throughout the century. Indeed, the challenging issue of adaptation and the evolution of the MMT is beyond the scope of this paper as recent literature suggest that populations can adapt to their local environment (Wu et al., 2024).

For the overall temperature effects, we first note that the loss of life expectancy for men is greater (in absolute value) than that for women. This difference is explained by men's higher sensitivity to cold temperatures compared to women (Figure 5), especially in the age groups between 65 and 84 years. In the three RCP scenarios presented, rising temperatures would limit the contribution of cold to mortality. This would conduct to a general reduction of temperature-related life expectancy loss. We also observe in the RCP8.5 scenario a gradual reduction of the gap between women and men from 2050 onward due to increasing temperatures.

At the same time, attention must be paid to the evolution of life expectancy loss related to extreme hot temperatures. These excesses lead to reductions in life expectancy loss that generally do not exceed 0.1 years, except in a few exceptional years, under the RCP2.6 scenario. This loss does not show any particular trend throughout the century. The RCP4.5 scenario shows a slight but
and difficulties in projecting model parameters. Additionally, cohort effect modeling should be justified by arguments supporting the existence of real cohort effects in the data (Hunt and Blake, 2021). This situation is not straightforward concerning French data. For example, Boumezoued et al. (2020) shows a significant reduction in isolated cohort effects on the French population by incorporating information on monthly birth counts. Another reason is the diagonal trends in the female population appear relatively diffuse and limited. Since the issue of cohort effects is beyond the scope of this paper, we choose to retain a model capturing only age and period effects. It is also worth noting that our chosen specification for handling temperature-attributable deaths could easily be adapted to another stochastic mortality model.

<figure><img src="image_8.png" /><figcaption>Figure 7: Estimated parameters ($\hat{A}_x$, $\hat{B}_x$, $\hat{K}_t$, $\hat{\alpha}_x^{(f)}$, $\hat{\beta}_x^{(f)}$, $\hat{\kappa}_t^{(f)}$, $\hat{\alpha}_x^{(m)}$, $\hat{\beta}_x^{(m)}$, $\hat{\kappa}_t^{(m)}$) of the Li-Lee model for the calibration period 1980-2019 and ages between 0-105 for the entire population of Metropolitan France (Common), females (Female), and males (Male). We distinguish between estimated parameters of model (5) represented by black solid lines, and those of model (26) in blue dashed lines.</figcaption></figure>

The joint dynamics of the parameters ($\hat{K}_t$, $\hat{\kappa}_t^{(f)}$, $\hat{\kappa}_t^{(m)}$) are determined using the time series model described in Section 4.4, employing respectively a RWD process and two AR(1) processes. The estimated parameters from this model are presented in Table 2, both for the original Li-Lee model (26) as well as the model with adjusted exposure to risk (5). Table 3 describes the correlation matrix of the residuals for these two sets of parameters.

To ensure compliance with the consistency constraint specified by the Lee-Li model, the dynamics