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..f6446f0a 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 @@ -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..d7e5d7db 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,18 @@ 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] + tax_by_firm = tax_by_sector[current_firm_sectors] # Demand-pull inflation demand_pull_inflation = np.zeros_like(prev_firm_prices) @@ -187,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 @@ -276,6 +287,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 +329,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 +341,16 @@ 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 @@ -371,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. @@ -383,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 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..a219aa6e 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,53 @@ 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 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: + 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") + input_em_ch4, + capital_em=self.firms.ts.current("capital_emissions") + capital_em_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, + ) + # 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. @@ -555,6 +617,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 +734,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 +752,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..880f7f01 --- /dev/null +++ b/macromodel/policy/output_based_price_system_can.py @@ -0,0 +1,250 @@ +"""Canada Output-Based Pricing System (OBPS) policy for the macroeconomic model. + +Computes a two-way carbon price signal for each regulated sector based on +emissions relative to an output-based benchmark: + + obps_cost[i] = (emissions[i] - limit[i]) * carbon_price[t] + +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. +""" + +import logging +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–2018 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. + 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. + 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 + 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 + + 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] + ) + + 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.baseline_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_by_year = dict(zip(self.df_rates["Date"], self.df_rates[self.country_name])) + + 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 signed cost of emissions relative to 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: 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)) + + 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 == 2018: + self.baseline_emission_intensity = np.divide( + self.reference_emission, + self.reference_production, + out=np.zeros_like(self.baseline_emission_intensity), + where=self.reference_production != 0, + ) + + if self.current_year < 2019: + 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 * current_price + + 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.baseline_emission_intensity[industry_idx] + + if self.current_year < 2023: + return production * B + + tightening_rate = row["tightening_rate"].values[0] + 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). + + 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.""" + self.current_t += 1 + self.current_year += 1 + + def reset(self) -> None: + """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/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, 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..3435395c --- /dev/null +++ b/tests/test_macromodel/unit/test_policy/test_obps_can.py @@ -0,0 +1,277 @@ +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.baseline_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) + + 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_2018(self): + obps = OutputBasedPriceSystemCAN(country_name="CAN", industries=INDUSTRIES, obps_data=_make_data()) + # 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]), + (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 = 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) + + 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) + + +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