derive an M matrix for IEF models

R/ExternalImportFactors.R

@@ -100,7 +100,7 @@ castImportFactors <- function(IFTable, model) {
 buildModelwithImportFactors <- function(model, configpaths = NULL) {
   # (see Palm et al. 2019)
 
-  logging::loginfo("Building Import A (A_m) accounting for ITA in Domestic FD.\n")
+  logging::loginfo("Building A_m (import requirements) accounting for international trade adjustment in domestic final demand.\n")
   # Re-derive import values in Use and final demand
   # _m denotes import-related structures
   model$UseTransactions_m <- model$UseTransactions - model$DomesticUseTransactions
@@ -126,5 +126,28 @@ buildModelwithImportFactors <- function(model, configpaths = NULL) {
 
   model$M_m <- M_m
 
+  model$M <- deriveMMatrix(model)
+
   return(model)
 }
+
+#' Derives an aggregate M matrix from M_d and M_m based on the ratio of commodity output to total imports.
+#' @param model, An EEIO model object with model specs and crosswalk table loaded
+#' @return An M matrix of flows x sector
+deriveMMatrix <- function(model) {
+  logging::loginfo("Deriving M matrix (total emissions and resource use per dollar) ...")
+  q <- model$q
+  import_code <- model$FinalDemandMeta[model$FinalDemandMeta$Group=="Import", "Code_Loc"]
+  # derive total imports (m) from the Use table
+  U_m <- model$U - model$U_d
+  # Exclude imports col when calculating total imports
+  m <- rowSums(U_m[model$Commodities$Code_Loc, !(colnames(U_m) %in% import_code)]) # drop VA
+
+  dr <- q / (q + m)
+  mr <- 1 - dr
+  # Derive M by taking the ratio of domestic vs imported goods
+  M <- model$M_d %*% diag(as.vector(dr)) + model$M_m %*% diag(as.vector(mr))
+  colnames(M) <- colnames(model$M_d)
+
+  return(M)
+}

R/ValidateModel.R

@@ -474,6 +474,19 @@ validateImportFactorsApproach <- function(model, demand = "Consumption"){
   cat("assuming model$M = model$M_m.\n")
   print(all.equal(LCIA_dm, LCIA))
 
+  cat("\nValidating that the derived M matrix has all values between M_d and M_m\n")
+
+  # Validate that M is between M_m and M_d
+  a <- signif(model$M_m, 6)
+  b <- signif(model$M, 6)
+  c <- signif(model$M_d, 6)
+  z <- ((b > a) & (b > c)) | ((b < a) & (b < c))
+  if(sum(z) == 0) {
+    print(TRUE)
+  } else {
+    comm <- colSums(z) > 0
+    print(paste0("Failures: ", names(comm)[comm == TRUE]))
+  }
 }
 
 #' Validate the calculation of household_emissions