Hubbry Logo
Berndt–Hall–Hall–Hausman algorithmBerndt–Hall–Hall–Hausman algorithmMain
Open search
Berndt–Hall–Hall–Hausman algorithm
Community hub
Berndt–Hall–Hall–Hausman algorithm
logo
7 pages, 0 posts
0 subscribers
Be the first to start a discussion here.
Be the first to start a discussion here.
Berndt–Hall–Hall–Hausman algorithm
Berndt–Hall–Hall–Hausman algorithm
Berndt–Hall–Hall–Hausman algorithm
View on Wikipedia
from Wikipedia

The Berndt–Hall–Hall–Hausman (BHHH) algorithm is a numerical optimization algorithm similar to the Newton–Raphson algorithm, but it replaces the observed negative Hessian matrix with the outer product of the gradient. This approximation is based on the information matrix equality and therefore only valid while maximizing a likelihood function.[1] The BHHH algorithm is named after the four originators: Ernst R. Berndt, Bronwyn Hall, Robert Hall, and Jerry Hausman.[2]

Usage

[edit]

If a nonlinear model is fitted to the data one often needs to estimate coefficients through optimization. A number of optimization algorithms have the following general structure. Suppose that the function to be optimized is Q(β). Then the algorithms are iterative, defining a sequence of approximations, βk given by

,

where is the parameter estimate at step k, and is a parameter (called step size) which partly determines the particular algorithm. For the BHHH algorithm λk is determined by calculations within a given iterative step, involving a line-search until a point βk+1 is found satisfying certain criteria. In addition, for the BHHH algorithm, Q has the form

and A is calculated using

In other cases, e.g. Newton–Raphson, can have other forms. The BHHH algorithm has the advantage that, if certain conditions apply, convergence of the iterative procedure is guaranteed.[citation needed]

See also

[edit]

References

[edit]

Further reading

[edit]
Revisions and contributorsEdit on WikipediaRead on Wikipedia

Berndt–Hall–Hall–Hausman algorithm

View on Grokipedia
from Grokipedia
The Berndt–Hall–Hall–Hausman (BHHH) algorithm is a numerical optimization technique designed for of parameters in nonlinear structural econometric models. Introduced in 1974 by economists Ernst R. Berndt, Bronwyn H. Hall, Robert E. Hall, and Jerry A. Hausman, it serves as an efficient hill-climbing method that approximates the of the log-likelihood function using the of the (score) vectors, thereby avoiding the need for explicit second derivatives and reducing computational demands. This approach leverages the equality, where the of the equals the negative expected Hessian, ensuring reliable convergence to the maximum under standard regularity conditions. The algorithm operates iteratively by computing the sample covariance matrix RR of the individual observation scores (gradients), which approximates the information matrix scaled by the sample size. Parameter updates follow the form βi+1=βi+λR−1g(βi)\beta^{i+1} = \beta^i + \lambda R^{-1} g(\beta^i), where g(βi)g(\beta^i) is the vector and λ\lambda is a step size selected to satisfy a convergence criterion, such as ensuring the likelihood increases monotonically. Upon convergence, the inverse of RR provides consistent estimates of the asymptotic variance- for inference purposes, including likelihood ratio tests distributed as chi-squared. Compared to full Newton-Raphson methods, BHHH is particularly advantageous for high-dimensional or nonlinear problems, such as simultaneous systems or production functions, as it requires only first derivatives and guarantees ascent toward the global maximum for concave log-likelihoods. In econometric practice, the BHHH algorithm has become a standard tool for estimating models involving nonlinearity in parameters and variables, including systems, frontiers, and dynamic specifications. Its implementation in software like RATS and has facilitated widespread adoption, with the method's efficiency making it suitable even for non-log-likelihood objectives, though standard errors may then require separate computation. The algorithm's development addressed key challenges in early computational , improving upon prior techniques like those of Eisenpress and Greenstadt by providing both theoretical guarantees and practical simplicity.

Background

Maximum Likelihood Estimation

Maximum likelihood estimation (MLE) is a fundamental statistical method used to estimate the parameters θ\theta of a probabilistic model by maximizing the L(θ)L(\theta), which represents the probability of observing the given data under the model. Equivalently, it maximizes the log-likelihood ln⁡L(θ)=∑i=1nln⁡f(yi∣θ)\ln L(\theta) = \sum_{i=1}^n \ln f(y_i \mid \theta), where f(yi∣θ)f(y_i \mid \theta) is the probability density (or mass) function for each observation yiy_i, assuming a sample of nn independent observations. This approach, pioneered by , selects the parameter values that make the observed data most probable, providing estimates that are consistent, asymptotically normal, and efficient under standard regularity conditions. Key assumptions underlying MLE include the independence of observations, such that the joint likelihood factors into the product of individual densities, and correct specification of the underlying model, meaning the chosen ff accurately reflects the data-generating process. Violations, such as dependence or misspecification, can lead to inconsistent estimates, though MLE remains robust in many cases when the model is reasonably well-specified. In , MLE plays a central role in estimating parameters of nonlinear structural models, where θ\theta captures underlying economic relationships, such as systems or simultaneous equations, that do not admit closed-form solutions under ordinary . These models often involve latent variables or non-normal errors, making MLE preferable for its use of full distributional information to achieve efficiency and enable inference via likelihood-based tests. For complex log-likelihood functions in such settings, numerical optimization techniques are typically required to locate the maximum. A representative example is the simple binary choice model, where an outcome yi∈{0,1}y_i \in \{0, 1\} follows a with success probability π\pi, and the log-likelihood is ln⁡L(π)=∑i=1n[(1−yi)ln⁡(1−π)+yiln⁡π]\ln L(\pi) = \sum_{i=1}^n \left[ (1 - y_i) \ln(1 - \pi) + y_i \ln \pi \right]. The MLE π^=1n∑i=1nyi\hat{\pi} = \frac{1}{n} \sum_{i=1}^n y_i matches the sample mean, illustrating how MLE recovers intuitive estimators even in discrete settings. In more general binary models, such as or , π\pi depends nonlinearly on covariates, underscoring MLE's applicability to structural econometric applications.

Optimization Challenges in Nonlinear Models

In nonlinear econometric models, the log-likelihood function, typically expressed as ln⁡Q(θ)=∑iQi(θ)\ln Q(\theta) = \sum_i Q_i(\theta) where Qi(θ)Q_i(\theta) represents the contribution from the ii-th , deviates from the of linear models, resulting in non-convexity and the potential for multiple local maxima. This non-quadratic structure arises because individual likelihood contributions Qi(θ)Q_i(\theta) are often highly nonlinear in the parameters θ\theta, leading to objective functions that may exhibit saddle points or flat regions, which complicate the identification of the global maximum. Such increases the risk of convergence to suboptimal solutions, particularly in models involving discrete choices, limited dependent variables, or structural relationships in . Computing the exact Hessian matrix for these models imposes significant demands, as it requires second derivatives that scale poorly with model size and dimensionality, often rendering analytical evaluation infeasible for large datasets or complex specifications. optimization methods, such as steepest ascent, can mitigate some Hessian-related costs but frequently exhibit slow convergence due to the irregular of the likelihood surface in nonlinear settings. Prior to 1974, practitioners relied on grid searches or rudimentary approximate techniques, which lacked guarantees of convergence and were computationally intensive even for modest models, often necessitating manual intervention or restrictive assumptions to ensure stability. In simultaneous equation models, endogeneity—where explanatory variables are correlated with error terms—necessitates to achieve efficiency, yet this amplifies optimization complexity by entangling parameter identification across equations and introducing further nonlinearity through dependencies. The resulting likelihood functions are particularly prone to ill-conditioning and sensitivity to initial conditions, exacerbating the challenges of and computational expense in empirical applications like demand systems or macroeconometric simulations.

History

Development Context

In the 1970s, the field of underwent a notable shift toward the adoption of nonlinear models to better represent intricate economic phenomena, including production functions and demand systems. These models accounted for the inherent nonlinearities in parameters and variables that linear approximations could not capture, driving the need for robust (MLE) techniques to derive parameter estimates. This evolution reflected a growing emphasis on theoretical accuracy in modeling economic behavior, where traditional linear frameworks proved inadequate for simultaneous relationships in economic systems. Existing optimization methods, such as the Newton-Raphson algorithm, posed significant challenges for these nonlinear estimations, as they demanded computationally expensive evaluations of second and even third derivatives, often leading to unreliable convergence in practical applications. The computational burden was particularly prohibitive given the limited processing power available at the time, prompting the search for alternatives that utilized only first derivatives—such as gradients—while providing guarantees of convergence to a local maximum. This gap highlighted the limitations of prior approaches in handling the complexity of nonlinear MLE without excessive resource demands. The motivation for improved methods was especially strong in structural econometric models characterized by simultaneity, where variables influenced each other across equations, as in systems of . In such contexts, full information maximum likelihood (FIML) estimation was favored over limited information techniques, like two-stage , because it exploited all available a priori restrictions system-wide to achieve consistent and asymptotically efficient parameter estimates. This preference underscored the theoretical superiority of FIML for capturing interdependencies, though its implementation required overcoming substantial numerical hurdles. These developments were bolstered by contemporaneous advances in computer technology, which expanded the feasibility of estimating nonlinear models on a larger scale and paved the way for innovations in . The era's improving hardware and software capabilities enabled econometricians to move beyond simple linear regressions toward more sophisticated analyses, fostering a wave of methodological progress that addressed the demands of increasingly data-rich and theoretically complex economic inquiries.

Publication and Authors

The Berndt–Hall–Hall–Hausman algorithm was developed by four economists: Ernst R. Berndt, an econometrician and Louis E. Seley Emeritus of Applied at the ; H. Hall, a expert and Emerita of at the ; E. Hall, a macroeconomist and Robert and Carole McNeil Senior Fellow at Stanford University's ; and Jerry A. Hausman, a labor economist and John and Jennie S. MacDonald of at MIT. At the time of publication, Berndt was at the , B. H. Hall at , and both R. E. Hall and Hausman at MIT. Note that the two Halls— and —are distinct individuals with no direct relation. The algorithm was introduced in the seminal paper "Estimation and Inference in Nonlinear Structural Models" by Berndt, B. H. Hall, R. E. Hall, and Hausman, published in the Annals of Economic and Social Measurement, Volume 3, Number 4 (October 1974), pages 653–665. In the paper, the authors describe the algorithm as requiring less computation than prior methods, such as , by eliminating the need for third derivatives and full evaluations in each iteration. They further claim that it guarantees convergence to a local maximum under standard regularity conditions for the , offering improvements over earlier approaches like those of Eisenpress and Greenstadt (1966) and Chow (1973). The work has had substantial impact, influencing the implementation of maximum likelihood routines in econometric software packages, such as Stata's technique(bhhh) option.

Mathematical Formulation

Objective and Gradient

The Berndt–Hall–Hall–Hausman (BHHH) algorithm addresses the maximization of the log-likelihood function for parameter estimation in nonlinear statistical models. The objective function is defined as the log-likelihood ln⁡Q(θ)=∑i=1NQi(θ)\ln Q(\theta) = \sum_{i=1}^N Q_i(\theta), where Qi(θ)=ln⁡f(yi,xi∣θ)Q_i(\theta) = \ln f(y_i, x_i \mid \theta) represents the log-density for the ii-th observation, with ff denoting the joint conditional on the θ\theta, covariates xix_i, and outcomes yiy_i. The first-order condition for maximization is given by the gradient vector, or score function, g(θ)=∂ln⁡Q(θ)∂θ=∑i=1N∂Qi(θ)∂θg(\theta) = \frac{\partial \ln Q(\theta)}{\partial \theta} = \sum_{i=1}^N \frac{\partial Q_i(\theta)}{\partial \theta}. At the maximum likelihood estimator θ^\hat{\theta}, the score satisfies g(θ^)=0g(\hat{\theta}) = 0, and under standard regularity conditions, E[g(θ)]=0E[g(\theta)] = 0 at the true parameter value, ensuring that the parameters align with the data's probabilistic structure under the model assumptions. A key property underpinning the algorithm is the information matrix equality, which states that under standard regularity conditions, the expected negative of the log-likelihood equals the expected of the score: −E[∂2ln⁡Q∂θ∂θ′]=E[g(θ)g(θ)′]=I(θ)-E\left[\frac{\partial^2 \ln Q}{\partial \theta \partial \theta'}\right] = E[g(\theta) g(\theta)'] = I(\theta), where I(θ)I(\theta) is the matrix. This equality connects the variance of the score to the curvature of the objective, facilitating efficient numerical optimization in high-dimensional parameter spaces. In practice, the g(θ)g(\theta) indicates the direction of ascent toward the likelihood maximum, with the algorithm exploiting the additive structure of individual contributions ∂Qi(θ)∂θ\frac{\partial Q_i(\theta)}{\partial \theta} to compute updates, particularly advantageous for models where evaluating the full Hessian is computationally intensive.

Hessian Approximation

In , the H(θ)=∂2ln⁡Q(θ)∂θ∂θ′=∑i=1N∂2Qi(θ)∂θ∂θ′H(\theta) = \frac{\partial^2 \ln Q(\theta)}{\partial \theta \partial \theta'} = \sum_{i=1}^N \frac{\partial^2 Q_i(\theta)}{\partial \theta \partial \theta'} provides second-order information essential for Newton-type optimization algorithms, but its direct computation requires evaluating costly second derivatives for each observation ii. The core innovation of the Berndt–Hall–Hall–Hausman (BHHH) algorithm lies in approximating the negative inverse Hessian using the outer product of gradients (OPG), derived from the information matrix equality. Specifically, the approximation A(θ)≈−H(θ)−1=[∑i=1N(∂Qi(θ)∂θ)(∂Qi(θ)∂θ)′]−1A(\theta) \approx -H(\theta)^{-1} = \left[ \sum_{i=1}^N \left( \frac{\partial Q_i(\theta)}{\partial \theta} \right) \left( \frac{\partial Q_i(\theta)}{\partial \theta} \right)' \right]^{-1} leverages the fact that, under regularity conditions, the expected value of the outer product of individual score vectors equals the negative expected Hessian. This substitution avoids explicit second derivatives while preserving asymptotic equivalence to the full Newton method near the true parameter value. At iteration kk, the BHHH approximation takes the form Ak=[∑i(∂Qi∂β(βk))(∂Qi∂β(βk))′]−1,A_k = \left[ \sum_i \left( \frac{\partial Q_i}{\partial \beta} (\beta_k) \right) \left( \frac{\partial Q_i}{\partial \beta} (\beta_k) \right)' \right]^{-1}, where β\beta denotes the vector. This matrix AkA_k is positive definite provided the individual gradients are non-zero, ensuring descent directions in the optimization landscape, and its validity stems from asymptotic equality under the information matrix equality, which holds as the sample size increases and near the maximum likelihood .

Algorithm Procedure

Initialization

The initialization phase of the Berndt–Hall–Hall–Hausman (BHHH) algorithm is crucial for ensuring reliable convergence in of nonlinear models, as the choice of initial parameter vector θ0\theta_0 directly influences the algorithm's trajectory toward the optimum of the log-likelihood function. Poor initial values can cause the iterations to diverge or converge to a local maximum instead of the global one, potentially yielding biased or inefficient estimates. Selecting appropriate starting points mitigates these risks and enhances computational efficiency, particularly in complex econometric models where the likelihood surface may be multimodal. Common strategies for obtaining θ0\theta_0 leverage simpler estimation techniques to approximate the nonlinear model. For instance, ordinary least squares (OLS) estimates from a linearized version of the model provide a sensible starting point, as they offer consistent but potentially inefficient approximations. The method of moments, which matches sample moments to theoretical ones, or grid search over a parameter space can also generate feasible initial values when OLS is unsuitable. The initial θ0\theta_0 must satisfy basic feasibility criteria to enable the algorithm's first iteration: the log-likelihood ln⁡Q(θ0)\ln Q(\theta_0) should be finite, ensuring the objective is well-defined, and the score vector (gradient) must be computable without numerical issues such as or undefined functions.

Iterative Updates and Convergence

The iterative process of the Berndt–Hall–Hall–Hausman (BHHH) algorithm begins with an initial parameter vector θ0\theta_0 and proceeds through successive updates to maximize the log-likelihood function ln⁡Q(θ)\ln Q(\theta). At each iteration kk, the score vector () g(θk)=∑i=1nsi(θk)g(\theta_k) = \sum_{i=1}^n s_i(\theta_k) and the outer product of gradients (OPG) approximation Ak=∑i=1nsi(θk)si(θk)TA_k = \sum_{i=1}^n s_i(\theta_k) s_i(\theta_k)^T are computed, where si(θk)s_i(\theta_k) denotes the contribution to the score from the ii-th . The parameter update follows the rule θk+1=θk+λkAk−1g(θk)\theta_{k+1} = \theta_k + \lambda_k A_k^{-1} g(\theta_k), where λk\lambda_k is the step size, typically initialized to 1 to approximate the scoring method for ascent. This step leverages the OPG as an estimate of the information matrix, ensuring the direction points toward increasing ln⁡Q\ln Q. If the update results in a decrease in ln⁡Q(θk+1)\ln Q(\theta_{k+1}) compared to ln⁡Q(θk)\ln Q(\theta_k), the step size is adjusted downward using a procedure until ascent is ensured or a minimum step size is reached. The iteration continues until convergence criteria are satisfied, including a small gradient norm ∣∣g(θk)∣∣<ϵ||g(\theta_k)|| < \epsilon (e.g., ϵ=10−5\epsilon = 10^{-5}), a small parameter change ∣∣θk+1−θk∣∣<δ||\theta_{k+1} - \theta_k|| < \delta (e.g., δ=10−6\delta = 10^{-6}), or a maximum number of iterations (e.g., 100–300). Additionally, the algorithm verifies the positive definiteness of AkA_k at each step to confirm a valid ascent direction; if AkA_k is not positive definite, the iteration may restart or terminate. These criteria balance computational efficiency with reliable convergence to a local maximum.

Properties

Advantages

The Berndt–Hall–Hall–Hausman (BHHH) algorithm provides notable computational efficiency for by relying solely on first-order derivatives to approximate the Hessian via the of gradients (OPG), avoiding the more demanding second- and third-order derivative evaluations required in Newton-type methods. This approach computes the OPG matrix in O(N p²) time, where N denotes the sample size and p the number of parameters, rendering it particularly practical for problems with moderate dimensions. Under standard regularity conditions—such as the log-likelihood being twice continuously differentiable, the existence of a unique maximum, and the being —the BHHH guarantees global convergence to a from any feasible starting value where the log-likelihood is defined, owing to its structure as a gradient method with appropriate for the step size λ_k. The of the OPG approximation further ensures that the search direction remains an ascent direction when λ_k is selected via or similar techniques, promoting reliable progress toward the optimum. The algorithm's suitability shines in large-sample settings, where the OPG consistently approximates the expected Hessian due to the asymptotic equality of the observed and expected information matrices, yielding robust and efficient parameter estimates without excessive computational overhead.

Limitations

The Berndt–Hall–Hall–Hausman (BHHH) algorithm's approximation of the Hessian matrix using the outer product of gradients (OPG) can be unreliable when the information matrix equality does not hold closely, such as in small samples or under model misspecification, leading to slower convergence and requiring more iterations to approach the maximum. This approximation performs poorly far from the likelihood maximum, where the OPG matrix deviates significantly from the true negative Hessian, often resulting in small parameter updates that increase the log-likelihood only marginally and prolong the optimization process. The algorithm is sensitive to the choice of step size in its updates, typically requiring a or to ensure monotonic improvement in the objective; without such adjustments, iterations may oscillate or diverge, particularly in nonquadratic likelihood surfaces. The BHHH procedure relies on first-order derivatives (scores) and assumes the structure of a log-likelihood function, making it strictly suited for and less effective for objectives involving penalties, robust losses, or non-likelihood criteria that violate the information equality. In high-dimensional problems with many parameters pp, the need to form and invert the p×pp \times p OPG matrix at each iteration incurs a computational cost of O(p3)O(p^3) for the inversion alone, which becomes prohibitive when pp is large, despite the overall relative to exact Hessian computation. The algorithm's dependence on accurate gradient evaluations further exacerbates issues in scenarios where is unstable or costly.

Comparisons

Newton–Raphson Method

The Newton–Raphson method is a second-order optimization employed to locate the maximum of a twice-differentiable objective function, such as the log-likelihood in (MLE). It generates a sequence of parameter estimates θk\theta_k through iterative updates given by θk+1=θk−H(θk)−1g(θk),\theta_{k+1} = \theta_k - H(\theta_k)^{-1} g(\theta_k), where g(θk)g(\theta_k) denotes the gradient vector (first derivatives) and H(θk)H(\theta_k) represents the (second derivatives) of the objective function evaluated at the current iterate θk\theta_k. A key strength of the method lies in its rapid convergence properties: under suitable conditions, such as the Hessian being positive definite and continuous near the optimum, it achieves quadratic convergence, whereby the error decreases quadratically, often requiring fewer s than alternatives for smooth, well-behaved problems. Despite these advantages, the method's drawbacks stem from the intensive computation of the exact Hessian, which involves evaluating p2p^2 second partial derivatives across NN observations at each , yielding an overall of O(Np2)O(N p^2) operations plus the O(p3)O(p^3) expense of matrix inversion for pp parameters. Additionally, the Hessian may fail to be positive definite, particularly distant from the solution, which can cause or and requires remedial steps like line searches or regularization. By comparison, the Berndt–Hall–Hall–Hausman (BHHH) algorithm addresses these issues by approximating the via the of gradients, maintaining the O(Np2)O(N p^2) scaling but eliminating the need for second derivatives, thereby reducing analytical burden and ensuring positive semi-definiteness while accepting slightly slower convergence in exchange for practicality in large-scale MLE settings.

Gauss–Newton Method

The Gauss–Newton method is an iterative algorithm designed for solving problems, where the objective is to minimize the sum of squared residuals, ∑i=1nri(θ)2\sum_{i=1}^n r_i(\theta)^2, with r(θ)\mathbf{r}(\theta) denoting the vector of residuals and θ\theta the parameter vector. In this approach, the of the objective function is approximated by J(θk)TJ(θk)\mathbf{J}(\theta_k)^T \mathbf{J}(\theta_k), where J(θk)\mathbf{J}(\theta_k) is the matrix of the residuals evaluated at the current iterate θk\theta_k. This approximation neglects the second-order derivatives of the residuals, simplifying computation while leveraging the structure of objectives. The parameter update is then given by θk+1=θk−(J(θk)TJ(θk))−1J(θk)Tr(θk),\theta_{k+1} = \theta_k - \left( \mathbf{J}(\theta_k)^T \mathbf{J}(\theta_k) \right)^{-1} \mathbf{J}(\theta_k)^T \mathbf{r}(\theta_k), which corresponds to a single step of linear least squares on the linearized residual model r(θk+δ)≈r(θk)+J(θk)δ\mathbf{r}(\theta_k + \delta) \approx \mathbf{r}(\theta_k) + \mathbf{J}(\theta_k) \delta. This formulation derives from applying Newton's method to the least squares objective and approximating the full Hessian, making it particularly efficient for problems where residuals represent model errors, such as in regression or curve fitting. In the context of (MLE), the Gauss–Newton method applies directly when the log-likelihood function takes a form, such as under Gaussian errors where ℓ(θ)=−12∑ri(θ)2+constant\ell(\theta) = -\frac{1}{2} \sum r_i(\theta)^2 + \text{constant}, transforming MLE into . For more general likelihoods, including those from distributions, the method can be extended by reparameterizing the score contributions, leading to an outer-product-of-gradients (OPG) approximation of the information matrix that mirrors the structure of the Berndt–Hall–Hall–Hausman (BHHH) algorithm. Specifically, under the information matrix equality, the expected outer product of the score vector equals the negative expected Hessian, allowing Gauss–Newton iterations to align with OPG-based updates for MLE in these cases. A key strength of the Gauss–Newton method is its avoidance of explicit second derivatives, relying solely on first-order information via the , which reduces computational cost compared to full Newton methods. Near the solution, where the is accurate and residuals are small, it exhibits quadratic convergence, often outperforming gradient-based methods in speed for well-behaved problems. However, far from the minimum, the approximation can fail if the neglected second-order terms are significant, leading to poor step directions or divergence, much like issues observed in the BHHH algorithm for MLE. This sensitivity to initial conditions and nonlinearity underscores the need for safeguards, such as or line searches, in practical implementations. The BHHH algorithm serves as the direct MLE analog to Gauss–Newton, replacing the of residuals with the of individual score vectors to approximate the matrix, rather than using residual . In interpretable as MLE, the two methods coincide, but BHHH generalizes to arbitrary likelihoods by exploiting the score's zero mean and variance properties. This distinction highlights Gauss–Newton's specialization to squared-error objectives, while BHHH extends the approximation logic to broader statistical estimation contexts.

Applications

Econometric Models

The Berndt–Hall–Hall–Hausman (BHHH) algorithm plays a key role in estimating parameters for simultaneous equation models through full information maximum likelihood (FIML), especially in systems addressing endogeneity, such as frameworks. By leveraging the outer product of gradients to approximate the information matrix, the algorithm efficiently optimizes the nonlinear , providing consistent and asymptotically efficient estimates while managing the computational challenges of correlated equations. This approach is particularly valuable in structural econometric models where ordinary fails due to simultaneity bias, as demonstrated in early applications to nonlinear systems. In models with limited dependent variables, the BHHH algorithm supports for binary outcomes in and specifications, where the introduces nonlinearity. It extends effectively to models, optimizing complex likelihoods that account for choice correlations and alternative-specific effects, such as in analysis of consumer behavior or labor market decisions. The algorithm's reliance on first derivatives alone makes it suitable for these settings, avoiding the need for second derivatives that can be cumbersome in high-dimensional choice sets. For time series models, the BHHH algorithm facilitates parameter estimation in ARMA processes with nonlinear components and GARCH models capturing conditional heteroskedasticity in financial returns. In GARCH specifications, it maximizes the by iteratively updating parameters based on score vectors, accommodating without requiring full Hessian computations, which enhances in volatile data environments. This application has been instrumental in asset returns and , where precise volatility estimates are critical. A notable example of the BHHH algorithm's application appears in the estimation of parameters within nonlinear structural models using , as illustrated in the seminal 1974 work. Here, optimizes the likelihood for systems incorporating firm-level heterogeneity and structural relationships, such as input-output dynamics, yielding reliable inferences on and despite the models' . This structural approach highlights the algorithm's utility in bridging theoretical with empirical .

Software Implementations

The Berndt–Hall–Hall–Hausman (BHHH) algorithm is implemented in several major statistical software packages, particularly those focused on econometrics and maximum likelihood estimation (MLE). In Stata, the algorithm is available through the ml command via the technique(bhhh) option, which allows users to specify BHHH as the optimization method for likelihood maximization. This implementation supports hybrid approaches, such as performing an initial set of iterations with BHHH before switching to the Newton–Raphson method, for example using technique(bhhh 10 nr 1000) to run 10 BHHH iterations followed by 1,000 Newton–Raphson steps. Stata also provides controls for step size adjustments and maximum iterations to enhance convergence reliability in complex models. In , the BHHH algorithm is supported in the maxLik package, which provides the maxBHHH function as a dedicated optimizer for MLE problems, approximating the Hessian via the of gradients. This function integrates with R's broader and allows customization of iteration limits and step sizes, often combined with other methods like BFGS for improved performance in econometric applications. While not part of the base optim function, maxBHHH serves as a reliable variant for likelihood-based estimation in packages like maxLik. For Python users, the BHHH algorithm lacks built-in support in core libraries like statsmodels or scipy.optimize, which primarily offer quasi-Newton methods such as BFGS. However, open-source implementations are available, including a pure-Python version on GitHub that handles unconstrained optimization and is derived from a MATLAB routine by Fedor Iskhakov. This repository provides a lightweight, customizable code base suitable for integration into custom MLE workflows, with options for iteration controls and gradient-based updates. Similar open-source MATLAB implementations exist for academic and research purposes, often shared via repositories for nonlinear optimization tasks. The BHHH algorithm is a common option in econometric software for MLE due to its robustness in handling misspecified starting values and noisy gradients, though Newton–Raphson remains the default in tools like and . These implementations facilitate its use in fields like econometric modeling, where reliable convergence is essential for parameter estimation.

References

  1. https://www.nber.org/system/files/chapters/c10206/c10206.pdf
  2. https://estima.com/webhelp/topics/bhhh.html
  3. https://pages.stern.nyu.edu/~wgreene/DiscreteChoice/Readings/Greene-Chapter-16.pdf
  4. https://discdown.org/microeconometrics/maximum-likelihood-estimation-1.html
  5. https://www.nber.org/system/files/chapters/c10205/c10205.pdf
  6. https://www.sciencedirect.com/science/article/pii/S1573441283010107
  7. https://mitsloan.mit.edu/faculty/directory/ernst-rudolf-berndt
  8. https://econ.berkeley.edu/profile/bronwyn-h-hall
  9. https://www.hoover.org/profiles/robert-e-hall
  10. https://economics.mit.edu/people/faculty/jerry-hausman
  11. https://www.nber.org/books-and-chapters/annals-economic-and-social-measurement-volume-3-number-4/estimation-and-inference-nonlinear-structural-models
  12. https://www.stata.com/manuals/rmaximize.pdf
  13. https://www.jstatsoft.org/article/view/v027i07/259
  14. https://eml.berkeley.edu/~bhhall/papers/BerndtHallHallHausman74.pdf
  15. https://www.sciencedirect.com/science/article/pii/S2666691X20300440
  16. https://ageconsearch.umn.edu/record/90849/files/software2.pdf
  17. https://cdn.aptech.com/www/uploads/2012/11/ml5.0.pdf
  18. https://eml.berkeley.edu/choice2/ch8.pdf
  19. https://s2.smu.edu/tfomby/eco6352/Notes/MLE_Lecture%206.pdf
  20. https://www.cirrelt.ca/documentstravail/cirrelt-2014-64.pdf
  21. https://bookdown.org/rdpeng/advstatcomp/newtons-method.html
  22. https://cs.nyu.edu/overton/NumericalComputing/newton.pdf
  23. https://www.geeksforgeeks.org/machine-learning/newtons-method-in-machine-learning/
  24. http://www.seas.ucla.edu/~vandenbe/236C/lectures/gn.pdf
  25. https://epubs.siam.org/doi/10.1137/050624935
  26. https://www2.isye.gatech.edu/~jeffwu/isye8813/mle.pdf
  27. https://arxiv.org/pdf/math/0610402
  28. https://dspace.mit.edu/bitstream/handle/1721.1/63493/efficientestimat00haus.pdf?sequence=1&origin=publication_detail
  29. https://www.sciencedirect.com/science/article/pii/0304407688900310
  30. https://www.jstor.org/stable/1391677
  31. https://faculty.washington.edu/ezivot/research/practicalgarchfinal.pdf
  32. https://www.jstor.org/stable/2285098
  33. https://www.stata.com/manuals/rml.pdf
  34. https://cran.r-project.org/web/packages/maxLik/maxLik.pdf
  35. https://www.rdocumentation.org/packages/maxLik/versions/0.5-6/topics/maxBHHH
  36. https://cran.r-project.org/web/packages/maxLik/vignettes/using-maxlik.pdf
  37. https://github.com/segsell/bhhh
  38. https://www.mathworks.com/matlabcentral/answers/6255-bhhh-method-for-nonlinear-regression
  39. http://kumlai.free.fr/RESEARCH/THESE/TEXTE/STATA/MLE_in_Stata.pdf
  40. https://www.jstatsoft.org/article/view/v041i11/495
Add your contribution
Related Hubs
User Avatar
No comments yet.