step5_damage_function.py

"""
Step 5: Damage Function

This module replicates ComputeDamageFunction.R to calculate damage functions
for different global temperature increases.

Original R code from ComputeDamageFunction.R:
# Code to compute "damage function" - i.e. projected global damages by 2100 at different values of mean temperature increase relative to pre-industrial.  
#   These values range from 0.8C above preindustrial (current warming above preindustrial) to 6C above preindustrial (=5.2C above current)

# To match to IAM data without have to interpolated between projected values, we evaluate over the mean temperatures that were evaluated by at least one of the three main IAMs reported in the Revesz et al 2014 data. 
#   This is computed for different regression models, and different underlying socioeconomic scenarios.

#   We do this for regression point estimates for all regression models and socioeconomics scenarios, and then for SSP5 and pooled model calculate for all bootstrapped regression runs to quantify full uncertainty (for figure 5d)

# Read in scenario data generated in ComputeMainProjections.R script
load(file="data/output/projectionOutput/popProjections.Rdata")
load(file="data/output/projectionOutput/growthProjections.Rdata")

# Read in projections of future temperature change, generated by the getTemperatureChange.R script
#   These are population-weighted country level projections averaged across all CMIP5 models
Tchg <- read.csv("data/input/CCprojections/CountryTempChange_RCP85.csv")
Tchg <- merge(popProjections[[1]][,1:3],Tchg,by.x="iso",by.y="GMI_CNTRY")
Tchg <- Tchg[Tchg$CNTRY_NAME%in%c("West Bank","Gaza Strip","Bouvet Island")==F,]  #these have the same iso code as israel and norway, respectively, which messes up the merge
Tconv <- Tchg$Tconv  #these are the "conversion factors":  i.e. what you multiply the global mean temp by to get the country-level change in temp.  again, this is based on RCP8.5 ensemble mean 

scens <- c("base","SSP"%&%1:5)  #vector defining the scenarios
sn = c(1,4,6)  #scenarios in the scens vector we are going to evaluate
yrs <- 2010:2099

# NOW COMPUTE DAMAGE FUNCTION ESTIMATES FOR FOUR MAIN HISTORICAL MODELS AND EACH SOCIOECONOMIC SCENARIO
#   Doing this for different amounts of assumed warming by 2100, as explained in text and SOM

iam <- read.csv("data/input/IAMdata/ProcessedKoppData.csv")
mods <- c("DICE","FUND","PAGE")
iam <- iam[iam$IAM%in%mods,]
inc <- sort(unique(iam$T[iam$T<=6 & iam$T>1]))  #temperatures under which some IAM was run
incs <- c(0.8,inc)  #adding in a zero increase relative to today (=0.8 increase relative to preindustrial)

####################################################################################################
# POOLED MODEL WITH NO LAGS
####################################################################################################

prj <- read.csv("data/output/bootstrap/bootstrap_noLag.csv")  #bootstrapped projections, all obs
np = dim(prj)[1]  #number of bootstrap replicates we ran. the first one is the one from the baseline model

tots = array(dim=c(length(incs),length(sn),2))  #array to hold total global GDP with and without climate change in 2099, for each scenario and temperature increase
dimnames(tots) <- list(incs,scens[sn],c("avgGDPcapCC","avgGDPcapNoCC"))

for (dt in 1:length(incs)) {  
  dtm <- (incs[dt] - 0.8)*Tconv  #calculate the country-specific temperature rise for the particular global mean increase, RELATIVE TO PRE-INDUSTRIAL so we can match IAMs.  
  # We assume an increase from pre-industrial to present day of 0.8C, so a 1C global increase to 2100 relative to pre-industrial means that we still have 1 - 0.8 = 0.2C left to go over the 21st century
  # we then convert this to country-specific estimates of warming, using our conversion factors from RCP8.5
  ccd <- dtm/length(yrs)  #rate of increase in temperature per year.
  
  # now loop over SSP scenarios and our own scenario. 
  for (scen in sn) {
    
    growthproj <- growthProjections[[scen]]
    popproj <- popProjections[[scen]]
    
    basegdp = popproj$gdpCap  #baseline GDP/cap
    temp <- popproj$meantemp #baseline temperature
    
    GDPcapCC = GDPcapNoCC = array(dim=c(dim(growthproj)[1],length(yrs)))  #array to fill with GDP/cap for each country
    dimnames(GDPcapCC) <- dimnames(GDPcapNoCC) <- list(growthproj[,1],yrs)
    GDPcapCC[,1] = GDPcapNoCC[,1] = basegdp  #initialize with baseline per cap GDP

      tt=1    # evaluating our regression point estimate 
      bg = prj$temp[tt]*temp + prj$temp2[tt]*temp*temp  #this finds the predicted growth level for each country's temperature for the particular bootstrap run
      for (i in 2:length(yrs)) {
        j = i - 1
        y = yrs[i]
        basegrowth <- growthproj[,which(names(growthproj)==y)]
        GDPcapNoCC[,i] = GDPcapNoCC[,j]*(1+basegrowth)  #last year's per cap GDP times this years growth rate, as projected by scenario
        newtemp = temp+j*ccd
        dg = prj$temp[tt]*newtemp + prj$temp2[tt]*newtemp*newtemp  #predicted growth under new temperature
        dg[newtemp>30] = prj$temp[tt]*30 + prj$temp2[tt]*30*30  #constrain response to response at 30C if temp goes above that.  this is so we are not projecting out of sample
        
        diff = dg - bg  #difference between predicted baseline growth and predicted growth under new temp
        GDPcapCC[,i] = GDPcapCC[,j]*(1+basegrowth + diff)  #last year's GDPcap (w/ climate change) times climate-adjusted growth rate for this year
      }
      wt = popproj[,which(names(popproj)==y)]  #population weights 
      tots[dt,which(sn==scen),1] <- round(weighted.mean(GDPcapCC[,90],wt),0)
      tots[dt,which(sn==scen),2] <- round(weighted.mean(GDPcapNoCC[,90],wt),0)
      GDPcapCC <- GDPcapNoCC <- NULL
  }  #end scenario loop
  }  # end temperature increase loop
  
save(tots,file="data/output/projectionOutput/DamageFunction_pooled.Rdata")

####################################################################################################
# RICH/POOR MODEL WITH NO LAGS
####################################################################################################

#   We have different response functions for countries that were above or below median income in the historical sample.  Countries move onto
#     the rich-country response function in the future if their income rises above the historical median income (and vice versa if it falls).
#     Depending where they are in the temperature distribution when this happens, the marginal effect of temperature on growth could change sign

prj <- read.csv("data/output/bootstrap/bootstrap_richpoor.csv")  #interacted model, all obs

tots = array(dim=c(length(incs),length(sn),2))  #array to hold total global GDP with and without climate change in 2099, for each scenario and temperature increase
dimnames(tots) <- list(incs,scens[sn],c("avgGDPcapCC","avgGDPcapNoCC"))

for (dt in 1:length(incs)) {  
  dtm <- (incs[dt] - 0.8)*Tconv   
  ccd <- dtm/length(yrs)  #rate of increase in temperature per year.
  
  for (scen in sn) {
    
    growthproj <- growthProjections[[scen]]
    popproj <- popProjections[[scen]]
    
    basegdp = popproj$gdpCap  #baseline GDP/cap
    medgdp <- median(basegdp)
    temp <- popproj$meantemp  #baseline temperature.  
    
    GDPcapCC = GDPcapNoCC = array(dim=c(dim(growthproj)[1],length(yrs)))  #array to fill with GDP/cap for each country
    GDPcapCC[,1] = GDPcapNoCC[,1] = basegdp  #initialize with baseline per cap GDP

  tt=1
  for (i in 2:length(yrs)) {
        j = i - 1
        y <- yrs[i]
        #baseline growth rate from SSP
        basegrowth <- growthproj[,which(names(growthproj)==y)]
        
        #was country poor last year in climate change world?
        poor <- GDPcapCC[,j]<=medgdp
        
        #calculate appropriate predicted growth rates in world without climate change
        bg = prj$temp[tt]*temp + prj$temp2[tt]*temp*temp  #predicted growth level for rich countries absent climate change
        bg[poor] = prj$temppoor[tt]*temp[poor] + prj$temp2poor[tt]*temp[poor]*temp[poor]  #predicted growth level for poor countries absent climate change
        
        GDPcapNoCC[,i] = GDPcapNoCC[,j]*(1+basegrowth)  #last year's per cap GDP times this years growth rate
        newtemp = temp+j*ccd
        dg = prj$temp[tt]*newtemp + prj$temp2[tt]*newtemp*newtemp  #predicted growth under new temperature for rich countries
        dg[poor] = prj$temppoor[tt]*newtemp[poor] + prj$temp2poor[tt]*newtemp[poor]*newtemp[poor]  #predicted growth for poor countries
        
        dg[newtemp>30 & poor==0] = prj$temp[tt]*30 + prj$temp2[tt]*30*30  #constrain response to response at 30C if temp goes above that for rich countries.  this is so we are not projecting out of sample
        dg[newtemp>30 & poor==1] = prj$temppoor[tt]*30 + prj$temp2poor[tt]*30*30  #constrain response to response at 30C for poor countries  
        diff = dg - bg  #difference between predicted baseline growth and predicted growth under new temp
        GDPcapCC[,i] = GDPcapCC[,j]*(1+basegrowth + diff)  #last year's GDPcap (w/ climate change) times climate adjusted growth rate for this year
      }
      wt = popproj[,which(names(popproj)==y)]  #population weights 
      tots[dt,which(sn==scen),1] <- round(weighted.mean(GDPcapCC[,90],wt),0)
      tots[dt,which(sn==scen),2] <- round(weighted.mean(GDPcapNoCC[,90],wt),0)
      GDPcapCC <- GDPcapNoCC <- NULL
    }  #end scenario loop
  }  # end temperature increase loop
  
save(tots,file="data/output/projectionOutput/DamageFunction_richpoor.Rdata")
"""

import pandas as pd
import numpy as np
import logging
from pathlib import Path
import pickle
import warnings
warnings.filterwarnings('ignore')

from config import *

# Set up logging
from config import setup_logging
logger = setup_logging()

class DamageFunction:
    """Class to handle damage function calculations."""
    
    def __init__(self):
        self.pop_projections = {}
        self.growth_projections = {}
        self.temperature_data = None
        self.bootstrap_data = {}
        self.iam_data = None
        
    def load_data(self):
        """Load all required data for damage function calculations."""
        logger.info("Loading data for damage function calculations...")
        
        # Load population and growth projections
        with open(OUTPUT_FILES['pop_projections'], 'rb') as f:
            self.pop_projections = pickle.load(f)
        
        with open(OUTPUT_FILES['growth_projections'], 'rb') as f:
            self.growth_projections = pickle.load(f)
        
        # Load temperature data
        self.temperature_data = pd.read_csv(OUTPUT_FILES['country_temp_change'], encoding='latin-1')
        
        # Load bootstrap data
        bootstrap_files = {
            'pooled_no_lag': OUTPUT_FILES['bootstrap_no_lag'],
            'rich_poor': OUTPUT_FILES['bootstrap_rich_poor']
        }
        
        for model, file_path in bootstrap_files.items():
            if file_path.exists():
                self.bootstrap_data[model] = pd.read_csv(file_path, encoding='latin-1')
                logger.info(f"Loaded bootstrap data for {model}: {len(self.bootstrap_data[model])} runs")
            else:
                logger.warning(f"Bootstrap file not found: {file_path}")
        
        # Load IAM data if available
        if INPUT_FILES['iam_data'].exists():
            self.iam_data = pd.read_csv(INPUT_FILES['iam_data'], encoding='latin-1')
            logger.info(f"Loaded IAM data: {len(self.iam_data)} rows")
        else:
            logger.warning("IAM data not found, using default temperature range")
        
        logger.info("Data loading completed")
    
    def get_temperature_increases(self):
        """Get temperature increases for damage function calculation."""
        logger.info("Setting up temperature increases for damage function...")
        
        if self.iam_data is not None:
            # Use IAM temperature scenarios
            iam_temp = self.iam_data[self.iam_data['IAM'].isin(['DICE', 'FUND', 'PAGE'])]
            temp_increases = sorted(iam_temp['T'].unique())
            temp_increases = [t for t in temp_increases if 1 < t <= 6]
        else:
            # Use default temperature range
            temp_increases = [1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0]
        
        # Add current warming level
        temp_increases = [0.8] + temp_increases
        
        logger.info(f"Temperature increases for damage function: {temp_increases}")
        return temp_increases
    
    def calculate_damage_function_pooled(self):
        """Calculate damage function for pooled model."""
        logger.info("Calculating damage function for pooled model...")
        
        if 'pooled_no_lag' not in self.bootstrap_data:
            logger.warning("Bootstrap data for pooled model not available")
            return
        
        bootstrap_coefs = self.bootstrap_data['pooled_no_lag']
        temp_increases = self.get_temperature_increases()
        
        # Scenarios to run
        scenarios = [0, 3, 5]  # baseline, SSP3, SSP5
        scenario_names = ['base', 'SSP3', 'SSP5']
        
        # Initialize results array
        n_temp = len(temp_increases)
        n_scenarios = len(scenarios)
        damage_results = np.zeros((n_temp, n_scenarios, 2))  # GDP with and without CC
        
        # Get conversion factors
        conversion_factors = self.temperature_data.set_index('GMI_CNTRY')['Tconv']
        
        for temp_idx, temp_increase in enumerate(temp_increases):
            logger.info(f"Processing temperature increase: {temp_increase}°C")
            
            # Calculate country-specific temperature changes
            # Assume 0.8°C current warming, so additional warming is temp_increase - 0.8
            additional_warming = temp_increase - 0.8
            years = PROJECTION_YEARS
            total_years = len(years) - 1
            
            for scenario_idx, scenario in enumerate(scenarios):
                if scenario not in self.pop_projections or scenario not in self.growth_projections:
                    continue
                
                # Get projections
                pop_proj = self.pop_projections[scenario]
                growth_proj = self.growth_projections[scenario]
                
                # Get baseline data
                base_gdp = pop_proj['gdpCap'].values
                base_temp = pop_proj['meantemp'].values
                countries = pop_proj['iso'].values
                
                # Use point estimate (first bootstrap run)
                coefs = bootstrap_coefs.iloc[0]
                temp_coef = coefs['temp']
                temp2_coef = coefs['temp2']
                
                # Calculate baseline growth level
                baseline_growth = temp_coef * base_temp + temp2_coef * base_temp ** 2
                
                # Project to 2099
                gdp_cap_cc = base_gdp.copy()
                gdp_cap_no_cc = base_gdp.copy()
                
                for year_idx in range(1, len(years)):
                    year = years[year_idx]
                    
                    # Get growth rate without climate change
                    growth_rate = growth_proj[year].values / 100
                    
                    # Project GDP without climate change
                    gdp_cap_no_cc = gdp_cap_no_cc * (1 + growth_rate)
                    
                    # Calculate new temperature
                    years_since_2010 = year - 2010
                    warming_per_year = additional_warming / total_years
                    
                    new_temp = base_temp + years_since_2010 * warming_per_year * np.array([
                        conversion_factors.get(iso, 1.0) for iso in countries
                    ])
                    
                    # Constrain temperature to 30°C maximum
                    new_temp = np.minimum(new_temp, MAX_TEMPERATURE)
                    
                    # Calculate growth with climate change
                    climate_growth = temp_coef * new_temp + temp2_coef * new_temp ** 2
                    
                    # Calculate climate impact
                    climate_impact = climate_growth - baseline_growth
                    
                    # Project GDP with climate change
                    gdp_cap_cc = gdp_cap_cc * (1 + growth_rate + climate_impact)
                
                # Calculate global weighted average GDP in 2099
                pop_weights = pop_proj[2099].values
                total_pop = pop_weights.sum()
                
                if total_pop > 0:
                    damage_results[temp_idx, scenario_idx, 0] = np.average(gdp_cap_cc, weights=pop_weights)
                    damage_results[temp_idx, scenario_idx, 1] = np.average(gdp_cap_no_cc, weights=pop_weights)
        
        # Save results
        self._save_damage_function(damage_results, temp_increases, scenario_names, 'pooled')
        
        logger.info("Damage function for pooled model completed")
        return damage_results
    
    def calculate_damage_function_rich_poor(self):
        """Calculate damage function for rich/poor model."""
        logger.info("Calculating damage function for rich/poor model...")
        
        if 'rich_poor' not in self.bootstrap_data:
            logger.warning("Bootstrap data for rich/poor model not available")
            return
        
        bootstrap_coefs = self.bootstrap_data['rich_poor']
        temp_increases = self.get_temperature_increases()
        
        # Scenarios to run
        scenarios = [0, 3, 5]  # baseline, SSP3, SSP5
        scenario_names = ['base', 'SSP3', 'SSP5']
        
        # Initialize results array
        n_temp = len(temp_increases)
        n_scenarios = len(scenarios)
        damage_results = np.zeros((n_temp, n_scenarios, 2))  # GDP with and without CC
        
        # Get conversion factors
        conversion_factors = self.temperature_data.set_index('GMI_CNTRY')['Tconv']
        
        for temp_idx, temp_increase in enumerate(temp_increases):
            logger.info(f"Processing temperature increase: {temp_increase}°C")
            
            # Calculate country-specific temperature changes
            additional_warming = temp_increase - 0.8
            years = PROJECTION_YEARS
            total_years = len(years) - 1
            
            for scenario_idx, scenario in enumerate(scenarios):
                if scenario not in self.pop_projections or scenario not in self.growth_projections:
                    continue
                
                # Get projections
                pop_proj = self.pop_projections[scenario]
                growth_proj = self.growth_projections[scenario]
                
                # Get baseline data
                base_gdp = pop_proj['gdpCap'].values
                base_temp = pop_proj['meantemp'].values
                countries = pop_proj['iso'].values
                median_gdp = np.median(base_gdp)
                
                # Use point estimate (first bootstrap run)
                coefs = bootstrap_coefs.iloc[0]
                temp_coef = coefs['temp']
                temp_poor_coef = coefs['temppoor']
                temp2_coef = coefs['temp2']
                temp2_poor_coef = coefs['temp2poor']
                
                # Project to 2099
                gdp_cap_cc = base_gdp.copy()
                gdp_cap_no_cc = base_gdp.copy()
                
                for year_idx in range(1, len(years)):
                    year = years[year_idx]
                    
                    # Get growth rate without climate change
                    growth_rate = growth_proj[year].values / 100
                    
                    # Project GDP without climate change
                    gdp_cap_no_cc = gdp_cap_no_cc * (1 + growth_rate)
                    
                    # Determine poor status based on current GDP
                    poor_status = gdp_cap_cc <= median_gdp
                    
                    # Calculate baseline growth level
                    baseline_growth = np.zeros(len(base_gdp))
                    baseline_growth[~poor_status] = (temp_coef * base_temp[~poor_status] + 
                                                   temp2_coef * base_temp[~poor_status] ** 2)
                    baseline_growth[poor_status] = (temp_poor_coef * base_temp[poor_status] + 
                                                  temp2_poor_coef * base_temp[poor_status] ** 2)
                    
                    # Calculate new temperature
                    years_since_2010 = year - 2010
                    warming_per_year = additional_warming / total_years
                    
                    new_temp = base_temp + years_since_2010 * warming_per_year * np.array([
                        conversion_factors.get(iso, 1.0) for iso in countries
                    ])
                    
                    # Constrain temperature to 30°C maximum
                    new_temp = np.minimum(new_temp, MAX_TEMPERATURE)
                    
                    # Calculate growth with climate change
                    climate_growth = np.zeros(len(base_gdp))
                    climate_growth[~poor_status] = (temp_coef * new_temp[~poor_status] + 
                                                  temp2_coef * new_temp[~poor_status] ** 2)
                    climate_growth[poor_status] = (temp_poor_coef * new_temp[poor_status] + 
                                                 temp2_poor_coef * new_temp[poor_status] ** 2)
                    
                    # Calculate climate impact
                    climate_impact = climate_growth - baseline_growth
                    
                    # Project GDP with climate change
                    gdp_cap_cc = gdp_cap_cc * (1 + growth_rate + climate_impact)
                
                # Calculate global weighted average GDP in 2099
                pop_weights = pop_proj[2099].values
                total_pop = pop_weights.sum()
                
                if total_pop > 0:
                    damage_results[temp_idx, scenario_idx, 0] = np.average(gdp_cap_cc, weights=pop_weights)
                    damage_results[temp_idx, scenario_idx, 1] = np.average(gdp_cap_no_cc, weights=pop_weights)
        
        # Save results
        self._save_damage_function(damage_results, temp_increases, scenario_names, 'richpoor')
        
        logger.info("Damage function for rich/poor model completed")
        return damage_results
    
    def _save_damage_function(self, damage_results, temp_increases, scenario_names, model_name):
        """Save damage function results."""
        logger.info(f"Saving damage function for {model_name}...")
        
        # Create output directory
        output_dir = OUTPUT_PATH / "projectionOutput"
        output_dir.mkdir(parents=True, exist_ok=True)
        
        # Create results dictionary
        results = {
            'damage_results': damage_results,
            'temp_increases': temp_increases,
            'scenario_names': scenario_names,
            'model_name': model_name
        }
        
        # Save as pickle file
        with open(output_dir / f"DamageFunction_{model_name}.pkl", 'wb') as f:
            pickle.dump(results, f)
        
        logger.info(f"Damage function saved for {model_name}")
    
    def calculate_damage_percentages(self):
        """Calculate damage as percentage of GDP."""
        logger.info("Calculating damage percentages...")
        
        # Load damage functions
        output_dir = OUTPUT_PATH / "projectionOutput"
        
        damage_percentages = {}
        
        for model_name in ['pooled', 'richpoor']:
            damage_file = output_dir / f"DamageFunction_{model_name}.pkl"
            
            if damage_file.exists():
                with open(damage_file, 'rb') as f:
                    damage_data = pickle.load(f)
                
                damage_results = damage_data['damage_results']
                temp_increases = damage_data['temp_increases']
                
                # Calculate damage percentages
                # Damage = (GDP_no_CC - GDP_with_CC) / GDP_no_CC * 100
                damage_pct = np.zeros_like(damage_results[:, :, 0])
                
                for i in range(damage_results.shape[0]):
                    for j in range(damage_results.shape[1]):
                        gdp_no_cc = damage_results[i, j, 1]
                        gdp_with_cc = damage_results[i, j, 0]
                        
                        if gdp_no_cc > 0:
                            damage_pct[i, j] = (gdp_no_cc - gdp_with_cc) / gdp_no_cc * 100
                
                damage_percentages[model_name] = {
                    'damage_pct': damage_pct,
                    'temp_increases': temp_increases
                }
                
                logger.info(f"Damage percentages calculated for {model_name}")
        
        return damage_percentages
    
    def run_all_damage_calculations(self):
        """Run all damage function calculations."""
        logger.info("Running all damage function calculations...")
        
        # Load data
        self.load_data()
        
        # Calculate damage functions
        self.calculate_damage_function_pooled()
        self.calculate_damage_function_rich_poor()
        
        # Calculate damage percentages
        damage_percentages = self.calculate_damage_percentages()
        
        logger.info("All damage function calculations completed")
        return damage_percentages

def run_step5():
    """Run Step 5: Damage Function."""
    logger.info("Starting Step 5: Damage Function")
    
    # Initialize
    processor = DamageFunction()
    
    # Run all damage calculations
    damage_percentages = processor.run_all_damage_calculations()
    
    logger.info("Step 5 completed successfully")
    return damage_percentages

if __name__ == "__main__":
    run_step5()