Correct zero inflated climate model simulations via vine copulas

The multivariate distribution functions of model and observed are estimated via vine copula estimation (see rvinecopulib::vine). This method provides the tool to model the multivariate distribution function of the climate data by Vine Copulas. The climate variables may be zero inflated. The quantiles are then mapped via (inverse) Rosenblatt transformation from the model to the observed distribution. The corrected data are mapped to the projection domain via delta mapping. The steps are equivalent to those in univariate bias correction via quantile delta mapping.

Usage

vbc(
  mp,
  mc,
  rc,
  var_names = colnames(rc),
  margins_controls = list(mult = NULL, xmin = NaN, xmax = NaN, bw = NA, deg = 2, type =
    "c"),
  time_mp = NA,
  ...
)

Arguments

mp

data.table OR list
Simulation data from a climate model during projection period. If this is a list, the algorithm expects each element to be a member of a model ensemble. Each list element is then a data table.

mc

data.table
Simulation data from a climate model during calibration period.

rc

data.table
Historical reference in calibration period.

var_names

character
Names of corrected climate data. Defaults to the column names of the observed data rc.

margins_controls

list
A list with arguments to be passed to kde1d::kde1d(). Currently, there can be

mult numeric vector of length one or d; all bandwidths for marginal kernel density estimation are multiplied with mult. Defaults to log(1 + d) where d is the number of climate variables.
xmin numeric vector of length d; see kde1d::kde1d().
xmax numeric vector of length d; see kde1d::kde1d().
bw numeric vector of length d; see kde1d::kde1d().
deg numeric vector of length one or d; kde1d::kde1d().
type character vector of length one or d; must be one of c, cont, continuous for continuous variables, one of d, disc, discrete for discrete integer variables, or one of zi, zinfl, zero-inflated for zero-inflated variables.

time_mp

numeric
Time vector of the projection period attached to the returned object as a column. Defaults to NA, which means no time vector is attached.

...

Arguments are passed to rvinecopulib::vinecop to specify the structure of vines and margins. Note that the ellipsis of observed and model data are specified with the same arguments.

Value

data.table OR list
The corrected projection period data in mp. Additionally the data frame contains the attributes vine_rc, kde_rc, vine_mp, and kde_mp which store the vine copula and kernel density estimation objects of the observed and model data. If mp is a list, the function returns a list of corrected data frames for each member of the model ensemble.

References

Czado, C. (2019). Analyzing dependent data with vine copulas. Lecture Notes in Statistics, Springer, 222.

Rosenblatt, M. (1952). Remarks on a multivariate transformation. The annals of mathematical statistics, 23(3), 470-472.

Cannon, A.J., S.R. Sobie, and T.Q. Murdock, (2015). Bias correction of simulated precipitation by quantile mapping: How well do methods preserve relative changes in quantiles and extremes? Journal of Climate, 28:6938-6959.

Examples

#' \dontrun{
#' # Start with a simple example:
#' # There is one timeline that should be corrected and projection and 
#' # calibaration period are both in 2010.
#' set.seed(1234L)
#' library(VBC)
#' library(data.table)
#' library(patchwork)
#' data("climate")
#' margins_controls <- list(xmin = c(NaN, 0, NaN, 0, 0),
#'                          type = c("c", "zi", "c", "zi", "c"))
#'                   
#' climate_2010 = lapply(climate, function(data) data[year(time) == 2010, ])
#' 
#' mp_vbc = vbc(climate_2010$mp[, -"time"], climate_2010$mp[, -"time"],
#'              climate_2010$rp[, -"time"], margins_controls = margins_controls,
#'              family_set = "tll", trunc_lvl = Inf)                 
#' 
#' plot_tails(climate_2010$mp, "pr", scale_d = 0.1, mult = 4, xmin = 0) 
#' plot_tails(climate_2010$rp, "pr", scale_d = 1, mult = 4, xmin = 0)
#' plot_tails(round(mp_vbc, 3), "pr", scale_d = 1, mult = 3, xmin = 0)
#' calc_wasserstein(climate_2010$mp[, "pr"], climate_2010$rp[, "pr"])
#' calc_wasserstein(climate_2010$rp[, "pr"], mp_vbc[, "pr"])
#' 
#' calc_wasserstein(climate_2010$rp[, -"time"], climate_2010$mp[, -"time"])
#' calc_wasserstein(climate_2010$rp[, -"time"], mp_vbc)
#' 
#' # An example with temporal subsetting using the method of fragments
#'  
#' temp_subs <- list(DJF = list(hours = seq(0, 21, 3), month = c(12, 1, 2)),
#'                   MAM = list(hours = seq(0, 21, 3), month = 3:5),
#'                   JJA = list(hours = seq(0, 21, 3), month = 6:8),
#'                   SON = list(hours = seq(0, 21, 3), month = 9:11))  
#'                                  
#' mp_vbc = vbc_tsub(climate$mp, climate$mc, climate$rc, t_subs = temp_subs,
#'                   margins_controls = margins_controls, family_set = "tll",
#'                   trunc_lvl = Inf)
#' class(mp_vbc)
#' attr(mp_vbc, "mvd")
#' 
#' measure = lapply(temp_subs, function(sub) {
#'    idx = which(hour(mp_vbc$time) %in% sub$hour &
#'                   month(mp_vbc$time) %in% sub$month)
#'    idxrp = which(hour(climate$rp$time) %in% sub$hour &
#'                      month(climate$rp$time) %in% sub$month)
#'    mci = calc_mci(climate$mp[idx, -"time"], mp_vbc[idx, -"time"],
#'                     time = climate$mp$time[idx])
#'    mci = attr(mci, "global_mci")
#'    
#'    uncorrected = calc_wasserstein(climate$rp[idxrp, -"time"],
#'                                     climate$mp[idx, -"time"])
#'    corrected = calc_wasserstein(climate$rp[idxrp, -"time"],
#'                                 mp_vbc[idx, -"time"])
#'    data.table(
#'        ModelCorrectionInconsistancy = mci,
#'        Wasserstein2_uncorrected = uncorrected[2],
#'        Wasserstein2_corrected = corrected[2],
#'        Improvement_Wasserstein2 = uncorrected[2] - corrected[2], 
#'        month = paste(sub$month, collapse = "-")
#'        )
#' })
#' rbindlist(measure)
#' }