R/calcDemandProjections.R

In pik-piam/moinput: Model Operations Input Data Library

#' @importFrom magclass setNames

calcDemandProjections<-function(
  scenario_name = "SSP2",
  pop_scen="pop_SSP2",
  gdp_scen="GDP_SSP2", 
  dem_regr_type="log_log_decl_time",
  ls_regr_type="u_shape",  
  calib_year_start="y2005",
  calib_year_end="y2100",
  calib_type="convergence",
  
  dem_aim=NULL, 
  dem_aim_type= "s",
  dem_aim_direction= "down",
  dem_aim_startyear="y2010",
  dem_aim_endyear="y2100",
  
  ls_aim=0.15, 
  ls_aim_type= "s",
  ls_aim_direction= "up",
  ls_aim_startyear="y2010",
  ls_aim_endyear="y2100",
  demand_input_country
) {
  

  #helpfunctions  
  
  demand_kcal_regression<-function(gdp_pc,type="") {
    if (!is.magpie(gdp_pc)) {stop("gdp_pc has to be a magpie object")}
    if (length(dim(unwrap(gdp_pc)))!=3) {stop("gdp_pc has to be a 3d object")}
    if (type=="log_log_lin_time") {
      regr_c<- 2.8251159
      regr_d<- 2.1313756*10^(-3)
      regr_e<- 0.1622219
      regr_f<- -3.12439*10^(-05)
      years<-new.magpie(cells_and_regions="GLO", years=getYears(gdp_pc),names="NULL",fill=getYears(gdp_pc,as.integer=TRUE))
      regr_a  <- exp(regr_c + regr_d * years)
      regr_b  <- regr_e + regr_f * years
      dummy<-(gdp_pc*0+1)*regr_b
      out<-regr_a * gdp_pc ^ (regr_b)
      #  exp(2.8251159+2.1313756*10^(-3)*2100) * 60000 ^ (0.1622219 + 3.12439*10^(-05) * 2100)
      
    } else if (type=="log_log_decl_time") {
      
      regr_am1 <- 387.473708
      regr_am2 <- 9.774726
      regr_am3 <- 933.888522
      regr_bm1 <- 0.008445119
      regr_bm2 <- -0.755692561
      regr_bm3 <- 0.08940805
      
      years<-new.magpie(cells_and_regions="GLO", years=getYears(gdp_pc),names="NULL",fill=getYears(gdp_pc,as.integer=TRUE))-1960
      regr_a  <- regr_am3 + (regr_am1 * years)/(years+regr_am2)
      regr_b  <- (regr_bm3 + (regr_bm1 * years)/(years+regr_bm2))
      out <- regr_a * gdp_pc ^ regr_b
      #   exp(2.8251159+2.1313756*10^(-3)*2100) * 60000 ^ (0.1622219 + 3.12439*10^(-05) * 2100)
      
    } else if (type=="log_log_no_time") {
      
      regr_a<- exp(7.074079)
      regr_b<- 0.099321
      out <-  regr_a * gdp_pc ^ regr_b
    } else {stop("unknown regression type")}
    return(out)
  }
  
  demand_ls_regression<-function(gdp_pc,type="u_shape") {
    length(dim(unwrap(gdp_pc)))
    if (type=="u_shape") {
      regr_p1 <- 1.371507*10^(-2)
      regr_p2 <- -5.295249*10^(-6)
      regr_q1 <- -1.102410*10^(-4)
      regr_q2 <- 6.403996*10^(-8)
      years<-new.magpie(cells_and_regions="GLO", years=getYears(gdp_pc),names="NULL",fill=getYears(gdp_pc,as.integer=TRUE))
      out <- (regr_p1 + regr_p2 * years)*((gdp_pc)^(0.5))*exp(-(regr_q1 + regr_q2 * years) * gdp_pc)
      
    } else if (type=="mult_lin_regr") {
      regr_a <- -36.732779363
      regr_b <- 4.497483702
      regr_c <- 0.016039027
      regr_d <- -0.002077227
      years<-new.magpie(cells_and_regions="GLO", years=getYears(gdp_pc),names="NULL",fill=getYears(gdp_pc,as.integer=TRUE))
      out <-  exp(regr_a + regr_b * log(gdp_pc) + regr_c*years + regr_d * log(gdp_pc)*years)
    } else {stop("unknown regression type")}
    return(out)
  }  
  
  
  ### readin
  gdp<-setNames(demand_input_country[,,gdp_scen],NULL)
  pop<-setNames(demand_input_country[,,pop_scen],NULL)
  gdp_pc <- gdp/pop
  kcal_pc_calib_to = setNames(demand_input_country[,,"kcal_fao"],NULL) / setNames(demand_input_country[,,"pop_hist"],NULL)
  ls_calib_to = setNames(demand_input_country[,,"livst_kcal_fao"],NULL) / setNames(demand_input_country[,,"kcal_fao"],NULL)
  #fish_pc_calib_to = setNames(demand_input_country[,,"fish_kcal_fao"],NULL)  / setNames(demand_input_country[,,"pop_hist"],NULL)
  material_pc_calib_to = setNames(demand_input_country[,,"material"],NULL)  / setNames(demand_input_country[,,"pop_hist"],NULL)
  
  ###checks
  
  ### dem projections
  # Apply regression on country data <- function
  
  kcal_pc_regr <- demand_kcal_regression(gdp_pc,type=dem_regr_type)
  ls_regr   <- demand_ls_regression(gdp_pc,type=ls_regr_type)
  
  # Calibrate country data
  
  
  kcal_pc_calib <- calibrate_it(
    origin=kcal_pc_regr, 
    cal_to=setNames(kcal_pc_calib_to,NULL), 
    cal_type=calib_type,
    cal_year=calib_year_start, 
    end_year=calib_year_end, 
    report_calibration_factors=FALSE)
  
  ls_calib <- calibrate_it(
    origin=ls_regr, 
    cal_to=setNames(ls_calib_to,NULL), 
    cal_type=calib_type,
    cal_year=calib_year_start, 
    end_year=calib_year_end, 
    report_calibration_factors=FALSE)
  
  
  # dem_aim_convergence
  
  if(!is.null(dem_aim)) {
    
    kcal_pc <- convergence(
      origin=      kcal_pc_calib,
      aim=         dem_aim,
      type=        dem_aim_type,
      start_year=  dem_aim_startyear,
      end_year=    dem_aim_endyear,   
      direction=   dem_aim_direction
    ) 
    
  } else {
    
    kcal_pc <- kcal_pc_calib
    
  } 
  
  if(!is.null(ls_aim)) {
    
    ls <- convergence(
      origin=      ls_calib,
      aim=         ls_aim,
      type=        ls_aim_type,
      start_year=  ls_aim_startyear,
      end_year=    ls_aim_endyear,   
      direction=   ls_aim_direction
    ) 
    
  } else {
    
    ls <- ls_calib
    
  } 
  
  ### Undernourishment
  
  
  hunger_shr <- 2674.855 * 0.997916997^kcal_pc / 100 
  hunger_shr[hunger_shr>1]<-1
  hunger <- hunger_shr * pop
  
  ### calculate total demand
  
  dem <- kcal_pc * pop
  l   <- kcal_pc * pop * ls  
  
  # material demand projections
  
  mat <- material_pc_calib_to * dem/setYears(dem[,"y1995",],NULL)
  
  # Fish is assumed to stay constant (Fish is part of livestock products, but not represented by MAGPIE)
  
  #fish<-dem
  #fish[,,] <- NA
  #fish[,,] <- setYears((fish_calib_to * pop)[,calib_year_start,],NULL)
  
  ### Add History before calibration point 
  
  if (scenario_name=="history"){
    years<-1:length(getYears(demand_input_country))
  } else {
    years<-1:(which(getYears(demand_input_country)==calib_year_start)-1)
  }
  
  gdp[,years,]<-setNames(demand_input_country[,years,"gdp_hist"],NULL)
  pop[,years,]<-setNames(demand_input_country[,years,"pop_hist"],NULL)
  kcal_pc[,years,]<-setNames(demand_input_country[,years,"kcal_fao"],NULL) / setNames(demand_input_country[,years,"pop_hist"],NULL) 
  #fish[,years,]<-(fish_calib_to * pop )[,years,]
  dem[,years,]<-kcal_pc[,years,]*pop[,years,]    
  l[,years,]<-setNames(demand_input_country[,years,"livst_kcal_fao"],NULL)
  ls[,years,]<-l[,years,]/dem[,years,]
  gdp_pc[,years,]<-gdp[,years,]/pop[,years,]  
  waste_shr=((kcal_pc>(2200/0.85))*(kcal_pc-2200)+(kcal_pc<=(2200/0.85))*(kcal_pc*0.15))/kcal_pc
  waste=((kcal_pc>(2200/0.85))*(kcal_pc-2200)+(kcal_pc<=(2200/0.85))*(kcal_pc*0.15))*pop
  
  # waste shr: at least 15% and all above 2200 kcal
  
  ### Aggregate to region  
  
  ### change units
  # kcal per day --> PJ per year
  #      1       --> 4.184 * 365 / 1000000
  dem = dem * 4.184*365/1000000
  l   = l * 4.184*365/1000000
  #fish = fish * 4.184*365/1000000
  waste = waste * 4.184*365/1000000  
  
  ### prepare output
  
  combine_it<-function(add_x,to_y,varname,scenname) {
    dimnames(add_x)[[3]]<-paste(varname,scenname,sep=".")
    out<-mbind(add_x,to_y)
    return(out)
  }
  
  country<-NULL
  
  country<-combine_it(dem,country,"dem",scenario_name)
  country<-combine_it(ls,country,"ls",scenario_name)
  country<-combine_it(pop,country,"pop",scenario_name)
  country<-combine_it(gdp,country,"gdp",scenario_name)
  country<-combine_it(gdp_pc,country,"gdp_pc",scenario_name)
  country<-combine_it(kcal_pc,country,"kcal_pc",scenario_name)
  country<-combine_it(l,country,"l",scenario_name)
  country<-combine_it(mat,country,"mat",scenario_name)
  #  country<-combine_it(fish,country,"fish",scenario_name)  
  country<-combine_it(waste,country,"waste",scenario_name)    
  country<-combine_it(waste_shr,country,"waste_shr",scenario_name)  
  country<-combine_it(hunger,country,"hunger",scenario_name) 
  country<-combine_it(hunger_shr,country,"hunger_shr",scenario_name) 
  
  country_weight<-country
  country_weight[,,c("dem","pop","gdp","l","mat","waste","hunger")]<-NA
  country_weight[,,c("gdp_pc","kcal_pc","hunger_shr")]<-setNames(country_weight[,,c("pop")],NULL)
  country_weight[,,c("ls","waste_shr")]  <-setNames(country_weight[,,c("dem")],NULL)
  
  return(list(x=country,weight=country_weight))
}