required_packages <- c("jsonlite", "ggplot2")
missing_packages <- required_packages[!vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)]
if (length(missing_packages) > 0) {
stop("Missing package(s): ", paste(missing_packages, collapse = ", "), call. = FALSE)
}Step 6 - Emissions
2 Step 6: AEA Production Emissions
2.1 Objective
Attach production-based CO2 coefficients to simulated sector outputs and compare baseline vs transition emissions paths.
2.2 Why this step matters
Ecological macro policy needs activity, structure, and emissions accounting in one framework. This step connects macro-structural dynamics to emissions indicators with explicit system boundaries.
2.3 What changed from Step 5
- Added emissions accounting layer to existing macro-IO simulations.
- Added sector-intensity alignment between IO sectors and AEA data.
- Kept macro dynamics unchanged to isolate emissions-accounting effects.
2.4 Equations
Economic question: how much of emissions change comes from activity scale, structural composition, and intensity assumptions?
Production-based emissions accounting: - Intensity: s_i = CO2_i / x_i - Period emissions: CO2_t = sum_i(s_i * x_i,t) - Reporting conversion: CO2_Mt,t = CO2_t / 1e9
Simple decomposition lens for scenario comparison: - Activity effect: emissions change from total output scale. - Structural effect: emissions change from sector composition at similar scale. - Intensity effect: emissions change from changing s_i (e.g., exogenous intensity decline).
Boundary clarification: - Implemented: production-based emissions (territorial production boundary). - Not implemented: consumption-based footprint requiring imported embodied intensities (MRIO logic).
Symbol glossary: - x_i,t: sector output path from macro-IO simulation. - s_i: direct production intensity by sector. - CO2_t: total production emissions in period t. - intensity_decline: exogenous decarbonization assumption in optional play.
Model status in Step 6: - Exogenous: intensity assumptions and optional decline scenario. - Endogenous: output-driven emissions path under fixed boundary definition. - Calibrated: base-year intensity alignment from AEA + IO data.
Reference note: production-vs-consumption boundary distinction is essential in MRIO-informed ecological macro analysis.
2.5 TFM (Step 6)
| Flow | Households | Government | Production | RoW | Sum |
|---|---|---|---|---|---|
Income Y |
+Y |
0 |
-Y |
0 |
0 |
Taxes T |
-T |
+T |
0 |
0 |
0 |
Consumption C |
-C |
0 |
+C |
0 |
0 |
Government demand G |
0 |
-G |
+G |
0 |
0 |
Exports EX |
0 |
0 |
+EX |
-EX |
0 |
Imports IM |
0 |
0 |
-IM |
+IM |
0 |
Financial stock changes Delta_V, Delta_B |
-Delta_V |
+Delta_B |
0 |
0 |
0 |
Change in net foreign assets Delta_NFA |
0 |
0 |
-Delta_NFA |
+Delta_NFA |
0 |
| Column sum | 0 | 0 | 0 | 0 | 0 |
From column sums: - Delta_V = YD - C, Delta_B = G - T, and Delta_V = Delta_B. - Delta_NFA = EX - IM = TB.
2.6 BSM (Step 6)
| Position | Households | Government | Production | RoW | Sum |
|---|---|---|---|---|---|
Wealth/debt stock V (=B) |
+V |
-B |
0 |
0 |
0 |
Net foreign assets NFA |
0 |
0 |
+NFA |
-NFA |
0 |
Net worth NW |
+NW_h |
+NW_g |
+NW_p |
+NW_row |
0 |
Net-worth identities: - NW_h = V - NW_g = -B - NW_p = NFA (production column carries domestic external position in this compact setup) - NW_row = -NFA - with V = B, NW_h + NW_g + NW_p + NW_row = 0.
Interpretation for Step 6: - Emissions accounting is added on top of the same institutional stock-flow structure as Step 4. - Production remains an explicit institutional column, and RoW mirrors the foreign-asset position.
Declare emissions configuration (first appearance in tutorial).
This configuration is separate from Steps 1-4, because AEA coverage constraints differ. It sets country/year/table choices for emission-intensity alignment and loading.
emissions_config <- make_emissions_config()2.7 Helper Logic Used in This Step (Detailed, not full code listing)
normalize_nace_code() (in R/helpers/emissions_load.R) standardizes sector codes before matching: 1. Convert to uppercase and normalize separators (- to _, remove dots). 2. Harmonize variant formats so IO sector codes and AEA codes can be compared one-to-one. 3. Return cleaned codes used in all later alignment checks.
align_iot_to_aea() ensures the IO system and emissions system use the same sector universe: 1. Get available AEA NACE codes for the selected emissions configuration. 2. Normalize IO sector codes and keep only sectors present in AEA. 3. Subset Z, A, and all sector vectors consistently to the retained sector set. 4. Recompute derived objects (L, F0, Y0, G0) and renormalize demand shares. 5. Return a reduced but internally consistent calibration object for emissions runs.
load_aea_emissions() builds production-emissions intensities from Eurostat AEA: 1. Enforce industry-by-industry mode (AEA is industry-coded). 2. Query AEA with chosen pollutant/year/country and parse values by NACE code. 3. Align parsed emissions to calibrated sector order. 4. Fail fast if any required sector is missing emissions data. 5. Compute intensity kg CO2 / M EUR as emissions divided by base output x0.
These helpers keep Step 6 simulation logic clean: the simulation loop only needs aligned sector outputs and a matched intensity vector.
2.8 Algorithm (pseudo-code)
- Load and align IO sectors with AEA emissions classifications.
- Run macro-IO simulation while storing sector outputs each year.
- Multiply sector outputs by intensities to compute production emissions.
- Compare baseline/transition and optional intensity-decline cases.
Define Step 6 simulation wrapper with stored sector output.
Algorithm: run the SIM+IOT+RoW macro loop while storing sectoral output vectors x_t each period. Later chunks multiply those x_t paths by AEA intensities to compute production CO2. This is the only new simulation code in Step 6; loading/alignment helpers are reused from earlier helper files.
simulate_row_with_x <- function(calib,
T = 20,
c_y = 0.82,
c_v = 0.03,
g_G = 0.01,
ex_growth = 0,
import_leakage_scale = 1,
transition_speed = 0.0015) {
years <- calib$base_year + seq_len(T) - 1
agg <- data.frame(
year = years,
GDP = NA_real_, TAX = NA_real_, YD = NA_real_, C = NA_real_, G = NA_real_,
V = NA_real_, DEF = NA_real_, B = NA_real_, EX = NA_real_, IM = NA_real_, TB = NA_real_,
betaC_brown = NA_real_, betaC_green = NA_real_
)
x_store <- matrix(NA_real_, nrow = T, ncol = calib$n)
V_prev <- calib$V0
B_prev <- calib$B0
G_prev <- calib$G0
x_prev <- as.numeric(calib$L %*% (calib$beta_C * sum(calib$I_i0) + calib$beta_G * calib$G0 + calib$I_i0 + calib$EX_i0))
beta_C <- calib$beta_C
for (tt in seq_len(T)) {
G_t <- if (tt == 1) calib$G0 else G_prev * (1 + g_G)
EX_i_t <- calib$EX_i0 * (1 + ex_growth)^(tt - 1)
if (is.finite(calib$idx$idx_green) && is.finite(calib$idx$idx_brown)) {
signal <- (x_prev[calib$idx$idx_green] - x_prev[calib$idx$idx_brown]) / pmax(sum(abs(x_prev)), 1)
shift <- max(0, transition_speed * (1 + signal))
shift <- min(shift, max(beta_C[calib$idx$idx_brown] - 1e-8, 0))
beta_C[calib$idx$idx_brown] <- beta_C[calib$idx$idx_brown] - shift
beta_C[calib$idx$idx_green] <- beta_C[calib$idx$idx_green] + shift
}
m_eff <- pmin(pmax(calib$m_i * import_leakage_scale, 0), 0.99)
S <- 1 - m_eff
k <- as.numeric(t(calib$va_coeff) %*% (calib$L %*% (S * beta_C)))
b <- as.numeric(t(calib$va_coeff) %*% (calib$L %*% (S * (calib$beta_G * G_t + calib$I_i0 + EX_i_t))))
den <- 1 - c_y * (1 - calib$tau_y) * k
C_t <- max((c_y * (1 - calib$tau_y) * b + c_v * V_prev) / den, 0)
fd_dom_t <- beta_C * C_t + calib$beta_G * G_t + calib$I_i0 + EX_i_t
im_i_t <- m_eff * fd_dom_t
fd_net_t <- pmax(fd_dom_t - im_i_t, 0)
x_t <- pmax(as.numeric(calib$L %*% fd_net_t), 0)
Y_t <- sum(calib$va_coeff * x_t)
TAX_t <- calib$tau_y * Y_t
YD_t <- Y_t - TAX_t
V_t <- V_prev + YD_t - C_t
DEF_t <- G_t - TAX_t
B_t <- B_prev + DEF_t
EX_t <- sum(EX_i_t)
IM_t <- sum(im_i_t)
TB_t <- EX_t - IM_t
agg[tt, c("GDP", "TAX", "YD", "C", "G", "V", "DEF", "B", "EX", "IM", "TB", "betaC_brown", "betaC_green")] <- c(
Y_t, TAX_t, YD_t, C_t, G_t, V_t, DEF_t, B_t, EX_t, IM_t, TB_t,
if (is.finite(calib$idx$idx_brown)) beta_C[calib$idx$idx_brown] else NA_real_,
if (is.finite(calib$idx$idx_green)) beta_C[calib$idx$idx_green] else NA_real_
)
x_store[tt, ] <- x_t
x_prev <- x_t
V_prev <- V_t
B_prev <- B_t
G_prev <- G_t
}
list(aggregate = agg, x = x_store)
}Prepare Step 6 calibration and baseline/transition runs.
Load IO and wealth with emissions configuration. Align IO sectors to AEA sector availability, then load direct CO2 intensities. Build baseline and transition-ready calibration once.
iot_emis <- download_or_load_iot(emissions_config)
if (isTRUE(emissions_config$enforce_co2_consistency)) iot_emis <- align_iot_to_aea(iot_emis, emissions_config)
wealth_emis <- load_wealth_init(emissions_config)
calib5 <- build_sim_iot_row_from_data(iot_emis, wealth_emis)
aea_intensity <- load_aea_emissions(emissions_config, iot_emis$sector_codes, iot_emis$x0)Core play: compare baseline and transition CO2.
Run baseline and transition simulations with stored x_t. Attach production intensities to both runs and compare resulting CO2 paths.
transition_speed_emis_play <- 0.0 # @exercise[id=step6_transition;kind=core;question_expr=0.0;prompt="Compare baseline and transition CO2 paths";hint="Use 0 for baseline and positive speed for transition"]
step6_baseline_sim <- simulate_row_with_x(calib5, T = 20, transition_speed = 0)
step6_transition_sim <- simulate_row_with_x(calib5, T = 20, transition_speed = transition_speed_emis_play)
step6_base <- attach_production_emissions(step6_baseline_sim, aea_intensity, iot_emis$sector_codes)
step6_trans <- attach_production_emissions(step6_transition_sim, aea_intensity, iot_emis$sector_codes)
step6_compare <- rbind(
data.frame(step6_base, scenario = paste0("base (s=", 0, ")")),
data.frame(step6_trans, scenario = paste0("play (s=", transition_speed_emis_play, ")"))
)Plot Step 6 core comparison.
Plot production CO2 for baseline and transition on the same axes. This is the main emissions comparison figure.
y_co2 <- range(step6_compare$CO2_Mt, na.rm = TRUE)
if (!all(is.finite(y_co2))) stop("Non-finite CO2 values in Step 6 core plot.", call. = FALSE)
if (y_co2[1] == y_co2[2]) y_co2 <- y_co2 + c(-1, 1) * max(1e-6, 0.05 * abs(y_co2[1]))
ggplot2::ggplot(step6_compare, ggplot2::aes(year, CO2_Mt, color = scenario)) +
ggplot2::geom_line(linewidth = 1) +
ggplot2::coord_cartesian(ylim = y_co2) +
ggplot2::labs(title = paste0("Step 6 Core: production CO2 (base s=0, play s=", transition_speed_emis_play, ")"), y = "Mt CO2") +
ggplot2::theme_minimal()
Decompose final-year change into activity and structural components (intensity held fixed).
x_base <- step6_baseline_sim$x[nrow(step6_baseline_sim$x), ]
x_trans <- step6_transition_sim$x[nrow(step6_transition_sim$x), ]
coefs <- setNames(aea_intensity$intensity_kg_per_meur, aea_intensity$sector_code)
idx <- match(iot_emis$sector_codes, names(coefs))
s_vec <- coefs[idx]
s_vec[!is.finite(s_vec)] <- 0
co2_base <- sum(s_vec * x_base)
co2_trans <- sum(s_vec * x_trans)
scale <- sum(x_trans) / pmax(sum(x_base), 1e-9)
x_activity_only <- x_base * scale
co2_activity_only <- sum(s_vec * x_activity_only)
activity_effect <- co2_activity_only - co2_base
structure_effect <- co2_trans - co2_activity_only
data.frame(component = c('activity_effect_kg', 'structure_effect_kg', 'total_change_kg'),
value = c(activity_effect, structure_effect, co2_trans - co2_base)) component value
1 activity_effect_kg 0
2 structure_effect_kg 0
3 total_change_kg 02.9 Insight (Step 6 core)
- Transition can change emissions through output scale and composition even before changing intensities.
- Reporting decomposition avoids confusing structural change with exogenous decarbonization assumptions.
Optional play: sensitivity to intensity decline.
Apply an exogenous decline factor to intensities. Recompute CO2 for the transition run to separate activity effects from intensity effects.
intensity_decline_base <- 0
intensity_decline_play <- 0.0 # @exercise[id=step6_intensity;kind=optional;question_expr=0.0;prompt="Apply annual intensity decline and compare CO2 path";hint="This is exogenous decarbonization of direct intensity"]
aea_intensity_base <- aea_intensity
aea_intensity_base$intensity_kg_per_meur <- aea_intensity_base$intensity_kg_per_meur * (1 - intensity_decline_base)
aea_intensity_play <- aea_intensity
aea_intensity_play$intensity_kg_per_meur <- aea_intensity_play$intensity_kg_per_meur * (1 - intensity_decline_play)
step6_base_opt <- attach_production_emissions(step6_transition_sim, aea_intensity_base, iot_emis$sector_codes)
step6_play <- attach_production_emissions(step6_transition_sim, aea_intensity_play, iot_emis$sector_codes)
step6_opt_compare <- rbind(
data.frame(step6_base_opt, scenario = paste0("base (d=", intensity_decline_base, ")")),
data.frame(step6_play, scenario = paste0("play (d=", intensity_decline_play, ")"))
)
y_co2_opt <- range(step6_opt_compare$CO2_Mt, na.rm = TRUE)
if (!all(is.finite(y_co2_opt))) stop("Non-finite CO2 values in Step 6 optional plot.", call. = FALSE)
if (y_co2_opt[1] == y_co2_opt[2]) y_co2_opt <- y_co2_opt + c(-1, 1) * max(1e-6, 0.05 * abs(y_co2_opt[1]))
ggplot2::ggplot(step6_opt_compare, ggplot2::aes(year, CO2_Mt, color = scenario)) +
ggplot2::geom_line(linewidth = 1) +
ggplot2::coord_cartesian(ylim = y_co2_opt) +
ggplot2::labs(title = paste0("Step 6 Optional: intensity decline (base d=", intensity_decline_base, ", play d=", intensity_decline_play, ")"), y = "Mt CO2") +
ggplot2::theme_minimal()
2.10 Insight (Step 6 optional intensity play)
- Intensity decline can offset activity-driven emissions growth, but does not remove scale effects automatically.
- If GDP/output rises strongly, rebound-like patterns can appear despite falling intensity assumptions.
2.11 What can go wrong / interpretation caveat
- Alignment mismatches between IO and AEA classifications can bias intensity estimates if unchecked.
- This is production accounting only; footprint conclusions require MRIO-style imported-emissions modeling.
2.12 Consumption-emissions caveat (equations/text only)
A consumption footprint would require domestic and imported embodied intensities, not only domestic production intensities.
3 Final Snapshot
Report final-year key aggregates from Step 6.
Extract the final year and display key macro/CO2 indicators by scenario.
last_year <- max(step6_compare$year)
summary_final <- step6_compare[step6_compare$year == last_year, c("scenario", "GDP", "V", "B", "TB", "CO2_Mt")]
names(summary_final)[names(summary_final) == "GDP"] <- "GDP_mEUR"
names(summary_final)[names(summary_final) == "V"] <- "V_mEUR"
names(summary_final)[names(summary_final) == "B"] <- "B_mEUR"
names(summary_final)[names(summary_final) == "TB"] <- "TB_mEUR"
summary_final <- summary_final[order(summary_final$scenario), ]
rownames(summary_final) <- NULL
summary_final scenario GDP_mEUR V_mEUR B_mEUR TB_mEUR CO2_Mt
1 base (s=0) 720818.9 2002408 -998732.3 415082.4 122.9031
2 play (s=0) 720818.9 2002408 -998732.3 415082.4 122.90314 Limitations
- Prices are fixed.
- No banking, portfolio, or inventories.
- SIM+IOT+RoW is pedagogical; MRIO is required for full footprint attribution.