ENH: combine demand- and supply-side market-specific computations

- Reduced Market initialization overhead.
- Re-use probabilities and their derivatives for supply-side calculations.

docs/api.rst

@@ -359,6 +359,7 @@ When errors occur, they will either be displayed as warnings or raised as except
    exceptions.NonpositiveSyntheticCostsError
    exceptions.InvalidParameterCovariancesError
    exceptions.InvalidMomentCovariancesError
+   exceptions.GenericNumericalError
    exceptions.DeltaNumericalError
    exceptions.CostsNumericalError
    exceptions.MicroMomentsNumericalError
@@ -372,6 +373,7 @@ When errors occur, they will either be displayed as warnings or raised as except
    exceptions.SyntheticCostsNumericalError
    exceptions.SyntheticMicroMomentsNumericalError
    exceptions.EquilibriumRealizationNumericalError
+   exceptions.JacobianRealizationNumericalError
    exceptions.PostEstimationNumericalError
    exceptions.AbsorptionError
    exceptions.ClippedSharesError

pyblp/exceptions.py

@@ -54,6 +54,10 @@ class InvalidMomentCovariancesError(Error):
     """Failed to compute a weighting matrix because of invalid estimated covariances of GMM moments."""
 
 
+class GenericNumericalError(NumericalError):
+    """Encountered a numerical error."""
+
+
 class DeltaNumericalError(NumericalError):
     r"""Encountered a numerical error when computing :math:`\delta`.
 
@@ -161,16 +165,10 @@ class EquilibriumRealizationNumericalError(NumericalError):
     """Encountered a numerical error when solving for a realization of equilibrium prices and shares."""
 
 
-class XiByThetaJacobianRealizationNumericalError(NumericalError):
+class JacobianRealizationNumericalError(NumericalError):
     r"""Encountered a numerical error when computing a realization of the Jacobian (holding :math:`\beta` fixed) of
-    :math:`\xi` (equivalently, of :math:`\delta`) with respect to :math:`\theta`.
-
-    """
-
-
-class OmegaByThetaJacobianRealizationNumericalError(NumericalError):
-    r"""Encountered a numerical error when computing a realization of the Jacobian (holding :math:`\gamma` fixed) of
-    :math:`\omega` (equivalently, of transformed marginal costs) with respect to :math:`\theta`.
+    :math:`\xi` (equivalently, of :math:`\delta`) or :math:`\omega` (equivalently, of transformed marginal costs)
+    with respect to :math:`\theta`.
 
     """
 

pyblp/markets/problem_market.py

@@ -14,16 +14,21 @@
 class ProblemMarket(Market):
     """A market underlying the BLP problem."""
 
-    def solve_demand(
-            self, delta: Array, last_delta: Array, moments: Moments, iteration: Iteration, fp_type: str,
-            shares_bounds: Bounds, compute_jacobians: bool, compute_micro_covariances: bool) -> (
-            Tuple[Array, Array, Array, Array, Array, Array, Array, Array, SolverStats, List[Error]]):
+    @NumericalErrorHandler(exceptions.GenericNumericalError)
+    def solve(
+            self, delta: Array, last_delta: Array, last_tilde_costs: Array, moments: Moments, iteration: Iteration,
+            fp_type: str, shares_bounds: Bounds, costs_bounds: Bounds, compute_jacobians: bool,
+            compute_micro_covariances: bool) -> (
+            Tuple[
+                Array, Array, Array, Array, Array, Array, Array, Array, SolverStats, Array, Array, Array, List[Error]
+            ]):
         """Compute the mean utility for this market that equates market shares to observed values by solving a fixed
         point problem. Then, if compute_jacobians is True, compute the Jacobian (holding beta fixed) of xi
-        (equivalently, of delta) with respect to theta. Finally, compute any contributions to micro moment values, their
-        Jacobian with respect to theta (if compute_jacobians is True), and, if compute_micro_covariances is True, their
-        covariances. Replace null elements in delta with their last values before computing micro moment contributions
-        and Jacobians.
+        (equivalently, of delta) with respect to theta. Compute any contributions to micro moment values, their Jacobian
+        with respect to theta (if compute_jacobians is True), and, if compute_micro_covariances is True, their
+        covariance contributions. If there is a supply side, compute transformed marginal costs for this market. Then,
+        if compute_jacobian is True, compute the Jacobian (holding gamma fixed) of omega (equivalently, of transformed
+        marginal costs) with respect to theta.
         """
         errors: List[Error] = []
 
@@ -36,16 +41,41 @@ def solve_demand(
         bad_delta_index = ~np.isfinite(delta)
         valid_delta[bad_delta_index] = last_delta[bad_delta_index]
 
-        # compute the Jacobian (keep derivatives of probabilities with respect to theta if there are micro moments)
+        # if needed, pre-compute probabilities
         probabilities = conditionals = None
+        if compute_jacobians or moments.MM > 0 or self.K3 > 0:
+            probabilities, conditionals = self.compute_probabilities(valid_delta)
+
+        # if needed, pre-compute derivatives of probabilities with respect to parameters
         probabilities_tangent_mapping: Dict[int, Array] = {}
+        conditionals_tangent_mapping: Dict[int, Array] = {}
+        if compute_jacobians:
+            assert probabilities is not None
+            for p, parameter in enumerate(self.parameters.unfixed):
+                probabilities_tangent, conditionals_tangent = self.compute_probabilities_by_parameter_tangent(
+                    parameter, probabilities, conditionals, valid_delta
+                )
+                probabilities_tangent_mapping[p] = probabilities_tangent
+                if self.K3 > 0:
+                    conditionals_tangent_mapping[p] = conditionals_tangent
+
+        # compute the Jacobian of xi (equivalently, of delta) with respect to theta
         xi_jacobian = np.full((self.J, self.parameters.P), np.nan, options.dtype)
         if compute_jacobians:
-            xi_jacobian, probabilities, conditionals, probabilities_tangent_mapping, xi_jacobian_errors = (
-                self.safely_compute_xi_by_theta_jacobian(valid_delta, keep_probabilities_tangents=moments.MM > 0)
+            assert probabilities is not None
+            xi_jacobian, xi_jacobian_errors = self.safely_compute_xi_by_theta_jacobian(
+                probabilities, conditionals, probabilities_tangent_mapping
             )
             errors.extend(xi_jacobian_errors)
 
+        # if needed, pre-compute derivatives of probabilities with respect to xi
+        probabilities_tensor = conditionals_tensor = None
+        if compute_jacobians and (moments.MM > 0 or self.K3 > 0):
+            assert probabilities is not None
+            probabilities_tensor, conditionals_tensor = self.compute_probabilities_by_xi_tensor(
+                probabilities, conditionals, compute_conditionals_tensor=self.K3 > 0
+            )
+
         # compute contributions to micro moments, their Jacobian, and their covariances
         if moments.MM == 0:
             micro_numerator = np.zeros((moments.MM, 0), options.dtype)
@@ -54,50 +84,53 @@ def solve_demand(
             micro_denominator_jacobian = np.full((moments.MM, self.parameters.P), np.nan, options.dtype)
             micro_covariances_numerator = np.full((moments.MM, moments.MM), np.nan, options.dtype)
         else:
+            assert probabilities is not None
             (
                 micro_numerator, micro_denominator, micro_numerator_jacobian, micro_denominator_jacobian,
                 micro_covariances_numerator, micro_errors
             ) = (
                 self.safely_compute_micro_contributions(
-                    moments, valid_delta, probabilities, conditionals, probabilities_tangent_mapping, xi_jacobian,
-                    compute_jacobians, compute_micro_covariances
+                    moments, valid_delta, probabilities, probabilities_tangent_mapping, probabilities_tensor,
+                    xi_jacobian, compute_jacobians, compute_micro_covariances
                 )
             )
             errors.extend(micro_errors)
 
-        return (
-            delta, xi_jacobian, micro_numerator, micro_denominator, micro_numerator_jacobian,
-            micro_denominator_jacobian, micro_covariances_numerator, clipped_shares, stats, errors
-        )
-
-    def solve_supply(
-            self, last_tilde_costs: Array, xi_jacobian: Array, costs_bounds: Bounds, compute_jacobian: bool) -> (
-            Tuple[Array, Array, Array, List[Error]]):
-        """Compute transformed marginal costs for this market. Then, if compute_jacobian is True, compute the Jacobian
-        (holding gamma fixed) of omega (equivalently, of transformed marginal costs) with respect to theta. Replace null
-        elements in transformed marginal costs with their last values before computing their Jacobian.
-        """
-        errors: List[Error] = []
-
-        # compute transformed marginal costs
-        tilde_costs, clipped_costs, tilde_costs_errors = self.safely_compute_tilde_costs(costs_bounds)
-        errors.extend(tilde_costs_errors)
-
-        # replace invalid transformed marginal costs with their last computed values
-        valid_tilde_costs = tilde_costs.copy()
-        bad_tilde_costs_index = ~np.isfinite(tilde_costs)
-        valid_tilde_costs[bad_tilde_costs_index] = last_tilde_costs[bad_tilde_costs_index]
+        # only compute supply side contributions if there is a supply side
+        if self.K3 == 0:
+            tilde_costs = np.full((self.J, 0), np.nan, options.dtype)
+            omega_jacobian = np.full((self.J, self.parameters.P), np.nan, options.dtype)
+            clipped_costs = np.zeros((self.J, 1), np.bool_)
+        else:
+            assert probabilities is not None
 
-        # compute the Jacobian, which is zero for clipped marginal costs
-        omega_jacobian = np.full((self.J, self.parameters.P), np.nan, options.dtype)
-        if compute_jacobian:
-            omega_jacobian, omega_jacobian_errors = self.safely_compute_omega_by_theta_jacobian(
-                valid_tilde_costs, xi_jacobian
+            # compute transformed marginal costs
+            tilde_costs, clipped_costs, tilde_costs_errors = self.safely_compute_tilde_costs(
+                probabilities, conditionals, costs_bounds
             )
-            errors.extend(omega_jacobian_errors)
-            omega_jacobian[clipped_costs.flat] = 0
+            errors.extend(tilde_costs_errors)
+
+            # replace invalid transformed marginal costs with their last computed values
+            valid_tilde_costs = tilde_costs.copy()
+            bad_tilde_costs_index = ~np.isfinite(tilde_costs)
+            valid_tilde_costs[bad_tilde_costs_index] = last_tilde_costs[bad_tilde_costs_index]
+
+            # compute the Jacobian, which is zero for clipped marginal costs
+            omega_jacobian = np.full((self.J, self.parameters.P), np.nan, options.dtype)
+            if compute_jacobians:
+                assert probabilities_tensor is not None
+                omega_jacobian, omega_jacobian_errors = self.safely_compute_omega_by_theta_jacobian(
+                    valid_tilde_costs, probabilities, conditionals, probabilities_tangent_mapping,
+                    conditionals_tangent_mapping, probabilities_tensor, conditionals_tensor, xi_jacobian
+                )
+                errors.extend(omega_jacobian_errors)
+                omega_jacobian[clipped_costs.flat] = 0
 
-        return tilde_costs, omega_jacobian, clipped_costs, errors
+        return (
+            delta, xi_jacobian, micro_numerator, micro_denominator, micro_numerator_jacobian,
+            micro_denominator_jacobian, micro_covariances_numerator, clipped_shares, stats, tilde_costs, omega_jacobian,
+            clipped_costs, errors
+        )
 
     @NumericalErrorHandler(exceptions.DeltaNumericalError)
     def safely_compute_delta(
@@ -113,26 +146,27 @@ def safely_compute_delta(
 
     @NumericalErrorHandler(exceptions.XiByThetaJacobianNumericalError)
     def safely_compute_xi_by_theta_jacobian(
-            self, delta: Array, keep_probabilities_tangents: bool) -> (
-            Tuple[Array, Array, Optional[Array], Dict[int, Array], List[Error]]):
+            self, probabilities: Array, conditionals: Optional[Array],
+            probabilities_tangent_mapping: Dict[int, Array]) -> Tuple[Array, List[Error]]:
         """Compute the Jacobian (holding beta fixed) of xi (equivalently, of delta) with respect to theta, handling any
         numerical errors.
         """
-        return self.compute_xi_by_theta_jacobian(delta, keep_probabilities_tangents)
+        return self.compute_xi_by_theta_jacobian(probabilities, conditionals, probabilities_tangent_mapping)
 
     @NumericalErrorHandler(exceptions.MicroMomentsNumericalError)
     def safely_compute_micro_contributions(
-            self, moments: Moments, delta: Array, probabilities: Optional[Array], conditionals: Optional[Array],
-            probabilities_tangent_mapping: Dict[int, Array], xi_jacobian: Array, compute_jacobians: bool,
-            compute_covariances: bool) -> Tuple[Array, Array, Array, Array, Array, List[Error]]:
+            self, moments: Moments, delta: Array, probabilities: Optional[Array],
+            probabilities_tangent_mapping: Dict[int, Array], probabilities_tensor: Optional[Array], xi_jacobian: Array,
+            compute_jacobians: bool, compute_covariances: bool) -> (
+            Tuple[Array, Array, Array, Array, Array, List[Error]]):
         """Compute micro moment value contributions, handling any numerical errors."""
         errors: List[Error] = []
         (
             micro_numerator, micro_denominator, micro_numerator_jacobian, micro_denominator_jacobian,
             micro_covariances_numerator
         ) = (
             self.compute_micro_contributions(
-                moments, delta, probabilities, conditionals, probabilities_tangent_mapping, xi_jacobian,
+                moments, delta, probabilities, probabilities_tangent_mapping, probabilities_tensor, xi_jacobian,
                 compute_jacobians, compute_covariances
             )
         )
@@ -142,12 +176,14 @@ def safely_compute_micro_contributions(
         )
 
     @NumericalErrorHandler(exceptions.CostsNumericalError)
-    def safely_compute_tilde_costs(self, costs_bounds: Bounds) -> Tuple[Array, Array, List[Error]]:
+    def safely_compute_tilde_costs(
+            self, probabilities: Array, conditionals: Optional[Array], costs_bounds: Bounds) -> (
+            Tuple[Array, Array, List[Error]]):
         """Compute transformed marginal costs, handling any numerical errors."""
         errors: List[Error] = []
 
         # compute marginal costs
-        eta, eta_errors = self.compute_eta()
+        eta, eta_errors = self.compute_eta(probabilities=probabilities, conditionals=conditionals)
         errors.extend(eta_errors)
         costs = self.products.prices - eta
 
@@ -170,8 +206,14 @@ def safely_compute_tilde_costs(self, costs_bounds: Bounds) -> Tuple[Array, Array
 
     @NumericalErrorHandler(exceptions.OmegaByThetaJacobianNumericalError)
     def safely_compute_omega_by_theta_jacobian(
-            self, tilde_costs: Array, xi_jacobian: Array) -> Tuple[Array, List[Error]]:
+            self, tilde_costs: Array, probabilities: Array, conditionals: Optional[Array],
+            probabilities_tangent_mapping: Dict[int, Array], conditionals_tangent_mapping: Dict[int, Optional[Array]],
+            probabilities_tensor: Array, conditionals_tensor: Optional[Array], xi_jacobian: Array) -> (
+            Tuple[Array, List[Error]]):
         """Compute the Jacobian (holding gamma fixed) of omega (equivalently, of transformed marginal costs) with
         respect to theta, handling any numerical errors.
         """
-        return self.compute_omega_by_theta_jacobian(tilde_costs, xi_jacobian)
+        return self.compute_omega_by_theta_jacobian(
+            tilde_costs, probabilities, conditionals, probabilities_tangent_mapping, conditionals_tangent_mapping,
+            probabilities_tensor, conditionals_tensor, xi_jacobian
+        )

pyblp/markets/results_market.py

@@ -1,6 +1,6 @@
 """Market-level structuring of BLP results."""
 
-from typing import Any, List, Optional, Tuple
+from typing import Any, Dict, List, Optional, Tuple
 
 import numpy as np
 
@@ -39,21 +39,45 @@ def safely_solve_equilibrium_realization(
         shares = self.compute_probabilities(delta, mu)[0] @ self.agents.weights
         return prices, shares, delta, stats, errors
 
-    @NumericalErrorHandler(exceptions.XiByThetaJacobianRealizationNumericalError)
-    def safely_compute_xi_by_theta_jacobian_realization(self) -> Tuple[Array, List[Error]]:
+    @NumericalErrorHandler(exceptions.JacobianRealizationNumericalError)
+    def safely_compute_jacobian_realizations(self, tilde_costs: Array) -> Tuple[Array, Array, List[Error]]:
         """Compute the Jacobian (holding beta fixed) of xi (equivalently, of delta) with respect to theta for a
-        realization of the market, handling any numerical errors.
+        realization of the market, handling any numerical errors. If there is a supply side, do the same for omega.
         """
-        jacobian, _, _, _, errors = self.compute_xi_by_theta_jacobian()
-        return jacobian, errors
+        errors: List[Error] = []
 
-    @NumericalErrorHandler(exceptions.OmegaByThetaJacobianRealizationNumericalError)
-    def safely_compute_omega_by_theta_jacobian_realization(
-            self, tilde_costs: Array, xi_jacobian: Array) -> Tuple[Array, List[Error]]:
-        """Compute the Jacobian (holding gamma fixed) of omega (equivalently, of transformed marginal costs) with
-        respect to theta for a realization of the market, handling any numerical errors.
-        """
-        return self.compute_omega_by_theta_jacobian(tilde_costs, xi_jacobian)
+        # if needed, pre-compute probabilities and their derivatives with respect to parameters
+        probabilities, conditionals = self.compute_probabilities()
+        probabilities_tangent_mapping: Dict[int, Array] = {}
+        conditionals_tangent_mapping: Dict[int, Array] = {}
+        for p, parameter in enumerate(self.parameters.unfixed):
+            probabilities_tangent, conditionals_tangent = self.compute_probabilities_by_parameter_tangent(
+                parameter, probabilities, conditionals
+            )
+            probabilities_tangent_mapping[p] = probabilities_tangent
+            if self.K3 > 0:
+                conditionals_tangent_mapping[p] = conditionals_tangent
+
+        # compute the demand-side Jacobian
+        xi_jacobian, xi_jacobian_errors = self.compute_xi_by_theta_jacobian(
+            probabilities, conditionals, probabilities_tangent_mapping
+        )
+        errors.extend(xi_jacobian_errors)
+
+        # compute the supply-side Jacobian if there is a supply side
+        if self.K3 == 0:
+            omega_jacobian = np.full((self.J, self.parameters.P), np.nan, options.dtype)
+        else:
+            probabilities_tensor, conditionals_tensor = self.compute_probabilities_by_xi_tensor(
+                probabilities, conditionals
+            )
+            omega_jacobian, omega_jacobian_errors = self.compute_omega_by_theta_jacobian(
+                tilde_costs, probabilities, conditionals, probabilities_tangent_mapping, conditionals_tangent_mapping,
+                probabilities_tensor, conditionals_tensor, xi_jacobian
+            )
+            errors.extend(omega_jacobian_errors)
+
+        return xi_jacobian, omega_jacobian, errors
 
     @NumericalErrorHandler(exceptions.DeltaNumericalError)
     def safely_compute_delta(

pyblp/markets/simulation_market.py

@@ -107,7 +107,8 @@ def safely_compute_tilde_costs(self, delta: Array) -> Tuple[Array, List[Error]]:
         errors: List[Error] = []
 
         # compute marginal costs
-        eta, eta_errors = self.compute_eta(delta=delta)
+        probabilities, conditionals = self.compute_probabilities(delta)
+        eta, eta_errors = self.compute_eta(probabilities=probabilities, conditionals=conditionals)
         errors.extend(eta_errors)
         costs = self.products.prices - eta