diff --git a/.github/workflows/build_macrocalib.yaml b/.github/workflows/build_macrocalib.yaml
index 35017ac1..3f71c196 100644
--- a/.github/workflows/build_macrocalib.yaml
+++ b/.github/workflows/build_macrocalib.yaml
@@ -2,9 +2,11 @@ name: Build and test calibration package
 
 on:
   pull_request:
-    # Run on PRs to any branch
+    branches:
+      - main
   push:
-    branches: [ main ]
+    branches:
+      - main
 
 jobs:
   build:
diff --git a/.github/workflows/build_macrodata.yml b/.github/workflows/build_macrodata.yml
index acd93596..0e2af7af 100644
--- a/.github/workflows/build_macrodata.yml
+++ b/.github/workflows/build_macrodata.yml
@@ -2,7 +2,8 @@ name: Build and test data package
 
 on:
   pull_request:
-    # Run on PRs to any branch
+    branches:
+      - main
   push:
     branches:
       - main
diff --git a/.github/workflows/build_macromodel.yml b/.github/workflows/build_macromodel.yml
index cbfc3e2f..7d57ee20 100644
--- a/.github/workflows/build_macromodel.yml
+++ b/.github/workflows/build_macromodel.yml
@@ -2,7 +2,8 @@ name: Build and test model package
 
 on:
   pull_request:
-    # Run on PRs to any branch
+    branches:
+      - main
   push:
     branches:
       - main
diff --git a/macromodel/agents/firms/firms.py b/macromodel/agents/firms/firms.py
index 40a98cc8..4eff8e3e 100644
--- a/macromodel/agents/firms/firms.py
+++ b/macromodel/agents/firms/firms.py
@@ -495,24 +495,29 @@ def set_targets(
             estimated_inflation: Expected inflation rate
             current_good_prices (np.ndarray): Industry-level average prices
         """
-        self.ts.limiting_intermediate_inputs.append(
-            self.functions["production"].compute_limiting_intermediate_inputs_stock(
-                intermediate_inputs_productivity_matrix=self.get_effective_intermediate_coefficients(),
-                intermediate_inputs_stock=self.ts.current("intermediate_inputs_stock"),
-                intermediate_inputs_utilisation_rate=self.intermediate_inputs_utilisation_rate,
-                goods_criticality_matrix=self.goods_criticality_matrix,
-                substitution_bundle_matrix=self.substitution_bundles,
-            )
+        # Calculate limiting inputs and apply TFP multiplier.
+        # TFP scales effective capacity from inputs, so limits should also scale with TFP.
+        # This ensures target-setting accounts for TFP-enhanced production capacity.
+        # Without this, production targets are capped at initial levels even when TFP grows.
+        tfp_multiplier = self.states["tfp_multiplier"]
+
+        limiting_intermediate = self.functions["production"].compute_limiting_intermediate_inputs_stock(
+            intermediate_inputs_productivity_matrix=self.get_effective_intermediate_coefficients(),
+            intermediate_inputs_stock=self.ts.current("intermediate_inputs_stock"),
+            intermediate_inputs_utilisation_rate=self.intermediate_inputs_utilisation_rate,
+            goods_criticality_matrix=self.goods_criticality_matrix,
+            substitution_bundle_matrix=self.substitution_bundles,
         )
-        self.ts.limiting_capital_inputs.append(
-            self.functions["production"].compute_limiting_capital_inputs_stock(
-                capital_inputs_productivity_matrix=self.get_effective_capital_coefficients(),
-                capital_inputs_stock=self.ts.current("capital_inputs_stock"),
-                capital_inputs_utilisation_rate=self.capital_inputs_utilisation_rate,
-                goods_criticality_matrix=self.goods_criticality_matrix,
-                substitution_bundle_matrix=self.substitution_bundles,
-            )
+        self.ts.limiting_intermediate_inputs.append(limiting_intermediate * tfp_multiplier)
+
+        limiting_capital = self.functions["production"].compute_limiting_capital_inputs_stock(
+            capital_inputs_productivity_matrix=self.get_effective_capital_coefficients(),
+            capital_inputs_stock=self.ts.current("capital_inputs_stock"),
+            capital_inputs_utilisation_rate=self.capital_inputs_utilisation_rate,
+            goods_criticality_matrix=self.goods_criticality_matrix,
+            substitution_bundle_matrix=self.substitution_bundles,
         )
+        self.ts.limiting_capital_inputs.append(limiting_capital * tfp_multiplier)
         self.ts.target_production.append(
             self.compute_target_production(
                 bank_overdraft_rate_on_firm_deposits=bank_overdraft_rate_on_firm_deposits,
@@ -590,6 +595,7 @@ def plan_productivity_investment(
             input_usage=self.ts.current("used_intermediate_inputs"),
             current_tech_multipliers=self.states["intermediate_tech_multipliers"],
             substitution_bundle_matrix=self.substitution_bundles,
+            firm_industries=self.states["Industry"],
         )
 
         return total_investment, tfp_investment, technical_investment
@@ -1267,18 +1273,60 @@ def prepare_buying_goods(
 
         # Target capital inputs
         if assume_zero_growth:
-            self.ts.target_capital_inputs.append(self.ts.initial("target_capital_inputs"))
+            target_capital_inputs = self.ts.initial("target_capital_inputs")
         else:
-            self.ts.target_capital_inputs.append(
-                self.functions["target_capital_inputs"].compute_target_capital_inputs(
-                    unconstrained_target_capital_inputs=self.ts.current("unconstrained_target_capital_inputs"),
-                    target_long_term_credit=self.ts.current("target_long_term_credit"),
-                    received_long_term_credit=self.ts.current("received_long_term_credit"),
-                    previous_good_prices=previous_good_prices,
-                    expected_inflation=expected_inflation,
-                )
+            target_capital_inputs = self.functions["target_capital_inputs"].compute_target_capital_inputs(
+                unconstrained_target_capital_inputs=self.ts.current("unconstrained_target_capital_inputs"),
+                target_long_term_credit=self.ts.current("target_long_term_credit"),
+                received_long_term_credit=self.ts.current("received_long_term_credit"),
+                previous_good_prices=previous_good_prices,
+                expected_inflation=expected_inflation,
             )
 
+        # [TFP_DEBUG] Add planned productivity investment as additional capital demand.
+        # This makes productivity investment a forward-looking target rather than a residual.
+        if len(self.ts.planned_productivity_investment) > 0:
+            planned_investment = self.ts.current("planned_productivity_investment")
+            if planned_investment is not None and np.any(planned_investment > 0):
+                import logging
+
+                logger = logging.getLogger(__name__)
+                logger.debug(
+                    f"[TFP_DEBUG] prepare_buying_goods: Adding planned_productivity_investment "
+                    f"(sum={np.sum(planned_investment):.2e}) to capital demand"
+                )
+
+                target_capital_inputs = target_capital_inputs.copy()
+                expected_prices = (1 + expected_inflation) * previous_good_prices
+                safe_prices = np.maximum(expected_prices, 1e-12)
+
+                # Distribute productivity investment across industries based on existing capital weights
+                target_totals = target_capital_inputs.sum(axis=1, keepdims=True)
+                weights = np.divide(
+                    target_capital_inputs,
+                    target_totals,
+                    out=np.zeros_like(target_capital_inputs),
+                    where=target_totals > 0,
+                )
+                # For firms with no existing capital targets, distribute evenly
+                if np.any(target_totals == 0):
+                    weights = np.where(
+                        target_totals == 0,
+                        1.0 / self.n_industries,
+                        weights,
+                    )
+
+                # Convert monetary investment to real quantities at expected prices
+                extra_capital_qty = (planned_investment[:, None] * weights) / safe_prices
+                target_capital_inputs += extra_capital_qty
+
+                logger.debug(
+                    f"[TFP_DEBUG] prepare_buying_goods: Added {np.sum(extra_capital_qty):.2e} "
+                    f"real units of capital from productivity investment"
+                )
+
+        self.ts.target_capital_inputs.append(target_capital_inputs)
+
         # Setting total real amount of goods to buy
         self.set_goods_to_buy(self.ts.current("target_intermediate_inputs") + self.ts.current("target_capital_inputs"))
 
diff --git a/macromodel/agents/firms/func/production.py b/macromodel/agents/firms/func/production.py
index 28fcb4d6..10b7390a 100644
--- a/macromodel/agents/firms/func/production.py
+++ b/macromodel/agents/firms/func/production.py
@@ -50,21 +50,19 @@ def compute_production(
         Returns:
             np.ndarray: Feasible production levels by firm
         """
+        # Limiting stock is already TFP-scaled in set_targets()
         limiting_stock = self.compute_limiting_stock(
             current_limiting_intermediate_inputs,
             current_limiting_capital_inputs,
         )
 
-        # Apply TFP multiplier if provided
+        # Apply TFP multiplier to labour only (limiting stock is pre-scaled)
         if tfp_multiplier is not None:
-            # TFP scales effective capacity from inputs
             effective_labour = current_labour_inputs * tfp_multiplier
-            effective_limiting_stock = limiting_stock * tfp_multiplier
         else:
             effective_labour = current_labour_inputs
-            effective_limiting_stock = limiting_stock
 
-        return np.amin([desired_production, effective_labour, effective_limiting_stock], axis=0)
+        return np.amin([desired_production, effective_labour, limiting_stock], axis=0)
 
     @abstractmethod
     def compute_limiting_intermediate_inputs_stock(
diff --git a/macromodel/agents/firms/func/productivity_growth.py b/macromodel/agents/firms/func/productivity_growth.py
index dee7ed69..382e5a62 100644
--- a/macromodel/agents/firms/func/productivity_growth.py
+++ b/macromodel/agents/firms/func/productivity_growth.py
@@ -1,7 +1,10 @@
+import logging
 from abc import ABC, abstractmethod
 
 import numpy as np
 
+logger = logging.getLogger(__name__)
+
 
 class ProductivityGrowth(ABC):
     """Abstract base class for computing Total Factor Productivity (TFP) growth.
@@ -134,6 +137,19 @@ def compute_tfp_growth(
         # Base growth applies to all firms
         tfp_growth = np.full_like(current_tfp, base_growth_rate)
 
+        # Handle empty arrays early
+        if len(current_tfp) == 0:
+            return tfp_growth
+
+        # [TFP_DEBUG] Log input parameters
+        logger.debug(
+            "[TFP_DEBUG] SimpleTFPGrowth.compute_tfp_growth called: "
+            f"base_growth_rate={base_growth_rate}, investment_elasticity={investment_elasticity}, "
+            f"investment_effectiveness={self.investment_effectiveness}, "
+            f"current_tfp(mean)={np.mean(current_tfp):.6f}, production(sum)={np.sum(production):.2e}, "
+            f"productivity_investment(sum)={np.sum(productivity_investment):.2e}"
+        )
+
         # Add investment-driven growth where production > 0
         positive_production = production > 0
         if np.any(positive_production):
@@ -146,6 +162,12 @@ def compute_tfp_growth(
             positive_investment = productivity_investment > 0
             valid_firms = positive_production & positive_investment
 
+            # [TFP_DEBUG] Log firm counts
+            logger.debug(
+                f"[TFP_DEBUG] Firms with positive_production={np.sum(positive_production)}, "
+                f"positive_investment={np.sum(positive_investment)}, valid_firms={np.sum(valid_firms)}"
+            )
+
             if np.any(valid_firms):
                 investment_intensity[valid_firms] = productivity_investment[valid_firms] / production[valid_firms]
 
@@ -156,6 +178,18 @@ def compute_tfp_growth(
 
                 tfp_growth += investment_contribution
 
+                # [TFP_DEBUG] Log investment contribution
+                logger.debug(
+                    f"[TFP_DEBUG] investment_intensity(mean of valid)={np.mean(investment_intensity[valid_firms]):.6f}, "
+                    f"investment_contribution(mean)={np.mean(investment_contribution):.6f}, "
+                    f"tfp_growth(mean)={np.mean(tfp_growth):.6f}"
+                )
+
+        # [TFP_DEBUG] Log final result
+        logger.debug(
+            f"[TFP_DEBUG] Final tfp_growth: mean={np.mean(tfp_growth):.6f}, min={np.min(tfp_growth):.6f}, max={np.max(tfp_growth):.6f}"
+        )
+
         return tfp_growth
 
 
diff --git a/macromodel/agents/firms/func/productivity_investment_planner.py b/macromodel/agents/firms/func/productivity_investment_planner.py
index e29c50d3..8c275b21 100644
--- a/macromodel/agents/firms/func/productivity_investment_planner.py
+++ b/macromodel/agents/firms/func/productivity_investment_planner.py
@@ -100,6 +100,7 @@ def plan_productivity_investment(
         input_usage: np.ndarray,
         current_tech_multipliers: np.ndarray,
         substitution_bundle_matrix: np.ndarray,
+        firm_industries: Optional[np.ndarray] = None,
     ) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
         """Plan productivity investment for each firm.
 
@@ -109,14 +110,16 @@ def plan_productivity_investment(
 
         Args:
             current_tfp (np.ndarray): Current TFP multipliers [n_firms]
-            current_production (np.ndarray): Current production levels [n_firms]
-            current_unit_costs (np.ndarray): Current unit costs of production [n_firms]
-            available_cash (np.ndarray): Cash available for investment [n_firms]
-            current_prices (np.ndarray): Current market prices [n_industries]
+            current_production (np.ndarray): Current production levels [n_firms] (physical units)
+            current_unit_costs (np.ndarray): Current unit costs of production [n_firms] ($/unit)
+            available_cash (np.ndarray): Cash available for investment [n_firms] ($)
+            current_prices (np.ndarray): Current market prices [n_industries] ($/unit)
             n_industries (int): Number of industries
             input_usage (np.ndarray): Input usage by firms [n_firms x n_industries]
             current_tech_multipliers (np.ndarray): Current technical multipliers [n_firms x n_industries]
             substitution_bundle_matrix (np.ndarray): Bundle matrix [n_industries x n_bundles]
+            firm_industries (np.ndarray | None): Industry index for each firm [n_firms], used to
+                get firm-specific prices from industry prices
 
         Returns:
             tuple[np.ndarray, np.ndarray, np.ndarray]: Tuple containing:
@@ -184,6 +187,64 @@ def compute_hurdle_adjusted_value(
 
         return cost_reduction_per_period * hurdle_discount_factor
 
+    def is_investment_profitable(
+        self,
+        productivity_investment: np.ndarray,
+        nominal_production: np.ndarray,
+        current_unit_costs: np.ndarray,
+    ) -> np.ndarray:
+        """Check if investment is profitable using the correct hurdle rate formula.
+
+        Following investment_decision.md, investment is profitable if:
+            (I/V)^(1-α) ≤ φ * (1+r_h)/r_h
+
+        where V = nominal production (price × physical production).
+
+        This is derived from the condition that the present value of cost savings
+        must exceed the investment cost. The derivation:
+        1. TFP growth: g = φ * (I/V)^α
+        2. Cost savings per period: g * V (where V is in monetary terms)
+        3. PV of savings: g * V * (1+r_h)/r_h
+        4. Profitable if: PV > I
+        5. Which rearranges to: (I/V)^(1-α) ≤ φ * (1+r_h)/r_h
+
+        Note: Using nominal production V (price × quantity) ensures dimensional
+        correctness. Both I and V are in monetary units ($), so I/V is dimensionless.
+
+        Args:
+            productivity_investment (np.ndarray): Investment amounts ($)
+            nominal_production (np.ndarray): Nominal production = price × quantity ($)
+            current_unit_costs (np.ndarray): Unit costs ($/unit), used to check
+                if cost savings are possible
+
+        Returns:
+            np.ndarray: Boolean array indicating profitability for each firm
+        """
+        # Compute investment intensity I/V (dimensionless: $ / $)
+        investment_intensity = np.divide(
+            productivity_investment,
+            nominal_production,
+            out=np.zeros_like(productivity_investment),
+            where=nominal_production > 0,
+        )
+
+        # Left-hand side: (I/V)^(1-α)
+        # Need to handle the case where investment_intensity is 0
+        lhs = np.power(
+            investment_intensity,
+            1 - self.investment_elasticity,
+            out=np.zeros_like(investment_intensity),
+            where=investment_intensity > 0,
+        )
+
+        # Right-hand side: φ * (1+r_h)/r_h (dimensionless)
+        rhs = self.investment_effectiveness * (1 + self.hurdle_rate) / self.hurdle_rate
+
+        # Investment is profitable if LHS ≤ RHS
+        # Also return True for zero investment (trivially profitable)
+        # If unit costs are 0, no investment is profitable (no cost savings possible)
+        return ((lhs <= rhs) & (current_unit_costs > 0)) | (productivity_investment <= 0)
+
     def compute_investment_budget(
         self,
         available_cash: np.ndarray,
@@ -461,6 +522,7 @@ def plan_productivity_investment(
         input_usage: np.ndarray,
         current_tech_multipliers: np.ndarray,
         substitution_bundle_matrix: Optional[np.ndarray],
+        firm_industries: Optional[np.ndarray] = None,
     ) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
         """Return zero productivity investment for all firms.
 
@@ -474,6 +536,7 @@ def plan_productivity_investment(
             input_usage (np.ndarray): Used intermediate inputs [n_firms x n_industries]
             current_tech_multipliers (np.ndarray): Technical coefficient multipliers [n_firms x n_industries]
             substitution_bundle_matrix (np.ndarray | None): Substitution bundle matrix
+            firm_industries (np.ndarray | None): Industry index for each firm
 
         Returns:
             tuple: Zero investments (total, TFP, technical by input)
@@ -546,6 +609,7 @@ def plan_productivity_investment(
         input_usage: np.ndarray,
         current_tech_multipliers: np.ndarray,
         substitution_bundle_matrix: Optional[np.ndarray],
+        firm_industries: Optional[np.ndarray] = None,
     ) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
         """Plan productivity investment using simple rules.
 
@@ -554,29 +618,39 @@ def plan_productivity_investment(
 
         Args:
             current_tfp (np.ndarray): Current TFP multipliers
-            current_production (np.ndarray): Current production levels
-            current_unit_costs (np.ndarray): Current unit costs of production
-            available_cash (np.ndarray): Cash available for investment
-            current_prices (np.ndarray): Current market prices by industry
+            current_production (np.ndarray): Current production levels (physical units)
+            current_unit_costs (np.ndarray): Current unit costs of production ($/unit)
+            available_cash (np.ndarray): Cash available for investment ($)
+            current_prices (np.ndarray): Current market prices by industry ($/unit)
             n_industries (int): Number of industries
             input_usage (np.ndarray): Used intermediate inputs [n_firms x n_industries]
             current_tech_multipliers (np.ndarray): Technical coefficient multipliers [n_firms x n_industries]
             substitution_bundle_matrix (np.ndarray | None): Substitution bundle matrix
+            firm_industries (np.ndarray | None): Industry index for each firm
 
         Returns:
             tuple: (total_investment, tfp_investment, technical_investment_by_input)
         """
-        # Compute available budget
-        budget = self.compute_investment_budget(available_cash, current_production)
+        # Compute firm-specific prices from industry prices
+        if firm_industries is not None:
+            firm_prices = current_prices[firm_industries]
+        else:
+            # Fallback: assume prices are 1 (unit costs will be used for profitability)
+            firm_prices = np.ones_like(current_production)
+
+        # Compute nominal production (price × physical production) for dimensional correctness
+        # The formula I/Y should use monetary values for both I and Y
+        nominal_production = firm_prices * current_production
+
+        # Compute available budget using nominal production
+        budget = self.compute_investment_budget(available_cash, nominal_production)
 
         # Candidate investment (fraction of budget)
         candidate_investment = self.investment_propensity * budget
 
-        # Compute hurdle-adjusted present value of cost savings
-        hurdle_value = self.compute_hurdle_adjusted_value(candidate_investment, current_production, current_unit_costs)
-
-        # Investment is profitable if hurdle-adjusted value exceeds investment cost
-        profitable = hurdle_value > candidate_investment
+        # Check if investment is profitable using the correct formula from investment_decision.md
+        # The condition is: (I/V)^(1-α) ≤ φ * (1+r_h)/r_h, where V = nominal production
+        profitable = self.is_investment_profitable(candidate_investment, nominal_production, current_unit_costs)
 
         # Only invest where profitable
         total_investment = np.where(profitable, candidate_investment, 0)
@@ -661,6 +735,7 @@ def plan_productivity_investment(
         input_usage: np.ndarray,
         current_tech_multipliers: np.ndarray,
         substitution_bundle_matrix: Optional[np.ndarray],
+        firm_industries: Optional[np.ndarray] = None,
     ) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
         """Plan productivity investment by optimizing expected returns.
 
@@ -669,14 +744,15 @@ def plan_productivity_investment(
 
         Args:
             current_tfp (np.ndarray): Current TFP multipliers
-            current_production (np.ndarray): Current production levels
-            current_unit_costs (np.ndarray): Current unit costs of production
-            available_cash (np.ndarray): Cash available for investment
-            current_prices (np.ndarray): Current market prices by industry
+            current_production (np.ndarray): Current production levels (physical units)
+            current_unit_costs (np.ndarray): Current unit costs of production ($/unit)
+            available_cash (np.ndarray): Cash available for investment ($)
+            current_prices (np.ndarray): Current market prices by industry ($/unit)
             n_industries (int): Number of industries
             input_usage (np.ndarray): Used intermediate inputs [n_firms x n_industries]
             current_tech_multipliers (np.ndarray): Technical coefficient multipliers [n_firms x n_industries]
             substitution_bundle_matrix (np.ndarray | None): Substitution bundle matrix
+            firm_industries (np.ndarray | None): Industry index for each firm
 
         Returns:
             tuple: (total_investment, tfp_investment, technical_investment_by_input)
@@ -684,41 +760,37 @@ def plan_productivity_investment(
         n_firms = len(current_production)
         total_investment = np.zeros(n_firms)
 
-        # Compute constraints
-        budget = self.compute_investment_budget(available_cash, current_production)
-
-        # For each firm, find optimal investment level
-        for i in range(n_firms):
-            if budget[i] <= 0 or current_production[i] <= 0 or current_unit_costs[i] <= 0:
-                continue
-
-            # Search over possible investment levels
-            # Start with zero investment as baseline (NPV = 0)
-            best_npv = 0
-            best_investment = 0
-
-            investment_levels = np.linspace(0, budget[i], self.search_steps)
+        # Compute firm-specific prices from industry prices
+        if firm_industries is not None:
+            firm_prices = current_prices[firm_industries]
+        else:
+            firm_prices = np.ones(n_firms)
 
-            for investment in investment_levels:
-                if investment == 0:
-                    continue
+        # Compute nominal production for dimensional correctness
+        nominal_production = firm_prices * current_production
 
-                # Compute hurdle-adjusted value for this investment level
-                inv_array = np.array([investment])
-                prod_array = np.array([current_production[i]])
-                costs_array = np.array([current_unit_costs[i]])
+        # Compute constraints using nominal production
+        budget = self.compute_investment_budget(available_cash, nominal_production)
 
-                hurdle_value = self.compute_hurdle_adjusted_value(inv_array, prod_array, costs_array)[0]
+        # For each firm, find optimal investment level using the correct formula
+        # From investment_decision.md: investment is profitable if (I/V)^(1-α) ≤ φ * (1+r_h)/r_h
+        # where V = nominal production (price × physical production)
+        for i in range(n_firms):
+            if budget[i] <= 0 or nominal_production[i] <= 0:
+                continue
 
-                # NPV: hurdle-adjusted value minus investment cost
-                npv = hurdle_value - investment
+            # Compute the maximum profitable investment intensity
+            # From (I/V)^(1-α) ≤ φ * (1+r_h)/r_h, solving for I/V:
+            # I/V ≤ (φ * (1+r_h)/r_h)^(1/(1-α))
+            rhs = self.investment_effectiveness[i] * (1 + self.hurdle_rate[i]) / self.hurdle_rate[i]
+            exponent = 1.0 / (1 - self.investment_elasticity[i])
+            max_profitable_intensity = np.power(rhs, exponent)
 
-                # Only invest if NPV is positive and better than current best
-                if npv > best_npv:
-                    best_npv = npv
-                    best_investment = investment
+            # Maximum profitable investment (now dimensionally correct: $ = ratio × $)
+            max_profitable_investment = max_profitable_intensity * nominal_production[i]
 
-            total_investment[i] = best_investment
+            # Investment is capped by budget and profitability
+            total_investment[i] = min(budget[i], max_profitable_investment)
 
         if substitution_bundle_matrix is not None:
             # Allocate between TFP and technical using bundle-aware logic
diff --git a/macromodel/util/get_histogram.py b/macromodel/util/get_histogram.py
index 879771bd..d0a2fb32 100644
--- a/macromodel/util/get_histogram.py
+++ b/macromodel/util/get_histogram.py
@@ -57,12 +57,53 @@ def get_histogram(values: np.ndarray, scale: Optional[int], bins: int = 40, norm
             values = (values - np.min(values)) / diff
         else:
             values = values - np.min(values)
+
+    # [TFP_DEBUG] Apply scaling and handle edge cases
     if scale is None:
-        hist, bin_edges = np.histogram(values, bins=bins)
+        scaled_values = values
     else:
-        hist, bin_edges = np.histogram(values / scale, bins=bins)
+        scaled_values = values / scale
+
+    # Filter out non-finite values
+    scaled_values = scaled_values[np.isfinite(scaled_values)]
+    if len(scaled_values) == 0:
+        return np.full((2, bins + 1), np.nan)
+
+    min_value = np.min(scaled_values)
+    max_value = np.max(scaled_values)
+    if not np.isfinite(min_value) or not np.isfinite(max_value):
+        return np.full((2, bins + 1), np.nan)
+
+    # Handle edge cases: zero-width or subnormal ranges
+    data_range = max_value - min_value
+    hist_range = None
+    if data_range <= 0:
+        # All values are identical - create a small range around them
+        low = np.nextafter(min_value, -np.inf)
+        high = np.nextafter(max_value, np.inf)
+        if low == high:
+            pad = max(np.finfo(float).tiny * bins, 1e-12)
+            low = min_value - pad
+            high = max_value + pad
+        hist_range = (low, high)
+    elif data_range / bins < np.finfo(float).tiny:
+        # Range is too small for the number of bins
+        pad = max(np.finfo(float).tiny * bins, 1e-12)
+        hist_range = (min_value - pad, max_value + pad)
+
+    try:
+        hist, bin_edges = np.histogram(scaled_values, bins=bins, range=hist_range)
+    except ValueError:
+        # Fallback to single bin if histogram fails
+        hist, bin_edges = np.histogram(scaled_values, bins=1, range=hist_range)
+        # Pad to expected size
+        hist = np.concatenate([hist, np.zeros(bins - 1)])
+        bin_edges = np.linspace(bin_edges[0], bin_edges[-1], bins + 1)
+
     hist = hist.astype(float)
-    hist /= hist.sum()
+    hist_sum = hist.sum()
+    if hist_sum > 0:
+        hist /= hist_sum
     return np.array([np.concatenate((hist, [np.nan])), bin_edges])
 
 
diff --git a/tests/test_macromodel/unit/test_agents/test_firms/func/test_production.py b/tests/test_macromodel/unit/test_agents/test_firms/func/test_production.py
index 6eb8be9c..3f091647 100644
--- a/tests/test_macromodel/unit/test_agents/test_firms/func/test_production.py
+++ b/tests/test_macromodel/unit/test_agents/test_firms/func/test_production.py
@@ -48,27 +48,40 @@ def test__compute_production_with_tfp_unity(self):
         )
 
     def test__compute_production_with_tfp_boost(self):
-        """Test that TFP>1.0 increases effective capacity."""
+        """Test that TFP>1.0 scales labour in compute_production().
+
+        Note: Limiting stock (intermediate + capital inputs) is pre-scaled by TFP
+        in set_targets(), so compute_production() only scales labour here.
+        This test passes pre-scaled limiting inputs to reflect real usage.
+        """
         tfp_boost = np.array([1.5, 1.2])  # 50% and 20% productivity boost
+        # Simulate pre-scaled limiting inputs (as done in set_targets())
+        # Raw limiting intermediate = [9, 11], raw limiting capital = [6, 8]
+        pre_scaled_limiting_intermediate = np.array([9.0, 11.0]) * tfp_boost  # [13.5, 13.2]
+        pre_scaled_limiting_capital = np.array([6.0, 8.0]) * tfp_boost  # [9.0, 9.6]
         result = PureLeontief().compute_production(
             desired_production=np.array([10.0, 10.0]),
             current_labour_inputs=np.array([9.0, 11.0]),
-            current_limiting_intermediate_inputs=np.array([9.0, 11.0]),
-            current_limiting_capital_inputs=np.array([6.0, 8.0]),
+            current_limiting_intermediate_inputs=pre_scaled_limiting_intermediate,
+            current_limiting_capital_inputs=pre_scaled_limiting_capital,
             tfp_multiplier=tfp_boost,
         )
-        # TFP scales inputs: min([10, 9*1.5, 9*1.5, 6*1.5]) = min([10, 13.5, 13.5, 9]) = 9
-        # TFP scales inputs: min([10, 11*1.2, 11*1.2, 8*1.2]) = min([10, 13.2, 13.2, 9.6]) = 9.6
+        # Labour is scaled in compute_production: 9*1.5=13.5, 11*1.2=13.2
+        # Limiting stock is already pre-scaled: min(intermediate, capital) = [9.0, 9.6]
+        # Final: min([10, 13.5, 9.0], [10, 13.2, 9.6]) = [9.0, 9.6]
         assert np.allclose(result, np.array([9.0, 9.6]))
 
     def test__compute_production_tfp_respects_target(self):
         """Test that TFP doesn't allow production above target."""
         tfp_high = np.array([2.0, 2.0])  # Double productivity
+        # Simulate pre-scaled limiting inputs (as done in set_targets())
+        pre_scaled_limiting_intermediate = np.array([9.0, 11.0]) * tfp_high  # [18.0, 22.0]
+        pre_scaled_limiting_capital = np.array([6.0, 8.0]) * tfp_high  # [12.0, 16.0]
         result = PureLeontief().compute_production(
             desired_production=np.array([5.0, 7.0]),  # Low targets
             current_labour_inputs=np.array([9.0, 11.0]),
-            current_limiting_intermediate_inputs=np.array([9.0, 11.0]),
-            current_limiting_capital_inputs=np.array([6.0, 8.0]),
+            current_limiting_intermediate_inputs=pre_scaled_limiting_intermediate,
+            current_limiting_capital_inputs=pre_scaled_limiting_capital,
             tfp_multiplier=tfp_high,
         )
         # Despite high TFP, production is limited by desired_production
diff --git a/tests/test_macromodel/unit/test_agents/test_firms/func/test_productivity_investment_planner.py b/tests/test_macromodel/unit/test_agents/test_firms/func/test_productivity_investment_planner.py
index 22c7b5a3..31bdc353 100644
--- a/tests/test_macromodel/unit/test_agents/test_firms/func/test_productivity_investment_planner.py
+++ b/tests/test_macromodel/unit/test_agents/test_firms/func/test_productivity_investment_planner.py
@@ -758,3 +758,172 @@ def test_parameter_validation_wrong_length(self):
                 n_firms=3,
                 hurdle_rate=[0.10, 0.15],  # Only 2 values for 3 firms
             )
+
+
+class TestHurdleRateSensitivity:
+    """Test that hurdle rate properly affects investment profitability.
+
+    This addresses bug #60: hurdle_rate insensitivity.
+    The fix uses the correct formula from investment_decision.md:
+    (I/V)^(1-α) ≤ φ * (1+r_h)/r_h, where V = nominal production
+    """
+
+    def test_high_hurdle_rate_rejects_investment(self):
+        """Test that a high hurdle rate rejects investment that lower rate accepts."""
+        n_firms = 1
+
+        # Use identical conditions for both planners except hurdle rate
+        common_params = {
+            "n_firms": n_firms,
+            "max_investment_fraction": 0.5,
+            "investment_effectiveness": 0.1,
+            "investment_elasticity": 0.3,
+            "investment_propensity": 1.0,
+        }
+
+        low_hurdle_planner = SimpleProductivityInvestmentPlanner(hurdle_rate=0.05, **common_params)
+        high_hurdle_planner = SimpleProductivityInvestmentPlanner(hurdle_rate=0.5, **common_params)
+
+        # Setup test scenario
+        current_tfp = np.array([1.0])
+        current_production = np.array([1000.0])  # Physical units
+        current_unit_costs = np.array([10.0])  # $/unit
+        available_cash = np.array([5000.0])  # $
+
+        test_params = {
+            "current_prices": np.ones(18),
+            "n_industries": 18,
+            "input_usage": np.ones((n_firms, 18)),
+            "current_tech_multipliers": np.ones((n_firms, 18)),
+            "substitution_bundle_matrix": None,
+            "firm_industries": np.array([0]),
+        }
+
+        low_inv, _, _ = low_hurdle_planner.plan_productivity_investment(
+            current_tfp=current_tfp,
+            current_production=current_production,
+            current_unit_costs=current_unit_costs,
+            available_cash=available_cash,
+            **test_params,
+        )
+
+        high_inv, _, _ = high_hurdle_planner.plan_productivity_investment(
+            current_tfp=current_tfp,
+            current_production=current_production,
+            current_unit_costs=current_unit_costs,
+            available_cash=available_cash,
+            **test_params,
+        )
+
+        # Low hurdle rate should allow investment, high should reject
+        assert low_inv[0] > 0, "Low hurdle rate should allow investment"
+        assert high_inv[0] < low_inv[0], "High hurdle rate should reduce/reject investment"
+
+    def test_hurdle_rate_threshold_formula(self):
+        """Test that the hurdle rate formula correctly computes profitability threshold.
+
+        The formula: (I/V)^(1-α) ≤ φ * (1+r_h)/r_h
+        Solving for I/V: I/V ≤ (φ * (1+r_h)/r_h)^(1/(1-α))
+        """
+        n_firms = 1
+        planner = SimpleProductivityInvestmentPlanner(
+            n_firms=n_firms,
+            hurdle_rate=0.15,
+            investment_effectiveness=0.1,
+            investment_elasticity=0.3,
+        )
+
+        # Expected max profitable intensity
+        # rhs = 0.1 * (1 + 0.15) / 0.15 = 0.1 * 7.667 = 0.7667
+        # max_intensity = 0.7667^(1/(1-0.3)) = 0.7667^1.4286 ≈ 0.68
+        rhs = 0.1 * (1 + 0.15) / 0.15
+        expected_max_intensity = np.power(rhs, 1 / (1 - 0.3))
+
+        # Test investment just below threshold should be profitable
+        nominal_production = np.array([1000.0])  # $
+        unit_costs = np.array([10.0])  # > 0 so savings possible
+        investment_below = np.array([expected_max_intensity * nominal_production[0] * 0.9])
+        investment_above = np.array([expected_max_intensity * nominal_production[0] * 1.1])
+
+        profitable_below = planner.is_investment_profitable(investment_below, nominal_production, unit_costs)
+        profitable_above = planner.is_investment_profitable(investment_above, nominal_production, unit_costs)
+
+        assert profitable_below[0], "Investment below threshold should be profitable"
+        assert not profitable_above[0], "Investment above threshold should not be profitable"
+
+    def test_zero_unit_costs_prevents_investment(self):
+        """Test that zero unit costs prevent investment (no cost savings possible)."""
+        n_firms = 1
+        planner = SimpleProductivityInvestmentPlanner(
+            n_firms=n_firms,
+            hurdle_rate=0.05,  # Very low hurdle rate
+            investment_propensity=1.0,
+        )
+
+        current_tfp = np.array([1.0])
+        current_production = np.array([1000.0])
+        current_unit_costs = np.array([0.0])  # Zero unit costs
+        available_cash = np.array([5000.0])
+
+        test_params = {
+            "current_prices": np.ones(18),
+            "n_industries": 18,
+            "input_usage": np.ones((n_firms, 18)),
+            "current_tech_multipliers": np.ones((n_firms, 18)),
+            "substitution_bundle_matrix": None,
+            "firm_industries": np.array([0]),
+        }
+
+        total_inv, _, _ = planner.plan_productivity_investment(
+            current_tfp=current_tfp,
+            current_production=current_production,
+            current_unit_costs=current_unit_costs,
+            available_cash=available_cash,
+            **test_params,
+        )
+
+        # Zero unit costs means no cost savings possible, so no investment
+        assert total_inv[0] == 0, "Zero unit costs should prevent investment"
+
+    def test_dimensional_correctness_with_prices(self):
+        """Test that using firm prices correctly computes nominal production."""
+        n_firms = 2
+        planner = SimpleProductivityInvestmentPlanner(
+            n_firms=n_firms,
+            hurdle_rate=0.15,
+            investment_propensity=1.0,
+        )
+
+        current_tfp = np.array([1.0, 1.0])
+        current_production = np.array([1000.0, 1000.0])  # Same physical production
+        current_unit_costs = np.array([10.0, 10.0])
+        available_cash = np.array([50000.0, 50000.0])
+
+        # Different prices for each firm's industry
+        current_prices = np.full(18, 1.0)
+        current_prices[0] = 1.0  # Firm 0's industry
+        current_prices[1] = 10.0  # Firm 1's industry - higher price
+
+        test_params = {
+            "current_prices": current_prices,
+            "n_industries": 18,
+            "input_usage": np.ones((n_firms, 18)),
+            "current_tech_multipliers": np.ones((n_firms, 18)),
+            "substitution_bundle_matrix": None,
+            "firm_industries": np.array([0, 1]),  # Different industries
+        }
+
+        total_inv, _, _ = planner.plan_productivity_investment(
+            current_tfp=current_tfp,
+            current_production=current_production,
+            current_unit_costs=current_unit_costs,
+            available_cash=available_cash,
+            **test_params,
+        )
+
+        # Both should produce valid investments (no NaN/inf)
+        assert np.all(np.isfinite(total_inv)), "Investment should be finite"
+        assert np.all(total_inv >= 0), "Investment should be non-negative"
+        # Higher-priced firm has higher nominal production, so can invest more
+        # (budget is capped by max_investment_fraction * nominal_production)
+        assert total_inv[1] >= total_inv[0], "Higher-priced firm should have higher/equal investment budget"
diff --git a/tests/test_macromodel/unit/test_simulation.py b/tests/test_macromodel/unit/test_simulation.py
index cc83a789..8f6427f6 100644
--- a/tests/test_macromodel/unit/test_simulation.py
+++ b/tests/test_macromodel/unit/test_simulation.py
@@ -696,6 +696,106 @@ def test_technical_only_investment_allocation(datawrapper, seed=42):
     ), "At least some effective coefficients should be strictly better than base coefficients"
 
 
+def test_tfp_growth_produces_real_gdp_growth(can_disagg_datawrapper):
+    """Test that TFP growth with productivity investment produces real GDP growth.
+
+    This test verifies the fix for the issue where TFP growth was configured but
+    real GDP remained flat. The fix ensures that planned_productivity_investment
+    is translated into actual capital demand in prepare_buying_goods.
+
+    Addresses concerns raised by Luke and Deven:
+    - Luke: "GDP does appear to grow (nominally)... but real GDP (GDP/CPI) is flat"
+    - Deven: Fixed by adding planned_productivity_investment to target_capital_inputs
+
+    This test runs a 20-quarter (5-year) simulation and verifies:
+    1. TFP multiplier increases over time
+    2. Real GDP grows (not just nominal)
+    3. Productivity investment is planned and executed
+    """
+    n_industries = can_disagg_datawrapper.n_industries
+
+    configuration = SimulationConfiguration(
+        seed=42,
+        t_max=20,
+        country_configurations={"CAN": CountryConfiguration.n_industry_default(n_industries=n_industries)},
+    )
+
+    # Configure TFP growth with investment
+    configuration.country_configurations["CAN"].firms.functions.productivity_growth.name = "SimpleTFPGrowth"
+    configuration.country_configurations["CAN"].firms.parameters.tfp_base_growth_rate = 0.002  # 0.2% per quarter
+    configuration.country_configurations["CAN"].firms.parameters.tfp_investment_elasticity = 0.5
+    configuration.country_configurations["CAN"].firms.functions.productivity_growth.parameters = {
+        "investment_effectiveness": 0.5,
+    }
+
+    # Configure productivity investment planner with favorable parameters
+    configuration.country_configurations["CAN"].firms.functions.productivity_investment_planner.name = (
+        "SimpleProductivityInvestmentPlanner"
+    )
+    configuration.country_configurations["CAN"].firms.functions.productivity_investment_planner.parameters = {
+        "n_firms": n_industries,
+        "hurdle_rate": 0.03,
+        "investment_effectiveness": 0.2,
+        "investment_elasticity": 0.5,
+        "max_investment_fraction": 0.2,
+        "technical_diminishing_returns": 0.07,
+    }
+
+    # Use stable prices (no price feedback) to isolate TFP effect
+    configuration.country_configurations["CAN"].firms.functions.prices.parameters["price_setting_speed_dp"] = 0.0
+    configuration.country_configurations["CAN"].firms.functions.prices.parameters["price_setting_speed_cp"] = 0.0
+
+    # Create simulation
+    simulation = Simulation.from_datawrapper(datawrapper=can_disagg_datawrapper, simulation_configuration=configuration)
+
+    country = simulation.countries["CAN"]
+
+    # Record initial values
+    initial_tfp = country.firms.states["tfp_multiplier"].mean()
+    initial_gdp = country.economy.gdp_output()[-1]
+    initial_cpi = country.economy.total_cpi_inflation()[-1]
+    initial_real_gdp = initial_gdp / initial_cpi
+
+    # Run simulation
+    n_periods = 20
+    for _ in range(n_periods):
+        simulation.iterate()
+
+    # Record final values
+    final_tfp = country.firms.states["tfp_multiplier"].mean()
+    final_gdp = country.economy.gdp_output()[-1]
+    final_cpi = country.economy.total_cpi_inflation()[-1]
+    final_real_gdp = final_gdp / final_cpi
+
+    # Calculate growth rates
+    tfp_growth = (final_tfp / initial_tfp - 1) * 100
+    real_gdp_growth = (final_real_gdp / initial_real_gdp - 1) * 100
+    cpi_change = (final_cpi / initial_cpi - 1) * 100
+
+    # Assertions
+    # 1. TFP should have grown significantly
+    assert final_tfp > initial_tfp, f"TFP should increase: {initial_tfp:.4f} -> {final_tfp:.4f}"
+    assert tfp_growth > 10, f"TFP growth should be > 10%, got {tfp_growth:.2f}%"
+
+    # 2. Real GDP should have grown (this was the original issue - it was flat)
+    assert (
+        final_real_gdp > initial_real_gdp
+    ), f"Real GDP should increase: {initial_real_gdp:.2e} -> {final_real_gdp:.2e}"
+    assert real_gdp_growth > 5, f"Real GDP growth should be > 5%, got {real_gdp_growth:.2f}%"
+
+    # 3. CPI should remain relatively stable (with price_setting_speed = 0)
+    assert abs(cpi_change) < 20, f"CPI change should be < 20%, got {cpi_change:.2f}%"
+
+    # 4. Productivity investment should have occurred
+    if len(country.firms.ts.executed_productivity_investment) > 0:
+        total_executed_investment = sum(inv.sum() for inv in country.firms.ts.executed_productivity_investment)
+        assert total_executed_investment > 0, "There should be positive executed productivity investment"
+
+    if len(country.firms.ts.planned_productivity_investment) > 0:
+        total_planned_investment = sum(inv.sum() for inv in country.firms.ts.planned_productivity_investment)
+        assert total_planned_investment > 0, "There should be positive planned productivity investment"
+
+
 def test_prehooks(datawrapper):
     """Test that pre-hooks execute correctly before each iteration."""
     # Track hook calls