From 9c68dfaacaf4284a6356554680013d1d7a5f4e7e Mon Sep 17 00:00:00 2001 From: reetiksahu Date: Tue, 19 May 2026 08:39:27 -0700 Subject: [PATCH 1/9] Installing firm level carbon taxes --- macro_data/processing/synthetic_country.py | 4 + macro_data/readers/default_readers.py | 34 ++- macro_data/readers/policy_data/__init__.py | 3 + .../readers/policy_data/obps_can_reader.py | 85 ++++++++ macromodel/agents/firms/firms.py | 4 + macromodel/agents/firms/func/prices.py | 21 +- .../configurations/country_configuration.py | 1 + macromodel/country/country.py | 54 ++++- macromodel/policy/__init__.py | 3 + .../policy/output_based_price_system_can.py | 194 ++++++++++++++++++ macromodel/simulation.py | 5 + 11 files changed, 402 insertions(+), 6 deletions(-) create mode 100644 macro_data/readers/policy_data/__init__.py create mode 100644 macro_data/readers/policy_data/obps_can_reader.py create mode 100644 macromodel/policy/__init__.py create mode 100644 macromodel/policy/output_based_price_system_can.py diff --git a/macro_data/processing/synthetic_country.py b/macro_data/processing/synthetic_country.py index 6d665c6d..75a377ad 100644 --- a/macro_data/processing/synthetic_country.py +++ b/macro_data/processing/synthetic_country.py @@ -112,6 +112,7 @@ from macro_data.readers.emissions.emissions_reader import CH4EmissionsDataCAN, EmissionsData from macro_data.readers.exo_prices.exo_prices_reader import SectorExoPrices from macro_data.readers.exogenous_data import ExogenousCountryData +from macro_data.readers.policy_data.obps_can_reader import OBPSCANData @dataclass @@ -173,6 +174,7 @@ class SyntheticCountry: emission_factors: EmissionsData emission_fractions: Optional[EmissionFractions] = None firm_exo_prices: Optional[SectorExoPrices] = None + obps_data: Optional[OBPSCANData] = None emission_factors_ch4: Optional[CH4EmissionsDataCAN] = None historical_emissions_df: Optional[pd.DataFrame] = None @@ -378,6 +380,7 @@ def eu_synthetic_country( firm_exo_prices=( SectorExoPrices.from_reader(readers.exo_prices) if readers.exo_prices is not None else None ), + obps_data=readers.obps_can.obps_data if readers.obps_can is not None else None, ) @classmethod @@ -594,6 +597,7 @@ def proxied_synthetic_country( firm_exo_prices=( SectorExoPrices.from_reader(readers.exo_prices) if readers.exo_prices is not None else None ), + obps_data=readers.obps_can.obps_data if readers.obps_can is not None else None, ) @classmethod diff --git a/macro_data/readers/default_readers.py b/macro_data/readers/default_readers.py index 69485972..2977337b 100644 --- a/macro_data/readers/default_readers.py +++ b/macro_data/readers/default_readers.py @@ -50,6 +50,7 @@ ) from macro_data.readers.io_tables.icio_reader import ICIOReader, split_gfcf_column from macro_data.readers.io_tables.industries import AGGREGATED_INDUSTRIES +from macro_data.readers.policy_data.obps_can_reader import OBPSCANReader from macro_data.readers.population_data.compustat_banks_reader import ( CompustatBanksReader, ) @@ -61,6 +62,21 @@ from macro_data.readers.util.prune_util import DataFilterWarning +@dataclass +class OBPSPaths: + """File paths for the Canada Output-Based Pricing System data. + + Attributes: + rates_path: CSV of carbon price rates by year and jurisdiction. + policy_path: CSV of per-industry reduction factors and tightening rates. + policy_elec_path: Optional CSV for electricity-specific tightening rates. + """ + + rates_path: Path + policy_path: Path + policy_elec_path: Optional[Path] = None + + @dataclass class DataPaths: """Manages file paths for all data sources used in the model. @@ -113,6 +129,7 @@ class DataPaths: emissions_fraction_path: Optional[Path] = None firm_prices_path: Optional[Path] = None ch4_emissions_path: Optional[Path] = None + obps_path: Optional[OBPSPaths] = None @classmethod def default_paths(cls, raw_data_path: Path, icio_years: Iterable[int]): @@ -150,6 +167,11 @@ def default_paths(cls, raw_data_path: Path, icio_years: Iterable[int]): ch4_emissions_path=raw_data_path / "emission_factors" / "EN-GHG_EconSectByGas-CA_Emissions_2014_2023_v4.csv", + obps_path=OBPSPaths( + rates_path=raw_data_path / "policy" / "output_based_price_system_rates.csv", + policy_path=raw_data_path / "policy" / "output_based_price_system_policy_values_disagg.csv", + policy_elec_path=raw_data_path / "policy" / "output_based_price_system_policy_values_elec.csv", + ), ) # @classmethod @@ -206,6 +228,7 @@ class DataReaders: emission_fractions: Optional[EmissionsFractionReader] = None exo_prices: Optional[SectorExoPricesReader] = None ch4_emissions: Optional[CH4EmissionsReaderCAN] = None + obps_can: Optional[OBPSCANReader] = None regions_dict: Optional[dict[Country, list[Region]]] = None @classmethod @@ -302,7 +325,7 @@ def get_investment_year(year: int, country_names_: Optional[list[Country | Regio if simulation_year != 2014: raise ValueError("Only 2014 is supported for this reader.") - disagg_path = raw_data_path / "icio" / "icio_can_2014_disagg.csv" + disagg_path = raw_data_path / "icio" / "sectoral_disagg_CAN_2014_v2.csv" df = pd.read_csv(disagg_path, header=[0, 1], index_col=[0, 1]) icio[simulation_year].iot = df industries = df.loc["ROW"].index.unique() @@ -481,6 +504,14 @@ def get_investment_year(year: int, country_names_: Optional[list[Country | Regio if datapaths.ch4_emissions_path is not None and datapaths.ch4_emissions_path.exists(): ch4_emissions = CH4EmissionsReaderCAN.read_data(datapaths.ch4_emissions_path) + obps_can = None + if datapaths.obps_path is not None: + obps_can = OBPSCANReader.read_from_raw_data( + rates_path=datapaths.obps_path.rates_path, + policy_path=datapaths.obps_path.policy_path, + policy_elec_path=datapaths.obps_path.policy_elec_path, + ) + return cls( icio=icio, wiod_sea=wiod_sea, @@ -500,6 +531,7 @@ def get_investment_year(year: int, country_names_: Optional[list[Country | Regio emission_fractions=emission_fractions, exo_prices=exo_prices, ch4_emissions=ch4_emissions, + obps_can=obps_can, regions_dict=regions_dict, ) diff --git a/macro_data/readers/policy_data/__init__.py b/macro_data/readers/policy_data/__init__.py new file mode 100644 index 00000000..bcf61f94 --- /dev/null +++ b/macro_data/readers/policy_data/__init__.py @@ -0,0 +1,3 @@ +from macro_data.readers.policy_data.obps_can_reader import OBPSCANData, OBPSCANReader + +__all__ = ["OBPSCANData", "OBPSCANReader"] diff --git a/macro_data/readers/policy_data/obps_can_reader.py b/macro_data/readers/policy_data/obps_can_reader.py new file mode 100644 index 00000000..c49496b8 --- /dev/null +++ b/macro_data/readers/policy_data/obps_can_reader.py @@ -0,0 +1,85 @@ +"""Reader and container for Canada Output-Based Pricing System (OBPS) policy data.""" + +from __future__ import annotations + +from dataclasses import dataclass +from pathlib import Path +from typing import Optional + +import pandas as pd + + +@dataclass +class OBPSCANData: + """Container for Canada OBPS policy tables loaded from CSV files. + + Attributes: + df_rates: Carbon price trajectory. Rows are years; columns are + 'Date' plus one column per province/jurisdiction code. + df_policy: Per-industry emission standards. Must have columns + 'Industry', 'reduction_factor', and 'tightening_rate'. + df_policy_elec: Optional electricity-specific tightening rates. + Same structure as df_policy but for electricity sectors only. + """ + + df_rates: pd.DataFrame + df_policy: pd.DataFrame + df_policy_elec: Optional[pd.DataFrame] = None + + +@dataclass +class OBPSCANReader: + """Reader for Canada OBPS CSV data files. + + CSV formats + ----------- + rates_path: + Date, CAN, CAN_AB, CAN_BC, ... + 2014, 15.0, 10.0, 12.0, ... + 2015, 20.0, 15.0, 15.0, ... + + policy_path: + Industry, reduction_factor, tightening_rate + B05a, 0.80, 0.02 + C19, 0.75, 0.02 + + policy_elec_path (optional): + same structure as policy_path but for electricity sectors. + """ + + obps_data: Optional[OBPSCANData] = None + + @classmethod + def read_from_raw_data( + cls, + rates_path: Path | str, + policy_path: Path | str, + policy_elec_path: Optional[Path | str] = None, + ) -> OBPSCANReader: + """Load OBPS data from CSV files. + + Args: + rates_path: Path to carbon price rates CSV. + policy_path: Path to per-industry policy values CSV. + policy_elec_path: Optional path to electricity-specific CSV. + + Returns: + OBPSCANReader with loaded OBPSCANData, or obps_data=None if any + required file is absent. + """ + rates_path = Path(rates_path) + policy_path = Path(policy_path) + + if not rates_path.exists() or not policy_path.exists(): + return cls(obps_data=None) + + df_rates = pd.read_csv(rates_path) + df_policy = pd.read_csv(policy_path) + + df_policy_elec = None + if policy_elec_path is not None: + policy_elec_path = Path(policy_elec_path) + if policy_elec_path.exists(): + df_policy_elec = pd.read_csv(policy_elec_path) + + return cls(obps_data=OBPSCANData(df_rates=df_rates, df_policy=df_policy, df_policy_elec=df_policy_elec)) diff --git a/macromodel/agents/firms/firms.py b/macromodel/agents/firms/firms.py index 53114c2a..4eb2184c 100644 --- a/macromodel/agents/firms/firms.py +++ b/macromodel/agents/firms/firms.py @@ -1011,6 +1011,7 @@ def compute_price( current_estimated_ppi_inflation: np.ndarray, previous_average_good_prices: np.ndarray, ppi_during: np.ndarray, + extra_marginal_taxes: Optional[np.ndarray] = None, ) -> np.ndarray: """Set prices for each firm's output. @@ -1026,6 +1027,8 @@ def compute_price( current_estimated_ppi_inflation (np.ndarray): Expected PPI inflation previous_average_good_prices (np.ndarray): Previous period prices ppi_during (np.ndarray): Producer price indices + extra_marginal_taxes (np.ndarray, optional): Per-sector marginal + tax added to sector average prices. Defaults to None. Returns: np.ndarray: New prices for each firm @@ -1051,6 +1054,7 @@ def compute_price( ), ppi_during=ppi_during, current_time=len(self.ts.historic("price")), + extra_marginal_taxes=extra_marginal_taxes, ) def compute_unconstrained_demand_for_intermediate_inputs( diff --git a/macromodel/agents/firms/func/prices.py b/macromodel/agents/firms/func/prices.py index 4de3470a..3a6ef229 100644 --- a/macromodel/agents/firms/func/prices.py +++ b/macromodel/agents/firms/func/prices.py @@ -1,4 +1,5 @@ from abc import ABC, abstractmethod +from typing import Optional import numpy as np from scipy.interpolate import interp1d @@ -74,6 +75,7 @@ def compute_price( prev_unit_costs: np.ndarray, ppi_during: np.ndarray, current_time: int, + extra_marginal_taxes: Optional[np.ndarray] = None, ) -> np.ndarray: """Calculate prices for each firm based on market conditions. @@ -139,11 +141,12 @@ def compute_price( current_time: int, min_inflation: float = -0.1, max_inflation: float = 0.1, + extra_marginal_taxes: Optional[np.ndarray] = None, ) -> np.ndarray: """Calculate prices using the default multi-factor strategy. The method: - 1. Maps sector average prices to firms + 1. Maps sector average prices to firms (plus any OBPS marginal tax) 2. Calculates demand-pull inflation based on market position 3. Calculates cost-push inflation from unit costs 4. Combines all factors with random noise @@ -158,11 +161,15 @@ def compute_price( Defaults to -0.1 (-10%). max_inflation (float, optional): Upper bound on inflation rates. Defaults to 0.1 (10%). + extra_marginal_taxes (np.ndarray, optional): Per-sector marginal + tax (e.g. OBPS) added to sector average prices seen by firms. + Shape (n_industries,). Defaults to None. Returns: np.ndarray: Updated prices by firm, guaranteed to be positive """ - average_price_by_firm = prev_average_good_prices[current_firm_sectors] + tax_by_sector = extra_marginal_taxes if extra_marginal_taxes is not None else np.zeros_like(prev_average_good_prices) + average_price_by_firm = (prev_average_good_prices + tax_by_sector)[current_firm_sectors] # Demand-pull inflation demand_pull_inflation = np.zeros_like(prev_firm_prices) @@ -276,6 +283,7 @@ def compute_price( current_time: int, min_inflation: float = -0.1, max_inflation: float = 0.1, + extra_marginal_taxes: Optional[np.ndarray] = None, ) -> np.ndarray: """Compute prices, overriding listed sectors with exogenous sector paths. @@ -317,6 +325,7 @@ def compute_price( current_time=current_time, min_inflation=min_inflation, max_inflation=max_inflation, + extra_marginal_taxes=extra_marginal_taxes, ) if self.firm_exo_prices is None or len(self.overriden_industries) == 0: @@ -328,12 +337,18 @@ def compute_price( else prev_average_good_prices ) + tax_by_firm = ( + extra_marginal_taxes[current_firm_sectors] + if extra_marginal_taxes is not None + else np.zeros_like(price) + ) + for industry_name in self.firm_exo_prices.prices.columns: if industry_name not in self.overriden_industries: continue ratio = self._normalised_price(industry_name, current_quarter=current_time) for idx in self._indices_for(industry_name): - price[idx] = base_prices[idx] * ratio + price[idx] = base_prices[idx] * ratio + tax_by_firm[idx] return price diff --git a/macromodel/configurations/country_configuration.py b/macromodel/configurations/country_configuration.py index 5a26b80d..dac6be69 100644 --- a/macromodel/configurations/country_configuration.py +++ b/macromodel/configurations/country_configuration.py @@ -39,6 +39,7 @@ class CountryConfiguration(BaseModel): assume_zero_noise: bool = False use_emission_multiplier: bool = False CH4_production_emissions_only: bool = False + use_obps_reg: bool = False @classmethod def n_industry_default( diff --git a/macromodel/country/country.py b/macromodel/country/country.py index e747f178..a98460f4 100644 --- a/macromodel/country/country.py +++ b/macromodel/country/country.py @@ -55,6 +55,7 @@ from macromodel.markets.credit_market.credit_market import CreditMarket from macromodel.markets.housing_market.housing_market import HousingMarket from macromodel.markets.labour_market.labour_market import LabourMarket +from macromodel.policy.output_based_price_system_can import OutputBasedPriceSystemCAN from macromodel.rest_of_the_world import RestOfTheWorld from macromodel.util.get_histogram import get_histogram @@ -143,6 +144,7 @@ def __init__( emitting_indices: Optional[np.ndarray] = None, emission_factors_lcu_ch4: Optional[np.ndarray] = None, emitting_indices_ch4: Optional[np.ndarray] = None, + obps: Optional[OutputBasedPriceSystemCAN] = None, ): """Initialize a new country economy. @@ -213,6 +215,10 @@ def __init__( self.emitting_indices_ch4 = emitting_indices_ch4 self.use_emission_multiplier = self.configuration.use_emission_multiplier + self.obps = obps + self.use_obps_reg = self.configuration.use_obps_reg + self.extra_marginal_taxes_firm = np.zeros(self.firms.n_industries) + @classmethod def from_pickled_country( cls, @@ -400,6 +406,14 @@ def from_pickled_country( scale=scale, ) + obps = None + if add_emissions and country_configuration.use_obps_reg and synthetic_country.obps_data is not None: + obps = OutputBasedPriceSystemCAN( + country_name=country_name, + industries=list(industries), + obps_data=synthetic_country.obps_data, + ) + return cls( country_name=country_name, scale=scale, @@ -425,6 +439,7 @@ def from_pickled_country( emitting_indices=emitting_indices, emission_factors_lcu_ch4=emission_factors_lcu_ch4, emitting_indices_ch4=emitting_indices_ch4, + obps=obps, ) def reset(self, configuration: CountryConfiguration) -> None: @@ -535,6 +550,35 @@ def target_setting_phase(self) -> None: ) ) + def update_extra_taxes(self, record_obps_reference: bool = True) -> None: + """Compute extra marginal taxes for firms from active policy instruments. + + Currently supports the Output-Based Pricing System (OBPS). The sectoral + tax cost is divided by production to obtain a per-unit marginal cost that + is added to the sector average price seen by firms during price-setting + and input-demand calculations. + + Args: + record_obps_reference: If True, accumulate 2017–2019 reference + emission data (should be True only in the planning phase). + """ + self.extra_marginal_taxes_firm = np.zeros(self.firms.n_industries) + + if self.use_obps_reg and self.obps is not None: + sectoral_tax = self.obps.compute_obps( + use_obps_reg=self.use_obps_reg, + record_obps_reference=record_obps_reference, + production=self.firms.ts.current("production"), + input_em=self.firms.ts.current("inputs_emissions")+self.firms.ts.current("inputs_emissions_ch4"), + capital_em=self.firms.ts.current("capital_emissions")+self.firms.ts.current("capital_emissions_ch4"), + ) + self.extra_marginal_taxes_firm = np.divide( + sectoral_tax, + self.firms.ts.current("production"), + out=np.zeros_like(sectoral_tax), + where=self.firms.ts.current("production") != 0, + ) + def clear_labour_market(self) -> None: """Execute labor market clearing. @@ -555,6 +599,9 @@ def update_planning_metrics(self) -> None: Computes expected profits, asset values, benefits, and other metrics used by agents in their planning decisions. """ + if self.add_emissions: + self.update_extra_taxes(record_obps_reference=True) + # Firms estimate profits self.firms.ts.expected_profits.append( self.firms.compute_estimated_profits( @@ -669,13 +716,15 @@ def update_planning_metrics(self) -> None: current_estimated_ppi_inflation=self.economy.ts.current("estimated_ppi_inflation")[0], previous_average_good_prices=self.economy.ts.current("good_prices"), ppi_during=self.exogenous.national_accounts_during["PPI (Value)"].values.flatten(), + extra_marginal_taxes=self.extra_marginal_taxes_firm, ) ) # Firm demand for goods self.firms.ts.unconstrained_target_intermediate_inputs.append( self.firms.compute_unconstrained_demand_for_intermediate_inputs( - good_prices=self.economy.ts.current("good_prices") + good_prices=self.economy.ts.current("good_prices"), + extra_taxes=self.extra_marginal_taxes_firm, ) ) self.firms.ts.unconstrained_target_intermediate_inputs_costs.append( @@ -685,7 +734,8 @@ def update_planning_metrics(self) -> None: ) self.firms.ts.unconstrained_target_capital_inputs.append( self.firms.compute_unconstrained_demand_for_capital_inputs( - good_prices=self.economy.ts.current("good_prices") + good_prices=self.economy.ts.current("good_prices"), + extra_taxes=self.extra_marginal_taxes_firm, ) ) self.firms.ts.unconstrained_target_capital_inputs_costs.append( diff --git a/macromodel/policy/__init__.py b/macromodel/policy/__init__.py new file mode 100644 index 00000000..8c256d56 --- /dev/null +++ b/macromodel/policy/__init__.py @@ -0,0 +1,3 @@ +from macromodel.policy.output_based_price_system_can import OutputBasedPriceSystemCAN + +__all__ = ["OutputBasedPriceSystemCAN"] diff --git a/macromodel/policy/output_based_price_system_can.py b/macromodel/policy/output_based_price_system_can.py new file mode 100644 index 00000000..43785468 --- /dev/null +++ b/macromodel/policy/output_based_price_system_can.py @@ -0,0 +1,194 @@ +"""Canada Output-Based Pricing System (OBPS) policy for the macroeconomic model. + +Calculates the tax that firms pay on emissions that exceed a sector-specific +prescribed limit. The tax is computed as: + + obps_cost[i] = max(0, emissions[i] - limit[i]) * carbon_price[t] + +Dividing by production gives a per-unit marginal cost passed to firms as +extra_marginal_taxes_firm in the country's target-setting phase. +""" + +from dataclasses import dataclass + +import numpy as np +import pandas as pd + +from macro_data.readers.policy_data.obps_can_reader import OBPSCANData + + +@dataclass +class OutputBasedPriceSystemCAN: + """Canada Output-Based Pricing System policy class. + + Calculates the tax that regulated firms pay on emissions above a + prescribed output-weighted limit. Reference emission intensities are + recorded during 2017–2019 and used to set sector-specific limits from + 2019 onwards. + + Attributes: + country_name: Jurisdiction column to use from the rates CSV. + industries: Ordered list of all model industry names. + regulated_industries: Industries subject to OBPS regulation. + regulated_indices: Array indices into the industries list for each + regulated industry. + df_policy: Per-industry reduction factors and tightening rates. + df_policy_elec: Electricity-specific tightening rates (optional). + df_rates: Carbon price schedule by year and jurisdiction. + reference_emission_intensity: Baseline emission intensity recorded + over the 2017–2019 reference period. + reference_emission: Cumulative emissions during the reference period. + reference_production: Cumulative production during the reference period. + emission_limit: Current period allowable emissions per industry. + price: Annual carbon price trajectory (indexed from 2014). + current_t: Current timestep index (incremented by update()). + current_year: Current calendar year. + """ + + country_name: str + industries: list[str] + regulated_industries: list[str] + regulated_indices: np.ndarray + df_policy: pd.DataFrame + df_policy_elec: pd.DataFrame + df_rates: pd.DataFrame + reference_emission_intensity: np.ndarray + reference_emission: np.ndarray + reference_production: np.ndarray + emission_limit: np.ndarray + price: np.ndarray + current_t: int = 0 + current_year: int = 2014 + + def __init__(self, country_name: str, industries: list[str], obps_data: OBPSCANData): + """Initialise the OBPS from a loaded OBPSCANData container. + + Args: + country_name: Jurisdiction code matching a column in obps_data.df_rates. + industries: Ordered list of all model industry names. + obps_data: Loaded OBPS CSV data. + """ + self.country_name = country_name + self.industries = industries + + # Industries regulated under OBPS (federal schedule) + self.regulated_industries = [ + "B05a", "B05b", "B05c", + "B07", "B09", + "C10T12", "C16", "C17", "C19", "C20", "C21", "C22", "C23", + "C24a", "C24b", "C29", "C30", + "D01b", "D01c", + ] + self.regulated_indices = np.array( + [list(industries).index(ind) for ind in self.regulated_industries if ind in industries] + ) + + n = len(industries) + self.reference_emission_intensity = np.zeros(n) + self.reference_emission = np.zeros(n) + self.reference_production = np.zeros(n) + self.emission_limit = np.zeros(n) + + self.df_policy = obps_data.df_policy + self.df_policy_elec = obps_data.df_policy_elec if obps_data.df_policy_elec is not None else pd.DataFrame() + self.df_rates = obps_data.df_rates + + self.price = np.zeros(len(self.df_rates)) + df_sub = self.df_rates[["Date", self.country_name]] + for t in range(len(self.df_rates)): + df_row = df_sub[df_sub["Date"] == t + 2014] + self.price[t] = df_row[self.country_name].values[0] + + def compute_obps( + self, + use_obps_reg: bool, + record_obps_reference: bool, + production: np.ndarray, + input_em: np.ndarray, + capital_em: np.ndarray, + ) -> np.ndarray: + """Compute per-sector OBPS tax cost. + + Records reference emission intensities during 2017–2019. From 2019 + onwards, computes the cost of emissions above the prescribed limit + at the current carbon price. + + Args: + use_obps_reg: If False, returns a zero array. + record_obps_reference: If True, accumulate reference period data. + production: Current-period production per industry. + input_em: Input-related CO₂e emissions per industry. + capital_em: Capital-related CO₂e emissions per industry. + + Returns: + np.ndarray: OBPS tax cost (dollars) per industry; zero for + unregulated industries or years before 2019. + """ + if not use_obps_reg: + return np.zeros(len(self.industries)) + + if record_obps_reference and self.current_year in (2017, 2018, 2019): + self.reference_emission += input_em + capital_em + self.reference_production += production + + if self.current_year == 2019: + self.reference_emission_intensity = np.divide( + self.reference_emission, + self.reference_production, + out=np.zeros_like(self.reference_emission), + where=self.reference_production != 0, + ) + + if self.current_year < 2019: + return np.zeros(len(self.industries)) + + obps_cost = np.zeros(len(self.industries)) + for i in self.regulated_indices: + if production[i] > 0: + limit = self.get_limit(i, production[i]) + self.emission_limit[i] = limit + difference = (input_em[i] + capital_em[i]) - limit + obps_cost[i] = difference * self.price[min(self.current_t, len(self.price) - 1)] + + return obps_cost + + def get_limit(self, industry_idx: int, production: float) -> float: + """Calculate the prescribed emission allowance for an industry. + + Uses the pre-2023 formula (reduction factor only) or the post-2023 + formula (reduction factor × tightening adjustment). + + Args: + industry_idx: Index into self.industries. + production: Current period output. + + Returns: + float: Allowable emissions in tCO₂e. + """ + industry_name = self.industries[industry_idx] + row = self.df_policy[self.df_policy["Industry"] == industry_name] + if row.empty: + return 0.0 + + reduction_factor = row["reduction_factor"].values[0] + B = reduction_factor * self.reference_emission_intensity[industry_idx] + + if self.current_year < 2023: + return production * B + + tightening_rate = row["tightening_rate"].values[0] + return production * (B - B * tightening_rate * (self.current_year - 2022)) + + def get_price(self) -> float: + """Return the current period carbon price ($/tCO₂e).""" + return self.price[self.current_t] + + def update(self) -> None: + """Advance the timestep by one annual period.""" + self.current_t += 1 + self.current_year += 1 + + def reset(self) -> None: + """Reset time variables to the initial year.""" + self.current_t = 0 + self.current_year = 2014 diff --git a/macromodel/simulation.py b/macromodel/simulation.py index fa8a74b4..40861c1c 100644 --- a/macromodel/simulation.py +++ b/macromodel/simulation.py @@ -302,6 +302,11 @@ def iterate(self, t: int = 0): # self.exchange_rates.set_current_exchange_rates(current_year=self.timestep.year) + for country in self.countries.values(): + if country.obps is not None: + while country.obps.current_year < self.timestep.year: + country.obps.update() + for ind, country in enumerate(self.countries.values()): exchange_rate = self.exchange_rates.get_current_exchange_rates_from_usd_to_lcu( country_name=country.country_name, From d2792ad5bad25a5c6055ed1806dad7b5a63679ac Mon Sep 17 00:00:00 2001 From: reetiksahu Date: Wed, 20 May 2026 13:36:17 -0700 Subject: [PATCH 2/9] style: apply ruff formatting to carbon tax files Co-Authored-By: Claude Sonnet 4.6 --- macro_data/readers/default_readers.py | 2 +- macromodel/agents/firms/func/prices.py | 8 +++---- macromodel/country/country.py | 4 ++-- .../policy/output_based_price_system_can.py | 24 +++++++++++++++---- 4 files changed, 26 insertions(+), 12 deletions(-) diff --git a/macro_data/readers/default_readers.py b/macro_data/readers/default_readers.py index 2977337b..f6446f0a 100644 --- a/macro_data/readers/default_readers.py +++ b/macro_data/readers/default_readers.py @@ -325,7 +325,7 @@ def get_investment_year(year: int, country_names_: Optional[list[Country | Regio if simulation_year != 2014: raise ValueError("Only 2014 is supported for this reader.") - disagg_path = raw_data_path / "icio" / "sectoral_disagg_CAN_2014_v2.csv" + disagg_path = raw_data_path / "icio" / "icio_can_2014_disagg.csv" df = pd.read_csv(disagg_path, header=[0, 1], index_col=[0, 1]) icio[simulation_year].iot = df industries = df.loc["ROW"].index.unique() diff --git a/macromodel/agents/firms/func/prices.py b/macromodel/agents/firms/func/prices.py index 3a6ef229..db95fd70 100644 --- a/macromodel/agents/firms/func/prices.py +++ b/macromodel/agents/firms/func/prices.py @@ -168,7 +168,9 @@ def compute_price( Returns: np.ndarray: Updated prices by firm, guaranteed to be positive """ - tax_by_sector = extra_marginal_taxes if extra_marginal_taxes is not None else np.zeros_like(prev_average_good_prices) + tax_by_sector = ( + extra_marginal_taxes if extra_marginal_taxes is not None else np.zeros_like(prev_average_good_prices) + ) average_price_by_firm = (prev_average_good_prices + tax_by_sector)[current_firm_sectors] # Demand-pull inflation @@ -338,9 +340,7 @@ def compute_price( ) tax_by_firm = ( - extra_marginal_taxes[current_firm_sectors] - if extra_marginal_taxes is not None - else np.zeros_like(price) + extra_marginal_taxes[current_firm_sectors] if extra_marginal_taxes is not None else np.zeros_like(price) ) for industry_name in self.firm_exo_prices.prices.columns: diff --git a/macromodel/country/country.py b/macromodel/country/country.py index a98460f4..7dad018d 100644 --- a/macromodel/country/country.py +++ b/macromodel/country/country.py @@ -569,8 +569,8 @@ def update_extra_taxes(self, record_obps_reference: bool = True) -> None: use_obps_reg=self.use_obps_reg, record_obps_reference=record_obps_reference, production=self.firms.ts.current("production"), - input_em=self.firms.ts.current("inputs_emissions")+self.firms.ts.current("inputs_emissions_ch4"), - capital_em=self.firms.ts.current("capital_emissions")+self.firms.ts.current("capital_emissions_ch4"), + input_em=self.firms.ts.current("inputs_emissions") + self.firms.ts.current("inputs_emissions_ch4"), + capital_em=self.firms.ts.current("capital_emissions") + self.firms.ts.current("capital_emissions_ch4"), ) self.extra_marginal_taxes_firm = np.divide( sectoral_tax, diff --git a/macromodel/policy/output_based_price_system_can.py b/macromodel/policy/output_based_price_system_can.py index 43785468..c8cf827e 100644 --- a/macromodel/policy/output_based_price_system_can.py +++ b/macromodel/policy/output_based_price_system_can.py @@ -73,11 +73,25 @@ def __init__(self, country_name: str, industries: list[str], obps_data: OBPSCAND # Industries regulated under OBPS (federal schedule) self.regulated_industries = [ - "B05a", "B05b", "B05c", - "B07", "B09", - "C10T12", "C16", "C17", "C19", "C20", "C21", "C22", "C23", - "C24a", "C24b", "C29", "C30", - "D01b", "D01c", + "B05a", + "B05b", + "B05c", + "B07", + "B09", + "C10T12", + "C16", + "C17", + "C19", + "C20", + "C21", + "C22", + "C23", + "C24a", + "C24b", + "C29", + "C30", + "D01b", + "D01c", ] self.regulated_indices = np.array( [list(industries).index(ind) for ind in self.regulated_industries if ind in industries] From 1acb72fd77b378a4b2f8821e4566236a39c61ff3 Mon Sep 17 00:00:00 2001 From: reetiksahu Date: Wed, 20 May 2026 13:46:19 -0700 Subject: [PATCH 3/9] fix: correct OBPS docstrings to reflect two-way cost formula MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Module and compute_obps docstrings claimed max(0, emissions - limit) but the implementation computes the signed difference with no clamp — sectors below their benchmark receive a rebate (negative cost). Co-Authored-By: Claude Sonnet 4.6 --- .../policy/output_based_price_system_can.py | 21 +++++++++++-------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/macromodel/policy/output_based_price_system_can.py b/macromodel/policy/output_based_price_system_can.py index c8cf827e..e6288494 100644 --- a/macromodel/policy/output_based_price_system_can.py +++ b/macromodel/policy/output_based_price_system_can.py @@ -1,12 +1,14 @@ """Canada Output-Based Pricing System (OBPS) policy for the macroeconomic model. -Calculates the tax that firms pay on emissions that exceed a sector-specific -prescribed limit. The tax is computed as: +Computes a two-way carbon price signal for each regulated sector based on +emissions relative to an output-based benchmark: - obps_cost[i] = max(0, emissions[i] - limit[i]) * carbon_price[t] + obps_cost[i] = (emissions[i] - limit[i]) * carbon_price[t] -Dividing by production gives a per-unit marginal cost passed to firms as -extra_marginal_taxes_firm in the country's target-setting phase. +Sectors emitting above their benchmark incur a positive cost; sectors below +receive a negative cost (rebate). Dividing by production gives a per-unit +marginal tax passed to firms as extra_marginal_taxes_firm in the country's +target-setting phase. """ from dataclasses import dataclass @@ -124,8 +126,8 @@ def compute_obps( """Compute per-sector OBPS tax cost. Records reference emission intensities during 2017–2019. From 2019 - onwards, computes the cost of emissions above the prescribed limit - at the current carbon price. + onwards, computes the signed cost of emissions relative to the + prescribed limit at the current carbon price. Args: use_obps_reg: If False, returns a zero array. @@ -135,8 +137,9 @@ def compute_obps( capital_em: Capital-related CO₂e emissions per industry. Returns: - np.ndarray: OBPS tax cost (dollars) per industry; zero for - unregulated industries or years before 2019. + np.ndarray: Signed OBPS cost (dollars) per industry — positive + when emissions exceed the limit, negative when below it; + zero for unregulated industries or years before 2019. """ if not use_obps_reg: return np.zeros(len(self.industries)) From ff2077875dfd45ca15028bc52c638d2250ae25e9 Mon Sep 17 00:00:00 2001 From: reetiksahu Date: Wed, 20 May 2026 22:14:26 -0700 Subject: [PATCH 4/9] Update tests for carbon tax --- .../test_readers/test_policy_data/__init__.py | 0 .../test_policy_data/test_obps_can_reader.py | 59 ++++++ .../unit/test_policy/__init__.py | 0 .../unit/test_policy/test_obps_can.py | 181 ++++++++++++++++++ 4 files changed, 240 insertions(+) create mode 100644 tests/test_macro_data/unit/test_readers/test_policy_data/__init__.py create mode 100644 tests/test_macro_data/unit/test_readers/test_policy_data/test_obps_can_reader.py create mode 100644 tests/test_macromodel/unit/test_policy/__init__.py create mode 100644 tests/test_macromodel/unit/test_policy/test_obps_can.py diff --git a/tests/test_macro_data/unit/test_readers/test_policy_data/__init__.py b/tests/test_macro_data/unit/test_readers/test_policy_data/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/tests/test_macro_data/unit/test_readers/test_policy_data/test_obps_can_reader.py b/tests/test_macro_data/unit/test_readers/test_policy_data/test_obps_can_reader.py new file mode 100644 index 00000000..b8499e62 --- /dev/null +++ b/tests/test_macro_data/unit/test_readers/test_policy_data/test_obps_can_reader.py @@ -0,0 +1,59 @@ +import pytest + +from macro_data.readers.policy_data.obps_can_reader import OBPSCANReader + +RATES_CSV = "Date,CAN,CAN_BC\n2014,15.0,12.0\n2019,40.0,30.0\n2030,170.0,150.0\n" +POLICY_CSV = "Industry,reduction_factor,tightening_rate\nC24a,0.80,0.02\nC19,0.75,0.02\n" +ELEC_CSV = "Industry,reduction_factor,tightening_rate\nD01b,0.70,0.01\n" + + +class TestOBPSCANReader: + def test_missing_rates_file_returns_none(self, tmp_path): + policy = tmp_path / "policy.csv" + policy.write_text(POLICY_CSV) + reader = OBPSCANReader.read_from_raw_data(tmp_path / "missing.csv", policy) + assert reader.obps_data is None + + def test_missing_policy_file_returns_none(self, tmp_path): + rates = tmp_path / "rates.csv" + rates.write_text(RATES_CSV) + reader = OBPSCANReader.read_from_raw_data(rates, tmp_path / "missing.csv") + assert reader.obps_data is None + + def test_loads_rates_and_policy(self, tmp_path): + rates = tmp_path / "rates.csv" + rates.write_text(RATES_CSV) + policy = tmp_path / "policy.csv" + policy.write_text(POLICY_CSV) + + reader = OBPSCANReader.read_from_raw_data(rates, policy) + + assert reader.obps_data is not None + assert list(reader.obps_data.df_rates.columns) == ["Date", "CAN", "CAN_BC"] + assert reader.obps_data.df_rates.loc[0, "CAN"] == pytest.approx(15.0) + assert list(reader.obps_data.df_policy["Industry"]) == ["C24a", "C19"] + assert reader.obps_data.df_policy.loc[0, "reduction_factor"] == pytest.approx(0.80) + + def test_optional_elec_file_loaded_when_present(self, tmp_path): + rates = tmp_path / "rates.csv" + rates.write_text(RATES_CSV) + policy = tmp_path / "policy.csv" + policy.write_text(POLICY_CSV) + elec = tmp_path / "elec.csv" + elec.write_text(ELEC_CSV) + + reader = OBPSCANReader.read_from_raw_data(rates, policy, policy_elec_path=elec) + + assert reader.obps_data.df_policy_elec is not None + assert list(reader.obps_data.df_policy_elec["Industry"]) == ["D01b"] + + def test_missing_optional_elec_file_gives_none(self, tmp_path): + rates = tmp_path / "rates.csv" + rates.write_text(RATES_CSV) + policy = tmp_path / "policy.csv" + policy.write_text(POLICY_CSV) + + reader = OBPSCANReader.read_from_raw_data(rates, policy, policy_elec_path=tmp_path / "missing.csv") + + assert reader.obps_data is not None + assert reader.obps_data.df_policy_elec is None diff --git a/tests/test_macromodel/unit/test_policy/__init__.py b/tests/test_macromodel/unit/test_policy/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/tests/test_macromodel/unit/test_policy/test_obps_can.py b/tests/test_macromodel/unit/test_policy/test_obps_can.py new file mode 100644 index 00000000..1321680d --- /dev/null +++ b/tests/test_macromodel/unit/test_policy/test_obps_can.py @@ -0,0 +1,181 @@ +import numpy as np +import pandas as pd +import pytest + +from macro_data.readers.policy_data.obps_can_reader import OBPSCANData +from macromodel.policy.output_based_price_system_can import OutputBasedPriceSystemCAN + +# Two industries: C24a is in the regulated list, A01 is not. +INDUSTRIES = ["C24a", "A01"] +CARBON_PRICE = 50.0 +REDUCTION_FACTOR = 0.8 +TIGHTENING_RATE = 0.02 +# reference intensity set directly on each fixture: 0.5 tCO2/unit +REFERENCE_INTENSITY = 0.5 + + +def _make_data( + carbon_price: float = CARBON_PRICE, + reduction_factor: float = REDUCTION_FACTOR, + tightening_rate: float = TIGHTENING_RATE, +) -> OBPSCANData: + df_rates = pd.DataFrame({"Date": list(range(2014, 2052)), "CAN": [carbon_price] * 38}) + df_policy = pd.DataFrame( + { + "Industry": ["C24a"], + "reduction_factor": [reduction_factor], + "tightening_rate": [tightening_rate], + } + ) + return OBPSCANData(df_rates=df_rates, df_policy=df_policy) + + +def _make_obps(year: int = 2020) -> OutputBasedPriceSystemCAN: + obps = OutputBasedPriceSystemCAN(country_name="CAN", industries=INDUSTRIES, obps_data=_make_data()) + obps.current_year = year + obps.current_t = year - 2014 + # pre-set reference intensity so tests are independent of the accumulation logic + obps.reference_emission_intensity[0] = REFERENCE_INTENSITY + return obps + + +def _call(obps, input_em, capital_em=None, production=None, record_reference=False): + if production is None: + production = np.array([100.0, 50.0]) + if capital_em is None: + capital_em = np.zeros(len(INDUSTRIES)) + return obps.compute_obps( + use_obps_reg=True, + record_obps_reference=record_reference, + production=production, + input_em=np.array(input_em), + capital_em=capital_em, + ) + + +class TestComputeOBPSDisabledAndEarlyYears: + def test_disabled_returns_all_zeros(self): + obps = _make_obps(year=2020) + cost = obps.compute_obps( + use_obps_reg=False, + record_obps_reference=False, + production=np.array([100.0, 50.0]), + input_em=np.array([999.0, 999.0]), + capital_em=np.zeros(2), + ) + assert np.all(cost == 0.0) + + def test_before_2019_returns_all_zeros(self): + obps = _make_obps(year=2017) + cost = _call(obps, input_em=[99.0, 99.0]) + assert np.all(cost == 0.0) + + def test_2018_returns_all_zeros(self): + obps = _make_obps(year=2018) + cost = _call(obps, input_em=[99.0, 99.0]) + assert np.all(cost == 0.0) + + +class TestComputeOBPSCosts: + # limit = production * reduction_factor * reference_intensity + # = 100 * 0.8 * 0.5 = 40.0 + + def test_positive_cost_when_above_limit(self): + # emissions = 60 > limit 40 → cost = (60-40)*50 = 1000 + obps = _make_obps(year=2020) + cost = _call(obps, input_em=[60.0, 0.0]) + assert cost[0] == pytest.approx(1000.0) + + def test_negative_cost_rebate_when_below_limit(self): + # emissions = 30 < limit 40 → cost = (30-40)*50 = -500 + obps = _make_obps(year=2020) + cost = _call(obps, input_em=[30.0, 0.0]) + assert cost[0] == pytest.approx(-500.0) + + def test_zero_cost_at_exact_limit(self): + # emissions == limit → cost = 0 + obps = _make_obps(year=2020) + cost = _call(obps, input_em=[40.0, 0.0]) + assert cost[0] == pytest.approx(0.0) + + def test_unregulated_industry_always_zero(self): + # A01 (index 1) is not in regulated_industries + obps = _make_obps(year=2020) + cost = _call(obps, input_em=[0.0, 999.0]) + assert cost[1] == pytest.approx(0.0) + + def test_zero_production_skipped(self): + obps = _make_obps(year=2020) + cost = _call(obps, input_em=[999.0, 0.0], production=np.array([0.0, 50.0])) + assert cost[0] == pytest.approx(0.0) + + def test_input_and_capital_emissions_summed(self): + # input_em=25, capital_em=20 → total=45, limit=40 → cost=(45-40)*50=250 + obps = _make_obps(year=2020) + cost = _call(obps, input_em=[25.0, 0.0], capital_em=np.array([20.0, 0.0])) + assert cost[0] == pytest.approx(250.0) + + def test_emission_limit_attribute_updated(self): + obps = _make_obps(year=2020) + _call(obps, input_em=[60.0, 0.0]) + assert obps.emission_limit[0] == pytest.approx(40.0) + + +class TestGetLimit: + def test_pre_2023_limit(self): + # limit = production * reduction_factor * reference_intensity = 100*0.8*0.5 = 40 + obps = _make_obps(year=2021) + limit = obps.get_limit(0, production=100.0) + assert limit == pytest.approx(40.0) + + def test_post_2022_tightening(self): + # 2024: limit = 100 * (0.4 - 0.4*0.02*(2024-2022)) = 100 * (0.4 - 0.016) = 38.4 + obps = _make_obps(year=2024) + limit = obps.get_limit(0, production=100.0) + assert limit == pytest.approx(38.4) + + def test_unknown_industry_returns_zero(self): + obps = _make_obps(year=2020) + # index 1 is A01, which has no row in df_policy + limit = obps.get_limit(1, production=100.0) + assert limit == pytest.approx(0.0) + + +class TestReferenceAccumulation: + def test_reference_intensity_computed_from_2017_to_2019(self): + obps = OutputBasedPriceSystemCAN( + country_name="CAN", industries=INDUSTRIES, obps_data=_make_data() + ) + # Simulate 2017, 2018, 2019 accumulation + for year, em, prod in [ + (2017, [30.0, 0.0], [100.0, 50.0]), + (2018, [30.0, 0.0], [100.0, 50.0]), + (2019, [30.0, 0.0], [100.0, 50.0]), + ]: + obps.current_year = year + obps.current_t = year - 2014 + obps.compute_obps( + use_obps_reg=True, + record_obps_reference=True, + production=np.array(prod), + input_em=np.array(em), + capital_em=np.zeros(2), + ) + # intensity = total_em / total_prod = 90 / 300 = 0.3 + assert obps.reference_emission_intensity[0] == pytest.approx(0.3) + assert obps.reference_emission_intensity[1] == pytest.approx(0.0) + + def test_reference_not_accumulated_when_flag_false(self): + obps = OutputBasedPriceSystemCAN( + country_name="CAN", industries=INDUSTRIES, obps_data=_make_data() + ) + obps.current_year = 2017 + obps.current_t = 3 + obps.compute_obps( + use_obps_reg=True, + record_obps_reference=False, + production=np.array([100.0, 50.0]), + input_em=np.array([99.0, 0.0]), + capital_em=np.zeros(2), + ) + assert np.all(obps.reference_emission == 0.0) From d2a3dba846dc1a745dc0f2cb4619727e2fc471fc Mon Sep 17 00:00:00 2001 From: reetiksahu Date: Wed, 20 May 2026 22:15:00 -0700 Subject: [PATCH 5/9] test: add unit tests for OBPS policy and reader Covers OutputBasedPriceSystemCAN (positive cost, negative rebate, tightening rate, reference accumulation, disabled/pre-2019 guards) and OBPSCANReader (CSV loading, missing files, optional elec file). Co-Authored-By: Claude Sonnet 4.6 --- tests/test_macromodel/unit/test_policy/test_obps_can.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/tests/test_macromodel/unit/test_policy/test_obps_can.py b/tests/test_macromodel/unit/test_policy/test_obps_can.py index 1321680d..a54f95db 100644 --- a/tests/test_macromodel/unit/test_policy/test_obps_can.py +++ b/tests/test_macromodel/unit/test_policy/test_obps_can.py @@ -143,9 +143,7 @@ def test_unknown_industry_returns_zero(self): class TestReferenceAccumulation: def test_reference_intensity_computed_from_2017_to_2019(self): - obps = OutputBasedPriceSystemCAN( - country_name="CAN", industries=INDUSTRIES, obps_data=_make_data() - ) + obps = OutputBasedPriceSystemCAN(country_name="CAN", industries=INDUSTRIES, obps_data=_make_data()) # Simulate 2017, 2018, 2019 accumulation for year, em, prod in [ (2017, [30.0, 0.0], [100.0, 50.0]), @@ -166,9 +164,7 @@ def test_reference_intensity_computed_from_2017_to_2019(self): assert obps.reference_emission_intensity[1] == pytest.approx(0.0) def test_reference_not_accumulated_when_flag_false(self): - obps = OutputBasedPriceSystemCAN( - country_name="CAN", industries=INDUSTRIES, obps_data=_make_data() - ) + obps = OutputBasedPriceSystemCAN(country_name="CAN", industries=INDUSTRIES, obps_data=_make_data()) obps.current_year = 2017 obps.current_t = 3 obps.compute_obps( From 87b65d6117f3e3f89aa697c7e3731ea66d0a93b7 Mon Sep 17 00:00:00 2001 From: reetiksahu Date: Fri, 22 May 2026 10:58:46 -0700 Subject: [PATCH 6/9] fix: address PR #94 review comments on OBPS implementation - Warn when use_obps_reg=True but no OBPS object is set (country.py) - Floor extra_marginal_taxes_firm at -good_prices so a negative rebate cannot push the effective sector price below zero (country.py) - Rename reference_emission_intensity -> baseline_emission_intensity to make clear the value is fixed from the 2017-2019 reference period - Add initial_year field so reset() and price-loading no longer hardcode 2014; current_year remains the mutable simulation state - Warn on init when regulated industries are absent from the model's industry set, with details on which codes were skipped - Update tests to reflect the baseline_emission_intensity rename --- macromodel/country/country.py | 11 ++++++ .../policy/output_based_price_system_can.py | 39 ++++++++++++++----- .../unit/test_policy/test_obps_can.py | 6 +-- 3 files changed, 44 insertions(+), 12 deletions(-) diff --git a/macromodel/country/country.py b/macromodel/country/country.py index 7dad018d..959a3a3e 100644 --- a/macromodel/country/country.py +++ b/macromodel/country/country.py @@ -564,6 +564,12 @@ def update_extra_taxes(self, record_obps_reference: bool = True) -> None: """ self.extra_marginal_taxes_firm = np.zeros(self.firms.n_industries) + if self.use_obps_reg and self.obps is None: + logging.warning( + "use_obps_reg is True for %s but no OBPS object is set — OBPS has no effect.", + self.country_name, + ) + if self.use_obps_reg and self.obps is not None: sectoral_tax = self.obps.compute_obps( use_obps_reg=self.use_obps_reg, @@ -578,6 +584,11 @@ def update_extra_taxes(self, record_obps_reference: bool = True) -> None: out=np.zeros_like(sectoral_tax), where=self.firms.ts.current("production") != 0, ) + # A negative rebate cannot bring the effective sector price below zero. + self.extra_marginal_taxes_firm = np.maximum( + -self.economy.ts.current("good_prices"), + self.extra_marginal_taxes_firm, + ) def clear_labour_market(self) -> None: """Execute labor market clearing. diff --git a/macromodel/policy/output_based_price_system_can.py b/macromodel/policy/output_based_price_system_can.py index e6288494..80fb0c9f 100644 --- a/macromodel/policy/output_based_price_system_can.py +++ b/macromodel/policy/output_based_price_system_can.py @@ -11,6 +11,7 @@ target-setting phase. """ +import logging from dataclasses import dataclass import numpy as np @@ -37,8 +38,8 @@ class OutputBasedPriceSystemCAN: df_policy: Per-industry reduction factors and tightening rates. df_policy_elec: Electricity-specific tightening rates (optional). df_rates: Carbon price schedule by year and jurisdiction. - reference_emission_intensity: Baseline emission intensity recorded - over the 2017–2019 reference period. + baseline_emission_intensity: Baseline emission intensity recorded + over the 2017–2019 reference period and fixed thereafter. reference_emission: Cumulative emissions during the reference period. reference_production: Cumulative production during the reference period. emission_limit: Current period allowable emissions per industry. @@ -54,11 +55,12 @@ class OutputBasedPriceSystemCAN: df_policy: pd.DataFrame df_policy_elec: pd.DataFrame df_rates: pd.DataFrame - reference_emission_intensity: np.ndarray + baseline_emission_intensity: np.ndarray reference_emission: np.ndarray reference_production: np.ndarray emission_limit: np.ndarray price: np.ndarray + initial_year: int = 2014 current_t: int = 0 current_year: int = 2014 @@ -99,8 +101,27 @@ def __init__(self, country_name: str, industries: list[str], obps_data: OBPSCAND [list(industries).index(ind) for ind in self.regulated_industries if ind in industries] ) + skipped = [ind for ind in self.regulated_industries if ind not in industries] + if skipped: + logging.warning( + "OBPS (%s): %d regulated industrie(s) not found in model and will be skipped: %s", + country_name, + len(skipped), + skipped, + ) + if len(self.regulated_indices) == 0: + logging.warning( + "OBPS (%s): none of the regulated industries appear in the model — OBPS will have no effect. " + "Check that industry codes match.", + country_name, + ) + + self.initial_year = 2014 + self.current_t = 0 + self.current_year = self.initial_year + n = len(industries) - self.reference_emission_intensity = np.zeros(n) + self.baseline_emission_intensity = np.zeros(n) self.reference_emission = np.zeros(n) self.reference_production = np.zeros(n) self.emission_limit = np.zeros(n) @@ -112,7 +133,7 @@ def __init__(self, country_name: str, industries: list[str], obps_data: OBPSCAND self.price = np.zeros(len(self.df_rates)) df_sub = self.df_rates[["Date", self.country_name]] for t in range(len(self.df_rates)): - df_row = df_sub[df_sub["Date"] == t + 2014] + df_row = df_sub[df_sub["Date"] == t + self.initial_year] self.price[t] = df_row[self.country_name].values[0] def compute_obps( @@ -149,10 +170,10 @@ def compute_obps( self.reference_production += production if self.current_year == 2019: - self.reference_emission_intensity = np.divide( + self.baseline_emission_intensity = np.divide( self.reference_emission, self.reference_production, - out=np.zeros_like(self.reference_emission), + out=np.zeros_like(self.baseline_emission_intensity), where=self.reference_production != 0, ) @@ -188,7 +209,7 @@ def get_limit(self, industry_idx: int, production: float) -> float: return 0.0 reduction_factor = row["reduction_factor"].values[0] - B = reduction_factor * self.reference_emission_intensity[industry_idx] + B = reduction_factor * self.baseline_emission_intensity[industry_idx] if self.current_year < 2023: return production * B @@ -208,4 +229,4 @@ def update(self) -> None: def reset(self) -> None: """Reset time variables to the initial year.""" self.current_t = 0 - self.current_year = 2014 + self.current_year = self.initial_year diff --git a/tests/test_macromodel/unit/test_policy/test_obps_can.py b/tests/test_macromodel/unit/test_policy/test_obps_can.py index a54f95db..62ea02d0 100644 --- a/tests/test_macromodel/unit/test_policy/test_obps_can.py +++ b/tests/test_macromodel/unit/test_policy/test_obps_can.py @@ -35,7 +35,7 @@ def _make_obps(year: int = 2020) -> OutputBasedPriceSystemCAN: obps.current_year = year obps.current_t = year - 2014 # pre-set reference intensity so tests are independent of the accumulation logic - obps.reference_emission_intensity[0] = REFERENCE_INTENSITY + obps.baseline_emission_intensity[0] = REFERENCE_INTENSITY return obps @@ -160,8 +160,8 @@ def test_reference_intensity_computed_from_2017_to_2019(self): capital_em=np.zeros(2), ) # intensity = total_em / total_prod = 90 / 300 = 0.3 - assert obps.reference_emission_intensity[0] == pytest.approx(0.3) - assert obps.reference_emission_intensity[1] == pytest.approx(0.0) + assert obps.baseline_emission_intensity[0] == pytest.approx(0.3) + assert obps.baseline_emission_intensity[1] == pytest.approx(0.0) def test_reference_not_accumulated_when_flag_false(self): obps = OutputBasedPriceSystemCAN(country_name="CAN", industries=INDUSTRIES, obps_data=_make_data()) From 87ac33ba446a5c210b0df17a29402c20d3f521f0 Mon Sep 17 00:00:00 2001 From: reetiksahu Date: Fri, 22 May 2026 15:54:06 -0700 Subject: [PATCH 7/9] fix: floor OBPS emission limit at zero to prevent negative limits The post-2022 tightening formula (B - B * tightening_rate * (self.current_year - 2022) becomes negativewhere B = reduction_factor * self.baseline_emission_intensity[industry_idx]once tightening accumulates past the baseline, causing the allowable emission limit to go below zero (~2042 onward for a 2% annual rate). Clamp the returned limit to max(0, ...) so regulated sectors never receive a pathologically negative allowance. Add a test covering the far-future tightening case. --- macromodel/policy/output_based_price_system_can.py | 2 +- tests/test_macromodel/unit/test_policy/test_obps_can.py | 6 ++++++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/macromodel/policy/output_based_price_system_can.py b/macromodel/policy/output_based_price_system_can.py index 80fb0c9f..5b9e4354 100644 --- a/macromodel/policy/output_based_price_system_can.py +++ b/macromodel/policy/output_based_price_system_can.py @@ -215,7 +215,7 @@ def get_limit(self, industry_idx: int, production: float) -> float: return production * B tightening_rate = row["tightening_rate"].values[0] - return production * (B - B * tightening_rate * (self.current_year - 2022)) + return max(0.0, production * (B - B * tightening_rate * (self.current_year - 2022))) def get_price(self) -> float: """Return the current period carbon price ($/tCO₂e).""" diff --git a/tests/test_macromodel/unit/test_policy/test_obps_can.py b/tests/test_macromodel/unit/test_policy/test_obps_can.py index 62ea02d0..2f6c7753 100644 --- a/tests/test_macromodel/unit/test_policy/test_obps_can.py +++ b/tests/test_macromodel/unit/test_policy/test_obps_can.py @@ -140,6 +140,12 @@ def test_unknown_industry_returns_zero(self): limit = obps.get_limit(1, production=100.0) assert limit == pytest.approx(0.0) + def test_limit_floored_at_zero_when_tightening_exceeds_baseline(self): + # 2022 + 50 years of tightening at 2%/yr → multiplier goes negative; floor clamps to 0 + obps = _make_obps(year=2072) + limit = obps.get_limit(0, production=100.0) + assert limit == pytest.approx(0.0) + class TestReferenceAccumulation: def test_reference_intensity_computed_from_2017_to_2019(self): From ea56671934d68ab50a511df64c7ed795fff603f4 Mon Sep 17 00:00:00 2001 From: reetiksahu Date: Tue, 23 Jun 2026 09:22:50 -0700 Subject: [PATCH 8/9] fix: address PR #94 review comments on OBPS implementation This commit systematically addresses all critical and major review comments: **Critical Issues Fixed:** 1. Sparse-year price schedule crash: Replace index-based price array with year-based dictionary lookup to handle non-contiguous years (2014, 2019, 2030). Use forward-filling for intermediate years. 2. Year-boundary double-counting: Move reference accumulation window from 2017-2019 to 2017-2018, with limit computation in 2019. Prevents 2019 from being counted in both baseline and levy periods. 3. Missing reset state: Implement reset() method to clear all accumulators (reference_emission, reference_production, baseline_emission_intensity, emission_limit) for multi-run workflows. **Major Issues Fixed:** 4. CH4 field nullability: Add None guards in country.update_extra_taxes() when inputs_emissions_ch4 and capital_emissions_ch4 are not configured. Substitute zeros to prevent array addition errors. 5. Cost-push sign ambiguity: Include tax in unit costs calculation so positive OBPS levy correctly raises prices instead of suppressing them. Tax now shifts both numerator and denominator equally. 6. Price lookup bounds: Consolidated into new get_price() method that handles sparse schedules robustly via sorted year lookup. **Testing:** - Add comprehensive sparse-year integration tests matching real Canada OBPS schedule (2014, 2019, 2030) - Test forward-fill behavior for intermediate years - Test reset functionality for multi-run scenarios - Update existing tests to reflect 2017-2018 reference period - All 25 OBPS tests pass --- macromodel/agents/firms/func/prices.py | 8 +- macromodel/country/country.py | 11 +- .../policy/output_based_price_system_can.py | 42 +++++--- .../unit/test_policy/test_obps_can.py | 100 +++++++++++++++++- 4 files changed, 141 insertions(+), 20 deletions(-) diff --git a/macromodel/agents/firms/func/prices.py b/macromodel/agents/firms/func/prices.py index db95fd70..41c3f235 100644 --- a/macromodel/agents/firms/func/prices.py +++ b/macromodel/agents/firms/func/prices.py @@ -172,6 +172,7 @@ def compute_price( extra_marginal_taxes if extra_marginal_taxes is not None else np.zeros_like(prev_average_good_prices) ) average_price_by_firm = (prev_average_good_prices + tax_by_sector)[current_firm_sectors] + tax_by_firm = tax_by_sector[current_firm_sectors] # Demand-pull inflation demand_pull_inflation = np.zeros_like(prev_firm_prices) @@ -196,12 +197,13 @@ def compute_price( ) demand_pull_inflation = np.maximum(min_inflation, np.minimum(max_inflation, demand_pull_inflation)) - # Cost-push inflation + # Cost-push inflation: include the tax in unit costs so positive tax raises prices + total_unit_costs = curr_unit_costs + tax_by_firm cost_push_inflation = ( np.divide( - curr_unit_costs, + total_unit_costs, average_price_by_firm, - out=np.ones_like(curr_unit_costs), + out=np.ones_like(total_unit_costs), where=average_price_by_firm != 0.0, ) - 1.0 diff --git a/macromodel/country/country.py b/macromodel/country/country.py index 959a3a3e..a219aa6e 100644 --- a/macromodel/country/country.py +++ b/macromodel/country/country.py @@ -571,12 +571,19 @@ def update_extra_taxes(self, record_obps_reference: bool = True) -> None: ) if self.use_obps_reg and self.obps is not None: + input_em_ch4 = self.firms.ts.current("inputs_emissions_ch4") + if input_em_ch4 is None: + input_em_ch4 = np.zeros_like(self.firms.ts.current("inputs_emissions")) + capital_em_ch4 = self.firms.ts.current("capital_emissions_ch4") + if capital_em_ch4 is None: + capital_em_ch4 = np.zeros_like(self.firms.ts.current("capital_emissions")) + sectoral_tax = self.obps.compute_obps( use_obps_reg=self.use_obps_reg, record_obps_reference=record_obps_reference, production=self.firms.ts.current("production"), - input_em=self.firms.ts.current("inputs_emissions") + self.firms.ts.current("inputs_emissions_ch4"), - capital_em=self.firms.ts.current("capital_emissions") + self.firms.ts.current("capital_emissions_ch4"), + input_em=self.firms.ts.current("inputs_emissions") + input_em_ch4, + capital_em=self.firms.ts.current("capital_emissions") + capital_em_ch4, ) self.extra_marginal_taxes_firm = np.divide( sectoral_tax, diff --git a/macromodel/policy/output_based_price_system_can.py b/macromodel/policy/output_based_price_system_can.py index 5b9e4354..880f7f01 100644 --- a/macromodel/policy/output_based_price_system_can.py +++ b/macromodel/policy/output_based_price_system_can.py @@ -26,7 +26,7 @@ class OutputBasedPriceSystemCAN: Calculates the tax that regulated firms pay on emissions above a prescribed output-weighted limit. Reference emission intensities are - recorded during 2017–2019 and used to set sector-specific limits from + recorded during 2017–2018 and used to set sector-specific limits from 2019 onwards. Attributes: @@ -130,11 +130,7 @@ def __init__(self, country_name: str, industries: list[str], obps_data: OBPSCAND self.df_policy_elec = obps_data.df_policy_elec if obps_data.df_policy_elec is not None else pd.DataFrame() self.df_rates = obps_data.df_rates - self.price = np.zeros(len(self.df_rates)) - df_sub = self.df_rates[["Date", self.country_name]] - for t in range(len(self.df_rates)): - df_row = df_sub[df_sub["Date"] == t + self.initial_year] - self.price[t] = df_row[self.country_name].values[0] + self.price_by_year = dict(zip(self.df_rates["Date"], self.df_rates[self.country_name])) def compute_obps( self, @@ -165,11 +161,11 @@ def compute_obps( if not use_obps_reg: return np.zeros(len(self.industries)) - if record_obps_reference and self.current_year in (2017, 2018, 2019): + if record_obps_reference and self.current_year in (2017, 2018): self.reference_emission += input_em + capital_em self.reference_production += production - if self.current_year == 2019: + if self.current_year == 2018: self.baseline_emission_intensity = np.divide( self.reference_emission, self.reference_production, @@ -181,12 +177,13 @@ def compute_obps( return np.zeros(len(self.industries)) obps_cost = np.zeros(len(self.industries)) + current_price = self.get_price() for i in self.regulated_indices: if production[i] > 0: limit = self.get_limit(i, production[i]) self.emission_limit[i] = limit difference = (input_em[i] + capital_em[i]) - limit - obps_cost[i] = difference * self.price[min(self.current_t, len(self.price) - 1)] + obps_cost[i] = difference * current_price return obps_cost @@ -218,8 +215,22 @@ def get_limit(self, industry_idx: int, production: float) -> float: return max(0.0, production * (B - B * tightening_rate * (self.current_year - 2022))) def get_price(self) -> float: - """Return the current period carbon price ($/tCO₂e).""" - return self.price[self.current_t] + """Return the current period carbon price ($/tCO₂e). + + For years not explicitly in the schedule, uses forward-filling: returns + the price from the last milestone year less than or equal to current_year. + """ + year = self.current_year + available_years = sorted(self.price_by_year.keys()) + + if year in self.price_by_year: + return self.price_by_year[year] + + past_years = [y for y in available_years if y <= year] + if past_years: + return self.price_by_year[max(past_years)] + + return self.price_by_year[min(available_years)] def update(self) -> None: """Advance the timestep by one annual period.""" @@ -227,6 +238,13 @@ def update(self) -> None: self.current_year += 1 def reset(self) -> None: - """Reset time variables to the initial year.""" + """Reset all internal state for a fresh simulation run. + + Resets time tracking, accumulators, and computed baselines to initial values. + """ self.current_t = 0 self.current_year = self.initial_year + self.reference_emission = np.zeros(len(self.industries)) + self.reference_production = np.zeros(len(self.industries)) + self.baseline_emission_intensity = np.zeros(len(self.industries)) + self.emission_limit = np.zeros(len(self.industries)) diff --git a/tests/test_macromodel/unit/test_policy/test_obps_can.py b/tests/test_macromodel/unit/test_policy/test_obps_can.py index 2f6c7753..3435395c 100644 --- a/tests/test_macromodel/unit/test_policy/test_obps_can.py +++ b/tests/test_macromodel/unit/test_policy/test_obps_can.py @@ -148,9 +148,9 @@ def test_limit_floored_at_zero_when_tightening_exceeds_baseline(self): class TestReferenceAccumulation: - def test_reference_intensity_computed_from_2017_to_2019(self): + def test_reference_intensity_computed_from_2017_to_2018(self): obps = OutputBasedPriceSystemCAN(country_name="CAN", industries=INDUSTRIES, obps_data=_make_data()) - # Simulate 2017, 2018, 2019 accumulation + # Simulate 2017, 2018, 2019 accumulation (reference recorded in 2017-2018 only) for year, em, prod in [ (2017, [30.0, 0.0], [100.0, 50.0]), (2018, [30.0, 0.0], [100.0, 50.0]), @@ -165,7 +165,7 @@ def test_reference_intensity_computed_from_2017_to_2019(self): input_em=np.array(em), capital_em=np.zeros(2), ) - # intensity = total_em / total_prod = 90 / 300 = 0.3 + # intensity = total_em / total_prod = 60 / 200 = 0.3 (accumulated during 2017-2018) assert obps.baseline_emission_intensity[0] == pytest.approx(0.3) assert obps.baseline_emission_intensity[1] == pytest.approx(0.0) @@ -181,3 +181,97 @@ def test_reference_not_accumulated_when_flag_false(self): capital_em=np.zeros(2), ) assert np.all(obps.reference_emission == 0.0) + + +class TestSparseYearSchedule: + """Test OBPS with sparse (non-contiguous) year schedules matching real Canada data.""" + + def _make_sparse_data(self) -> OBPSCANData: + """Create a realistic sparse carbon price schedule: 2014, 2019, 2030.""" + df_rates = pd.DataFrame( + { + "Date": [2014, 2019, 2030], + "CAN": [15.0, 50.0, 170.0], + } + ) + df_policy = pd.DataFrame( + { + "Industry": ["C24a"], + "reduction_factor": [0.8], + "tightening_rate": [0.02], + } + ) + return OBPSCANData(df_rates=df_rates, df_policy=df_policy) + + def test_sparse_schedule_initializes_without_error(self): + """Sparse schedule should not crash during initialization.""" + obps = OutputBasedPriceSystemCAN(country_name="CAN", industries=INDUSTRIES, obps_data=self._make_sparse_data()) + obps.baseline_emission_intensity[0] = 0.5 + assert obps.get_price() == pytest.approx(15.0) + + def test_sparse_schedule_forward_fills_intermediate_years(self): + """Years between milestones should use the last milestone's price (forward-fill).""" + obps = OutputBasedPriceSystemCAN(country_name="CAN", industries=INDUSTRIES, obps_data=self._make_sparse_data()) + obps.baseline_emission_intensity[0] = 0.5 + + obps.current_year = 2014 + assert obps.get_price() == pytest.approx(15.0) + + obps.current_year = 2016 + assert obps.get_price() == pytest.approx(15.0) + + obps.current_year = 2019 + assert obps.get_price() == pytest.approx(50.0) + + obps.current_year = 2025 + assert obps.get_price() == pytest.approx(50.0) + + obps.current_year = 2030 + assert obps.get_price() == pytest.approx(170.0) + + obps.current_year = 2035 + assert obps.get_price() == pytest.approx(170.0) + + def test_sparse_schedule_obps_cost_with_milestone_prices(self): + """OBPS cost should use correct price at each milestone year.""" + obps = OutputBasedPriceSystemCAN(country_name="CAN", industries=INDUSTRIES, obps_data=self._make_sparse_data()) + obps.baseline_emission_intensity[0] = 0.5 + + # At 2019: price = 50, limit = 100*0.8*0.5 = 40, emissions = 60 → cost = 20 * 50 = 1000 + obps.current_year = 2019 + obps.current_t = 5 + cost = obps.compute_obps( + use_obps_reg=True, + record_obps_reference=False, + production=np.array([100.0, 50.0]), + input_em=np.array([60.0, 0.0]), + capital_em=np.zeros(2), + ) + assert cost[0] == pytest.approx(1000.0) + + # At 2030: price = 170, B = 0.4, limit = 100*(0.4 - 0.4*0.02*8) = 33.6, emissions = 60 + # → cost = 26.4 * 170 = 4488 + obps.current_year = 2030 + obps.current_t = 16 + cost = obps.compute_obps( + use_obps_reg=True, + record_obps_reference=False, + production=np.array([100.0, 50.0]), + input_em=np.array([60.0, 0.0]), + capital_em=np.zeros(2), + ) + assert cost[0] == pytest.approx(4488.0) + + def test_sparse_schedule_with_reset(self): + """Reset should clear accumulated state on sparse schedule.""" + obps = OutputBasedPriceSystemCAN(country_name="CAN", industries=INDUSTRIES, obps_data=self._make_sparse_data()) + + obps.current_year = 2030 + obps.baseline_emission_intensity[0] = 0.5 + obps.reference_emission[0] = 100.0 + + obps.reset() + + assert obps.current_year == 2014 + assert obps.baseline_emission_intensity[0] == 0.0 + assert obps.reference_emission[0] == 0.0 From dec3d60e9adb86bbc95ab70829c08174f3c5c3ff Mon Sep 17 00:00:00 2001 From: reetiksahu Date: Tue, 23 Jun 2026 09:27:27 -0700 Subject: [PATCH 9/9] fix: add extra_marginal_taxes parameter to ExogenousPriceSetter for interface consistency Ensures ExogenousPriceSetter.compute_price() has the same signature as other PriceSetter subclasses, preventing TypeError when OBPS combines with exogenous pricing setters. --- macromodel/agents/firms/func/prices.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/macromodel/agents/firms/func/prices.py b/macromodel/agents/firms/func/prices.py index 41c3f235..d7e5d7db 100644 --- a/macromodel/agents/firms/func/prices.py +++ b/macromodel/agents/firms/func/prices.py @@ -388,6 +388,7 @@ def compute_price( current_time: int, min_inflation: float = -0.1, max_inflation: float = 0.1, + extra_marginal_taxes: Optional[np.ndarray] = None, ) -> np.ndarray: """Set prices according to exogenous PPI path. @@ -400,6 +401,8 @@ def compute_price( current_time (int): Current period index min_inflation (float, optional): Unused. Defaults to -0.1. max_inflation (float, optional): Unused. Defaults to 0.1. + extra_marginal_taxes (np.ndarray, optional): Unused. Accepts for + interface consistency. Defaults to None. Returns: np.ndarray: Price level from exogenous PPI path