Seemingly Unrelated Regression

Definition

The Seemingly Unrelated Regression (SUR) model consists of $m$ regression equations

$y_j = X_j \beta_j + u_j, \quad j = 1, \ldots, m$

where $y_j$ ( $n \times 1$ ) are the dependent variables, $X_j$ ( $n \times p_j$ ) are design matrices (potentially different across equations), $\beta_j$ ( $p_j \times 1$ ) are coefficient vectors, and the error terms are cross-sectionally correlated: $E[u_i u_j'] = \omega_{ij} I$ with $\Omega = (\omega_{ij})$ an $m \times m$ positive definite matrix. The equations are "seemingly unrelated" because each has its own regressors, yet the cross-equation error covariance creates informational linkages.

Key Ideas

GLS efficiency gain: When errors are correlated across equations, generalized least squares (GLS) (Zellner's efficient estimator) dominates equation-by-equation ordinary least squares (OLS) by exploiting cross-equation covariance. The gain is largest when the $X_j$ matrices are very different across equations and $|\omega_{ij}|$ is large.
Identification: If all equations share the same regressors ( $X_1 = \cdots = X_m$ ), GLS and OLS coincide — the cross-equation correlation provides no additional information in that special case.
Bayesian approaches: Zellner (1971) derived conditional posteriors for $\beta | \Omega$ and $\Omega | \beta$ under diffuse priors — both are standard (Normal and Inverse-Wishart), enabling Gibbs sampling. However, in high dimensions ( $p_j$ large relative to $n$ ), the Gibbs sampler has severe convergence problems.
Triangular reparameterization: A recursive transformation (Zellner et al. 1988) rewrites the model in a triangular form where the new parameters $\{b_j, \sigma_j^2\}$ are independent across equations under a modified Jeffreys' prior. This enables Direct Monte Carlo (DMC) sampling without Markov Chain Monte Carlo (MCMC).

How It Works

GLS Estimator (Zellner 1962)

Stack the system as $y = X\beta + u$ with $E[uu'] = \Sigma_\varepsilon \otimes I_T$ . The infeasible GLS (Aitken) estimator is: $b^* = (X'(\Sigma_\varepsilon^{-1} \otimes I_T)X)^{-1} X'(\Sigma_\varepsilon^{-1} \otimes I_T)y$

Two-stage Aitken (feasible GLS): Stage 1 — run OLS equation by equation and form residuals $\hat{u}_\mu$ ; estimate $(T-l)\hat{s}_{\mu\mu'} = \hat{u}_\mu'\hat{u}_{\mu'}$ to get $\hat{S}=(\hat{s}_{\mu\mu'})$ . Stage 2 — $b = (X'\hat{S}^{-1}X)^{-1}X'\hat{S}^{-1}y$ . Under mild regularity, $b$ has the same asymptotic distribution as $b^*$ ; $V(b) = (X'\Sigma_\varepsilon^{-1}X)^{-1} + o(T^{-1})$ .

Efficiency conditions: SUR $\equiv$ OLS if and only if (i) $\sigma_{\mu\mu'}=0$ for all $\mu\neq\mu'$ (uncorrelated errors) or (ii) $X_1=\cdots=X_M$ (all equations share the same regressors). Maximum efficiency gain occurs when regressors are orthogonal across equations and $|\sigma_{\mu\mu'}|$ is large.

Efficiency gain formula (two equations, Zellner 1962 eq. 3.8): $|V(b_1^*)| = \frac{(1-\rho^2)^{l_1}}{\prod_{\mu=1}^{l_1}(1-\rho^2 r_\mu^2)} \cdot |\sigma^2(X_1'X_1)^{-1}|$ where $r_\mu$ are canonical correlations between the column spaces of $X_1$ and $X_2$ , and $\rho = \sigma_{12}/\sqrt{\sigma_{11}\sigma_{22}}$ . The ratio approaches $(1-\rho^2)$ as $r_\mu\to 0$ and equals 1 when $r_\mu=1$ for all $\mu$ .

Grunfeld (1958) application: GE and Westinghouse investment equations ( $T=20$ , 1935–1954; regressors: $C_{t-1}$ , $F_{t-1}$ , constant). Estimated cross-equation error correlation $\hat\rho=0.81$ ; SUR variance approximately 20% below equation-by-equation OLS. The aggregation bias $F$ -test $\tilde{F}=3.452 > F_{3,34}(0.95)=2.88$ rejects $\beta_1=\beta_2$ at 5%.

Aggregation Bias Testing (Zellner 1962)

When a macro equation is viewed as an aggregate of micro equations $y_\mu = X_\mu\beta_\mu + u_\mu$ , aggregation is unbiased only if all $\beta_\mu$ are equal.

Micro test (eq. 4.4): $H_0: \beta_1=\cdots=\beta_M$ . Compute $\tilde{F} = (b_0'X'Sy - b'X'Sy) / (b'X'Sy) \cdot (n-m)/q$ , where $q=(M-1)l$ is the number of restrictions. Asymptotically $\chi^2_q/q$ under $H_0$ .

Macro test (eq. 4.8–4.11): Regress aggregate $\bar{y}_t = \sum_\mu\omega_\mu y_{\mu t}$ on $\bar{x}_t = \sum_\mu\omega_\mu x_{\mu t}$ and a residual term $w(t)\bar{x}(t)$ capturing deviations. If the coefficient on the residual term is zero, the macro equation is unbiased.

Bayesian Analysis: Independent Heterogeneous Variances (Tiao-Zellner 1964)

A precursor to the full SUR Bayesian analysis. When two equations share the same slope vector $\beta$ but have independent unknown variances $\sigma_1^2$ and $\sigma_2^2$ (rather than an unrestricted $\Omega$ ), independent noninformative priors on $\log\sigma_1$ and $\log\sigma_2$ yield a posterior that is the product of two multivariate t-distributions — the "double-t" distribution:

$p(\beta|y_1,y_2) \propto \left\{1 + \frac{Q(\beta,\hat\beta_1,Z_1)}{\nu_1 s_1^2}\right\}^{-\frac{1}{2}(\nu_1+p)}\!\!\left\{1 + \frac{Q(\beta,\hat\beta_2,Z_2)}{\nu_2 s_2^2}\right\}^{-\frac{1}{2}(\nu_2+p)}$

The normalizing constant has no closed form for $p > 1$ ; Tiao-Zellner develop an asymptotic expansion in $\nu_1^{-1}$ and $\nu_2^{-1}$ to approximate it. The posterior mean converges to the precision-weighted average $\bar\beta = (M_1+M_2)^{-1}(M_1\hat\beta_1 + M_2\hat\beta_2)$ — the Bayesian counterpart of Theil's (1963) mixed estimator — as $\nu_1\to\infty$ with $\sigma_1^2$ known. See Tiao-Zellner (1964) and Bayesian Linear Regression.

Bayesian Analysis under Diffuse Prior (Zellner 1971)

Under Jeffreys' invariant prior $\pi_1(\beta, \Omega) \propto |\Omega|^{-(m+1)/2}$ , the conditional posteriors are: $\beta | \Omega, D \sim \mathcal{N}(\hat\beta_{\text{GLS}}, \hat\Omega_\beta), \qquad \Omega | \beta, D \sim \text{IW}(R, n)$ where $R_{ij} = (y_i - X_i\beta_i)'(y_j - X_j\beta_j)$ . The two blocks are coupled: $\hat\beta_{\text{GLS}}$ depends on $\Omega$ and $R$ depends on $\beta$ . Gibbs sampling iterates between them.

Predictive Density Approximations (Percy 1992)

Percy (1992) is the first systematic Bayesian treatment of prediction for the SUR model. The goal is the predictive density $f(y_{n+1}|X_{n+1},\mathbb{T})$ for a new observation given training data $\mathbb{T} = \{y_1,\ldots,y_n, X_1,\ldots,X_n\}$ .

Intractability result. Under Jeffreys' prior $f(\beta,\Phi)\propto|\Phi|^{-(p+1)/2}$ , integrating $\Phi$ from the joint predictive density yields an integral over $\beta$ of the form $\int_\beta |\sum_{i=1}^{n+1}(y_i-X_i\beta)(y_i-X_i\beta)^T|^{-(n+1)/2}\,d\beta$ , which Drèze (1977) showed is intractable unless $x_{i1}=\cdots=x_{ip}$ for all $i$ (the traditional multivariate regression case, which gives a multivariate Student- $t$ predictive density; Zellner-Chetty 1965).

Gibbs sampler. Three conditional distributions (eqs. 6–8) form a valid Gibbs cycle: $f(y_{n+1}|\beta,\Phi,\mathbb{T}) \equiv N_p[X_{n+1}\beta,\;\Phi^{-1}]$ $f(\Phi|\beta,y_{n+1},\mathbb{T}) \equiv W_p\!\left[n+1,\;\Bigl(\textstyle\sum_{i=1}^{n+1}(y_i-X_i\beta)(y_i-X_i\beta)^T\Bigr)^{-1}\right]$ $f(\beta|\Phi,y_{n+1},\mathbb{T}) \equiv N_q\!\left[\Bigl(\textstyle\sum X_i^T\Phi X_i\Bigr)^{-1}\sum X_i^T\Phi y_i,\;\Bigl(\sum X_i^T\Phi X_i\Bigr)^{-1}\right]$

After $t=100$ burn-in iterations, the predictive density is approximated by averaging $m=200$ normal densities at sampled $(\beta^{(j)},\Phi^{(j)})$ draws (eq. 9). Geman-Geman (1984) convergence applies in both the $L_1$ - and supremum norms.

First-order approximation. Integrate $\beta$ from the predictive density analytically, then substitute the modal Bayes estimate $\hat\Phi$ (maximiser of the log-integrated posterior, eq. 10, computed via Newton's method). The result is a closed-form multivariate normal (eq. 11): $f(y_{n+1}|X_{n+1},\mathbb{T}) \equiv N_p[\mu,\;\Sigma], \qquad \mu = X_{n+1}\hat{S}_{XX}^{-1}\hat{S}_{Xy}, \quad \Sigma = \hat\Phi^{-1} + X_{n+1}\hat{S}_{XX}^{-1}X_{n+1}^T$

Missing responses. If some $y_i$ components are unobserved, their conditional given the observed components is multivariate normal by standard partitioning. Adding this as a fourth Gibbs block extends the sampler; the first-order approximation extends by redefining $\hat\Phi^{(0)}$ on the observed data only (eq. 13).

Simulation. Four bivariate ( $p=2$ ) cases ( $n=10$ ): A–C have known exact Bayesian densities (degenerate sub-models); D is the true SUR case. Contour plots show both approximations closely matching the exact density for A–C, and the two approximations agreeing with each other for D.

Triangular Transformation and Direct Monte Carlo (Zellner et al. 1988; Zellner-Ando 2010)

Reparameterize sequentially: $y_1 = X_1\beta_1 + e_1, \qquad y_j = X_j\beta_j + \sum_{l=1}^{j-1}\rho_{jl}u_l + e_j, \quad j \geq 2$ where $e_j \perp e_{j'}$ for $j \neq j'$ with $\text{Var}(e_j) = \sigma_j^2 I$ . Let $Z_j = [X_j, u_1, \ldots, u_{j-1}]$ and $b_j = (\beta_j', \rho_{j1}, \ldots, \rho_{j,j-1})'$ . Under Jeffreys' prior $\pi_3(b, \Sigma) \propto \prod_j (\sigma_j^2)^{-1}$ , the conditionals are: $b_j | b_{j-1},\ldots,b_1, \sigma_j^2, D \sim \mathcal{N}\!\left(\hat{b}_j,\; \sigma_j^2 (Z_j'Z_j)^{-1}\right)$ $\sigma_j^2 | b_{j-1},\ldots,b_1, D \sim \text{IG}\!\left(\tfrac{\hat\lambda_j}{2},\, \tfrac{\hat\nu_j}{2}\right)$ with $\hat{b}_j = (Z_j'Z_j)^{-1}Z_j'y_j$ , $\hat\lambda_j = (y_j - Z_j\hat{b}_j)'(y_j - Z_j\hat{b}_j)$ , $\hat\nu_j = n - m - p_j + j + 1$ .

Key property: these conditionals are independent across equations, enabling direct sequential sampling — no iteration required. The Jacobian of the transformation from $\Omega$ to $\Sigma$ is $|J_m| = \prod_j (\sigma_j^2)^{m-j}$ , so $\pi_1$ and $\pi_3$ are not the same prior.

DMC Algorithm

Set $j = 1$ . Sample $\sigma_1^{2(k)}$ from $\text{IG}$ and $b_1^{(k)}$ from $\mathcal{N}$ , for $k = 1,\ldots,N$ .
Increment $j$ ; sample $\sigma_j^{2(k)}$ and $b_j^{(k)}$ from the conditional posteriors using $\{b_1^{(k)},\ldots,b_{j-1}^{(k)}\}$ .
Repeat until $j = m$ . Transform $\{b^{(k)}, \Sigma^{(k)}\}$ back to $\{\beta^{(k)}, \Omega^{(k)}\}$ via the recursive relation $\omega_j^2 = \sum_k \rho_{jk}^2 \omega_k^2 + \sigma_j^2$ (DMC2).

All $N$ draws are independent (inefficiency factor (INEF) = 1); no burn-in, no convergence checking, no proposal density needed.

Hierarchical SUR via MCMC (Chib-Greenberg 1995b)

Chib and Greenberg (1995b) embed the SUR in a three-level hierarchical prior and estimate it by MCMC, then extend to correlated-error and time-varying-parameter variants.

Hierarchical structure: $y_t = X_t \beta + \varepsilon_t, \quad \varepsilon_t \sim \mathcal{N}_p(0, \Omega)$ $\beta = A_0 \beta_0 + u, \quad u \sim \mathcal{N}_k(0, B_0)$ $\beta_0 = A_1 \mu + \eta, \quad \eta \sim \mathcal{N}_m(0, B_1)$

All full conditionals are standard (Normal and Wishart), so the Gibbs sampler converges by the Roberts-Smith (1993) sufficient conditions. The hierarchy accommodates cross-equation pooling (set $A_0 = (I_k, \ldots, I_k)'$ ) and model adequacy checking (replace $B_0 = \tau^2 I_k$ , draw $\tau^2$ as a scalar shrinkage parameter).

Vector autoregressive (VAR(1)) errors: When $\varepsilon_t = \Phi \varepsilon_{t-1} + u_t$ , the quasi-differenced data $(y_t^* = y_t - \Phi y_{t-1},\; X_t^* = X_t - \Phi X_{t-1})$ restore the standard Gibbs conditionals for $(\beta, \Omega)$ . The AR matrix $\Phi$ has a matrix-normal full conditional truncated to the stationary region and is sampled by accept-reject within Gibbs.

Vector moving-average (VMA(1)) errors: When $\varepsilon_t = u_t + \Theta u_{t-1}$ , the MA matrix $\Theta$ has a non-standard posterior (proportional to the innovations likelihood). A Metropolis-within-Gibbs step uses the candidate generating density $q(\theta | Y_n, \Omega) \propto |V|^{-1/2} \exp\!\left[-\tfrac{1}{2}(\theta - \hat\theta)' V^{-1} (\theta - \hat\theta)\right] \cdot I_V(\Theta)$ derived by linearising the innovations $u_t(\theta)$ around the nonlinear least squares (NLS) estimate $\bar\Theta$ (eq. 12). The acceptance rate is approximately 50% in simulations.

Time-varying-parameter (TVP)-SUR: Let $\beta_t = A_0 \theta_t + u_t$ with $\theta_t = \theta_{t-1} + \eta_t$ (random walk). The key efficiency result (eq. 16): sample all states $\{\theta_t\}$ jointly from $p(\theta_0, \theta_1, \ldots, \theta_n | \{\beta_t\}, \psi)$ via the factorisation (eq. 17) and the Kalman backward smoother. Given $\{\beta_t\}$ , the model is a linear Gaussian state-space system — the Kalman filter runs forward and then backward draws are: $\theta_t \mid \{\beta_t\}, \theta_{t+1}, \psi \sim \mathcal{N}_k\!\left(\hat\theta_t, R_t\right), \quad t < n$ where $\hat\theta_t = \hat\theta_{t|t} + M_t(\theta_{t+1} - \hat\theta_{t+1|t})$ , $R_t = R_{t|t} - M_t R_{t+1|t} M_t'$ , and $M_t = R_{t|t} R_{t+1|t}^{-1}$ . This is an early derivation of what Carter and Kohn (1994) call the forward filter backward sampler (FFBS) smoother. Sampling the joint path is one Gibbs block; sampling each $\theta_t$ individually (the naive approach) requires $n+1$ blocks and mixes orders of magnitude slower.

Partial Bayes factors (Dawid 1984): Use the first $n_0$ observations as a training sample to form a proper posterior for all competing models; then estimate partial marginal densities for the remaining data by Monte Carlo averaging (eq. 7). Avoids the computational complexity of full Bayes factors while remaining simulation-consistent.

Organisation for Economic Co-operation and Development (OECD) application: Five-country gross national product (GNP) growth (1960–1987), 40 parameters total. Strong evidence for the pooled shrinkage model (log partial Bayes factor $\approx 40$ ); posterior of the random-walk persistence $\phi$ concentrated near zero, confirming constant parameters.

Multi-Set SUR with Cross-Set Coefficient Pooling (Griffiths-Valenzuela 2002)

Griffiths and Valenzuela (2002) extend the standard single-set SUR to H groups of equations, where each group shares a common coefficient vector $\eta$ with all other groups while retaining its own set-specific parameters $\Theta_h$ and error covariance $\Omega_h$ . The model is:

$Y_h = Z_h\Theta_h + X_h\eta + e_h, \quad e_h \sim \mathcal{N}(0,\, \Omega_h \otimes I_{M_h}), \quad h = 1,\ldots,H$

Key structural distinction: each set has its own unrestricted $\Omega_h$ , unlike Baltagi (1995) who assumes a common $\Omega$ across sets. The $Z_h$ matrices contain set-varying regressors (mapped to $\Theta_h$ ); the $X_h$ matrices contain set-invariant regressors (mapped to $\eta$ ) and may differ across sets in dimension.

All-conjugate three-block Gibbs sampler:

Block 1: $\Omega_h \mid \text{rest} \sim \text{IW}$ — one draw per set from the Inverted-Wishart with scale $A_h$ = within-set residual outer product
Block 2: $\Theta_h \mid \text{rest} \sim \mathcal{N}(\hat\Theta_h,\, [Z_h'(\Omega_h^{-1}\otimes I)Z_h]^{-1})$ — within-set GLS draw after removing the common $\eta$ contribution
Block 3: $\eta \mid \text{rest} \sim \mathcal{N}(W^{-1}Q,\, W^{-1})$ where $W = \sum_h X_h'(\Omega_h^{-1}\otimes I_{M_h})X_h$ ; this pooled-GLS update is the methodological innovation: it aggregates cross-set information for $\eta$ in a single conjugate Normal draw

Applications:

Extended Linear Expenditure System (ELES) equivalence scales (Application 1): H=8 Australian demographic household types, $N_h$ =11 commodity equations per type, 5,532 households; $\eta$ = marginal budget shares (common), $\Theta_h$ = subsistence quantities (type-specific); 18,000 draws, 3,000 burn-in
Merino wool translog cost (Application 2): H=23 years, time-varying intercepts $\Theta_h$ , common technology parameters $\alpha_i$ , $\alpha_{ij}$ in $\eta$ ; 23,000 draws, 3,000 burn-in

See Griffiths-Valenzuela (2002) and Equivalence Scale.

MCMC Algorithms and Extensions (Griffiths 2001)

Griffiths (2001) provides a unified treatment of Bayesian inference for the M-equation SUR under the noninformative prior $f(\beta,\Sigma) \propto |\Sigma|^{-(M+1)/2}$ . The joint posterior is:

$f(\beta,\Sigma|y) \propto |\Sigma|^{-(T+M+1)/2} \exp\!\left\{-\tfrac{1}{2}\operatorname{tr}(A\Sigma^{-1})\right\}$

where $[A]_{ij} = (y_i - X_i\beta_i)'(y_j - X_j\beta_j)$ . The marginal $f(\beta|y) \propto |A|^{-T/2}$ is intractable, motivating MCMC.

Full conditionals:

$\beta|\Sigma \sim \mathcal{N}(\hat\beta_{\text{GLS}},\; [X'(\Sigma^{-1}\otimes I_T)X]^{-1})$ where $\hat\beta_{\text{GLS}}$ is the GLS estimator
$\Sigma|\beta \sim \text{IW}(T, A)$

Algorithm 1 — Joint Gibbs( $\beta,\Sigma$ ): Alternate draws from these two full conditionals. Straightforward but requires updating the full $KM\times KM$ covariance at each step.

Algorithm 2 — Equation-by-equation Gibbs: Marginalise $\Sigma$ analytically for each equation. The conditional $\beta_i|\beta_{-i}$ follows a multivariate t:

$\beta_i|\beta_{-i},y \sim t_{v_i}\!\left(\tilde\beta_i,\; \tilde{V}_i\right), \quad v_i = T - K_i$

where $\tilde\beta_i$ and $\tilde{V}_i$ depend on the cross-equation residuals (Griffiths eq. 19). This eliminates the need to draw $\Sigma$ explicitly and reduces block size from $KM$ to $K_i$ .

Algorithm 3 — Metropolis-Hastings random walk: Draw a candidate $\beta^* = \beta^{(t-1)} + \epsilon$ , $\epsilon\sim\mathcal{N}(0,cI)$ , and accept with probability $\min(1, f(\beta^*|y)/f(\beta^{(t-1)}|y))$ . The only algorithm applicable to nonlinear SUR, where $A$ is not quadratic in $\beta$ .

Extensions:

Inequality restrictions $g(\beta)\geq 0$ : multiply the prior by the indicator $I_S(\beta)$ ; for linear restrictions $R\beta=r$ substitute $\beta = Z\gamma + e$ to reduce dimension; for general inequality restrictions, $\beta|\Sigma$ follows a truncated multivariate normal (MVN) sampled by rejection.
Nonlinear SUR: only Algorithm 3 applies; the GLS and IW conditionals are no longer standard.
Missing data/Tobit: augment the Gibbs cycle with conditional normal draws for unobserved $y_i$ components.

Predictive density: $f(y^*|\beta,y)$ is multivariate t with $v^* = T-M+1$ df, mean $X^*\beta$ , and scale $A/(v^*-2)$ . In practice, average over MCMC draws of $\beta$ (contrast with Percy 1992, who samples $(y^*,\beta,\Sigma)$ jointly in Gibbs).

Applications: wheat yield (5 WA shires, mild inequality → Algorithm 2); translog cost for Merino woolgrowers (310 observations $\times$ 23 years, severe monotonicity constraints → Algorithm 3); nonlinear expenditure functions (1,834 Bangkok households → Algorithm 3).

See Griffiths (2001) and William E. Griffiths.

Sample Size Requirements (Griffiths-Skeels-Chotikapanich 2001)

The standard textbook condition — $T \geq \max(M, k_{\max} + 1)$ — is both incomplete and potentially seriously misleading. The correct requirements differ by estimator and depend on the rank structure of the full combined regressor matrix $X = [X_1,\ldots,X_M]$ .

Key parameters:

$\rho = \text{rank}([X_1,\ldots,X_M])$ — rank of the combined regressor matrix ( $\geq k_{\max}$ , $\leq T$ )
$\omega$ = dimension of the regressor subspace common to all equations ( $\omega > 0$ whenever all equations share regressors, e.g., an intercept)
$d$ = number of equations whose regressors do not form a basis for $X$
$\eta = \text{rank}(D)$ where $D = [P_{Z_1}y_1,\ldots,P_{Z_M}y_M]$ captures the projection of each $y_j$ onto the part of $X$ -space not spanned by $X_j$ ; with probability one $\eta = \min(d, \rho - \omega)$

Theorem 1 — Two-Stage Feasible GLS (FGLS) (necessary condition for $\hat\Sigma$ nonsingular with probability one):

$T \geq M + \rho - \eta \quad \Longleftrightarrow \quad T \geq \begin{cases} M + \omega & \text{if } \eta = \rho - \omega \\ M + \rho - d & \text{if } \eta = d \end{cases}$

The key insight: OLS residuals $\hat{e}_j = M_{X_j}y_j$ decompose as $\hat{e}_j = M_V y_j + P_{Z_j}y_j$ , so $\hat\Sigma = \hat{E}'\hat{E} = Y'M_V Y + D'D$ . The rank of $D$ (= $\eta$ ) relaxes the requirement relative to the multivariate regression case.

Theorem 2 — maximum likelihood (ML) and Bayesian (necessary condition for a bounded likelihood / proper posterior):

$T \geq M + \rho$

This equals the two-stage requirement only when $\eta = 0$ (all equations share the same regressors). When equations have distinct regressors ( $\eta > 0$ ), ML and Bayesian need strictly more observations. Under the noninformative prior $f(\beta,\Sigma) \propto |\Sigma|^{-(M+1)/2}$ , the marginal posterior $f(\beta|y) \propto |S|^{-T/2}$ is improper at $\beta$ 's where $S = E'E$ is singular — which necessarily happens when $T < M + \rho$ .

Why the gap? Two-stage minimizes a quadratic in $\beta$ ; ML/Bayes minimize $|S|$ , a polynomial of degree $2M$ in $\beta$ . The higher-dimensional criterion has larger minimal information requirements for M > 1.

Special cases:

Multivariate regression ( $X_1 = \cdots = X_M = X^*$ ): both estimators require $T \geq M + k^*$ (gap vanishes)
Fully distinct regressors (M=2, each equation has its own unique regressor): two-stage requires $T \geq 2$ ; ML/Bayes require $T \geq M + \rho$
M equations, each with intercept + 2 unique regressors, T=19: two-stage works up to M=18; ML/Bayes break down at $M \geq 7$ — confirmed numerically with SHAZAM and Bayesian Gibbs sampler

Practical warning: Standard software does not reliably detect the failure mode. SHAZAM sometimes found local maxima instead of detecting unboundedness. The Bayesian Gibbs sampler got stuck in a narrow nonsensical range or broke from singularities — without obvious error messages. Users must verify $T \geq M + \rho$ before ML/Bayesian SUR estimation.

Relation to moments: Existence of moments of the two-stage estimator (needed for valid standard inference) requires the stronger condition $T > M + \rho + 1$ (Srivastava-Raj 1979), which is more demanding than mere feasibility.

See Griffiths-Skeels-Chotikapanich (2001), William E. Griffiths, Christopher L. Skeels, Duangkamon Chotikapanich.

Robust Estimation: S-Estimators (Bilodeau-Duchesne 2000)

Zellner's GLS estimator has zero breakdown point (BP) — one outlier can move estimates arbitrarily. Bilodeau and Duchesne (2000) adapt S-estimators (Rousseeuw-Yohai 1984) to the full SUR system. The S-estimator solves:

$\min_{(\beta,\Sigma)} |\Sigma|, \quad \text{subject to} \quad \frac{1}{n}\sum_{i=1}^n \rho\!\left[(e_i'\Sigma^{-1}e_i)^{1/2}\right] = b$

where $b = E_{F_0}\rho(|r|)$ at the elliptical target $E_q(0,I)$ and the $\rho$ -function is the biweight. The estimating equations downweight observations by $u(d_i)=\rho'(d_i)/d_i$ in the GLS formula and in the scatter update — reducing to the maximum likelihood estimator (MLE) when $\rho(d)=d^2$ . Key properties: (i) affine equivariance (detects multivariate outliers, unlike Koenker-Portnoy M-estimators); (ii) BP up to 50%; (iii) $\sqrt{n}$ -consistent and asymptotically normal, enabling bootstrap standard error (SE). Asymptotic efficiency relative to MLE: at BP=40%, $\hat\beta$ efficiency $\approx 74\%$ , $\hat\Sigma$ efficiency $\approx 58\%$ .

Grunfeld diagnostic: GE+Westinghouse 1935–1954 ( $n=20$ ). Robust $\hat\rho=0.85$ vs. MLE's 0.77. Bivariate Mahalanobis distance plot of residuals identifies 1946 (WWII transition) and 1952 (GE jet-engine expansion) as joint multivariate outliers — invisible in univariate residual plots.

Event Study Application: Bayesian SUR with Equicorrelation (Brav 2000)

Brav (2000) applies SUR to long-horizon event studies where N firms (e.g., 113 or 1,521 initial public offerings (IPOs) grouped by industry) are followed over T months. The key features differ from the classical SUR setup:

Regression: $y_i = f_i + v_i$ (firm return = benchmark return + residual); the "regressor" is a single benchmark portfolio return. Stacked as $Y = F\iota + V$ .
Covariance structure: Rather than estimating a full $N \times N$ matrix (infeasible for large $N$ ), Brav imposes $\Sigma = SRS$ where $S$ is diagonal with firm-specific standard deviations (SDs) and $R$ is an equicorrelation matrix with a single common $\rho$ per industry. This reduces inference to $N + 1$ parameters per industry.
Prior: Lognormal Empirical Bayes prior on $\sigma_i$ : $\log(\sigma_i) \sim N(\bar s, \delta_\sigma)$ centered at the within-industry grand mean SD. This is a random-effects / hierarchical shrinkage prior that pools extreme firm estimates toward the industry average.
Non-conjugate sampling: The conditionals for $\rho$ and each $\sigma_i$ are non-standard (not closed-form Normal or IG). Both are sampled via Griddy-Gibbs: evaluate the conditional on a grid, approximate the inverse cumulative distribution function (CDF), and invert a Uniform draw.
Predictive use: The posterior draws are used not for point estimation but to simulate the predictive distribution of the sample mean buy-and-hold return — a posterior predictive check applied to abnormal performance testing.

Cross-sectional correlation $\rho \approx 2.5$ – $2.7\%$ within the computer/data services industry — small relative to prior studies but widens the predictive null density tails by 10–12 percentage points. See Brav (2000) and Long-Horizon Event Study.

Extension to Nonstationary Systems: DSUR (Mark-Ogaki-Sul 2003)

When the regressors $x_{it}$ are I(1) and $y_{it}$ cointegrates with $x_{it}$ , OLS is biased because $\Delta x_{it}$ is correlated with the equilibrium error $u^y_{it}$ (endogeneity). The fix — leads and lags of $\Delta x$ augmented into each equation — is the Dynamic OLS (DOLS) idea. Dynamic SUR (DSUR) applies this to the full SUR system:

Augment each equation $i$ with $\Delta x_{jt\pm p}$ for all $j = 1,\ldots,N$ (not just own regressors).
Apply Zellner GLS across the stacked augmented system using estimated $\hat\Omega$ .

The result inherits all properties of classical SUR (GLS efficiency gain from cross-equation correlation) while also correcting for nonstationarity-induced endogeneity. The cross-equation $\Delta x$ augmentation step is essential: DOLS stacked naively (System DOLS / SDOLS) without cross-equation $\Delta x$ leads/lags does not achieve asymptotic efficiency when errors are correlated. DSUR is asymptotically equivalent to the semiparametric MLE (Mark-Ogaki-Sul 2003, Proposition 4), and offers 46–73% mean squared error (MSE) reductions over equation-by-equation DOLS when cross-equation error correlation is high.

The key practical caveat: DSUR is infeasible for large $N$ because $\Omega$ has $N(N+1)/2$ free parameters. See Dynamic Seemingly Unrelated Regression for full details.

Why It Matters

SUR is the standard framework for systems of regressions in macro, finance, and actuarial science whenever equations share correlated disturbances (e.g., seemingly unrelated industry output equations, multi-country VARs with a small number of observations per country, or mortality models for multiple age groups).
Gibbs sampling for SUR is unreliable when $p_j$ is large: Zellner-Ando (2010) report 92/100 convergence failures at $p_j = 25$ , $n = 50$ , with INEF = 31.6 for off-diagonal $\Omega$ elements. DMC eliminates this problem entirely.
The triangular transformation generalizes: it underlies efficient inference in simultaneous equations, hierarchical SUR models, and Student- $t$ SUR extensions.

Open Questions

Extension to nonlinear SUR (different link functions per equation) requires equation-specific DMC adaptations; the triangular structure no longer yields standard conditionals.
With very large $m$ (many equations), the sequential structure of DMC is inherently serial — parallelization is limited.
Prior sensitivity: the non-equivalence of $\pi_1$ and $\pi_3$ means the two DMC variants make different prior assumptions; which is more appropriate depends on the context.