R/reportPrices.R
In pik-piam/remind: The REMIND R Package
Defines functions
reportPrices
Documented in
reportPrices
#' Read in GDX and calculate prices, used in convGDX2MIF.R for the reporting
#'
#' Read in price information from GDX file, information used in convGDX2MIF.R
#' for the reporting
#'
#'
#' @param gdx a GDX object as created by readGDX, or the path to a gdx
#' @param output a magpie object containing all needed variables generated by other report*.R functions
#' @param regionSubsetList a list containing regions to create report variables region
#' aggregations. If NULL (default value) only the global region aggregation "GLO" will
#' be created.
#' @return MAgPIE object - contains the price variables
#' @importFrom luscale speed_aggregate
#' @author David Klein
#' @seealso \code{\link{convGDX2MIF}}
#' @examples
#'
#' \dontrun{reportPrices(gdx)}
#'
#' @importFrom dplyr %>% case_when distinct filter_ inner_join tibble
#' @importFrom gdx readGDX
#' @importFrom luscale speed_aggregate
#' @importFrom magclass mbind getYears getRegions setNames dimSums new.magpie lowpass complete_magpie
#' @importFrom quitte df.2.named.vector getColValues
#' @importFrom readr read_csv
#' @export
reportPrices <- function(gdx,output=NULL,regionSubsetList=NULL) {
if(is.null(output)){
output <- reportPE(gdx,regionSubsetList)
output <- mbind(output,reportSE(gdx,regionSubsetList))
output <- mbind(output,reportFE(gdx,regionSubsetList))
output <- mbind(output,reportEmi(gdx,regionSubsetList))
output <- mbind(output,reportExtraction(gdx,regionSubsetList))
output <- mbind(output,reportMacroEconomy(gdx,regionSubsetList)[,getYears(output),])
}
#---- Functions
find_real_module <- function(module_set, module_name){
return(module_set[module_set$modules == module_name,2])
}
compute_agg_price_fe = function(mobj_prices,mobj_output,sector){
if (!sector %in% c("Buildings","Industry","Transport")) stop("argument 'sector' must be in c('Buildings','Industry','Transport')")
names_prices = grep(paste0("^Price\\|Final Energy\\|[A-z ]+\\|",sector," \\(US\\$2005\\/GJ\\)"),getNames(mobj_prices),value = T)
names_fe = gsub(paste0("^Price\\|Final Energy\\|([A-z ]+)\\|",sector," \\(US\\$2005\\/GJ\\)"),
paste0("FE|",sector,"\\|\\1 (EJ/yr)"), names_prices)
names_fe = gsub("Heating Oil","Liquids",names_fe)
names(names_fe) = names_prices
# price_total = sum(i,price_i*quant_i)/quant_total
mobj_prices = mbind(mobj_prices,
setNames(
dimSums(
setNames(
mobj_prices[,,names_prices],names_fe[names_prices])
* mobj_output[regs,time,names_fe]
,dim = 3)
/collapseNames(mobj_output[regs,time,paste0("FE|",sector," (EJ/yr)")]),
paste0("Price|Final Energy|",sector," (US$2005/GJ)")
))
# Moving Averages
mobj_prices = mbind(mobj_prices,
setNames( lowpass(
dimSums(
setNames(
mobj_prices[,,names_prices],names_fe[names_prices])
* mobj_output[regs,time,names_fe]
,dim = 3)
/collapseNames(mobj_output[regs,time,paste0("FE|",sector," (EJ/yr)")]), fix="both", altFilter=match(2010,time)),
paste0("Price|Final Energy|",sector,"|Moving Avg (US$2005/GJ)")
))
return(mobj_prices)
}
#---- End of Functions
####### get realisations #########
realisation <- readGDX(gdx, "module2realisation")
####### get power realisations #########
realisation <- readGDX(gdx, "module2realisation")
power_realisation <- if ("power" %in% realisation[, 1]) realisation[which(realisation[,1] == "power"),2]
####### conversion factors ##########
s_GWP_CH4 <- readGDX(gdx,c("sm_gwpCH4","s_gwpCH4","s_GWP_CH4"),format="first_found", react = "silent")
s_GWP_N2O <- readGDX(gdx,c("s_gwpN2O","s_GWP_N2O"),format="first_found", react = "silent")
TWa_2_EJ <- 31.536
tdptwyr2dpgj <- 31.71 #TerraDollar per TWyear to Dollar per GJ
p80_subset <- c("perm","good","peur","peoil","pegas","pecoal","pebiolc") #TODO: read in from gdx as sets trade
pebal_subset <- c("pebiolc","peoil","pegas","pecoal","peur")
####### read in needed data #########
## sets
fe_transp_fety35 <- readGDX(gdx,name="FE_Transp_fety35",types="sets",format="first_found", react = "silent")
se2fe <- readGDX(gdx,name="se2fe",types="sets",format="first_found", react = "silent")
sety <- readGDX(gdx,c("entySe","sety"),types="sets",format="first_found", react = "silent")
fety <- readGDX(gdx,name=c("entyFe","fety"),types="sets",format="first_found", react = "silent")
finenbal <- readGDX(gdx,name=c("fe2ppfEn","finenbal"),types="sets",format="first_found", react = "silent")
all_esty <- readGDX(gdx,name=c("all_esty"),types="sets",format="first_found", react = "silent")
if(!is.null(all_esty)){
fe2es <- readGDX(gdx,name=c("fe2es"),types="sets",format="first_found", react = "silent")
fe2es <- unique(fe2es[c("all_enty","all_esty")])
colnames(fe2es) <- gsub("all_esty","all_in", colnames(fe2es))
} else {
fe2es = NULL
}
ppfenfromes <- readGDX(gdx,name=c("ppfEnFromEs","ppfenfromes"),types="sets",format="first_found", react = "silent")
finenbalfehos <- finenbal %>% filter_(~all_enty == "fehos")
ppfKap <- readGDX(gdx,"ppfKap")
ppfen_stat <- readGDX(gdx,c("ppfen_stationary_dyn38","ppfen_stationary_dyn28","ppfen_stationary"),format="first_found", react = "silent")
if (length(ppfen_stat) == 0) ppfen_stat = NULL
ppfen_build <- readGDX(gdx,c("ppfen_buildings_dyn36","ppfen_buildings_dyn28","ppfen_buildings"),format="first_found", react = "silent")
ue_dyn36 <- readGDX(gdx,c("ue_dyn36"),format="first_found", react = "silent")
ppfen_build <- setdiff(ppfen_build, ue_dyn36)
ppfen_ind <- readGDX(gdx,c("ppfen_industry_dyn37","ppfen_industry_dyn28","ppfen_industry"),format="first_found", react = "silent")
ppfen_stat_build_ind <- c(ppfen_stat,ppfen_build,ppfen_ind)
edge_scen <- readGDX(gdx,"EDGE_scenario",format="first_found", react = "silent")
# Realisation of the different modules
if ( !is.null(realisation) & (!"CES_structure" %in% realisation[, 1] )){
stat_mod = find_real_module(realisation,"stationary")
tran_mod = find_real_module(realisation,"transport")
indu_mod = find_real_module(realisation,"industry")
buil_mod = find_real_module(realisation,"buildings")
} else {
if ( !is.null(ppfen_stat)){ # In case the set realisation did not exist, find out whether it was stationary or buildings-industry
stat_mod = "simple"
indu_mod = "off"
buil_mod = "off"
} else {
stat_mod = "off"
indu_mod = "fixed_shares"
buil_mod = "simple"
}
if ( !is.null(edge_scen)){ # In case the set realization did not exist, find out whether "edge_esm" or "complex" realization was chosen
tran_mod = "edge_esm"
} else {
tran_mod = "complex"
}
}
if (indu_mod %in% c('fixed_shares', 'subsectors')) {
fe2ppfen37 <- readGDX(gdx, 'fe2ppfen') %>%
filter(.data$all_in %in% ppfen_ind)
}
## parameter
shift_p <- readGDX(gdx,name="p30_pebiolc_pricshift",format="first_found")
mult_p <- readGDX(gdx,name="p30_pebiolc_pricmult",format="first_found")
pric_mag <- readGDX(gdx,name="p30_pebiolc_pricemag",format="first_found")
pric_emu_pre <- readGDX(gdx,name="p30_pebiolc_price_emu_preloop",format="first_found")
pric_emu_pre_shifted <- readGDX(gdx,name="p30_pebiolc_price_emu_preloop_shifted",format="first_found")
bio_tax_factor <- readGDX(gdx,name="p21_tau_bioenergy_tax",format="first_found")
pm_pvp <- readGDX(gdx,name=c("pm_pvp","p80_pvp"),format="first_found")[,,p80_subset]
pm_ts <- readGDX(gdx,name='pm_ts',format="first_found")
pm_dataemi <- readGDX(gdx,name=c("pm_emifac","pm_dataemi"),format="first_found",restore_zeros=FALSE)[,,c("pegas.seel.ngt.co2","pecoal.seel.pc.co2")]
pm_pvpRegi <- readGDX(gdx,name='pm_pvpRegi',format="first_found")[,,"perm"]
pm_taxCO2eq <- readGDX(gdx,name=c("pm_taxCO2eq","pm_tau_CO2_tax"),format="first_found")
pm_taxCO2eqPower <- readGDX(gdx,name=c("pm_taxCO2eqPower"),format="first_found")
pm_taxCO2eqSCC <- readGDX(gdx,name='pm_taxCO2eqSCC',format="first_found")
pm_delta_kap <- readGDX(gdx, c("pm_delta_kap", "p_delta_kap"), format = "first_found")
if(is.null(getNames(pm_delta_kap))) pm_delta_kap = setNames(pm_delta_kap,"kap")
## variables
pric_emu <- readGDX(gdx,name="vm_pebiolc_price",field="l",format="first_found")
prodFE <- readGDX(gdx,name='vm_prodFe',field="l",format="first_found",restore_zeros = FALSE)
prodSE <- readGDX(gdx,name='vm_prodSe',field="l",format="first_found",restore_zeros = FALSE)
demSe <- readGDX(gdx,name=c("vm_demSe","v_demSe"),types="variables",field="l",format="first_found",restore_zeros=FALSE)
cesIO <- readGDX(gdx,name='vm_cesIO',types="variables",field="l",format="first_found",restore_zeros=FALSE)
fe_taxCES <- readGDX(gdx, name=c('p21_tau_fe_tax','pm_tau_fe_tax'), format="first_found", react = F)[ppfen_stat_build_ind]
fe_subCES <- readGDX(gdx, name=c('p21_tau_fe_sub','pm_tau_fe_sub'), format="first_found", react = F)[ppfen_stat_build_ind]
if(is.null(fe_taxCES) & is.null(fe_subCES)){
fe_taxCES = readGDX(gdx, name=c('p21_tau_fe_tax_bit_st'))[,,ppfen_stat_build_ind]
fe_subCES = readGDX(gdx, name=c('p21_tau_fe_sub_bit_st'))[,,ppfen_stat_build_ind]
}
getSets(fe_taxCES) = gsub("all_enty","all_in", getSets(fe_taxCES))
getSets(fe_subCES) = gsub("all_enty","all_in", getSets(fe_subCES))
if (!is.null(all_esty)){
fe_taxES = readGDX(gdx, name=c("pm_tau_fe_tax_ES_st",'p21_tau_fe_tax_ES_st'),format = "first_found", react = "silent")
fe_subES = readGDX(gdx, name=c('pm_tau_fe_sub_ES_st','p21_tau_fe_sub_ES_st'),format = "first_found", react = "silent")
getSets(fe_taxES) = gsub("all_esty","all_in", getSets(fe_taxES))
getSets(fe_subES) = gsub("all_esty","all_in", getSets(fe_subES))
} else {
fe_taxES = NULL
fe_subES = NULL
}
## equations
pebal.m <- readGDX(gdx,name=c("q_balPe","qm_pebal"),types = "equations",field = "m",format = "first_found")[,,pebal_subset]
budget.m <- readGDX(gdx,name='qm_budget',types = "equations",field = "m",format = "first_found") # Alternative: calcPrice
sebal.m <- readGDX(gdx,name=c("q_balSe","q_sebal"),types="equations",field="m",format="first_found")
if (!is.null(power_realisation)) {
sebalSeel.m <- readGDX(gdx,name="q32_balSe",types="equations",field="m",format="first_found")
}
esm2macro.m <- readGDX(gdx,name='q_esm2macro',types="equations",field="m",format="first_found")
febal.m <- readGDX(gdx,name=c("q_balFe","q_febal"),types="equations",field="m",format="first_found")[,,fety]
balfinen.m <- readGDX(gdx,name=c("qm_balFeForCesAndEs","qm_balFeForCes","q_balFeForCes","q_balfinen"),types="equations",field="m",format="first_found", react = "silent")
cm_emiscen <- readGDX(gdx,name='cm_emiscen',format="first_found")
q37_limit_secondary_steel_share.m <- readGDX(
gdx, name = 'q37_limit_secondary_steel_share', types = 'equation',
field = 'm', react = 'silent')
#####################################
#choose the CES entries names for transport
name_trsp=c("fepet","ueLDVt","fedie","ueHDVt","feelt","ueelTt")
name_trsp=name_trsp[name_trsp%in%getNames(esm2macro.m)]
#####################################
####### calculate minimal temporal resolution #####
y <- Reduce(intersect,list(getYears(febal.m),getYears(prodFE),getYears(sebal.m),getYears(prodSE)))
febal.m <- febal.m[,y,]
balfinen.m <- balfinen.m[,y,]
budget.m <- budget.m[,y,]
sebal.m <- sebal.m[,y,]
if (!is.null(power_realisation)) {
sebalSeel.m <- sebalSeel.m[,y,]
}
pm_pvp <- pm_pvp[,y,]
pebal.m <- pebal.m[,y,]
shift_p <- shift_p[,y,]
mult_p <- mult_p[,y,]
pric_mag <- pric_mag[,y,]
pric_emu_pre <- pric_emu_pre[,y,]
pric_emu_pre_shifted <- pric_emu_pre_shifted[,y,]
pric_emu <- pric_emu[,y,]
bio_tax_factor <- bio_tax_factor[,y,]
esm2macro.m <- esm2macro.m[,y,]
pm_dataemi <- pm_dataemi[,y,]
pm_ts <- pm_ts[,y,]
pm_pvpRegi <- pm_pvpRegi[,y,]
pm_taxCO2eq <- pm_taxCO2eq[,y,]
if (!is.null(pm_taxCO2eqPower)) {
pm_taxCO2eqPower <- pm_taxCO2eqPower[,y,]
}
if (!is.null(pm_taxCO2eqSCC)) {
pm_taxCO2eqSCC <- pm_taxCO2eqSCC[,y,]
}
cesIO <- cesIO[,y,]
####### pm_pvp = EPS results in fantasy carbon prices
for(t in getYears(pm_pvp)){
if(pm_pvp[,t,"good"] == 0){
pm_pvp[,t,"good"] = 0.0001;
}
}
####### calculate reporting parameters ############
tmp <- NULL
# Calculate and append new variables to magpie object
# Costs and prices for purpose-grown 2nd gen. bioenergy
# - Shiftfactor
# - imported values from MAgPIE reporting
# - results from emulator based on MAgPIE demand (before main solve)
# - results from emulator based on MAgPIE demand multiplied with shiftfactor (before main solve)
# - results from emulator based on REMIND demand (after main solve)
# - results from emulator based on REMIND demand multiplied with shiftfactor (after main solve)
time <- getYears(budget.m,as.integer=TRUE)
regs <- getRegions(budget.m)
# report variables
tmp <- NULL
tmp <- mbind(tmp,setNames(pebal.m[,,"pebiolc"] / (budget.m+1e-10) * tdptwyr2dpgj, "Price|Biomass|Primary Level (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(lowpass(pebal.m[,,"pebiolc"] / (budget.m+1e-10), fix="both", altFilter=match(2010,time)) * tdptwyr2dpgj, "Price|Biomass|Primary Level|Moving Avg (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(shift_p, "Price|Biomass|Shiftfactor ()"))
tmp <- mbind(tmp,setNames(mult_p, "Price|Biomass|Multfactor ()"))
tmp <- mbind(tmp,setNames(pric_mag * tdptwyr2dpgj, "Price|Biomass|MAgPIE (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(pric_emu_pre * tdptwyr2dpgj, "Price|Biomass|Emulator presolve (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(pric_emu_pre_shifted * tdptwyr2dpgj, "Price|Biomass|Emulator presolve shifted (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(pric_emu * tdptwyr2dpgj, "Price|Biomass|Emulator shifted (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(pric_emu * bio_tax_factor * tdptwyr2dpgj, "Price|Biomass|Bioenergy tax (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(pebal.m[,,"peoil"]/(budget.m+1e-10) * tdptwyr2dpgj, "Price|Crude Oil|Primary Level (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(pebal.m[,,"pegas"]/(budget.m+1e-10) * tdptwyr2dpgj, "Price|Natural Gas|Primary Level (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(lowpass(pebal.m[,,"pegas"]/(budget.m+1e-10), fix="both", altFilter=match(2010,time)) * tdptwyr2dpgj, "Price|Natural Gas|Primary Level|Moving Avg (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(pebal.m[,,"pecoal"]/(budget.m+1e-10) * tdptwyr2dpgj, "Price|Coal|Primary Level (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(pebal.m[,,"peur"]/(budget.m+1e-10) * tdptwyr2dpgj, "Price|Uranium|Primary Level (US$2005/GJ)"))
# seondary energy prices
if (!is.null(power_realisation)) {
tmp <- mbind(tmp,setNames(sebalSeel.m[,,"seel"]/(budget.m+1e-10) * tdptwyr2dpgj, "Price|Secondary Energy|Electricity (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(lowpass(sebalSeel.m[,,"seel"]/(budget.m+1e-10), fix="both", altFilter=match(2010,time)) * tdptwyr2dpgj, "Price|Secondary Energy|Electricity|Moving Avg (US$2005/GJ)"))
} else {
tmp <- mbind(tmp,setNames(sebal.m[,,"seel"]/(budget.m+1e-10) * tdptwyr2dpgj, "Price|Secondary Energy|Electricity (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(lowpass(sebal.m[,,"seel"]/(budget.m+1e-10), fix="both", altFilter=match(2010,time)) * tdptwyr2dpgj, "Price|Secondary Energy|Electricity|Moving Avg (US$2005/GJ)"))
}
# for biomass price for secondary level is the same as primary level
tmp <- mbind(tmp,setNames(pebal.m[,,"pebiolc"]/(budget.m+1e-10) * tdptwyr2dpgj, "Price|Secondary Energy|Biomass (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(sebal.m[,,"seh2"]/(budget.m+1e-10) * tdptwyr2dpgj , "Price|Secondary Energy|Hydrogen (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(sebal.m[,,"sehe"]/(budget.m+1e-10) * tdptwyr2dpgj , "Price|Secondary Energy|Heat (US$2005/GJ)"))
# energy services
tmp <- mbind(tmp,setNames(abs(esm2macro.m[,,name_trsp[2]]/(budget.m+1e-10)) * tdptwyr2dpgj , "Price|Energy Service|Transport nonLDV (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(abs(esm2macro.m[,,name_trsp[1]]/(budget.m+1e-10)) * tdptwyr2dpgj , "Price|Energy Service|Transport LDV (US$2005/GJ)"))
#Final energy prices
#Transport (Taxes already included)
if (tran_mod == "complex"){
#Translate the marginal utility of the constraint into the marginal income (price).
#For transport liquids use average between diesel and petrol.
tmp <- mbind(tmp,setNames(abs(febal.m[,,"feelt"]/(budget.m+1e-10)) * tdptwyr2dpgj, "Price|Final Energy|Electricity|Transport (US$2005/GJ)"))
if("fegat" %in% getNames(febal.m)){
tmp <- mbind(tmp,setNames(abs(febal.m[,,"fegat"]/(budget.m+1e-10)) * tdptwyr2dpgj, "Price|Final Energy|Gases|Transport (US$2005/GJ)"))
}
tmp <- mbind(tmp,setNames(abs(lowpass(febal.m[,,"feelt"]/(budget.m+1e-10), fix="both", altFilter=match(2010,time))) * tdptwyr2dpgj, "Price|Final Energy|Electricity|Transport|Moving Avg (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(abs(febal.m[,,"feh2t"]/(budget.m+1e-10)) * tdptwyr2dpgj, "Price|Final Energy|Hydrogen|Transport (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(abs(febal.m[,,"fedie"]/(budget.m+1e-10)) * tdptwyr2dpgj, "Price|Final Energy|Liquids|Transport (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(abs(lowpass((febal.m[,,"fedie"]+febal.m[,,"fepet"])/(2*budget.m+1e-10), fix="both", altFilter=match(2010,time))) * tdptwyr2dpgj, "Price|Final Energy|Liquids|Transport|Moving Avg (US$2005/GJ)"))
} else if(tran_mod == "edge_esm"){
#Translate the marginal utility of the constraint into the marginal income (price)
#For transport liquids use average between diesel and petrol.
tmp <- mbind(tmp,setNames(abs(balfinen.m[,,"feelt"]/(budget.m+1e-10)) * tdptwyr2dpgj, "Price|Final Energy|Electricity|Transport (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(abs(balfinen.m[,,"fegat"]/(budget.m+1e-10)) * tdptwyr2dpgj, "Price|Final Energy|Gases|Transport (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(abs(lowpass(balfinen.m[,,"feelt"]/(budget.m+1e-10), fix="both", altFilter=match(2010,time))) * tdptwyr2dpgj, "Price|Final Energy|Electricity|Transport|Moving Avg (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(abs(balfinen.m[,,"feh2t"]/(budget.m+1e-10)) * tdptwyr2dpgj, "Price|Final Energy|Hydrogen|Transport (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(abs((balfinen.m[,,"fedie"]+balfinen.m[,,"fepet"])/(2*budget.m+1e-10)) * tdptwyr2dpgj, "Price|Final Energy|Liquids|Transport (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(abs(lowpass((balfinen.m[,,"fedie"]+balfinen.m[,,"fepet"])/(2*budget.m+1e-10), fix="both", altFilter=match(2010,time))) * tdptwyr2dpgj, "Price|Final Energy|Liquids|Transport|Moving Avg (US$2005/GJ)"))
}
tmp = compute_agg_price_fe(tmp,output,"Transport")
#Buildings and Industry (Taxes need to be added)
a = abs(mselect(balfinen.m[,y,][rbind(finenbal,fe2es)]))/abs((budget.m+1e-10)) # Translate the marginal utility of the constraint into the marginal income (price)
b = complete_magpie(a) # Due to a strange behaviour of magclass objects addition, we need to use complete_magpie to make the addition
prices_fe_bi = (b + mbind(fe_taxCES[,y,getColValues(finenbal,"all_in")],
fe_taxES[,y,getColValues(fe2es,"all_in")])
+ mbind(fe_subCES[,y,getColValues(finenbal,"all_in")],
fe_subES[,y,getColValues(fe2es,"all_in")])) * tdptwyr2dpgj # add the taxes and subsidies for the prices of buildings and industry. For transport, they are in the marginal of febal
prices_fe_bi = prices_fe_bi[,,getNames(a)]
if (stat_mod == "simple" ){
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"feels.feels"], "Price|Final Energy|Electricity|Stationary (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"fegas.fegas"], "Price|Final Energy|Gases (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"fehes.fehes"], "Price|Final Energy|Heat (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"fehos.fehos"], "Price|Final Energy|Heating Oil (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"feh2s.feh2s"], "Price|Final Energy|Hydrogen|Stationary (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"fesos.fesos"], "Price|Final Energy|Solids (US$2005/GJ)"))
}
if (buil_mod == "simple"){
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"feels.feelb"], "Price|Final Energy|Electricity|Buildings (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(lowpass(prices_fe_bi[,,"feels.feelb"], fix="both", altFilter=match(2010,time)) , "Price|Final Energy|Electricity|Buildings|Moving Avg (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"fegas.fegab"], "Price|Final Energy|Gases|Buildings (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"fehes.feheb"], "Price|Final Energy|Heat|Buildings (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"fehos.fehob"], "Price|Final Energy|Heating Oil|Buildings (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"feh2s.feh2b"], "Price|Final Energy|Hydrogen|Buildings (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"fesos.fesob"], "Price|Final Energy|Solids|Buildings (US$2005/GJ)"))
tmp = compute_agg_price_fe(tmp,output,"Buildings")
}
if (buil_mod %in% c("services_putty", "services_with_capital")){
# Prices across energy services should be identical, so we take the price from only one service (cooking and water heating)
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"feels.uecwelb"], "Price|Final Energy|Electricity|Buildings (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(lowpass(prices_fe_bi[,,"feels.uecwelb"], fix="both", altFilter=match(2010,time)) , "Price|Final Energy|Electricity|Buildings|Moving Avg (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"fegas.uecwgab"], "Price|Final Energy|Gases|Buildings (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"fehes.uecwheb"], "Price|Final Energy|Heat|Buildings (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"fehos.uecwhob"], "Price|Final Energy|Heating Oil|Buildings (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"feh2s.uecwh2b"], "Price|Final Energy|Hydrogen|Buildings (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"fesos.uecwsob"], "Price|Final Energy|Solids|Buildings (US$2005/GJ)"))
tmp = compute_agg_price_fe(tmp,output,"Buildings")
}
if (indu_mod %in% c('fixed_shares', 'subsectors')) {
industry_fe_mixer <- inner_join(
fe2ppfen37,
paste(
'all_enty, fe',
'fesos, Solids',
'fehos, Heating Oil',
'fegas, Gases',
'feh2s, Hydrogen',
'fehes, Heat',
'feels, Electricity',
sep = '\n') %>%
read_csv(),
'all_enty'
)
for (.fe in unique(industry_fe_mixer$fe)) {
fe2ppfen <- industry_fe_mixer %>%
filter(.data$fe == .fe) %>%
distinct(.data$all_enty, .keep_all = TRUE) %>%
select(-.data$fe) %>%
paste(collapse = '.')
tmp <- mbind(
tmp,
setNames(
prices_fe_bi[,,fe2ppfen],
paste0('Price|Final Energy|', .fe, '|Industry (US$2005/GJ)')
)
)
tmp <- mbind(
tmp,
setNames(
lowpass(prices_fe_bi[,,fe2ppfen] , fix="both", altFilter=match(2010,time)) ,
paste0('Price|Final Energy|', .fe, '|Industry|Moving Avg (US$2005/GJ)')
)
)
}
rm(industry_fe_mixer, fe2ppfen)
tmp = compute_agg_price_fe(tmp,output,"Industry")
}
if ("seliq" %in% sety) {
tmp <- mbind(tmp,setNames(sebal.m[,,"seliq"]/(budget.m+1e-10)*tdptwyr2dpgj, "Price|SE|seliq (US$2005/GJ)" ))
} else if ("sepet" %in% sety) {
tmp <- mbind(tmp,setNames(sebal.m[,,"sepet"]/(budget.m+1e-10)*tdptwyr2dpgj, "Price|SE|sepet (US$2005/GJ)" ))
tmp <- mbind(tmp,setNames(sebal.m[,,"sedie"]/(budget.m+1e-10)*tdptwyr2dpgj, "Price|SE|sedie (US$2005/GJ)" ))
} else {
tmp <- mbind(tmp,setNames(sebal.m[,,"seliqbio"]/(budget.m+1e-10)*tdptwyr2dpgj, "Price|SE|seliqbio (US$2005/GJ)" ))
tmp <- mbind(tmp,setNames(sebal.m[,,"seliqfos"]/(budget.m+1e-10)*tdptwyr2dpgj, "Price|SE|seliqfos (US$2005/GJ)" ))
}
tmp <- mbind(tmp,setNames((pebal.m[,,"pegas"]/pm_ts - pm_dataemi[,,"pegas.seel.ngt.co2"]*pm_pvpRegi) / (budget.m/pm_ts + 1e-10) * tdptwyr2dpgj, "Price|Natural Gas|Secondary Level (US$2005/GJ)"))
tmp <- mbind(tmp,setNames((pebal.m[,,"pecoal"]/pm_ts - pm_dataemi[,,"pecoal.seel.pc.co2"]*pm_pvpRegi) / (budget.m/pm_ts + 1e-10) * tdptwyr2dpgj, "Price|Coal|Secondary Level (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(abs(pm_pvpRegi / (pm_pvp[,,"good"] + 1e-10)) * 1000 * 12/44, "Price|Carbon (US$2005/t CO2)"))
tmp <- mbind(tmp,setNames(abs(pm_taxCO2eq) * 1000 * 12/44, "Price|Carbon|Guardrail (US$2005/t CO2)"))
if (is.null(regionSubsetList$EUR)) {
tmp <- mbind(tmp,setNames(pm_taxCO2eq * 1000 * 12/44, "Price|Carbon|EU-wide Regulation For All Sectors (US$2005/t CO2)"))
} else {
co2EUprice <- pm_taxCO2eq * 1000 * 12/44
co2EUprice[getRegions(pm_taxCO2eq)[!getRegions(pm_taxCO2eq) %in% regionSubsetList$EUR],,] <- 0
priceReg <- regionSubsetList$EUR[!regionSubsetList$EUR %in% c("DEU","FRA","EUC","EUW")][1] #select region to define EU price (excluding Germany and France)
for (r in regionSubsetList$EUR){
co2EUprice[r,,] <- as.vector(co2EUprice[priceReg,,])
}
tmp <- mbind(tmp,setNames(co2EUprice, "Price|Carbon|EU-wide Regulation For All Sectors (US$2005/t CO2)"))
}
if (!is.null(pm_taxCO2eqPower)) {
tmp <- mbind(tmp,setNames((pm_taxCO2eqPower + pm_taxCO2eq) * 1000 * 12/44, "Price|Carbon|Power (US$2005/t CO2)"))
}
if (!is.null(pm_taxCO2eqSCC)) {
tmp <- mbind(tmp,setNames(abs(pm_taxCO2eqSCC) * 1000 * 12/44, "Price|Carbon|SCC (US$2005/t CO2)"))
}
tmp <- mbind(tmp,setNames(tmp[,,"Price|Carbon (US$2005/t CO2)"] * s_GWP_N2O, "Price|N2O (US$2005/t N2O)"))
tmp <- mbind(tmp,setNames(tmp[,,"Price|Carbon (US$2005/t CO2)"] * s_GWP_CH4, "Price|CH4 (US$2005/t CH4)"))
# some precalculations
x1 <- mselect(prodFE,all_enty1=c("fepet","fedie"))
x1 <- dimSums(x1,dim=c(3.1,3.3)) # reducing to fepet, fedie dimensions only
x2 <- mselect(cesIO,all_in = finenbalfehos$all_in)
if ((length(magclass::getNames(x2)) == 1)){
x2 <- setNames(x2,NULL)
}
v1 <- dimSums(abs(mselect(febal.m[,y,],all_enty=c("fedie","fepet")))*x1,dim=3) +
dimSums(abs(mselect(balfinen.m[,y,],all_enty=c("fehos")))*x2,dim=3)
v2 <- dimSums(x1,dim=3) + dimSums(x2,dim=3)
dimnames(v1)[[3]] <- "fhos"
dimnames(v2)[[3]] <- "fhos"
tmp <- mbind(tmp,setNames(v1/v2/(abs(budget.m + 1e-10)) * tdptwyr2dpgj, "Price|Final Energy|Liquids (US$2005/GJ)"))
# some precalculations
if ("seliq" %in% sety) {
tmp <- mbind(tmp,setNames(sebal.m[,,"seliq"]/(budget.m+1e-10) * tdptwyr2dpgj, "Price|Secondary Energy|Liquids (US$2005/GJ)"))
} else if ("sepet" %in% sety) {
x1 <- dimSums(mselect(demSe,all_enty="sedie"),dim=3) * sebal.m[,,"sedie"] +
dimSums(mselect(demSe,all_enty="sepet"),dim=3) * sebal.m[,,"sepet"]
x2 <- dimSums(mselect(demSe,all_enty="sedie"),dim=3) + dimSums(mselect(demSe,all_enty="sepet"),dim=3)
tmp <- mbind(tmp,setNames(x1/x2/(abs(budget.m) + 1e-10) * tdptwyr2dpgj, "Price|Secondary Energy|Liquids (US$2005/GJ)"))
} else {
x1 <- dimSums( dimSums(mselect(demSe,all_enty="seliqbio"),dim=3) * abs(sebal.m[,,"seliqbio"]) +
dimSums(mselect(demSe,all_enty="seliqfos"),dim=3) * abs(sebal.m[,,"seliqfos"]) )
x2 <- dimSums(mselect(demSe,all_enty="seliqbio"),dim=3) + dimSums(mselect(demSe,all_enty="seliqfos"),dim=3)
tmp <- mbind(tmp,setNames(x1/x2/(abs(budget.m) + 1e-10) * tdptwyr2dpgj, "Price|Secondary Energy|Liquids (US$2005/GJ)"))
}
if ("seso" %in% sety) {
tmp <- mbind(tmp,setNames(sebal.m[,,"seso"]/(budget.m+1e-10) * tdptwyr2dpgj , "Price|Secondary Energy|Solids (US$2005/GJ)"))
} else {
x1 <- dimSums( dimSums(mselect(demSe,all_enty="sesobio"),dim=3) * abs(sebal.m[,,"sesobio"]) +
dimSums(mselect(demSe,all_enty="sesofos"),dim=3) * abs(sebal.m[,,"sesofos"]) )
x2 <- dimSums(mselect(demSe,all_enty="sesobio"),dim=3) + dimSums(mselect(demSe,all_enty="sesofos"),dim=3)
tmp <- mbind(tmp,setNames(x1/x2/(abs(budget.m) + 1e-10) * tdptwyr2dpgj, "Price|Secondary Energy|Solids (US$2005/GJ)"))
}
if ("sega" %in% sety) {
tmp <- mbind(tmp,setNames(sebal.m[,,"sega"]/(budget.m+1e-10) * tdptwyr2dpgj , "Price|Secondary Energy|Gases (US$2005/GJ)"))
} else {
x1 <- dimSums( dimSums(mselect(demSe,all_enty="segabio"),dim=3) * abs(sebal.m[,,"segabio"]) +
dimSums(mselect(demSe,all_enty="segafos"),dim=3) * abs(sebal.m[,,"segafos"]) )
x2 <- dimSums(mselect(demSe,all_enty="segabio"),dim=3) + dimSums(mselect(demSe,all_enty="segafos"),dim=3)
tmp <- mbind(tmp,setNames(x1/x2/(abs(budget.m) + 1e-10) * tdptwyr2dpgj, "Price|Secondary Energy|Gases (US$2005/GJ)"))
}
# some precalculations
if ("seliq" %in% sety) {
x1 <- sebal.m[,,"seliq"]/budget.m * tdptwyr2dpgj
x1[which(dimSums(mselect(prodSE,all_enty="pebiolc",all_enty1="seliq"),dim=3) < 1e-5)] <- 0
} else if ("sepet" %in% sety) {
x1 <- sebal.m[,,"sedie"]/budget.m * tdptwyr2dpgj
x1[which(dimSums(mselect(prodSE,all_enty="pebiolc",all_enty1="sedie"),dim=3) < 1e-5)] <- 0
} else {
x1 <- sebal.m[,,"seliqbio"]/(budget.m+1e-10)*tdptwyr2dpgj
}
tmp <- mbind(tmp,setNames(x1, "Price|Secondary Energy|Liquids|Biomass (US$2005/GJ)"))
## average prices
# some precalculations
x1 <- mselect(febal.m,all_enty=fe_transp_fety35)
x2 <- mselect(prodFE,all_enty1=fe_transp_fety35)
x2 <- dimSums(x2,dim=c(3.1,3.3)) # reducing to fepet, fedie, feh2t and feelt dimensions only
x1 <- x1[,,order(magclass::getNames(x1))]
x2 <- x2[,,order(magclass::getNames(x2))]
dimnames(x2)[[3]] <- magclass::getNames(x1)
tmp <- mbind(tmp,setNames(dimSums(abs(x1)*x2,dim=3)/dimSums(x2,dim=3)/(abs(budget.m) + 1e-10) * tdptwyr2dpgj, "Average price for transport energy|FE (US$2005/GJ)"))
# some precalculations
x1 <- mselect(febal.m,all_enty=c("fedie","fepet"))
x2 <- mselect(prodFE,all_enty1=c("fedie","fepet"))
x2 <- dimSums(x2,dim=c(3.1,3.3)) # reducing to fepet, fedie dimensions only
tmp <- mbind(tmp,setNames(dimSums(abs(x1)*x2,dim=3)/dimSums(x2,dim=3)/(abs(budget.m) + 1e-10) * tdptwyr2dpgj, "Average price for liquid transport energy|FE (US$2005/GJ)"))
# some precalculations
share_seel_t <- dimSums(mselect(demSe,all_enty="seel",all_enty1="feelt"),dim=3)/dimSums(mselect(demSe,all_enty="seel"),dim=3)
share_seh2_t <- dimSums(mselect(demSe,all_enty="seh2",all_enty1="feh2t"),dim=3)/dimSums(mselect(demSe,all_enty="seh2"),dim=3)
if ("seliq" %in% sety) {
share_seliq_t <- (dimSums(mselect(demSe,all_enty="seliq",all_enty1="fedie"),dim=3) + dimSums(mselect(demSe,all_enty="seliq",all_enty1="fepet"),dim=3))/dimSums(mselect(demSe,all_enty="seliq"),dim=3)
if (!is.null(power_realisation)) {
x1 <- share_seel_t * dimSums(mselect(demSe,all_enty="seel"),dim=3) * abs(sebalSeel.m[,,"seel"]) +
share_seh2_t * dimSums(mselect(demSe,all_enty="seh2"),dim=3) * abs(sebal.m[,,"seh2"]) +
share_seliq_t * dimSums(mselect(demSe,all_enty="seliq"),dim=3) * abs(sebal.m[,,"seliq"])
} else {
x1 <- share_seel_t * dimSums(mselect(demSe,all_enty="seel"),dim=3) * abs(sebal.m[,,"seel"]) +
share_seh2_t * dimSums(mselect(demSe,all_enty="seh2"),dim=3) * abs(sebal.m[,,"seh2"]) +
share_seliq_t * dimSums(mselect(demSe,all_enty="seliq"),dim=3) * abs(sebal.m[,,"seliq"])
}
x2 <- share_seel_t * dimSums(mselect(demSe,all_enty="seel"),dim=3) +
share_seh2_t * dimSums(mselect(demSe,all_enty="seh2"),dim=3) +
share_seliq_t * dimSums(mselect(demSe,all_enty="seliq"),dim=3)
} else if ("sepet" %in% sety) {
share_sedie_t <- dimSums(mselect(demSe,all_enty="sedie",all_enty1="fedie"),dim=3)/dimSums(mselect(demSe,all_enty="sedie"),dim=3)
if (!is.null(power_realisation)) {
x1 <- share_seel_t * dimSums(mselect(demSe,all_enty="seel"),dim=3) * abs(sebalSeel.m[,,"seel"]) +
share_seh2_t * dimSums(mselect(demSe,all_enty="seh2"),dim=3) * abs(sebal.m[,,"seh2"]) +
share_sedie_t * dimSums(mselect(demSe,all_enty="sedie"),dim=3) * abs(sebal.m[,,"sedie"]) +
dimSums(mselect(demSe,all_enty="sepet"),dim=3) * abs(sebal.m[,,"sepet"])
} else {
x1 <- share_seel_t * dimSums(mselect(demSe,all_enty="seel"),dim=3) * abs(sebal.m[,,"seel"]) +
share_seh2_t * dimSums(mselect(demSe,all_enty="seh2"),dim=3) * abs(sebal.m[,,"seh2"]) +
share_sedie_t * dimSums(mselect(demSe,all_enty="sedie"),dim=3) * abs(sebal.m[,,"sedie"]) +
dimSums(mselect(demSe,all_enty="sepet"),dim=3) * abs(sebal.m[,,"sepet"])
}
x2 <- share_seel_t * dimSums(mselect(demSe,all_enty="seel"),dim=3) +
share_seh2_t * dimSums(mselect(demSe,all_enty="seh2"),dim=3) +
share_sedie_t * dimSums(mselect(demSe,all_enty="sedie"),dim=3) +
dimSums(mselect(demSe,all_enty="sepet"),dim=3)
} else {
share_seliqbio_t <- (dimSums(mselect(demSe,all_enty="seliqbio",all_enty1="fedie"),dim=3) + dimSums(mselect(demSe,all_enty="seliqbio",all_enty1="fepet"),dim=3))/dimSums(mselect(demSe,all_enty=c("seliqbio","seliqfos")),dim=3)
share_seliqfos_t <- (dimSums(mselect(demSe,all_enty="seliqfos",all_enty1="fedie"),dim=3) + dimSums(mselect(demSe,all_enty="seliqfos",all_enty1="fepet"),dim=3))/dimSums(mselect(demSe,all_enty=c("seliqbio","seliqfos")),dim=3)
if (!is.null(power_realisation)) {
x1 <- share_seel_t * dimSums(mselect(demSe,all_enty="seel"),dim=3) * abs(sebalSeel.m[,,"seel"]) +
share_seh2_t * dimSums(mselect(demSe,all_enty="seh2"),dim=3) * abs(sebal.m[,,"seh2"]) +
share_seliqbio_t * dimSums(mselect(demSe,all_enty="seliqbio"),dim=3) * abs(sebal.m[,,"seliqbio"]) +
share_seliqfos_t * dimSums(mselect(demSe,all_enty="seliqfos"),dim=3) * abs(sebal.m[,,"seliqfos"])
} else {
x1 <- share_seel_t * dimSums(mselect(demSe,all_enty="seel"),dim=3) * abs(sebal.m[,,"seel"]) +
share_seh2_t * dimSums(mselect(demSe,all_enty="seh2"),dim=3) * abs(sebal.m[,,"seh2"]) +
share_seliqbio_t * dimSums(mselect(demSe,all_enty="seliqbio"),dim=3) * abs(sebal.m[,,"seliqbio"]) +
share_seliqfos_t * dimSums(mselect(demSe,all_enty="seliqfos"),dim=3) * abs(sebal.m[,,"seliqfos"])
}
x2 <- share_seel_t * dimSums(mselect(demSe,all_enty="seel"),dim=3) +
share_seh2_t * dimSums(mselect(demSe,all_enty="seh2"),dim=3) +
share_seliqbio_t * dimSums(mselect(demSe,all_enty="seliqbio"),dim=3) +
share_seliqfos_t * dimSums(mselect(demSe,all_enty="seliqbio"),dim=3)
}
tmp <- mbind(tmp,setNames(x1/x2/(abs(budget.m) + 1e-10) * tdptwyr2dpgj, "Average price for transport energy|SE (US$2005/GJ)"))
# some precalculations
if ("seliq" %in% sety) {
tmp <- mbind(tmp,setNames(sebal.m[,,"seliq"]/(budget.m+1e-10) * tdptwyr2dpgj, "Average price for liquid transport energy|SE (US$2005/GJ)"))
} else if ("sepet" %in% sety) {
x1 <- share_sedie_t * dimSums(mselect(demSe,all_enty="sedie"),dim=3) * abs(sebal.m[,,"sedie"]) +
dimSums(mselect(demSe,all_enty="sepet"),dim=3) * abs(sebal.m[,,"sepet"])
x2 <- share_sedie_t * dimSums(mselect(demSe,all_enty="sedie"),dim=3) +
dimSums(mselect(demSe,all_enty="sepet"),dim=3)
tmp <- mbind(tmp,setNames(x1/x2/(abs(budget.m) + 1e-10) * tdptwyr2dpgj, "Average price for liquid transport energy|SE (US$2005/GJ)"))
} else {
x1 <- share_seliqbio_t * dimSums(mselect(demSe,all_enty="seliqbio"),dim=3) * abs(sebal.m[,,"seliqbio"]) +
share_seliqfos_t * dimSums(mselect(demSe,all_enty="seliqfos"),dim=3) * abs(sebal.m[,,"seliqfos"])
x2 <- share_seliqbio_t * dimSums(mselect(demSe,all_enty="seliqbio"),dim=3) +
share_seliqfos_t * dimSums(mselect(demSe,all_enty="seliqfos"),dim=3)
tmp <- mbind(tmp,setNames(x1/x2/(abs(budget.m) + 1e-10) * tdptwyr2dpgj, "Average price for liquid transport energy|SE (US$2005/GJ)"))
}
# some precalculations
x1 <- abs(balfinen.m[finenbal] / (budget.m + 1e-10)) * 1000
dimnames(x1)[[3]] <- dimnames(cesIO[,,finenbal[,2]])[[3]]
x1 <- dimSums(x1 * cesIO[,,finenbal[,2]],dim=3)
x2 <- dimSums(abs(esm2macro.m[,y,ppfenfromes] / (budget.m + 1e-10)) * 1000 * cesIO[,,ppfenfromes],dim=3)
x3 <- dimSums(cesIO[finenbal],dim=3) / TWa_2_EJ
x4 <- dimSums(cesIO[,,ppfenfromes],dim=3) / TWa_2_EJ
tmp <- mbind(tmp,setNames((x1 + x2)/(x3 + x4), "Average price for energy|CES (US$2005/GJ)"))
# ---- compute CES prices ----
tmp <- mbind(
tmp,
calc_CES_marginals(gdx, NULL) %>%
select('regi', 't', 'pf', 'price') %>%
left_join(
pm_delta_kap %>%
as.data.frame() %>%
select('regi' = 'Region', 'pf' = 'Data1',
'pm_delta_kap' = 'Value') %>%
filter(0 != !!sym('pm_delta_kap')),
c('regi', 'pf')
) %>%
mutate(
!!sym('price') := case_when(
# subtract depriciation rate from capital prices
!!sym('pf') %in% !!sym('ppfKap') ~ (
!!sym('price') - !!sym('pm_delta_kap')),
# convert unit of energy prices
!!sym('pf') %in% !!sym('ppfen_stat_build_ind') ~ (
!!sym('price') * tdptwyr2dpgj),
# keep all other prices
TRUE ~ !!sym('price')
),
!!sym('pf') := paste0(
'Price|CES_input|', !!sym('pf'), ' ',
case_when(
!!sym('pf') %in% !!sym('ppfKap') ~ '(Net of depreciation)',
!!sym('pf') %in% !!sym('ppfen_stat_build_ind') ~ '(US$2005/GJ)',
TRUE ~ '()'))) %>%
select(-'pm_delta_kap') %>%
as.magpie())
# find weights for prices for global aggregation
int2ext_prices_ces <- tibble(
intensive = grep('Price\\|CES_input', getNames(tmp), value = TRUE)) %>%
mutate(!!sym('extensive') := sub('^Price\\|', '', !!sym('intensive')),
!!sym('extensive') := sub('\\(US\\$2005/GJ\\)', '(EJ/yr)',
!!sym('extensive')),
!!sym('extensive') := sub('\\(Net of depreciation\\)',
'(billion US$2005)', !!sym('extensive'))) %>%
filter(!!sym('extensive') %in% getNames(output)) %>%
df.2.named.vector()
# ---- mapping of weights for the variables for global aggregation ----
int2ext <- c(
"Price|Biomass|Primary Level (US$2005/GJ)" = "PE|+|Biomass (EJ/yr)",
"Price|Biomass|Primary Level|Moving Avg (US$2005/GJ)" = "PE|+|Biomass (EJ/yr)",
"Price|Biomass|Emulator presolve (US$2005/GJ)" = "Primary Energy Production|Biomass|Energy Crops MAgPIE (EJ/yr)",
"Price|Biomass|Emulator presolve shifted (US$2005/GJ)" = "Primary Energy Production|Biomass|Energy Crops MAgPIE (EJ/yr)",
"Price|Biomass|Emulator shifted (US$2005/GJ)" = "Primary Energy Production|Biomass|Energy Crops (EJ/yr)",
"Price|Biomass|MAgPIE (US$2005/GJ)" = "Primary Energy Production|Biomass|Energy Crops MAgPIE (EJ/yr)",
"Price|Biomass|Bioenergy tax (US$2005/GJ)" = "Primary Energy Production|Biomass|Energy Crops (EJ/yr)",
"Price|N2O (US$2005/t N2O)" = "Emi|N2O (kt N2O/yr)",
"Price|CH4 (US$2005/t CH4)" = "Emi|CH4 (Mt CH4/yr)",
"Price|Coal|Primary Level (US$2005/GJ)" = "PE|+|Coal (EJ/yr)",
"Price|Coal|Secondary Level (US$2005/GJ)" = "PE|+|Coal (EJ/yr)",
"Price|Crude Oil|Primary Level (US$2005/GJ)" = "PE|+|Oil (EJ/yr)",
"Price|Natural Gas|Primary Level (US$2005/GJ)" = "PE|+|Gas (EJ/yr)",
"Price|Natural Gas|Primary Level|Moving Avg (US$2005/GJ)" = "PE|+|Gas (EJ/yr)",
"Price|Natural Gas|Secondary Level (US$2005/GJ)" = "PE|+|Gas (EJ/yr)",
"Price|Final Energy|Liquids (US$2005/GJ)" = "FE|+|Liquids (EJ/yr)",
"Price|Final Energy|Electricity|Transport (US$2005/GJ)" = "FE|Transport|Electricity (EJ/yr)",
"Price|Final Energy|Electricity|Transport|Moving Avg (US$2005/GJ)" = "FE|Transport|Electricity (EJ/yr)",
"Price|Final Energy|Hydrogen|Transport (US$2005/GJ)" = "FE|Transport|Hydrogen (EJ/yr)",
"Price|Final Energy|Liquids|Transport (US$2005/GJ)" = "FE|Transport|Liquids (EJ/yr)",
"Price|Final Energy|Liquids|Transport|Moving Avg (US$2005/GJ)"= "FE|Transport|Liquids (EJ/yr)",
"Price|Final Energy|Transport (US$2005/GJ)" = "FE|Transport (EJ/yr)",
"Price|Final Energy|Transport|Moving Avg (US$2005/GJ)" = "FE|Transport (EJ/yr)",
"Price|Energy Service|Transport LDV (US$2005/GJ)" = "CES_input|fepet (EJ/yr)", ## TODO: check units
"Price|Energy Service|Transport nonLDV (US$2005/GJ)" = "CES_input|fedie (EJ/yr)", ## TODO: check units
"Average price for transport energy|FE (US$2005/GJ)" = "FE|Transport (EJ/yr)",
"Average price for transport energy|SE (US$2005/GJ)" = "FE|Transport (EJ/yr)",
"Average price for liquid transport energy|FE (US$2005/GJ)" = "FE|Transport|Liquids (EJ/yr)",
"Average price for liquid transport energy|SE (US$2005/GJ)" = "FE|Transport|Liquids (EJ/yr)",
"Price|Uranium|Primary Level (US$2005/GJ)" = "PE|+|Nuclear (EJ/yr)",
"Price|Secondary Energy|Electricity (US$2005/GJ)" = "SE|Electricity (EJ/yr)",
"Price|Secondary Energy|Electricity|Moving Avg (US$2005/GJ)" = "SE|Electricity (EJ/yr)",
"Price|Secondary Energy|Biomass (US$2005/GJ)" = "SE|Biomass (EJ/yr)",
"Price|Secondary Energy|Hydrogen (US$2005/GJ)" = "SE|Hydrogen (EJ/yr)",
"Price|Secondary Energy|Solids (US$2005/GJ)" = "SE|Solids (EJ/yr)",
"Price|Secondary Energy|Gases (US$2005/GJ)" = "SE|Gases (EJ/yr)",
"Price|Secondary Energy|Heat (US$2005/GJ)" = "SE|Heat (EJ/yr)",
"Price|Secondary Energy|Liquids|Biomass (US$2005/GJ)" = "SE|Liquids|Biomass (EJ/yr)"
)
int2ext <- c(int2ext, int2ext_prices_ces)
if ("seliq" %in% sety) {
int2ext <- c(int2ext,
"Price|Secondary Energy|Liquids (US$2005/GJ)" = "SE|Liquids (EJ/yr)"
)
} else if ("sepet" %in% sety) {
int2ext <- c(int2ext,
"Price|SE|sepet (US$2005/GJ)" = "SE|Liquids|sepet (EJ/yr)",
"Price|SE|sedie (US$2005/GJ)" = "SE|Liquids|sedie (EJ/yr)"
)
} else {
int2ext <- c(int2ext,
"Price|Secondary Energy|Liquids (US$2005/GJ)" = "SE|Liquids (EJ/yr)",
"Price|SE|seliqbio (US$2005/GJ)" = "SE|Liquids|Biomass (EJ/yr)" ,
"Price|SE|seliqfos (US$2005/GJ)" = "SE|Liquids|Fossil (EJ/yr)"
)
}
if (all(cm_emiscen == 9)) int2ext <- c(int2ext, c( "Price|Carbon (US$2005/t CO2)" = "Emi|CO2|FFaI|Emi to which CO2tax is applied (Mt CO2/yr)"))
if (stat_mod == "simple" ){
int2ext <- c(int2ext, c( "Price|Final Energy|Heating Oil (US$2005/GJ)" = "FE|Other Sector|Liquids (EJ/yr)"),
"Price|Final Energy|Gases (US$2005/GJ)" = "CES_input|fegas (EJ/yr)",
"Price|Final Energy|Solids (US$2005/GJ)" = "CES_input|fesos (EJ/yr)",
"Price|Final Energy|Heat (US$2005/GJ)" = "CES_input|fehes (EJ/yr)",
"Price|Final Energy|Electricity|Stationary (US$2005/GJ)" = "CES_input|feels (EJ/yr)",
"Price|Final Energy|Hydrogen|Stationary (US$2005/GJ)" = "CES_input|feh2s (EJ/yr)"
)
}
if (buil_mod %in% c("simple", "services_putty","services_with_capital")){
int2ext <- c(int2ext,c(
"Price|Final Energy|Heating Oil|Buildings (US$2005/GJ)" = "FE|Buildings|Liquids (EJ/yr)",
"Price|Final Energy|Gases|Buildings (US$2005/GJ)" = "FE|Buildings|Gases (EJ/yr)",
"Price|Final Energy|Solids|Buildings (US$2005/GJ)" = "FE|Buildings|Solids (EJ/yr)",
"Price|Final Energy|Heat|Buildings (US$2005/GJ)" = "FE|Buildings|Heat (EJ/yr)",
"Price|Final Energy|Electricity|Buildings (US$2005/GJ)" = "FE|Buildings|Electricity (EJ/yr)",
"Price|Final Energy|Electricity|Buildings|Moving Avg (US$2005/GJ)" = "FE|Buildings|Electricity (EJ/yr)",
"Price|Final Energy|Hydrogen|Buildings (US$2005/GJ)" = "FE|Buildings|Hydrogen (EJ/yr)",
"Price|Final Energy|Buildings (US$2005/GJ)" = "FE|Buildings (EJ/yr)",
"Price|Final Energy|Buildings|Moving Avg (US$2005/GJ)" = "FE|Buildings (EJ/yr)"
))
}
if (indu_mod %in% c('fixed_shares', 'subsectors')){
int2ext <- c(int2ext,c(
"Price|Final Energy|Heating Oil|Industry (US$2005/GJ)" = "FE|Industry|Liquids (EJ/yr)",
"Price|Final Energy|Gases|Industry (US$2005/GJ)" = "FE|Industry|Gases (EJ/yr)",
"Price|Final Energy|Solids|Industry (US$2005/GJ)" = "FE|Industry|Solids (EJ/yr)",
"Price|Final Energy|Heat|Industry (US$2005/GJ)" = "FE|Industry|Heat (EJ/yr)",
"Price|Final Energy|Electricity|Industry (US$2005/GJ)" = "FE|Industry|Electricity (EJ/yr)",
"Price|Final Energy|Hydrogen|Industry (US$2005/GJ)" = "FE|Industry|Hydrogen (EJ/yr)",
"Price|Final Energy|Industry (US$2005/GJ)" = "FE|Industry (EJ/yr)",
"Price|Final Energy|Industry|Moving Avg (US$2005/GJ)" = "FE|Industry (EJ/yr)"
))
}
if("fegat" %in% getNames(febal.m)){
int2ext <- c(int2ext, "Price|Final Energy|Gases|Transport (US$2005/GJ)" = "FE|Transport|Gases (EJ/yr)")
}
output[is.na(output)] <- 0 # substitute na by 0
# add global prices
map <- data.frame(region=getRegions(tmp),world="GLO",stringsAsFactors=FALSE)
tmp_GLO <- new.magpie("GLO",getYears(tmp),magclass::getNames(tmp),fill=0)
for (i2e in names(int2ext)){
tmp_GLO["GLO",,i2e] <- speed_aggregate(tmp[,,i2e],map,weight=output[map$region,,int2ext[i2e]])
for(t in getYears(tmp)){
if(all(output[map$region,t,int2ext[i2e]]==0)){
tmp_GLO["GLO",t,i2e] <- NA
}
}
}
tmp <- mbind(tmp,tmp_GLO)
# add other region aggregations
if (!is.null(regionSubsetList)){
tmp_RegAgg <- new.magpie(names(regionSubsetList),getYears(tmp),magclass::getNames(tmp),fill=0)
for(region in names(regionSubsetList)){
tmp_RegAgg_ie2 <- do.call("mbind",lapply(names(int2ext), function(i2e) {
map <- data.frame(region=regionSubsetList[[region]],parentRegion=region,stringsAsFactors=FALSE)
result <- speed_aggregate(tmp[regionSubsetList[[region]],,i2e],map,weight=output[regionSubsetList[[region]],,int2ext[i2e]])
getRegions(result) <- region
for(t in getYears(tmp)){
if(all(output[regionSubsetList[[region]],t,int2ext[i2e]]==0)){
result[region,t,i2e] <- NA
}
}
return(result)
}))
tmp_RegAgg[region,,names(int2ext)] <- tmp_RegAgg_ie2[region,,names(int2ext)]
}
tmp <- mbind(tmp,tmp_RegAgg)
}
tmp[is.na(tmp)] <- 0 # tmp is NA if weight is zero for all regions within the GLO or the specific region aggregation. Therefore, we replace all NAs with zeros.
# prices that are same for all regions, including GLO
glob_price <- new.magpie(getRegions(tmp),getYears(tmp),fill=0)
for(i in getRegions(tmp)) glob_price[i,,] <- pm_pvp[,,"peur"]*1000
tmp <- mbind(tmp,setNames(glob_price, "PVP1|Uranium (billionDpktU)"))
for(i in getRegions(tmp)) glob_price[i,,] <- pm_pvp[,,"peoil"]*1000
tmp <- mbind(tmp,setNames(glob_price, "PVP1|Oil (billionDpTWyr)"))
for(i in getRegions(tmp)) glob_price[i,,] <- pm_pvp[,,"pegas"]*1000
tmp <- mbind(tmp,setNames(glob_price, "PVP1|Gas (billionDpTWyr)"))
for(i in getRegions(tmp)) glob_price[i,,] <- pm_pvp[,,"pecoal"]*1000
tmp <- mbind(tmp,setNames(glob_price, "PVP1|Coal (billionDpTWyr)"))
for(i in getRegions(tmp)) glob_price[i,,] <- pm_pvp[,,"pebiolc"]*1000
tmp <- mbind(tmp,setNames(glob_price, "PVP1|Biomass (billionDpTWyr)"))
for(i in getRegions(tmp)) glob_price[i,,] <- pm_pvp[,,"good"]*1000
tmp <- mbind(tmp,setNames(glob_price, "PVP2|Good ()"))
for(i in getRegions(tmp)) glob_price[i,,] <- pm_pvp[,,"perm"]
tmp <- mbind(tmp,setNames(glob_price, "PVP3|Permit ()"))
for(i in getRegions(tmp)) glob_price[i,,] <- pm_pvp[,,"peur"]/pm_pvp[,,"good"] * tdptwyr2dpgj
tmp <- mbind(tmp,setNames(glob_price, "Price|Uranium|World Market (US$2005/GJ)"))
for(i in getRegions(tmp)) glob_price[i,,] <- pm_pvp[,,"peoil"]/pm_pvp[,,"good"] * tdptwyr2dpgj
tmp <- mbind(tmp,setNames(glob_price, "Price|Oil|World Market (US$2005/GJ)"))
for(i in getRegions(tmp)) glob_price[i,,] <- pm_pvp[,,"pegas"]/pm_pvp[,,"good"] * tdptwyr2dpgj
tmp <- mbind(tmp,setNames(glob_price, "Price|Gas|World Market (US$2005/GJ)"))
for(i in getRegions(tmp)) glob_price[i,,] <- pm_pvp[,,"pecoal"]/pm_pvp[,,"good"] * tdptwyr2dpgj
tmp <- mbind(tmp,setNames(glob_price, "Price|Coal|World Market (US$2005/GJ)"))
for(i in getRegions(tmp)) glob_price[i,,] <- pm_pvp[,,"pebiolc"]/pm_pvp[,,"good"] * tdptwyr2dpgj
tmp <- mbind(tmp,setNames(glob_price, "Price|Biomass|World Market (US$2005/GJ)"))
## special global prices
#AJS calc global carbon price as average over regional pm_pvpRegi's, weighted by total emissions.
regi_on_gdx <- unique(readGDX(gdx, name = "regi2iso")[,1])
tmp["GLO",,"Price|Carbon (US$2005/t CO2)"] <-
dimSums( pm_pvpRegi[regi_on_gdx,,"perm"] * output[regi_on_gdx,,"Emi|CO2|FFaI|Emi to which CO2tax is applied (Mt CO2/yr)"],dim=1 ) /
dimSums(output[regi_on_gdx,,"Emi|CO2|FFaI|Emi to which CO2tax is applied (Mt CO2/yr)"],dim=1) /
(pm_pvp[1,,"good"] + 1e-10) * 1000*12/44
# add other region aggregations carbon price as average over regional pm_pvpRegi's, weighted by total emissions.
if (!is.null(regionSubsetList)){
for(region in names(regionSubsetList)){
tmp[region,,"Price|Carbon (US$2005/t CO2)"] <- dimSums( pm_pvpRegi[regionSubsetList[[region]],,"perm"] * output[regionSubsetList[[region]],,"Emi|CO2|FFaI|Emi to which CO2tax is applied (Mt CO2/yr)"],dim=1 ) /
dimSums(output[regionSubsetList[[region]],,"Emi|CO2|FFaI|Emi to which CO2tax is applied (Mt CO2/yr)"],dim=1) /
(pm_pvp[1,,"good"] + 1e-10) * 1000*12/44;
}
}
## not meaningful global prices set to NA
tmp["GLO",,"Price|Biomass|Shiftfactor ()"] <- NA
if (!is.null(regionSubsetList$EUR))
tmp["EUR",,"Price|Carbon|EU-wide Regulation For All Sectors (US$2005/t CO2)"] <- as.vector(co2EUprice[priceReg,,])
## not meaningful region aggregation prices set to NA
if (!is.null(regionSubsetList)){
for(region in names(regionSubsetList)){
tmp[region,,"Price|Biomass|Shiftfactor ()"] <- NA
}
}
# ---- debug information for industry/subsectors ----
if ('subsectors' == indu_mod & !is.null(q37_limit_secondary_steel_share.m)) {
.x <- q37_limit_secondary_steel_share.m[,y,] / budget.m
tmp2 <- mbind(
# fake some GLO data
setNames(
mbind(.x, dimSums(.x * NA, dim = 1)),
'Debug|Industry|Secondary Steel Premium (US$2005)'),
mbind(
lapply(
list(
c('ue_industry', '', 'arbitrary unit', 1),
c('ue_cement', '|Cement', 't cement', 1e3),
c('ue_chemicals', '|Chemicals', 'arbitrary unit', 1),
c('ue_steel', '|Steel', 't Steel', 1e3),
c('ue_steel_primary', '|Steel|Primary', 't Steel', 1e3),
c('ue_steel_secondary', '|Steel|Secondary', 't Steel', 1e3),
c('ue_otherInd', '|other', 'arbitrary unit', 1)),
function(x) {
setNames(
( tmp[,,paste0('Price|CES_input|', x[1], ' ('), pmatch = 'left']
* as.numeric(x[4])
),
paste0('Debug|Price|Industry', x[2], ' (US$2005/', x[3], ')'))
})
)
)
tmp <- mbind(
tmp,
mbind(
tmp2,
calc_regionSubset_sums(tmp2, regionSubsetList)
)
)
}
return(tmp)
}
pik-piam/remind documentation built on Sept. 9, 2021, 1:09 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions reportPrices
Documented in reportPrices
#' Read in GDX and calculate prices, used in convGDX2MIF.R for the reporting
#'
#' Read in price information from GDX file, information used in convGDX2MIF.R
#' for the reporting
#'
#'
#' @param gdx a GDX object as created by readGDX, or the path to a gdx
#' @param output a magpie object containing all needed variables generated by other report*.R functions
#' @param regionSubsetList a list containing regions to create report variables region
#' aggregations. If NULL (default value) only the global region aggregation "GLO" will
#' be created.
#' @return MAgPIE object - contains the price variables
#' @importFrom luscale speed_aggregate
#' @author David Klein
#' @seealso \code{\link{convGDX2MIF}}
#' @examples
#'
#' \dontrun{reportPrices(gdx)}
#'
#' @importFrom dplyr %>% case_when distinct filter_ inner_join tibble
#' @importFrom gdx readGDX
#' @importFrom luscale speed_aggregate
#' @importFrom magclass mbind getYears getRegions setNames dimSums new.magpie lowpass complete_magpie
#' @importFrom quitte df.2.named.vector getColValues
#' @importFrom readr read_csv
#' @export
reportPrices <- function(gdx,output=NULL,regionSubsetList=NULL) {
if(is.null(output)){
output <- reportPE(gdx,regionSubsetList)
output <- mbind(output,reportSE(gdx,regionSubsetList))
output <- mbind(output,reportFE(gdx,regionSubsetList))
output <- mbind(output,reportEmi(gdx,regionSubsetList))
output <- mbind(output,reportExtraction(gdx,regionSubsetList))
output <- mbind(output,reportMacroEconomy(gdx,regionSubsetList)[,getYears(output),])
}
#---- Functions
find_real_module <- function(module_set, module_name){
return(module_set[module_set$modules == module_name,2])
}
compute_agg_price_fe = function(mobj_prices,mobj_output,sector){
if (!sector %in% c("Buildings","Industry","Transport")) stop("argument 'sector' must be in c('Buildings','Industry','Transport')")
names_prices = grep(paste0("^Price\\|Final Energy\\|[A-z ]+\\|",sector," \\(US\\$2005\\/GJ\\)"),getNames(mobj_prices),value = T)
names_fe = gsub(paste0("^Price\\|Final Energy\\|([A-z ]+)\\|",sector," \\(US\\$2005\\/GJ\\)"),
paste0("FE|",sector,"\\|\\1 (EJ/yr)"), names_prices)
names_fe = gsub("Heating Oil","Liquids",names_fe)
names(names_fe) = names_prices
# price_total = sum(i,price_i*quant_i)/quant_total
mobj_prices = mbind(mobj_prices,
setNames(
dimSums(
setNames(
mobj_prices[,,names_prices],names_fe[names_prices])
* mobj_output[regs,time,names_fe]
,dim = 3)
/collapseNames(mobj_output[regs,time,paste0("FE|",sector," (EJ/yr)")]),
paste0("Price|Final Energy|",sector," (US$2005/GJ)")
))
# Moving Averages
mobj_prices = mbind(mobj_prices,
setNames( lowpass(
dimSums(
setNames(
mobj_prices[,,names_prices],names_fe[names_prices])
* mobj_output[regs,time,names_fe]
,dim = 3)
/collapseNames(mobj_output[regs,time,paste0("FE|",sector," (EJ/yr)")]), fix="both", altFilter=match(2010,time)),
paste0("Price|Final Energy|",sector,"|Moving Avg (US$2005/GJ)")
))
return(mobj_prices)
}
#---- End of Functions
####### get realisations #########
realisation <- readGDX(gdx, "module2realisation")
####### get power realisations #########
realisation <- readGDX(gdx, "module2realisation")
power_realisation <- if ("power" %in% realisation[, 1]) realisation[which(realisation[,1] == "power"),2]
####### conversion factors ##########
s_GWP_CH4 <- readGDX(gdx,c("sm_gwpCH4","s_gwpCH4","s_GWP_CH4"),format="first_found", react = "silent")
s_GWP_N2O <- readGDX(gdx,c("s_gwpN2O","s_GWP_N2O"),format="first_found", react = "silent")
TWa_2_EJ <- 31.536
tdptwyr2dpgj <- 31.71 #TerraDollar per TWyear to Dollar per GJ
p80_subset <- c("perm","good","peur","peoil","pegas","pecoal","pebiolc") #TODO: read in from gdx as sets trade
pebal_subset <- c("pebiolc","peoil","pegas","pecoal","peur")
####### read in needed data #########
## sets
fe_transp_fety35 <- readGDX(gdx,name="FE_Transp_fety35",types="sets",format="first_found", react = "silent")
se2fe <- readGDX(gdx,name="se2fe",types="sets",format="first_found", react = "silent")
sety <- readGDX(gdx,c("entySe","sety"),types="sets",format="first_found", react = "silent")
fety <- readGDX(gdx,name=c("entyFe","fety"),types="sets",format="first_found", react = "silent")
finenbal <- readGDX(gdx,name=c("fe2ppfEn","finenbal"),types="sets",format="first_found", react = "silent")
all_esty <- readGDX(gdx,name=c("all_esty"),types="sets",format="first_found", react = "silent")
if(!is.null(all_esty)){
fe2es <- readGDX(gdx,name=c("fe2es"),types="sets",format="first_found", react = "silent")
fe2es <- unique(fe2es[c("all_enty","all_esty")])
colnames(fe2es) <- gsub("all_esty","all_in", colnames(fe2es))
} else {
fe2es = NULL
}
ppfenfromes <- readGDX(gdx,name=c("ppfEnFromEs","ppfenfromes"),types="sets",format="first_found", react = "silent")
finenbalfehos <- finenbal %>% filter_(~all_enty == "fehos")
ppfKap <- readGDX(gdx,"ppfKap")
ppfen_stat <- readGDX(gdx,c("ppfen_stationary_dyn38","ppfen_stationary_dyn28","ppfen_stationary"),format="first_found", react = "silent")
if (length(ppfen_stat) == 0) ppfen_stat = NULL
ppfen_build <- readGDX(gdx,c("ppfen_buildings_dyn36","ppfen_buildings_dyn28","ppfen_buildings"),format="first_found", react = "silent")
ue_dyn36 <- readGDX(gdx,c("ue_dyn36"),format="first_found", react = "silent")
ppfen_build <- setdiff(ppfen_build, ue_dyn36)
ppfen_ind <- readGDX(gdx,c("ppfen_industry_dyn37","ppfen_industry_dyn28","ppfen_industry"),format="first_found", react = "silent")
ppfen_stat_build_ind <- c(ppfen_stat,ppfen_build,ppfen_ind)
edge_scen <- readGDX(gdx,"EDGE_scenario",format="first_found", react = "silent")
# Realisation of the different modules
if ( !is.null(realisation) & (!"CES_structure" %in% realisation[, 1] )){
stat_mod = find_real_module(realisation,"stationary")
tran_mod = find_real_module(realisation,"transport")
indu_mod = find_real_module(realisation,"industry")
buil_mod = find_real_module(realisation,"buildings")
} else {
if ( !is.null(ppfen_stat)){ # In case the set realisation did not exist, find out whether it was stationary or buildings-industry
stat_mod = "simple"
indu_mod = "off"
buil_mod = "off"
} else {
stat_mod = "off"
indu_mod = "fixed_shares"
buil_mod = "simple"
}
if ( !is.null(edge_scen)){ # In case the set realization did not exist, find out whether "edge_esm" or "complex" realization was chosen
tran_mod = "edge_esm"
} else {
tran_mod = "complex"
}
}
if (indu_mod %in% c('fixed_shares', 'subsectors')) {
fe2ppfen37 <- readGDX(gdx, 'fe2ppfen') %>%
filter(.data$all_in %in% ppfen_ind)
}
## parameter
shift_p <- readGDX(gdx,name="p30_pebiolc_pricshift",format="first_found")
mult_p <- readGDX(gdx,name="p30_pebiolc_pricmult",format="first_found")
pric_mag <- readGDX(gdx,name="p30_pebiolc_pricemag",format="first_found")
pric_emu_pre <- readGDX(gdx,name="p30_pebiolc_price_emu_preloop",format="first_found")
pric_emu_pre_shifted <- readGDX(gdx,name="p30_pebiolc_price_emu_preloop_shifted",format="first_found")
bio_tax_factor <- readGDX(gdx,name="p21_tau_bioenergy_tax",format="first_found")
pm_pvp <- readGDX(gdx,name=c("pm_pvp","p80_pvp"),format="first_found")[,,p80_subset]
pm_ts <- readGDX(gdx,name='pm_ts',format="first_found")
pm_dataemi <- readGDX(gdx,name=c("pm_emifac","pm_dataemi"),format="first_found",restore_zeros=FALSE)[,,c("pegas.seel.ngt.co2","pecoal.seel.pc.co2")]
pm_pvpRegi <- readGDX(gdx,name='pm_pvpRegi',format="first_found")[,,"perm"]
pm_taxCO2eq <- readGDX(gdx,name=c("pm_taxCO2eq","pm_tau_CO2_tax"),format="first_found")
pm_taxCO2eqPower <- readGDX(gdx,name=c("pm_taxCO2eqPower"),format="first_found")
pm_taxCO2eqSCC <- readGDX(gdx,name='pm_taxCO2eqSCC',format="first_found")
pm_delta_kap <- readGDX(gdx, c("pm_delta_kap", "p_delta_kap"), format = "first_found")
if(is.null(getNames(pm_delta_kap))) pm_delta_kap = setNames(pm_delta_kap,"kap")
## variables
pric_emu <- readGDX(gdx,name="vm_pebiolc_price",field="l",format="first_found")
prodFE <- readGDX(gdx,name='vm_prodFe',field="l",format="first_found",restore_zeros = FALSE)
prodSE <- readGDX(gdx,name='vm_prodSe',field="l",format="first_found",restore_zeros = FALSE)
demSe <- readGDX(gdx,name=c("vm_demSe","v_demSe"),types="variables",field="l",format="first_found",restore_zeros=FALSE)
cesIO <- readGDX(gdx,name='vm_cesIO',types="variables",field="l",format="first_found",restore_zeros=FALSE)
fe_taxCES <- readGDX(gdx, name=c('p21_tau_fe_tax','pm_tau_fe_tax'), format="first_found", react = F)[ppfen_stat_build_ind]
fe_subCES <- readGDX(gdx, name=c('p21_tau_fe_sub','pm_tau_fe_sub'), format="first_found", react = F)[ppfen_stat_build_ind]
if(is.null(fe_taxCES) & is.null(fe_subCES)){
fe_taxCES = readGDX(gdx, name=c('p21_tau_fe_tax_bit_st'))[,,ppfen_stat_build_ind]
fe_subCES = readGDX(gdx, name=c('p21_tau_fe_sub_bit_st'))[,,ppfen_stat_build_ind]
}
getSets(fe_taxCES) = gsub("all_enty","all_in", getSets(fe_taxCES))
getSets(fe_subCES) = gsub("all_enty","all_in", getSets(fe_subCES))
if (!is.null(all_esty)){
fe_taxES = readGDX(gdx, name=c("pm_tau_fe_tax_ES_st",'p21_tau_fe_tax_ES_st'),format = "first_found", react = "silent")
fe_subES = readGDX(gdx, name=c('pm_tau_fe_sub_ES_st','p21_tau_fe_sub_ES_st'),format = "first_found", react = "silent")
getSets(fe_taxES) = gsub("all_esty","all_in", getSets(fe_taxES))
getSets(fe_subES) = gsub("all_esty","all_in", getSets(fe_subES))
} else {
fe_taxES = NULL
fe_subES = NULL
}
## equations
pebal.m <- readGDX(gdx,name=c("q_balPe","qm_pebal"),types = "equations",field = "m",format = "first_found")[,,pebal_subset]
budget.m <- readGDX(gdx,name='qm_budget',types = "equations",field = "m",format = "first_found") # Alternative: calcPrice
sebal.m <- readGDX(gdx,name=c("q_balSe","q_sebal"),types="equations",field="m",format="first_found")
if (!is.null(power_realisation)) {
sebalSeel.m <- readGDX(gdx,name="q32_balSe",types="equations",field="m",format="first_found")
}
esm2macro.m <- readGDX(gdx,name='q_esm2macro',types="equations",field="m",format="first_found")
febal.m <- readGDX(gdx,name=c("q_balFe","q_febal"),types="equations",field="m",format="first_found")[,,fety]
balfinen.m <- readGDX(gdx,name=c("qm_balFeForCesAndEs","qm_balFeForCes","q_balFeForCes","q_balfinen"),types="equations",field="m",format="first_found", react = "silent")
cm_emiscen <- readGDX(gdx,name='cm_emiscen',format="first_found")
q37_limit_secondary_steel_share.m <- readGDX(
gdx, name = 'q37_limit_secondary_steel_share', types = 'equation',
field = 'm', react = 'silent')
#####################################
#choose the CES entries names for transport
name_trsp=c("fepet","ueLDVt","fedie","ueHDVt","feelt","ueelTt")
name_trsp=name_trsp[name_trsp%in%getNames(esm2macro.m)]
#####################################
####### calculate minimal temporal resolution #####
y <- Reduce(intersect,list(getYears(febal.m),getYears(prodFE),getYears(sebal.m),getYears(prodSE)))
febal.m <- febal.m[,y,]
balfinen.m <- balfinen.m[,y,]
budget.m <- budget.m[,y,]
sebal.m <- sebal.m[,y,]
if (!is.null(power_realisation)) {
sebalSeel.m <- sebalSeel.m[,y,]
}
pm_pvp <- pm_pvp[,y,]
pebal.m <- pebal.m[,y,]
shift_p <- shift_p[,y,]
mult_p <- mult_p[,y,]
pric_mag <- pric_mag[,y,]
pric_emu_pre <- pric_emu_pre[,y,]
pric_emu_pre_shifted <- pric_emu_pre_shifted[,y,]
pric_emu <- pric_emu[,y,]
bio_tax_factor <- bio_tax_factor[,y,]
esm2macro.m <- esm2macro.m[,y,]
pm_dataemi <- pm_dataemi[,y,]
pm_ts <- pm_ts[,y,]
pm_pvpRegi <- pm_pvpRegi[,y,]
pm_taxCO2eq <- pm_taxCO2eq[,y,]
if (!is.null(pm_taxCO2eqPower)) {
pm_taxCO2eqPower <- pm_taxCO2eqPower[,y,]
}
if (!is.null(pm_taxCO2eqSCC)) {
pm_taxCO2eqSCC <- pm_taxCO2eqSCC[,y,]
}
cesIO <- cesIO[,y,]
####### pm_pvp = EPS results in fantasy carbon prices
for(t in getYears(pm_pvp)){
if(pm_pvp[,t,"good"] == 0){
pm_pvp[,t,"good"] = 0.0001;
}
}
####### calculate reporting parameters ############
tmp <- NULL
# Calculate and append new variables to magpie object
# Costs and prices for purpose-grown 2nd gen. bioenergy
# - Shiftfactor
# - imported values from MAgPIE reporting
# - results from emulator based on MAgPIE demand (before main solve)
# - results from emulator based on MAgPIE demand multiplied with shiftfactor (before main solve)
# - results from emulator based on REMIND demand (after main solve)
# - results from emulator based on REMIND demand multiplied with shiftfactor (after main solve)
time <- getYears(budget.m,as.integer=TRUE)
regs <- getRegions(budget.m)
# report variables
tmp <- NULL
tmp <- mbind(tmp,setNames(pebal.m[,,"pebiolc"] / (budget.m+1e-10) * tdptwyr2dpgj, "Price|Biomass|Primary Level (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(lowpass(pebal.m[,,"pebiolc"] / (budget.m+1e-10), fix="both", altFilter=match(2010,time)) * tdptwyr2dpgj, "Price|Biomass|Primary Level|Moving Avg (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(shift_p, "Price|Biomass|Shiftfactor ()"))
tmp <- mbind(tmp,setNames(mult_p, "Price|Biomass|Multfactor ()"))
tmp <- mbind(tmp,setNames(pric_mag * tdptwyr2dpgj, "Price|Biomass|MAgPIE (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(pric_emu_pre * tdptwyr2dpgj, "Price|Biomass|Emulator presolve (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(pric_emu_pre_shifted * tdptwyr2dpgj, "Price|Biomass|Emulator presolve shifted (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(pric_emu * tdptwyr2dpgj, "Price|Biomass|Emulator shifted (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(pric_emu * bio_tax_factor * tdptwyr2dpgj, "Price|Biomass|Bioenergy tax (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(pebal.m[,,"peoil"]/(budget.m+1e-10) * tdptwyr2dpgj, "Price|Crude Oil|Primary Level (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(pebal.m[,,"pegas"]/(budget.m+1e-10) * tdptwyr2dpgj, "Price|Natural Gas|Primary Level (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(lowpass(pebal.m[,,"pegas"]/(budget.m+1e-10), fix="both", altFilter=match(2010,time)) * tdptwyr2dpgj, "Price|Natural Gas|Primary Level|Moving Avg (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(pebal.m[,,"pecoal"]/(budget.m+1e-10) * tdptwyr2dpgj, "Price|Coal|Primary Level (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(pebal.m[,,"peur"]/(budget.m+1e-10) * tdptwyr2dpgj, "Price|Uranium|Primary Level (US$2005/GJ)"))
# seondary energy prices
if (!is.null(power_realisation)) {
tmp <- mbind(tmp,setNames(sebalSeel.m[,,"seel"]/(budget.m+1e-10) * tdptwyr2dpgj, "Price|Secondary Energy|Electricity (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(lowpass(sebalSeel.m[,,"seel"]/(budget.m+1e-10), fix="both", altFilter=match(2010,time)) * tdptwyr2dpgj, "Price|Secondary Energy|Electricity|Moving Avg (US$2005/GJ)"))
} else {
tmp <- mbind(tmp,setNames(sebal.m[,,"seel"]/(budget.m+1e-10) * tdptwyr2dpgj, "Price|Secondary Energy|Electricity (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(lowpass(sebal.m[,,"seel"]/(budget.m+1e-10), fix="both", altFilter=match(2010,time)) * tdptwyr2dpgj, "Price|Secondary Energy|Electricity|Moving Avg (US$2005/GJ)"))
}
# for biomass price for secondary level is the same as primary level
tmp <- mbind(tmp,setNames(pebal.m[,,"pebiolc"]/(budget.m+1e-10) * tdptwyr2dpgj, "Price|Secondary Energy|Biomass (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(sebal.m[,,"seh2"]/(budget.m+1e-10) * tdptwyr2dpgj , "Price|Secondary Energy|Hydrogen (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(sebal.m[,,"sehe"]/(budget.m+1e-10) * tdptwyr2dpgj , "Price|Secondary Energy|Heat (US$2005/GJ)"))
# energy services
tmp <- mbind(tmp,setNames(abs(esm2macro.m[,,name_trsp[2]]/(budget.m+1e-10)) * tdptwyr2dpgj , "Price|Energy Service|Transport nonLDV (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(abs(esm2macro.m[,,name_trsp[1]]/(budget.m+1e-10)) * tdptwyr2dpgj , "Price|Energy Service|Transport LDV (US$2005/GJ)"))
#Final energy prices
#Transport (Taxes already included)
if (tran_mod == "complex"){
#Translate the marginal utility of the constraint into the marginal income (price).
#For transport liquids use average between diesel and petrol.
tmp <- mbind(tmp,setNames(abs(febal.m[,,"feelt"]/(budget.m+1e-10)) * tdptwyr2dpgj, "Price|Final Energy|Electricity|Transport (US$2005/GJ)"))
if("fegat" %in% getNames(febal.m)){
tmp <- mbind(tmp,setNames(abs(febal.m[,,"fegat"]/(budget.m+1e-10)) * tdptwyr2dpgj, "Price|Final Energy|Gases|Transport (US$2005/GJ)"))
}
tmp <- mbind(tmp,setNames(abs(lowpass(febal.m[,,"feelt"]/(budget.m+1e-10), fix="both", altFilter=match(2010,time))) * tdptwyr2dpgj, "Price|Final Energy|Electricity|Transport|Moving Avg (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(abs(febal.m[,,"feh2t"]/(budget.m+1e-10)) * tdptwyr2dpgj, "Price|Final Energy|Hydrogen|Transport (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(abs(febal.m[,,"fedie"]/(budget.m+1e-10)) * tdptwyr2dpgj, "Price|Final Energy|Liquids|Transport (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(abs(lowpass((febal.m[,,"fedie"]+febal.m[,,"fepet"])/(2*budget.m+1e-10), fix="both", altFilter=match(2010,time))) * tdptwyr2dpgj, "Price|Final Energy|Liquids|Transport|Moving Avg (US$2005/GJ)"))
} else if(tran_mod == "edge_esm"){
#Translate the marginal utility of the constraint into the marginal income (price)
#For transport liquids use average between diesel and petrol.
tmp <- mbind(tmp,setNames(abs(balfinen.m[,,"feelt"]/(budget.m+1e-10)) * tdptwyr2dpgj, "Price|Final Energy|Electricity|Transport (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(abs(balfinen.m[,,"fegat"]/(budget.m+1e-10)) * tdptwyr2dpgj, "Price|Final Energy|Gases|Transport (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(abs(lowpass(balfinen.m[,,"feelt"]/(budget.m+1e-10), fix="both", altFilter=match(2010,time))) * tdptwyr2dpgj, "Price|Final Energy|Electricity|Transport|Moving Avg (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(abs(balfinen.m[,,"feh2t"]/(budget.m+1e-10)) * tdptwyr2dpgj, "Price|Final Energy|Hydrogen|Transport (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(abs((balfinen.m[,,"fedie"]+balfinen.m[,,"fepet"])/(2*budget.m+1e-10)) * tdptwyr2dpgj, "Price|Final Energy|Liquids|Transport (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(abs(lowpass((balfinen.m[,,"fedie"]+balfinen.m[,,"fepet"])/(2*budget.m+1e-10), fix="both", altFilter=match(2010,time))) * tdptwyr2dpgj, "Price|Final Energy|Liquids|Transport|Moving Avg (US$2005/GJ)"))
}
tmp = compute_agg_price_fe(tmp,output,"Transport")
#Buildings and Industry (Taxes need to be added)
a = abs(mselect(balfinen.m[,y,][rbind(finenbal,fe2es)]))/abs((budget.m+1e-10)) # Translate the marginal utility of the constraint into the marginal income (price)
b = complete_magpie(a) # Due to a strange behaviour of magclass objects addition, we need to use complete_magpie to make the addition
prices_fe_bi = (b + mbind(fe_taxCES[,y,getColValues(finenbal,"all_in")],
fe_taxES[,y,getColValues(fe2es,"all_in")])
+ mbind(fe_subCES[,y,getColValues(finenbal,"all_in")],
fe_subES[,y,getColValues(fe2es,"all_in")])) * tdptwyr2dpgj # add the taxes and subsidies for the prices of buildings and industry. For transport, they are in the marginal of febal
prices_fe_bi = prices_fe_bi[,,getNames(a)]
if (stat_mod == "simple" ){
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"feels.feels"], "Price|Final Energy|Electricity|Stationary (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"fegas.fegas"], "Price|Final Energy|Gases (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"fehes.fehes"], "Price|Final Energy|Heat (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"fehos.fehos"], "Price|Final Energy|Heating Oil (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"feh2s.feh2s"], "Price|Final Energy|Hydrogen|Stationary (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"fesos.fesos"], "Price|Final Energy|Solids (US$2005/GJ)"))
}
if (buil_mod == "simple"){
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"feels.feelb"], "Price|Final Energy|Electricity|Buildings (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(lowpass(prices_fe_bi[,,"feels.feelb"], fix="both", altFilter=match(2010,time)) , "Price|Final Energy|Electricity|Buildings|Moving Avg (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"fegas.fegab"], "Price|Final Energy|Gases|Buildings (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"fehes.feheb"], "Price|Final Energy|Heat|Buildings (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"fehos.fehob"], "Price|Final Energy|Heating Oil|Buildings (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"feh2s.feh2b"], "Price|Final Energy|Hydrogen|Buildings (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"fesos.fesob"], "Price|Final Energy|Solids|Buildings (US$2005/GJ)"))
tmp = compute_agg_price_fe(tmp,output,"Buildings")
}
if (buil_mod %in% c("services_putty", "services_with_capital")){
# Prices across energy services should be identical, so we take the price from only one service (cooking and water heating)
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"feels.uecwelb"], "Price|Final Energy|Electricity|Buildings (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(lowpass(prices_fe_bi[,,"feels.uecwelb"], fix="both", altFilter=match(2010,time)) , "Price|Final Energy|Electricity|Buildings|Moving Avg (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"fegas.uecwgab"], "Price|Final Energy|Gases|Buildings (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"fehes.uecwheb"], "Price|Final Energy|Heat|Buildings (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"fehos.uecwhob"], "Price|Final Energy|Heating Oil|Buildings (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"feh2s.uecwh2b"], "Price|Final Energy|Hydrogen|Buildings (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(prices_fe_bi[,,"fesos.uecwsob"], "Price|Final Energy|Solids|Buildings (US$2005/GJ)"))
tmp = compute_agg_price_fe(tmp,output,"Buildings")
}
if (indu_mod %in% c('fixed_shares', 'subsectors')) {
industry_fe_mixer <- inner_join(
fe2ppfen37,
paste(
'all_enty, fe',
'fesos, Solids',
'fehos, Heating Oil',
'fegas, Gases',
'feh2s, Hydrogen',
'fehes, Heat',
'feels, Electricity',
sep = '\n') %>%
read_csv(),
'all_enty'
)
for (.fe in unique(industry_fe_mixer$fe)) {
fe2ppfen <- industry_fe_mixer %>%
filter(.data$fe == .fe) %>%
distinct(.data$all_enty, .keep_all = TRUE) %>%
select(-.data$fe) %>%
paste(collapse = '.')
tmp <- mbind(
tmp,
setNames(
prices_fe_bi[,,fe2ppfen],
paste0('Price|Final Energy|', .fe, '|Industry (US$2005/GJ)')
)
)
tmp <- mbind(
tmp,
setNames(
lowpass(prices_fe_bi[,,fe2ppfen] , fix="both", altFilter=match(2010,time)) ,
paste0('Price|Final Energy|', .fe, '|Industry|Moving Avg (US$2005/GJ)')
)
)
}
rm(industry_fe_mixer, fe2ppfen)
tmp = compute_agg_price_fe(tmp,output,"Industry")
}
if ("seliq" %in% sety) {
tmp <- mbind(tmp,setNames(sebal.m[,,"seliq"]/(budget.m+1e-10)*tdptwyr2dpgj, "Price|SE|seliq (US$2005/GJ)" ))
} else if ("sepet" %in% sety) {
tmp <- mbind(tmp,setNames(sebal.m[,,"sepet"]/(budget.m+1e-10)*tdptwyr2dpgj, "Price|SE|sepet (US$2005/GJ)" ))
tmp <- mbind(tmp,setNames(sebal.m[,,"sedie"]/(budget.m+1e-10)*tdptwyr2dpgj, "Price|SE|sedie (US$2005/GJ)" ))
} else {
tmp <- mbind(tmp,setNames(sebal.m[,,"seliqbio"]/(budget.m+1e-10)*tdptwyr2dpgj, "Price|SE|seliqbio (US$2005/GJ)" ))
tmp <- mbind(tmp,setNames(sebal.m[,,"seliqfos"]/(budget.m+1e-10)*tdptwyr2dpgj, "Price|SE|seliqfos (US$2005/GJ)" ))
}
tmp <- mbind(tmp,setNames((pebal.m[,,"pegas"]/pm_ts - pm_dataemi[,,"pegas.seel.ngt.co2"]*pm_pvpRegi) / (budget.m/pm_ts + 1e-10) * tdptwyr2dpgj, "Price|Natural Gas|Secondary Level (US$2005/GJ)"))
tmp <- mbind(tmp,setNames((pebal.m[,,"pecoal"]/pm_ts - pm_dataemi[,,"pecoal.seel.pc.co2"]*pm_pvpRegi) / (budget.m/pm_ts + 1e-10) * tdptwyr2dpgj, "Price|Coal|Secondary Level (US$2005/GJ)"))
tmp <- mbind(tmp,setNames(abs(pm_pvpRegi / (pm_pvp[,,"good"] + 1e-10)) * 1000 * 12/44, "Price|Carbon (US$2005/t CO2)"))
tmp <- mbind(tmp,setNames(abs(pm_taxCO2eq) * 1000 * 12/44, "Price|Carbon|Guardrail (US$2005/t CO2)"))
if (is.null(regionSubsetList$EUR)) {
tmp <- mbind(tmp,setNames(pm_taxCO2eq * 1000 * 12/44, "Price|Carbon|EU-wide Regulation For All Sectors (US$2005/t CO2)"))
} else {
co2EUprice <- pm_taxCO2eq * 1000 * 12/44
co2EUprice[getRegions(pm_taxCO2eq)[!getRegions(pm_taxCO2eq) %in% regionSubsetList$EUR],,] <- 0
priceReg <- regionSubsetList$EUR[!regionSubsetList$EUR %in% c("DEU","FRA","EUC","EUW")][1] #select region to define EU price (excluding Germany and France)
for (r in regionSubsetList$EUR){
co2EUprice[r,,] <- as.vector(co2EUprice[priceReg,,])
}
tmp <- mbind(tmp,setNames(co2EUprice, "Price|Carbon|EU-wide Regulation For All Sectors (US$2005/t CO2)"))
}
if (!is.null(pm_taxCO2eqPower)) {
tmp <- mbind(tmp,setNames((pm_taxCO2eqPower + pm_taxCO2eq) * 1000 * 12/44, "Price|Carbon|Power (US$2005/t CO2)"))
}
if (!is.null(pm_taxCO2eqSCC)) {
tmp <- mbind(tmp,setNames(abs(pm_taxCO2eqSCC) * 1000 * 12/44, "Price|Carbon|SCC (US$2005/t CO2)"))
}
tmp <- mbind(tmp,setNames(tmp[,,"Price|Carbon (US$2005/t CO2)"] * s_GWP_N2O, "Price|N2O (US$2005/t N2O)"))
tmp <- mbind(tmp,setNames(tmp[,,"Price|Carbon (US$2005/t CO2)"] * s_GWP_CH4, "Price|CH4 (US$2005/t CH4)"))
# some precalculations
x1 <- mselect(prodFE,all_enty1=c("fepet","fedie"))
x1 <- dimSums(x1,dim=c(3.1,3.3)) # reducing to fepet, fedie dimensions only
x2 <- mselect(cesIO,all_in = finenbalfehos$all_in)
if ((length(magclass::getNames(x2)) == 1)){
x2 <- setNames(x2,NULL)
}
v1 <- dimSums(abs(mselect(febal.m[,y,],all_enty=c("fedie","fepet")))*x1,dim=3) +
dimSums(abs(mselect(balfinen.m[,y,],all_enty=c("fehos")))*x2,dim=3)
v2 <- dimSums(x1,dim=3) + dimSums(x2,dim=3)
dimnames(v1)[[3]] <- "fhos"
dimnames(v2)[[3]] <- "fhos"
tmp <- mbind(tmp,setNames(v1/v2/(abs(budget.m + 1e-10)) * tdptwyr2dpgj, "Price|Final Energy|Liquids (US$2005/GJ)"))
# some precalculations
if ("seliq" %in% sety) {
tmp <- mbind(tmp,setNames(sebal.m[,,"seliq"]/(budget.m+1e-10) * tdptwyr2dpgj, "Price|Secondary Energy|Liquids (US$2005/GJ)"))
} else if ("sepet" %in% sety) {
x1 <- dimSums(mselect(demSe,all_enty="sedie"),dim=3) * sebal.m[,,"sedie"] +
dimSums(mselect(demSe,all_enty="sepet"),dim=3) * sebal.m[,,"sepet"]
x2 <- dimSums(mselect(demSe,all_enty="sedie"),dim=3) + dimSums(mselect(demSe,all_enty="sepet"),dim=3)
tmp <- mbind(tmp,setNames(x1/x2/(abs(budget.m) + 1e-10) * tdptwyr2dpgj, "Price|Secondary Energy|Liquids (US$2005/GJ)"))
} else {
x1 <- dimSums( dimSums(mselect(demSe,all_enty="seliqbio"),dim=3) * abs(sebal.m[,,"seliqbio"]) +
dimSums(mselect(demSe,all_enty="seliqfos"),dim=3) * abs(sebal.m[,,"seliqfos"]) )
x2 <- dimSums(mselect(demSe,all_enty="seliqbio"),dim=3) + dimSums(mselect(demSe,all_enty="seliqfos"),dim=3)
tmp <- mbind(tmp,setNames(x1/x2/(abs(budget.m) + 1e-10) * tdptwyr2dpgj, "Price|Secondary Energy|Liquids (US$2005/GJ)"))
}
if ("seso" %in% sety) {
tmp <- mbind(tmp,setNames(sebal.m[,,"seso"]/(budget.m+1e-10) * tdptwyr2dpgj , "Price|Secondary Energy|Solids (US$2005/GJ)"))
} else {
x1 <- dimSums( dimSums(mselect(demSe,all_enty="sesobio"),dim=3) * abs(sebal.m[,,"sesobio"]) +
dimSums(mselect(demSe,all_enty="sesofos"),dim=3) * abs(sebal.m[,,"sesofos"]) )
x2 <- dimSums(mselect(demSe,all_enty="sesobio"),dim=3) + dimSums(mselect(demSe,all_enty="sesofos"),dim=3)
tmp <- mbind(tmp,setNames(x1/x2/(abs(budget.m) + 1e-10) * tdptwyr2dpgj, "Price|Secondary Energy|Solids (US$2005/GJ)"))
}
if ("sega" %in% sety) {
tmp <- mbind(tmp,setNames(sebal.m[,,"sega"]/(budget.m+1e-10) * tdptwyr2dpgj , "Price|Secondary Energy|Gases (US$2005/GJ)"))
} else {
x1 <- dimSums( dimSums(mselect(demSe,all_enty="segabio"),dim=3) * abs(sebal.m[,,"segabio"]) +
dimSums(mselect(demSe,all_enty="segafos"),dim=3) * abs(sebal.m[,,"segafos"]) )
x2 <- dimSums(mselect(demSe,all_enty="segabio"),dim=3) + dimSums(mselect(demSe,all_enty="segafos"),dim=3)
tmp <- mbind(tmp,setNames(x1/x2/(abs(budget.m) + 1e-10) * tdptwyr2dpgj, "Price|Secondary Energy|Gases (US$2005/GJ)"))
}
# some precalculations
if ("seliq" %in% sety) {
x1 <- sebal.m[,,"seliq"]/budget.m * tdptwyr2dpgj
x1[which(dimSums(mselect(prodSE,all_enty="pebiolc",all_enty1="seliq"),dim=3) < 1e-5)] <- 0
} else if ("sepet" %in% sety) {
x1 <- sebal.m[,,"sedie"]/budget.m * tdptwyr2dpgj
x1[which(dimSums(mselect(prodSE,all_enty="pebiolc",all_enty1="sedie"),dim=3) < 1e-5)] <- 0
} else {
x1 <- sebal.m[,,"seliqbio"]/(budget.m+1e-10)*tdptwyr2dpgj
}
tmp <- mbind(tmp,setNames(x1, "Price|Secondary Energy|Liquids|Biomass (US$2005/GJ)"))
## average prices
# some precalculations
x1 <- mselect(febal.m,all_enty=fe_transp_fety35)
x2 <- mselect(prodFE,all_enty1=fe_transp_fety35)
x2 <- dimSums(x2,dim=c(3.1,3.3)) # reducing to fepet, fedie, feh2t and feelt dimensions only
x1 <- x1[,,order(magclass::getNames(x1))]
x2 <- x2[,,order(magclass::getNames(x2))]
dimnames(x2)[[3]] <- magclass::getNames(x1)
tmp <- mbind(tmp,setNames(dimSums(abs(x1)*x2,dim=3)/dimSums(x2,dim=3)/(abs(budget.m) + 1e-10) * tdptwyr2dpgj, "Average price for transport energy|FE (US$2005/GJ)"))
# some precalculations
x1 <- mselect(febal.m,all_enty=c("fedie","fepet"))
x2 <- mselect(prodFE,all_enty1=c("fedie","fepet"))
x2 <- dimSums(x2,dim=c(3.1,3.3)) # reducing to fepet, fedie dimensions only
tmp <- mbind(tmp,setNames(dimSums(abs(x1)*x2,dim=3)/dimSums(x2,dim=3)/(abs(budget.m) + 1e-10) * tdptwyr2dpgj, "Average price for liquid transport energy|FE (US$2005/GJ)"))
# some precalculations
share_seel_t <- dimSums(mselect(demSe,all_enty="seel",all_enty1="feelt"),dim=3)/dimSums(mselect(demSe,all_enty="seel"),dim=3)
share_seh2_t <- dimSums(mselect(demSe,all_enty="seh2",all_enty1="feh2t"),dim=3)/dimSums(mselect(demSe,all_enty="seh2"),dim=3)
if ("seliq" %in% sety) {
share_seliq_t <- (dimSums(mselect(demSe,all_enty="seliq",all_enty1="fedie"),dim=3) + dimSums(mselect(demSe,all_enty="seliq",all_enty1="fepet"),dim=3))/dimSums(mselect(demSe,all_enty="seliq"),dim=3)
if (!is.null(power_realisation)) {
x1 <- share_seel_t * dimSums(mselect(demSe,all_enty="seel"),dim=3) * abs(sebalSeel.m[,,"seel"]) +
share_seh2_t * dimSums(mselect(demSe,all_enty="seh2"),dim=3) * abs(sebal.m[,,"seh2"]) +
share_seliq_t * dimSums(mselect(demSe,all_enty="seliq"),dim=3) * abs(sebal.m[,,"seliq"])
} else {
x1 <- share_seel_t * dimSums(mselect(demSe,all_enty="seel"),dim=3) * abs(sebal.m[,,"seel"]) +
share_seh2_t * dimSums(mselect(demSe,all_enty="seh2"),dim=3) * abs(sebal.m[,,"seh2"]) +
share_seliq_t * dimSums(mselect(demSe,all_enty="seliq"),dim=3) * abs(sebal.m[,,"seliq"])
}
x2 <- share_seel_t * dimSums(mselect(demSe,all_enty="seel"),dim=3) +
share_seh2_t * dimSums(mselect(demSe,all_enty="seh2"),dim=3) +
share_seliq_t * dimSums(mselect(demSe,all_enty="seliq"),dim=3)
} else if ("sepet" %in% sety) {
share_sedie_t <- dimSums(mselect(demSe,all_enty="sedie",all_enty1="fedie"),dim=3)/dimSums(mselect(demSe,all_enty="sedie"),dim=3)
if (!is.null(power_realisation)) {
x1 <- share_seel_t * dimSums(mselect(demSe,all_enty="seel"),dim=3) * abs(sebalSeel.m[,,"seel"]) +
share_seh2_t * dimSums(mselect(demSe,all_enty="seh2"),dim=3) * abs(sebal.m[,,"seh2"]) +
share_sedie_t * dimSums(mselect(demSe,all_enty="sedie"),dim=3) * abs(sebal.m[,,"sedie"]) +
dimSums(mselect(demSe,all_enty="sepet"),dim=3) * abs(sebal.m[,,"sepet"])
} else {
x1 <- share_seel_t * dimSums(mselect(demSe,all_enty="seel"),dim=3) * abs(sebal.m[,,"seel"]) +
share_seh2_t * dimSums(mselect(demSe,all_enty="seh2"),dim=3) * abs(sebal.m[,,"seh2"]) +
share_sedie_t * dimSums(mselect(demSe,all_enty="sedie"),dim=3) * abs(sebal.m[,,"sedie"]) +
dimSums(mselect(demSe,all_enty="sepet"),dim=3) * abs(sebal.m[,,"sepet"])
}
x2 <- share_seel_t * dimSums(mselect(demSe,all_enty="seel"),dim=3) +
share_seh2_t * dimSums(mselect(demSe,all_enty="seh2"),dim=3) +
share_sedie_t * dimSums(mselect(demSe,all_enty="sedie"),dim=3) +
dimSums(mselect(demSe,all_enty="sepet"),dim=3)
} else {
share_seliqbio_t <- (dimSums(mselect(demSe,all_enty="seliqbio",all_enty1="fedie"),dim=3) + dimSums(mselect(demSe,all_enty="seliqbio",all_enty1="fepet"),dim=3))/dimSums(mselect(demSe,all_enty=c("seliqbio","seliqfos")),dim=3)
share_seliqfos_t <- (dimSums(mselect(demSe,all_enty="seliqfos",all_enty1="fedie"),dim=3) + dimSums(mselect(demSe,all_enty="seliqfos",all_enty1="fepet"),dim=3))/dimSums(mselect(demSe,all_enty=c("seliqbio","seliqfos")),dim=3)
if (!is.null(power_realisation)) {
x1 <- share_seel_t * dimSums(mselect(demSe,all_enty="seel"),dim=3) * abs(sebalSeel.m[,,"seel"]) +
share_seh2_t * dimSums(mselect(demSe,all_enty="seh2"),dim=3) * abs(sebal.m[,,"seh2"]) +
share_seliqbio_t * dimSums(mselect(demSe,all_enty="seliqbio"),dim=3) * abs(sebal.m[,,"seliqbio"]) +
share_seliqfos_t * dimSums(mselect(demSe,all_enty="seliqfos"),dim=3) * abs(sebal.m[,,"seliqfos"])
} else {
x1 <- share_seel_t * dimSums(mselect(demSe,all_enty="seel"),dim=3) * abs(sebal.m[,,"seel"]) +
share_seh2_t * dimSums(mselect(demSe,all_enty="seh2"),dim=3) * abs(sebal.m[,,"seh2"]) +
share_seliqbio_t * dimSums(mselect(demSe,all_enty="seliqbio"),dim=3) * abs(sebal.m[,,"seliqbio"]) +
share_seliqfos_t * dimSums(mselect(demSe,all_enty="seliqfos"),dim=3) * abs(sebal.m[,,"seliqfos"])
}
x2 <- share_seel_t * dimSums(mselect(demSe,all_enty="seel"),dim=3) +
share_seh2_t * dimSums(mselect(demSe,all_enty="seh2"),dim=3) +
share_seliqbio_t * dimSums(mselect(demSe,all_enty="seliqbio"),dim=3) +
share_seliqfos_t * dimSums(mselect(demSe,all_enty="seliqbio"),dim=3)
}
tmp <- mbind(tmp,setNames(x1/x2/(abs(budget.m) + 1e-10) * tdptwyr2dpgj, "Average price for transport energy|SE (US$2005/GJ)"))
# some precalculations
if ("seliq" %in% sety) {
tmp <- mbind(tmp,setNames(sebal.m[,,"seliq"]/(budget.m+1e-10) * tdptwyr2dpgj, "Average price for liquid transport energy|SE (US$2005/GJ)"))
} else if ("sepet" %in% sety) {
x1 <- share_sedie_t * dimSums(mselect(demSe,all_enty="sedie"),dim=3) * abs(sebal.m[,,"sedie"]) +
dimSums(mselect(demSe,all_enty="sepet"),dim=3) * abs(sebal.m[,,"sepet"])
x2 <- share_sedie_t * dimSums(mselect(demSe,all_enty="sedie"),dim=3) +
dimSums(mselect(demSe,all_enty="sepet"),dim=3)
tmp <- mbind(tmp,setNames(x1/x2/(abs(budget.m) + 1e-10) * tdptwyr2dpgj, "Average price for liquid transport energy|SE (US$2005/GJ)"))
} else {
x1 <- share_seliqbio_t * dimSums(mselect(demSe,all_enty="seliqbio"),dim=3) * abs(sebal.m[,,"seliqbio"]) +
share_seliqfos_t * dimSums(mselect(demSe,all_enty="seliqfos"),dim=3) * abs(sebal.m[,,"seliqfos"])
x2 <- share_seliqbio_t * dimSums(mselect(demSe,all_enty="seliqbio"),dim=3) +
share_seliqfos_t * dimSums(mselect(demSe,all_enty="seliqfos"),dim=3)
tmp <- mbind(tmp,setNames(x1/x2/(abs(budget.m) + 1e-10) * tdptwyr2dpgj, "Average price for liquid transport energy|SE (US$2005/GJ)"))
}
# some precalculations
x1 <- abs(balfinen.m[finenbal] / (budget.m + 1e-10)) * 1000
dimnames(x1)[[3]] <- dimnames(cesIO[,,finenbal[,2]])[[3]]
x1 <- dimSums(x1 * cesIO[,,finenbal[,2]],dim=3)
x2 <- dimSums(abs(esm2macro.m[,y,ppfenfromes] / (budget.m + 1e-10)) * 1000 * cesIO[,,ppfenfromes],dim=3)
x3 <- dimSums(cesIO[finenbal],dim=3) / TWa_2_EJ
x4 <- dimSums(cesIO[,,ppfenfromes],dim=3) / TWa_2_EJ
tmp <- mbind(tmp,setNames((x1 + x2)/(x3 + x4), "Average price for energy|CES (US$2005/GJ)"))
# ---- compute CES prices ----
tmp <- mbind(
tmp,
calc_CES_marginals(gdx, NULL) %>%
select('regi', 't', 'pf', 'price') %>%
left_join(
pm_delta_kap %>%
as.data.frame() %>%
select('regi' = 'Region', 'pf' = 'Data1',
'pm_delta_kap' = 'Value') %>%
filter(0 != !!sym('pm_delta_kap')),
c('regi', 'pf')
) %>%
mutate(
!!sym('price') := case_when(
# subtract depriciation rate from capital prices
!!sym('pf') %in% !!sym('ppfKap') ~ (
!!sym('price') - !!sym('pm_delta_kap')),
# convert unit of energy prices
!!sym('pf') %in% !!sym('ppfen_stat_build_ind') ~ (
!!sym('price') * tdptwyr2dpgj),
# keep all other prices
TRUE ~ !!sym('price')
),
!!sym('pf') := paste0(
'Price|CES_input|', !!sym('pf'), ' ',
case_when(
!!sym('pf') %in% !!sym('ppfKap') ~ '(Net of depreciation)',
!!sym('pf') %in% !!sym('ppfen_stat_build_ind') ~ '(US$2005/GJ)',
TRUE ~ '()'))) %>%
select(-'pm_delta_kap') %>%
as.magpie())
# find weights for prices for global aggregation
int2ext_prices_ces <- tibble(
intensive = grep('Price\\|CES_input', getNames(tmp), value = TRUE)) %>%
mutate(!!sym('extensive') := sub('^Price\\|', '', !!sym('intensive')),
!!sym('extensive') := sub('\\(US\\$2005/GJ\\)', '(EJ/yr)',
!!sym('extensive')),
!!sym('extensive') := sub('\\(Net of depreciation\\)',
'(billion US$2005)', !!sym('extensive'))) %>%
filter(!!sym('extensive') %in% getNames(output)) %>%
df.2.named.vector()
# ---- mapping of weights for the variables for global aggregation ----
int2ext <- c(
"Price|Biomass|Primary Level (US$2005/GJ)" = "PE|+|Biomass (EJ/yr)",
"Price|Biomass|Primary Level|Moving Avg (US$2005/GJ)" = "PE|+|Biomass (EJ/yr)",
"Price|Biomass|Emulator presolve (US$2005/GJ)" = "Primary Energy Production|Biomass|Energy Crops MAgPIE (EJ/yr)",
"Price|Biomass|Emulator presolve shifted (US$2005/GJ)" = "Primary Energy Production|Biomass|Energy Crops MAgPIE (EJ/yr)",
"Price|Biomass|Emulator shifted (US$2005/GJ)" = "Primary Energy Production|Biomass|Energy Crops (EJ/yr)",
"Price|Biomass|MAgPIE (US$2005/GJ)" = "Primary Energy Production|Biomass|Energy Crops MAgPIE (EJ/yr)",
"Price|Biomass|Bioenergy tax (US$2005/GJ)" = "Primary Energy Production|Biomass|Energy Crops (EJ/yr)",
"Price|N2O (US$2005/t N2O)" = "Emi|N2O (kt N2O/yr)",
"Price|CH4 (US$2005/t CH4)" = "Emi|CH4 (Mt CH4/yr)",
"Price|Coal|Primary Level (US$2005/GJ)" = "PE|+|Coal (EJ/yr)",
"Price|Coal|Secondary Level (US$2005/GJ)" = "PE|+|Coal (EJ/yr)",
"Price|Crude Oil|Primary Level (US$2005/GJ)" = "PE|+|Oil (EJ/yr)",
"Price|Natural Gas|Primary Level (US$2005/GJ)" = "PE|+|Gas (EJ/yr)",
"Price|Natural Gas|Primary Level|Moving Avg (US$2005/GJ)" = "PE|+|Gas (EJ/yr)",
"Price|Natural Gas|Secondary Level (US$2005/GJ)" = "PE|+|Gas (EJ/yr)",
"Price|Final Energy|Liquids (US$2005/GJ)" = "FE|+|Liquids (EJ/yr)",
"Price|Final Energy|Electricity|Transport (US$2005/GJ)" = "FE|Transport|Electricity (EJ/yr)",
"Price|Final Energy|Electricity|Transport|Moving Avg (US$2005/GJ)" = "FE|Transport|Electricity (EJ/yr)",
"Price|Final Energy|Hydrogen|Transport (US$2005/GJ)" = "FE|Transport|Hydrogen (EJ/yr)",
"Price|Final Energy|Liquids|Transport (US$2005/GJ)" = "FE|Transport|Liquids (EJ/yr)",
"Price|Final Energy|Liquids|Transport|Moving Avg (US$2005/GJ)"= "FE|Transport|Liquids (EJ/yr)",
"Price|Final Energy|Transport (US$2005/GJ)" = "FE|Transport (EJ/yr)",
"Price|Final Energy|Transport|Moving Avg (US$2005/GJ)" = "FE|Transport (EJ/yr)",
"Price|Energy Service|Transport LDV (US$2005/GJ)" = "CES_input|fepet (EJ/yr)", ## TODO: check units
"Price|Energy Service|Transport nonLDV (US$2005/GJ)" = "CES_input|fedie (EJ/yr)", ## TODO: check units
"Average price for transport energy|FE (US$2005/GJ)" = "FE|Transport (EJ/yr)",
"Average price for transport energy|SE (US$2005/GJ)" = "FE|Transport (EJ/yr)",
"Average price for liquid transport energy|FE (US$2005/GJ)" = "FE|Transport|Liquids (EJ/yr)",
"Average price for liquid transport energy|SE (US$2005/GJ)" = "FE|Transport|Liquids (EJ/yr)",
"Price|Uranium|Primary Level (US$2005/GJ)" = "PE|+|Nuclear (EJ/yr)",
"Price|Secondary Energy|Electricity (US$2005/GJ)" = "SE|Electricity (EJ/yr)",
"Price|Secondary Energy|Electricity|Moving Avg (US$2005/GJ)" = "SE|Electricity (EJ/yr)",
"Price|Secondary Energy|Biomass (US$2005/GJ)" = "SE|Biomass (EJ/yr)",
"Price|Secondary Energy|Hydrogen (US$2005/GJ)" = "SE|Hydrogen (EJ/yr)",
"Price|Secondary Energy|Solids (US$2005/GJ)" = "SE|Solids (EJ/yr)",
"Price|Secondary Energy|Gases (US$2005/GJ)" = "SE|Gases (EJ/yr)",
"Price|Secondary Energy|Heat (US$2005/GJ)" = "SE|Heat (EJ/yr)",
"Price|Secondary Energy|Liquids|Biomass (US$2005/GJ)" = "SE|Liquids|Biomass (EJ/yr)"
)
int2ext <- c(int2ext, int2ext_prices_ces)
if ("seliq" %in% sety) {
int2ext <- c(int2ext,
"Price|Secondary Energy|Liquids (US$2005/GJ)" = "SE|Liquids (EJ/yr)"
)
} else if ("sepet" %in% sety) {
int2ext <- c(int2ext,
"Price|SE|sepet (US$2005/GJ)" = "SE|Liquids|sepet (EJ/yr)",
"Price|SE|sedie (US$2005/GJ)" = "SE|Liquids|sedie (EJ/yr)"
)
} else {
int2ext <- c(int2ext,
"Price|Secondary Energy|Liquids (US$2005/GJ)" = "SE|Liquids (EJ/yr)",
"Price|SE|seliqbio (US$2005/GJ)" = "SE|Liquids|Biomass (EJ/yr)" ,
"Price|SE|seliqfos (US$2005/GJ)" = "SE|Liquids|Fossil (EJ/yr)"
)
}
if (all(cm_emiscen == 9)) int2ext <- c(int2ext, c( "Price|Carbon (US$2005/t CO2)" = "Emi|CO2|FFaI|Emi to which CO2tax is applied (Mt CO2/yr)"))
if (stat_mod == "simple" ){
int2ext <- c(int2ext, c( "Price|Final Energy|Heating Oil (US$2005/GJ)" = "FE|Other Sector|Liquids (EJ/yr)"),
"Price|Final Energy|Gases (US$2005/GJ)" = "CES_input|fegas (EJ/yr)",
"Price|Final Energy|Solids (US$2005/GJ)" = "CES_input|fesos (EJ/yr)",
"Price|Final Energy|Heat (US$2005/GJ)" = "CES_input|fehes (EJ/yr)",
"Price|Final Energy|Electricity|Stationary (US$2005/GJ)" = "CES_input|feels (EJ/yr)",
"Price|Final Energy|Hydrogen|Stationary (US$2005/GJ)" = "CES_input|feh2s (EJ/yr)"
)
}
if (buil_mod %in% c("simple", "services_putty","services_with_capital")){
int2ext <- c(int2ext,c(
"Price|Final Energy|Heating Oil|Buildings (US$2005/GJ)" = "FE|Buildings|Liquids (EJ/yr)",
"Price|Final Energy|Gases|Buildings (US$2005/GJ)" = "FE|Buildings|Gases (EJ/yr)",
"Price|Final Energy|Solids|Buildings (US$2005/GJ)" = "FE|Buildings|Solids (EJ/yr)",
"Price|Final Energy|Heat|Buildings (US$2005/GJ)" = "FE|Buildings|Heat (EJ/yr)",
"Price|Final Energy|Electricity|Buildings (US$2005/GJ)" = "FE|Buildings|Electricity (EJ/yr)",
"Price|Final Energy|Electricity|Buildings|Moving Avg (US$2005/GJ)" = "FE|Buildings|Electricity (EJ/yr)",
"Price|Final Energy|Hydrogen|Buildings (US$2005/GJ)" = "FE|Buildings|Hydrogen (EJ/yr)",
"Price|Final Energy|Buildings (US$2005/GJ)" = "FE|Buildings (EJ/yr)",
"Price|Final Energy|Buildings|Moving Avg (US$2005/GJ)" = "FE|Buildings (EJ/yr)"
))
}
if (indu_mod %in% c('fixed_shares', 'subsectors')){
int2ext <- c(int2ext,c(
"Price|Final Energy|Heating Oil|Industry (US$2005/GJ)" = "FE|Industry|Liquids (EJ/yr)",
"Price|Final Energy|Gases|Industry (US$2005/GJ)" = "FE|Industry|Gases (EJ/yr)",
"Price|Final Energy|Solids|Industry (US$2005/GJ)" = "FE|Industry|Solids (EJ/yr)",
"Price|Final Energy|Heat|Industry (US$2005/GJ)" = "FE|Industry|Heat (EJ/yr)",
"Price|Final Energy|Electricity|Industry (US$2005/GJ)" = "FE|Industry|Electricity (EJ/yr)",
"Price|Final Energy|Hydrogen|Industry (US$2005/GJ)" = "FE|Industry|Hydrogen (EJ/yr)",
"Price|Final Energy|Industry (US$2005/GJ)" = "FE|Industry (EJ/yr)",
"Price|Final Energy|Industry|Moving Avg (US$2005/GJ)" = "FE|Industry (EJ/yr)"
))
}
if("fegat" %in% getNames(febal.m)){
int2ext <- c(int2ext, "Price|Final Energy|Gases|Transport (US$2005/GJ)" = "FE|Transport|Gases (EJ/yr)")
}
output[is.na(output)] <- 0 # substitute na by 0
# add global prices
map <- data.frame(region=getRegions(tmp),world="GLO",stringsAsFactors=FALSE)
tmp_GLO <- new.magpie("GLO",getYears(tmp),magclass::getNames(tmp),fill=0)
for (i2e in names(int2ext)){
tmp_GLO["GLO",,i2e] <- speed_aggregate(tmp[,,i2e],map,weight=output[map$region,,int2ext[i2e]])
for(t in getYears(tmp)){
if(all(output[map$region,t,int2ext[i2e]]==0)){
tmp_GLO["GLO",t,i2e] <- NA
}
}
}
tmp <- mbind(tmp,tmp_GLO)
# add other region aggregations
if (!is.null(regionSubsetList)){
tmp_RegAgg <- new.magpie(names(regionSubsetList),getYears(tmp),magclass::getNames(tmp),fill=0)
for(region in names(regionSubsetList)){
tmp_RegAgg_ie2 <- do.call("mbind",lapply(names(int2ext), function(i2e) {
map <- data.frame(region=regionSubsetList[[region]],parentRegion=region,stringsAsFactors=FALSE)
result <- speed_aggregate(tmp[regionSubsetList[[region]],,i2e],map,weight=output[regionSubsetList[[region]],,int2ext[i2e]])
getRegions(result) <- region
for(t in getYears(tmp)){
if(all(output[regionSubsetList[[region]],t,int2ext[i2e]]==0)){
result[region,t,i2e] <- NA
}
}
return(result)
}))
tmp_RegAgg[region,,names(int2ext)] <- tmp_RegAgg_ie2[region,,names(int2ext)]
}
tmp <- mbind(tmp,tmp_RegAgg)
}
tmp[is.na(tmp)] <- 0 # tmp is NA if weight is zero for all regions within the GLO or the specific region aggregation. Therefore, we replace all NAs with zeros.
# prices that are same for all regions, including GLO
glob_price <- new.magpie(getRegions(tmp),getYears(tmp),fill=0)
for(i in getRegions(tmp)) glob_price[i,,] <- pm_pvp[,,"peur"]*1000
tmp <- mbind(tmp,setNames(glob_price, "PVP1|Uranium (billionDpktU)"))
for(i in getRegions(tmp)) glob_price[i,,] <- pm_pvp[,,"peoil"]*1000
tmp <- mbind(tmp,setNames(glob_price, "PVP1|Oil (billionDpTWyr)"))
for(i in getRegions(tmp)) glob_price[i,,] <- pm_pvp[,,"pegas"]*1000
tmp <- mbind(tmp,setNames(glob_price, "PVP1|Gas (billionDpTWyr)"))
for(i in getRegions(tmp)) glob_price[i,,] <- pm_pvp[,,"pecoal"]*1000
tmp <- mbind(tmp,setNames(glob_price, "PVP1|Coal (billionDpTWyr)"))
for(i in getRegions(tmp)) glob_price[i,,] <- pm_pvp[,,"pebiolc"]*1000
tmp <- mbind(tmp,setNames(glob_price, "PVP1|Biomass (billionDpTWyr)"))
for(i in getRegions(tmp)) glob_price[i,,] <- pm_pvp[,,"good"]*1000
tmp <- mbind(tmp,setNames(glob_price, "PVP2|Good ()"))
for(i in getRegions(tmp)) glob_price[i,,] <- pm_pvp[,,"perm"]
tmp <- mbind(tmp,setNames(glob_price, "PVP3|Permit ()"))
for(i in getRegions(tmp)) glob_price[i,,] <- pm_pvp[,,"peur"]/pm_pvp[,,"good"] * tdptwyr2dpgj
tmp <- mbind(tmp,setNames(glob_price, "Price|Uranium|World Market (US$2005/GJ)"))
for(i in getRegions(tmp)) glob_price[i,,] <- pm_pvp[,,"peoil"]/pm_pvp[,,"good"] * tdptwyr2dpgj
tmp <- mbind(tmp,setNames(glob_price, "Price|Oil|World Market (US$2005/GJ)"))
for(i in getRegions(tmp)) glob_price[i,,] <- pm_pvp[,,"pegas"]/pm_pvp[,,"good"] * tdptwyr2dpgj
tmp <- mbind(tmp,setNames(glob_price, "Price|Gas|World Market (US$2005/GJ)"))
for(i in getRegions(tmp)) glob_price[i,,] <- pm_pvp[,,"pecoal"]/pm_pvp[,,"good"] * tdptwyr2dpgj
tmp <- mbind(tmp,setNames(glob_price, "Price|Coal|World Market (US$2005/GJ)"))
for(i in getRegions(tmp)) glob_price[i,,] <- pm_pvp[,,"pebiolc"]/pm_pvp[,,"good"] * tdptwyr2dpgj
tmp <- mbind(tmp,setNames(glob_price, "Price|Biomass|World Market (US$2005/GJ)"))
## special global prices
#AJS calc global carbon price as average over regional pm_pvpRegi's, weighted by total emissions.
regi_on_gdx <- unique(readGDX(gdx, name = "regi2iso")[,1])
tmp["GLO",,"Price|Carbon (US$2005/t CO2)"] <-
dimSums( pm_pvpRegi[regi_on_gdx,,"perm"] * output[regi_on_gdx,,"Emi|CO2|FFaI|Emi to which CO2tax is applied (Mt CO2/yr)"],dim=1 ) /
dimSums(output[regi_on_gdx,,"Emi|CO2|FFaI|Emi to which CO2tax is applied (Mt CO2/yr)"],dim=1) /
(pm_pvp[1,,"good"] + 1e-10) * 1000*12/44
# add other region aggregations carbon price as average over regional pm_pvpRegi's, weighted by total emissions.
if (!is.null(regionSubsetList)){
for(region in names(regionSubsetList)){
tmp[region,,"Price|Carbon (US$2005/t CO2)"] <- dimSums( pm_pvpRegi[regionSubsetList[[region]],,"perm"] * output[regionSubsetList[[region]],,"Emi|CO2|FFaI|Emi to which CO2tax is applied (Mt CO2/yr)"],dim=1 ) /
dimSums(output[regionSubsetList[[region]],,"Emi|CO2|FFaI|Emi to which CO2tax is applied (Mt CO2/yr)"],dim=1) /
(pm_pvp[1,,"good"] + 1e-10) * 1000*12/44;
}
}
## not meaningful global prices set to NA
tmp["GLO",,"Price|Biomass|Shiftfactor ()"] <- NA
if (!is.null(regionSubsetList$EUR))
tmp["EUR",,"Price|Carbon|EU-wide Regulation For All Sectors (US$2005/t CO2)"] <- as.vector(co2EUprice[priceReg,,])
## not meaningful region aggregation prices set to NA
if (!is.null(regionSubsetList)){
for(region in names(regionSubsetList)){
tmp[region,,"Price|Biomass|Shiftfactor ()"] <- NA
}
}
# ---- debug information for industry/subsectors ----
if ('subsectors' == indu_mod & !is.null(q37_limit_secondary_steel_share.m)) {
.x <- q37_limit_secondary_steel_share.m[,y,] / budget.m
tmp2 <- mbind(
# fake some GLO data
setNames(
mbind(.x, dimSums(.x * NA, dim = 1)),
'Debug|Industry|Secondary Steel Premium (US$2005)'),
mbind(
lapply(
list(
c('ue_industry', '', 'arbitrary unit', 1),
c('ue_cement', '|Cement', 't cement', 1e3),
c('ue_chemicals', '|Chemicals', 'arbitrary unit', 1),
c('ue_steel', '|Steel', 't Steel', 1e3),
c('ue_steel_primary', '|Steel|Primary', 't Steel', 1e3),
c('ue_steel_secondary', '|Steel|Secondary', 't Steel', 1e3),
c('ue_otherInd', '|other', 'arbitrary unit', 1)),
function(x) {
setNames(
( tmp[,,paste0('Price|CES_input|', x[1], ' ('), pmatch = 'left']
* as.numeric(x[4])
),
paste0('Debug|Price|Industry', x[2], ' (US$2005/', x[3], ')'))
})
)
)
tmp <- mbind(
tmp,
mbind(
tmp2,
calc_regionSubset_sums(tmp2, regionSubsetList)
)
)
}
return(tmp)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.