#' Calculate CES Marginals
#'
#' Calculate marginals on the REMIND CES function and combine them to prices.
#'
#' Marginals are calculated analytically
#' \deqn{\frac{\partial V_i}{\partial V_o} = \xi_i (\theta_i \delta_i)^{\rho_o}
#' {V_o}^{1 - \rho_o} {V_i}^{\rho_o - 1}}
#' and prices by recursively applying the chain rule
#' \deqn{\pi_i = \frac{\partial V_i}{\partial V_o} \pi_o
#' \quad \forall (i,o) \in CES}
#'
#' @md
#' @param gdxName Vector of paths to `.gdx` files.
#' @param id If several `.gdx` files are read, an id column is appended to the
#' result; either `file`, with the paths of the originating `.gdx` files,
#' or `scenario`, with the content of `c_expname`.
#'
#' @return A `data frame` with columns `pf` (production factor), `t`, `regi`,
#' `marginal`, `price`, and `file` (path to originating `.gdx` file).
#'
#' @importFrom quitte read.gdx
#' @importFrom dplyr %>% left_join filter sym select rename mutate pull
#' @importFrom tidyr pivot_wider drop_na
#' @importFrom gdxrrw gdxInfo
#' @importFrom rlang is_empty
#' @export
calc_CES_marginals <- function(gdxName, id = 'file') {
if (all(!is.null(id), !id %in% c('file', 'scenario'))) {
warning('id must be either "file" or "scenario". Defaulting to "file".')
id <- 'file'
}
gdxName <- path.expand(gdxName)
.calc_CES_marginals <- function(gdxName, id) {
# ---- read required items from gdx ----
pm_cesdata <- read.gdx(gdxName, 'pm_cesdata', factors = FALSE,
colNames = c('t', 'regi', 'pf', 'param', 'value'))
vm_effGr <- read.gdx(gdxName, 'vm_effGr', factors = FALSE,
colNames = c('t', 'regi', 'pf', 'effGr'))
vm_cesIO <- read.gdx(gdxName, 'vm_cesIO', factors = FALSE,
colNames = c('t', 'regi', 'pf', 'value'))
cesOut2cesIn <- read.gdx(gdxName, 'cesOut2cesIn', factors = FALSE,
colNames = c('pf.out', 'pf.in'))
# ---- calculate marginals ----
marginals <- cesOut2cesIn %>%
left_join(
pm_cesdata %>%
filter(!!sym('param') %in% c('xi', 'eff')) %>%
pivot_wider(names_from = 'param') %>%
drop_na(),
c('pf.in' = 'pf')
) %>%
left_join(
pm_cesdata %>%
filter('rho' == !!sym('param')) %>%
select(-'param', 'rho' = 'value'),
c('t', 'regi', 'pf.out' = 'pf')
) %>%
left_join(
vm_effGr,
c('t', 'regi', 'pf.in' = 'pf')
) %>%
left_join(
vm_cesIO %>%
rename('value.in' = 'value'),
c('t', 'regi', 'pf.in' = 'pf')
) %>%
left_join(
vm_cesIO %>%
rename('value.out' = 'value'),
c('t', 'regi', 'pf.out' = 'pf')
) %>%
mutate(
# ^ !!sym() doesn't work, so use the explicit function call
!!sym('marginal') := !!sym('xi')
* (!!sym('eff') * !!sym('effGr')) ^ (!!sym('rho'))
* `^`(!!sym('value.out'), 1 - !!sym('rho'))
* `^`(!!sym('value.in'), !!sym('rho') - 1))
# ---- calculate prices recursively using the chain rule ----
CES_root <- setdiff(cesOut2cesIn$pf.out, cesOut2cesIn$pf.in)
prices <- marginals %>%
filter(!!sym('pf.out') %in% CES_root) %>%
select('pf' = 'pf.in', 't', 'regi', 'price' = 'marginal')
CES_root <- cesOut2cesIn %>%
filter(!!sym('pf.out') %in% CES_root) %>%
pull('pf.in')
while (!is_empty(CES_root)) {
prices <- bind_rows(
prices,
marginals %>%
filter(!!sym('pf.out') %in% CES_root) %>%
select('pf' = 'pf.in', 't', 'regi', 'price' = 'marginal') %>%
left_join(cesOut2cesIn, c('pf' = 'pf.in')) %>%
left_join(
prices %>%
rename('price.out' = 'price'),
c('t', 'regi', 'pf.out' = 'pf')
) %>%
mutate(!!sym('price') := !!sym('price') * !!sym('price.out')) %>%
select('pf', 't', 'regi', 'price')
)
CES_root <- cesOut2cesIn %>%
filter(!!sym('pf.out') %in% CES_root) %>%
pull('pf.in')
}
# ---- bind marginals and prices together ----
r <- bind_rows(
marginals %>%
select('pf' = 'pf.in', 't', 'regi', 'value' = 'marginal') %>%
mutate(!!sym('name') := 'marginal'),
prices %>%
rename('value' = 'price') %>%
mutate(!!sym('name') := 'price')
) %>%
pivot_wider()
if (id) {
r <- r %>%
mutate(!!sym('scenario') := read.gdx(gdxName, 'c_expname',
colNames = 'c_expname',
factors = FALSE) %>%
pull('c_expname'))
}
return(r)
}
# ---- turn gdxInfo errors on gdx files into warnings and remove paths ----
for (p in gdxName) {
gdxName <- tryCatch(
expr = {
gdxInfo(p, dump = FALSE)
gdxName
},
error = function(e) {
warning(e)
return(setdiff(gdxName, p))
}
)
}
# ---- bind results for all valid input files together ----
r <- bind_rows(
lapply(gdxName, function(gdxName) {
.calc_CES_marginals(gdxName, id = all(!is.null(id), id == 'scenario')) %>%
mutate(file = gdxName)
})
)
if (any(is.null(id), 'file' != id)) {
r <- r %>%
select(-file)
}
return(r)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.