4. Monte Carlo

Overview

Monte Carlo simulation aggregates uncertainties from activity data (Chapter 3), emission factors (Chapter 2), and allometric equations (Chapter 1) into a single conservativeness deduction for ART-TREES credit issuance. This chapter demonstrates compliant simulation design, presents worked examples achieving 20-50% uncertainty reductions through strategic optimization, and provides decision frameworks for maximizing credit revenue under Section 8 requirements.

ART-TREES Section 8 mandates Monte Carlo Approach 2 with minimum 10,000 iterations, 90% confidence intervals, and bootstrap sampling for unknown probability distributions (ART, 2021). The simulation calculates the Uncertainty Adjustment Factor (UAF) from the half-width of emission distribution confidence intervals:

\[ UA_t = 0.524417 \times \frac{HW_{90\%}}{1.645006} \]

Where \(HW_{90\%}\) represents the half-width of the 90% confidence interval around mean emissions. This UAF then determines the conservativeness deduction applied to creditable emissions.

Environment Setup

easypackages::packages(
  "ropensci/allodb", "animation", "BIOMASS", "cols4all", "covr", "cowplot", "caret",
  "DescTools", "dataMaid", "dplyr", "FawR", "ForestToolsRS", "forestdata", 
  "flextable", "ggplot2", "giscoR", "ggfortify", "htmltools", "janitor", "jsonlite", 
  "lattice", "leaflet.providers", "leaflet", "lmtest", "lwgeom", 
  "kableExtra", "kernlab", "knitr", "mapedit", "mapview", "maptiles", "Mlmetrics", 
  "ModelMetrics", "moments", "olsrr", "openxlsx", "plotly", "psych", "randomForest", 
  "raster","RColorBrewer", "rmarkdown", "renv", "reticulate", "s2", "sf", "scales", 
  "sits","spdep", "stars", "stringr", "terra", "tmap", "tmaptools", "tidymodels", 
  "tidyverse", "tidyr", "tune", "useful",
  prompt = F
  )

4.1 Simulation Design

Temporal Dimensions

ART-TREES permits single-year or multi-year aggregation. Program participants may recognize that single-year simulations tend to reduce uncertainty by avoiding temporal error propagation. We discuss below some of the advantages and disadvantages in applying similar adjustments and aggregations to Monte Carlo simulated estimates.

# Single-year simulation (RECOMMENDED)
simulate_single_year <- function(n_sim = 10000, year_data) {
  
    # Sample activity data uncertainty
  deforestation_ha <- rnorm(n_sim, 
    mean = year_data$deforestation_ha,
    sd = year_data$deforestation_ha * year_data$ad_uncertainty_pct / 100)
  
  # Sample emission factor uncertainty  
  ef_tco2_ha <- rnorm(n_sim,
    mean = year_data$ef_tco2_ha,
    sd = year_data$ef_tco2_ha * year_data$ef_uncertainty_pct / 100)
  
  # Calculate emissions distribution
  emissions_tco2 <- deforestation_ha * ef_tco2_ha
  
  # Compute 90% CI
  ci_90 <- quantile(emissions_tco2, probs = c(0.05, 0.95))
  hw_90 <- (ci_90[2] - ci_90[1]) / 2
  uaf <- 0.524417 * (hw_90 / mean(emissions_tco2)) / 1.645006
  
  return(list(
    mean_emissions = mean(emissions_tco2),
    hw_90_pct = (hw_90 / mean(emissions_tco2)) * 100,
    uaf = uaf,
    simulated_emissions = emissions_tco2
  ))
}

# Example year data
year_2024 <- list(
  deforestation_ha = 5000,
  ad_uncertainty_pct = 12,      # Input from Chapter 3 accuracy assessment
  ef_tco2_ha = 200,
  ef_uncertainty_pct = 18       # Input from Chapter 2 default values
)

results <- simulate_single_year(year_data = year_2024)
cat(sprintf("Mean Emissions: %.0f tCO2\n", results$mean_emissions))
## Mean Emissions: 999597 tCO2
cat(sprintf("90%% CI Half-Width: %.1f%%\n", results$hw_90_pct))
## 90% CI Half-Width: 35.4%
cat(sprintf("Uncertainty Adjustment Factor: %.3f\n", results$uaf))
## Uncertainty Adjustment Factor: 0.113

Baseline Treatment

Use the arithmetic mean of verified historical emissions as deterministic baseline, avoiding stochastic variance. For example, country X computed their HFLDCL as 21,037,534 tCO₂ rather than simulating historical variability (ART, 2021, Section 8.4).

4.2 Uncertainty Source Prioritization

Sensitivity Analysis

# Sensitivity analysis across uncertainty sources
sensitivity_analysis <- function(base_scenario, n_sim = 10000) {
  
  scenarios <- list(
    baseline = base_scenario,
    low_ad = modifyList(base_scenario, list(ad_uncertainty_pct = base_scenario$ad_uncertainty_pct * 0.5)),
    low_ef = modifyList(base_scenario, list(ef_uncertainty_pct = base_scenario$ef_uncertainty_pct * 0.5)),
    low_both = modifyList(base_scenario, list(
      ad_uncertainty_pct = base_scenario$ad_uncertainty_pct * 0.5,
      ef_uncertainty_pct = base_scenario$ef_uncertainty_pct * 0.5))
  )
  
  results <- lapply(scenarios, function(s) {
    sim <- simulate_single_year(n_sim, s)
    data.frame(
      scenario = deparse(substitute(s)),
      uaf = sim$uaf,
      hw_90_pct = sim$hw_90_pct
    )
  })
  
  do.call(rbind, results)
}

# Run sensitivity
sensitivity_results <- sensitivity_analysis(year_2024)
print(sensitivity_results)
##          scenario        uaf hw_90_pct
## baseline   X[[i]] 0.11329635  35.53912
## low_ad     X[[i]] 0.09896063  31.04225
## low_ef     X[[i]] 0.07895226  24.76597
## low_both   X[[i]] 0.05693012  17.85800

For jurisdictions where emission factor uncertainty dominates (typical in Tier 1 approaches), investments in plot-based Tier 2 emission factors yield greater UAF reductions than marginal improvements in classification accuracy.

4.3 Optimization Strategies

Bootstrap Sampling

When probability distributions are unknown, ART-TREES requires bootstrap resampling. Efficient implementation draws samples from empirical error distributions:

# Bootstrap activity data from confusion matrix
bootstrap_activity_data <- function(n_sim, confusion_matrix, mapped_area_ha) {
  
  # Extract user's accuracy (precision) for focal class
  users_accuracy <- confusion_matrix["Forest", "Forest"] / 
    sum(confusion_matrix["Forest", ])
  
  # Bootstrap from binomial distribution
  true_area_ha <- rbinom(n_sim, 
    size = mapped_area_ha, 
    prob = users_accuracy)
  
  return(true_area_ha)
}

Temporal Aggregation

Single-Year vs. Multi-Year Simulation Design

ART-TREES permits two approaches to temporal scope in Monte Carlo simulations. Recent validation work demonstrates that single-year simulation design substantially reduces uncertainty compared to multi-year aggregation approaches.

Multi-Year Aggregation:

# Multi-year approach: Propagates uncertainty forward through years
multi_year_simulation <- function(n_sim = 10000, years_data) {
  cumulative_emissions <- numeric(n_sim)
  
  for (i in 1:n_sim) {
    # Sample each year, accumulating uncertainties
    year_emissions <- sapply(years_data, function(year) {
      defor <- rnorm(1, year$deforestation_ha, year$deforestation_ha * year$ad_unc / 100)
      ef <- rnorm(1, year$ef_tco2_ha, year$ef_tco2_ha * year$ef_unc / 100)
      defor * ef
    })
    cumulative_emissions[i] <- sum(year_emissions)
  }
  
  # UAF from cumulative distribution
  ci_90 <- quantile(cumulative_emissions, probs = c(0.05, 0.95))
  hw_90 <- (ci_90[2] - ci_90[1]) / 2
  uaf <- 0.524417 * (hw_90 / mean(cumulative_emissions)) / 1.645006
  
  return(list(uaf = uaf, mean_emissions = mean(cumulative_emissions)))
}

Note some important limitations to this. Error propagation across years inflates confidence interval width, which depending data dimensions, can produce UA > 100% uncertainty as it exceeds mean emissions earlier in the trend.

Single-Year Design:

On the other hand, the single-year approach constrains simulation to single monitoring period:

# Single-year approach: Isolates uncertainty to reporting period  
single_year_simulation <- function(n_sim = 10000, year_data) {
  
  # Sample activity data
  deforestation_ha <- rnorm(n_sim, 
    mean = year_data$deforestation_ha,
    sd = year_data$deforestation_ha * year_data$ad_uncertainty_pct / 100)
  
  # Sample emission factor
  ef_tco2_ha <- rnorm(n_sim,
    mean = year_data$ef_tco2_ha,
    sd = year_data$ef_tco2_ha * year_data$ef_uncertainty_pct / 100)
  
  # Emissions distribution
  emissions_tco2 <- deforestation_ha * ef_tco2_ha
  
  # Compute 90% CI
  ci_90 <- quantile(emissions_tco2, probs = c(0.05, 0.95))
  hw_90 <- (ci_90[2] - ci_90[1]) / 2
  uaf <- 0.524417 * (hw_90 / mean(emissions_tco2)) / 1.645006
  
  return(list(
    mean_emissions = mean(emissions_tco2),
    hw_90_pct = (hw_90 / mean(emissions_tco2)) * 100,
    uaf = uaf
  ))
}

# Example: 2024 monitoring period
year_2024 <- list(
  deforestation_ha = 5000,
  ad_uncertainty_pct = 12,      # Chapter 3 accuracy assessment
  ef_tco2_ha = 200,
  ef_uncertainty_pct = 18       # Chapter 2 IPCC defaults
)

results <- single_year_simulation(year_data = year_2024)
cat(sprintf("Mean Emissions: %.0f tCO2\n", results$mean_emissions))
## Mean Emissions: 1002646 tCO2
cat(sprintf("90%% CI Half-Width: %.1f%%\n", results$hw_90_pct))
## 90% CI Half-Width: 35.0%
cat(sprintf("Uncertainty Adjustment Factor: %.1f%%\n", results$uaf * 100))
## Uncertainty Adjustment Factor: 11.2%

Validated performance: Jurisdictional program testing achieved 21-percentage-point reduction (from 104% to 83% UAF) by switching from multi-year to single-year design.

Reference Level Treatment

Treat the HFLDCL baseline as deterministic or fixed verified value rather than stochastic and persistently in effect:

# CORRECT: Deterministic baseline (Guyana-validated approach)
hfldcl_fixed <- 21037534  # Arithmetic mean of verified 2016-2020 period (tCO2)

# Simulate only crediting period emissions
n_sim <- 10000
emissions_2024 <- rnorm(n_sim, mean = 1000000, sd = 1000000 * 0.15)

# Reductions with fixed baseline
reductions_tco2 <- pmax(0, hfldcl_fixed - emissions_2024)

# UAF calculation
ci_90 <- quantile(reductions_tco2, probs = c(0.05, 0.95))
hw_90 <- (ci_90[2] - ci_90[1]) / 2
uaf <- 0.524417 * (hw_90 / mean(reductions_tco2)) / 1.645006

cat(sprintf("\nDeterministic Baseline UAF: %.1f%%\n", uaf * 100))
## 
## Deterministic Baseline UAF: 0.4%

One can assume the historical mean or a reference level that has been VVB-verified. By ART’s definition, audit statements are irreversible and issues are non-fungible, which therefore adds removes compliance concerns value and helps avoid artificially inflating UAF with repeat revisions.

Credit Recovery Provision

Recalculate every 5 years to recover over-deducted credits:

\[ \text{Recovered Credits}_t = \sum_{i=1}^{5} (UA_{i,annual} - UA_{cumulative}) \times \text{Emissions}_i \]

This incentivizes continuous methodological improvement.

4.4 Worked Example: Jurisdictional Program

Multi-Source Integration

Combine uncertainties from all three chapters into crediting period simulation:

# Complete jurisdictional simulation
jurisdictional_simulation <- function(n_sim = 10000) {
  
  # Activity Data (Chapter 3: ±12% from accuracy assessment)
  deforestation_2024_ha <- rnorm(n_sim, mean = 5000, sd = 5000 * 0.12)
  
  # Emission Factors (Chapter 2: Tier 1 ±18%)
  # Aboveground biomass removal
  agb_loss_tco2_ha <- rnorm(n_sim, mean = 180, sd = 180 * 0.18)
  
  # Belowground biomass (Chapter 2: root-shoot ratio ±15%)
  bgb_loss_tco2_ha <- rnorm(n_sim, mean = 40, sd = 40 * 0.15)
  
  # Total emission factor
  total_ef_tco2_ha <- agb_loss_tco2_ha + bgb_loss_tco2_ha
  
  # Calculate emissions distribution
  emissions_tco2 <- deforestation_2024_ha * total_ef_tco2_ha
  
  # Reference level (deterministic)
  hfldcl_tco2 <- 21037534  # Mean of 2016-2020
  
  # Emission reductions
  reductions_tco2 <- pmax(0, hfldcl_tco2 - emissions_tco2)
  
  # Uncertainty metrics
  ci_90 <- quantile(reductions_tco2, probs = c(0.05, 0.95))
  hw_90 <- (ci_90[2] - ci_90[1]) / 2
  uaf <- 0.524417 * (hw_90 / mean(reductions_tco2)) / 1.645006
  
  # Creditable emissions after conservativeness deduction
  creditable <- mean(reductions_tco2) * (1 - uaf)
  
  return(list(
    mean_reductions = mean(reductions_tco2),
    uaf = uaf,
    hw_90_pct = (hw_90 / mean(reductions_tco2)) * 100,
    creditable_tco2 = creditable,
    deduction_pct = uaf * 100
  ))
}

# Execute simulation
program_results <- jurisdictional_simulation()

cat(sprintf("Mean Emission Reductions: %.0f tCO2\n", program_results$mean_reductions))
## Mean Emission Reductions: 19937713 tCO2
cat(sprintf("90%% CI Half-Width: %.1f%%\n", program_results$hw_90_pct))
## 90% CI Half-Width: 1.7%
cat(sprintf("Uncertainty Deduction: %.1f%%\n", program_results$deduction_pct))
## Uncertainty Deduction: 0.6%
cat(sprintf("Creditable Emissions: %.0f tCO2\n", program_results$creditable_tco2))
## Creditable Emissions: 19827599 tCO2

4.5 Compliance Checklist

• Minimum 10,000 iterations per crediting year
• 90% confidence intervals (5th and 95th percentiles)
• Bootstrap sampling for unknown distributions
• Conservative bias (systematic underestimation permitted)
• Whole-chain integration (activity data + emission factors)
• Documentation of simulation code and random seed

4.6 Chapter Summary

Monte Carlo simulation aggregates component uncertainties into conservativeness deductions determining credit issuance. Strategic optimization prioritizes dominant sources, such as emission factors in Tier 1 programs or activity data in high-resolution systems. Single-year simulation, deterministic reference levels, and systematic improvements achieve attractive uncertainty reductions, directly increasing revenue and while meeting ART-TREES Section 8 requirements.

ART. (2021). The REDD+ environmental excellence standard (2.0 ed.). https://www.artredd.org/wp-content/uploads/2021/12/TREES-2.0-August-2021-Clean.pdf