Merge pull request #9 from deleip/usdaUpdate

Usda update

.buildlibrary

@@ -1,8 +1,9 @@
-ValidationKey: '1200840'
+ValidationKey: '1422201'
 AutocreateReadme: yes
 AcceptedWarnings:
 - 'Warning: package ''.*'' was built under R version'
 - 'Warning: namespace ''.*'' is not available and has been replaced'
 AcceptedNotes: .* includes the non-default packages.*
 allowLinterWarnings: no
 enforceVersionUpdate: no
+skipCoverage: no

.github/workflows/check.yaml

@@ -63,6 +63,6 @@ jobs:
         shell: Rscript {0}
         run: |
           nonDummyTests <- setdiff(list.files("./tests/testthat/"), c("test-dummy.R", "_snaps"))
-          if(length(nonDummyTests) > 0) covr::codecov(quiet = FALSE)
+          if(length(nonDummyTests) > 0 && !lucode2:::loadBuildLibraryConfig()[["skipCoverage"]]) covr::codecov(quiet = FALSE)
         env:
           NOT_CRAN: "true"

CITATION.cff

@@ -2,8 +2,8 @@ cff-version: 1.2.0
 message: If you use this software, please cite it using the metadata from this file.
 type: software
 title: 'mrfactors: MADRaT based package on factor inputs'
-version: 0.6.0
-date-released: '2024-10-18'
+version: 0.7.1
+date-released: '2024-11-04'
 abstract: This package provides functions for MAgPIE input data on factor inputs to
   agricultural production (with a focus on capital and labor).
 authors:

DESCRIPTION

@@ -1,8 +1,8 @@
 Type: Package
 Package: mrfactors
 Title: MADRaT based package on factor inputs
-Version: 0.6.0
-Date: 2024-10-18
+Version: 0.7.1
+Date: 2024-11-04
 Authors@R:
     c(person("Debbora", "Leip", , "leip@pik-potsdam.de", role = c("aut", "cre")),
       person("Edna", "Molina Bacca", role = "aut"))

R/calcAgCapLabourShare.R

@@ -3,8 +3,12 @@
 #' @description This function calculates historical capital shares (Capital + Labour)
 #' of the factor requirements using USDA
 #'
+#' @param fillWithRegression boolean: should missing values be filled based on a regression between capital share
+#' and GDPpcPPP (calibrated to countries)
+#' @param projection either FALSE or SSP on which projections should be based. Only relevant if fillWithRegression is
+#' TRUE.
 #' @return MAgPIE object
-#' @author Edna J. Molina Bacca
+#' @author Edna J. Molina Bacca, Debbora Leip
 #' @seealso [calcAgCapLabourShare()]
 #'
 #' @examples
@@ -14,38 +18,91 @@
 #' @importFrom madrat calcOutput
 #' @importFrom magclass collapseDim collapseNames dimSums getYears
 #'
-calcAgCapLabourShare <- function() {
-  # Reads requirements shares for capital and labour based on USDA
-  fraction <- collapseNames(calcOutput("FractionInputsUSDA", aggregate = FALSE))[, , c("Capital", "Labor")]
+calcAgCapLabourShare <- function(fillWithRegression = TRUE, projection = FALSE) {
+  # raw USDA cost shares -- labor-capital ratio is the same for kcr and kli, as both are shared inputs
+  # We use kcr here, but use full factor costs (crops + livst) as aggregation weight below
+  shares <- calcOutput("FractionInputsUSDA", products = "kcr", aggregate = FALSE,
+                       keepConstantExtrapolation = FALSE, interpolate = TRUE)
 
   # Determines fraction between capital and labour
-  fractionCapital <- fraction[, , "Capital"] / (dimSums(fraction, dim = 3.1))
+  fractionCapital <- shares[, , "Capital"] / (shares[, , "Labor", drop = TRUE] + shares[, , "Capital", drop = TRUE])
+  fractionCapital[is.na(fractionCapital)] <- 0
 
-  # In case of values different to finite makes them 0
-  fractionCapital[!is.finite(fractionCapital)] <- 0
+  # drop mexico case study due to unreasonably high capital share
+  mapping <- toolGetMapping("caseStudiesUSDATFP.csv", where = "mrfactors", type = "regional")
+  countries <- mapping$ISO[mapping$CaseStudiesUsed == "MEX"]
+  fractionCapital[countries, , ] <- 0
+
+  # fill gaps with regression
+  if (isTRUE(fillWithRegression)) {
+
+    regCoeff <- calcOutput("RegFactorShare", aggregate = FALSE)
+
+    # calculate GDPpc [USD2017MER] for regression
+    gdp <- calcOutput("GDPpc", naming = "scenario", unit = "constant 2017 Int$PPP", aggregate = FALSE)
+    if (!isFALSE(projection)) {
+      years <- setdiff(getItems(gdp, dim = 2), paste0("y", seq(2105, 2150, 5)))
+      gdp <- gdp[, years, projection]
+    } else {
+      gdp <- gdp[, getItems(fractionCapital, dim = 2), "SSP2"]
+    }
+
+    # add years with GDP data to hourlyCosts object
+    fractionCapital <- magpiesort(add_columns(fractionCapital, dim = 2, fill = 0,
+                                              addnm = setdiff(getItems(gdp, dim = 2),
+                                                              getItems(fractionCapital, dim = 2))))
+    years <- getYears(fractionCapital, as.integer = TRUE)
+
+    .fillTimeseries <- function(var, country, regressor) {
+      ctryEst <- regCoeff[, , "intercept", drop = TRUE] + regCoeff[, , "slope", drop = TRUE] * regressor[ctry, , ]
+      y <- where(var[ctry, , ] != 0)$true$years
+      if (length(y) == 0) {
+        var[ctry, , ] <- ctryEst
+      } else {
+        y <- as.integer(str_split(y, "y", simplify = TRUE)[, 2])
+        yPast <- years[years < min(y)]
+        yFuture <- years[years > max(y)]
+        yGaps <- setdiff(years[(years >= min(y)) & (years <= max(y))], y)
+
+        if (length(yGaps) > 0) {
+          var[ctry, sort(c(y, yGaps)), ] <- time_interpolate(dataset = var[ctry, y, ],
+                                                             interpolated_year = yGaps,
+                                                             integrate_interpolated_years = TRUE)
+        }
+
+        if (length(yPast) > 0) {
+          calibPast <- var[ctry, min(y), ] - ctryEst[, min(y), ]
+          var[ctry, yPast, ] <- ctryEst[, yPast, ] + calibPast
+        }
+
+        if (length(yFuture) > 0) {
+          calibFuture <- var[ctry, max(y), ] - ctryEst[, max(y), ]
+          var[ctry, yFuture, ] <- ctryEst[, yFuture, ] + calibFuture
+        }
+      }
+      return(var)
+    }
+
+    logGdp <- log(gdp, base = 10)
+
+    for (ctry in getItems(fractionCapital, dim = 1)) {
+      fractionCapital <- .fillTimeseries(var = fractionCapital, country = ctry, regressor = logGdp)
+    }
+  }
 
   # factor costs as weight
   factorCostsCrops <- dimSums(calcOutput("FactorCostsCrops", aggregate = FALSE), dim = 3.1)
   factorCostsLivst <- dimSums(calcOutput("FactorCostsLivst", aggregate = FALSE), dim = 3.1)
   weight <- factorCostsCrops + factorCostsLivst
 
-  # add missing years to weight
+  # keep weight constant for missing years in the past and projection
   missingYears <- setdiff(getYears(fractionCapital, as.integer = TRUE), getYears(weight, as.integer = TRUE))
-  if (any(missingYears > min(getYears(weight, as.integer = TRUE)))) {
-    stop("Need to fix weight for new years not covered in factor costs dataset")
-  }
-  minYear <- min(getYears(weight, as.integer = TRUE))
-  weight <- magpiesort(add_columns(weight, dim = 2, addnm = paste0("y", missingYears)))
-  weight[, missingYears, ] <- weight[, minYear, ]
-
-  # Give 0 weigh to countries with unexpectedly high capital shares
-  weight[c("BLZ", "CRI", "DOM", "HND", "JAM", "MEX", "NIC", "PAN", "SLV"), , ] <- 0
-
-  years <- intersect(getYears(weight), getYears(fractionCapital))
+  weight <- time_interpolate(dataset = weight, interpolated_year = missingYears,
+                             integrate_interpolated_years = TRUE, extrapolation_type = "constant")
+  weight <- weight[, getItems(fractionCapital, dim = 2), ]
 
-  weight <- weight[, years, ]
-  weight[fractionCapital[, years, ] == 0] <- 0
-  x <- setNames(fractionCapital[, years, ], NULL)
+  weight[fractionCapital == 0] <- 0
+  x <- setNames(fractionCapital, NULL)
 
   return(list(x = x,
               weight = weight,

R/calcFacReq.R

@@ -1,6 +1,6 @@
 #' @title calcFacReq
 #'
-#' @description This function calculates factor requirement costs (05USDMER/tDM) on regional level,
+#' @description This function calculates factor requirement costs (constant 2017 US$MER/tDM) on regional level,
 #' based FAO databases
 #'
 #' @param splitSectors if TRUE requirements for labor and capital will be reported separately
@@ -18,7 +18,7 @@
 #'
 calcFacReq <- function(splitSectors = FALSE) {
 
-  # get factor requirements in US$PPP05
+  # get factor requirements in 2017 US$MER/tDM
   facReq <- calcOutput("FactorIntensity", aggregate = FALSE)
   facReq <- facReq[, where(facReq != 0)$true$years, ]
 

R/calcFacReqGLO.R

@@ -18,7 +18,7 @@ calcFacReqGLO <- function() {
 
   .Deprecated(paste0("calcOutput", "(\"FacReq\", aggregate = \"GLO\", years = 2005)"))
 
-  # Reads factor requirements in USD05ppp
+  # Reads factor requirements in constant 2017 US$MER
   factors <- dimSums(calcOutput("FactorIntensity", aggregate = FALSE, years = 2005), dim = 3.1)
 
   # Reads production as weight

R/calcFactorCosts.R

@@ -0,0 +1,36 @@
+#' @title calcFactorCosts
+#' @description combines factor costs for crop and livestock production in mio. US$MER05
+#' @param datasource only source available is "USDA" (calculates factor costs by applying factor cost share from USDA
+#' to VoP from FAO)
+#' @param otherLivst boolean: should FAO livestock categories that can't be matched to MAgPIE categories (i.e. beeswax,
+#' wool, silkworms, and honey) be reported as "livst_other"?
+#' @param subsectors boolean: should output be aggregated or reported by crop and livestock subsectors?
+#' @param unit output currency unit based on the toolConvertGDP function from the  GDPuc library
+#' @return List of magpie objects with results on country level, weight on country level, unit and description.
+#' @author Debbora Leip
+#' @examples
+#' \dontrun{
+#' calcOutput("FactorCosts")
+#' }
+#' @importFrom magclass setNames dimSums time_interpolate
+
+calcFactorCosts <- function(datasource = "USDA", otherLivst = FALSE, subsectors = FALSE,
+                            unit = "constant 2017 US$MER") {
+
+  facCostsCrops <- calcOutput("FactorCostsCrops", datasource = datasource, unit = unit, aggregate = FALSE)
+  facCostsLivst <- calcOutput("FactorCostsLivst", datasource = datasource, otherLivst = otherLivst,
+                              unit = unit, aggregate = FALSE)
+
+  factorCosts <- mbind(facCostsCrops, facCostsLivst)
+
+  if (isFALSE(subsectors)) {
+    factorCosts <- dimSums(factorCosts, dim = 3)
+  }
+
+  units <- paste0("mio ", unit)
+
+  return(list(x = factorCosts,
+              weight = NULL,
+              unit = units,
+              description = "Factor costs for crop and livestock production"))
+}

R/calcFactorCostsCrops.R

@@ -14,7 +14,7 @@
 calcFactorCostsCrops <- function(datasource = "USDA", unit = "constant 2017 US$MER") {
 
   if (datasource == "USDA") {
-    # Value of Production for livestock in US$MER2017 (including FAO livst categories not mapped to MAgPIE categories)
+    # Value of Production for livestock in US$MER2017
     vopCrops <- calcOutput("VoPcrops", fillGaps = TRUE, aggregate = FALSE, unit = "constant 2017 US$MER")
 
     # no VoP data before 1991, data for 2019 incomplete
@@ -36,12 +36,6 @@ calcFactorCostsCrops <- function(datasource = "USDA", unit = "constant 2017 US$M
     shares <- toolFillWithRegionAvg(shares[, y, ], valueToReplace = 0, weight = weight,
                                     regionmapping = h12, verbose = FALSE, warningThreshold = 1)
 
-    # for REF in 1990 no country has a value, so toolFillWithRegionAvg assigns NA. Use values from 1995 instead:
-    if ("y1990" %in% y) { # subsidy data starts only in 2005
-      ref <- h12$CountryCode[h12$RegionCode == "REF"]
-      shares[ref, 1990, ]  <- shares[ref, 1995, ]
-    }
-
     # interpolate between the five-year-steps
     shares <- time_interpolate(shares,
                                interpolated_year = setdiff(years, getYears(shares, as.integer = TRUE)),

R/calcFactorCostsLivst.R

@@ -4,6 +4,7 @@
 #' to VoP from FAO)
 #' @param otherLivst boolean: should FAO livestock categories that can't be matched to MAgPIE categories (i.e. beeswax,
 #' wool, silkworms, and honey) be reported as "livst_other"?
+#' @param inclFish boolean: should fish factor costs be included?
 #' @param unit output currency unit based on the toolConvertGDP function from the  GDPuc library
 #' @return List of magpie objects with results on country level, weight on country level, unit and description.
 #' @author Debbora Leip
@@ -13,18 +14,26 @@
 #' }
 #' @importFrom magclass setNames dimSums time_interpolate
 
-calcFactorCostsLivst <- function(datasource = "USDA", otherLivst = FALSE, unit = "constant 2017 US$MER") {
+calcFactorCostsLivst <- function(datasource = "USDA", otherLivst = FALSE, inclFish = FALSE,
+                                 unit = "constant 2017 US$MER") {
 
   if (datasource == "USDA") {
-    # Value of Production for livestock in US$MER2005 (including FAO livst categories not mapped to MAgPIE categories)
+    # Value of Production for livestock in US$MER2017 (including FAO livst categories not mapped to MAgPIE categories)
     vopLivst <- calcOutput("VoPlivst", fillGaps = TRUE, aggregate = FALSE, other = otherLivst,
                            unit = "constant 2017 US$MER")
 
     # no VoP data before 1991, data for 2019 incomplete
     years <- setdiff(getYears(vopLivst, as.integer = TRUE), c(1960:1990, 2019))
 
+    # VoP of fish (reduces years) -- will use livestock cost shares also for fish
+    if (isTRUE(inclFish)) {
+      vopFish <- calcOutput("VoPAFF", aggregate = FALSE)[, , "Fisheries"]
+      years <- intersect(years, getYears(vopFish, as.integer = TRUE))
+      vopLivst <- mbind(vopLivst[, years, ], vopFish[, years, ])
+    }
+
     # USDA labor cost shares
-    shares <- calcOutput("FractionInputsUSDA", products = "kli", aggregate = FALSE)
+    shares <- calcOutput("FractionInputsUSDA", products = "kli", aggregate = FALSE, keepConstantExtrapolation = TRUE)
     shares <- dimSums(shares[, , c("Labor", "Capital")], dim = 3)
 
     # closest 5-year step before and after start of VoP data needed for interpolation of shares
@@ -39,12 +48,6 @@ calcFactorCostsLivst <- function(datasource = "USDA", otherLivst = FALSE, unit =
     shares <- toolFillWithRegionAvg(shares[, y, ], valueToReplace = 0, weight = weight,
                                     regionmapping = h12, verbose = FALSE, warningThreshold = 1)
 
-    # for REF in 1990 no country has a value, so toolFillWithRegionAvg assigns NA. Use values from 1995 instead:
-    if ("y1990" %in% y) { # subsidy data starts only in 2005
-      ref <- h12$CountryCode[h12$RegionCode == "REF"]
-      shares[ref, 1990, ]  <- shares[ref, 1995, ]
-    }
-
     # interpolate between the five-year-steps
     shares <- time_interpolate(shares,
                                interpolated_year = setdiff(years, getYears(shares, as.integer = TRUE)),

R/calcFactorIntensity.R

@@ -1,6 +1,6 @@
 #' @title calcFactorIntensity
-#' @description Calculates factor intensity for labour and/or capital from USDA (Inputs share)
-#' and FAO (Value of Production).
+#' @description Calculates factor intensity in crop production for labour and/or capital from USDA (Inputs share)
+#' and FAO (Value of Production)in constant 2017 US$MER per ton.
 #' Capital intensity and requirements can also be calculated from FAO's CapitalStock database.
 #'
 #'
@@ -142,6 +142,8 @@ calcFactorIntensity <- function(output = "intensities", method = "USDA", unit =
                         replace_NAs = "no_conversion")
   }
 
+  weight <- weight + 1e-20
+
   return(list(x = x,
               weight = weight,
               mixed_aggregation = NULL,