Modeling Covariance Matrices via Partial Autocorrelations

2.1 Reparameterization in terms of Partial Autocorrelations

The notion of PACF is known to be indispensable in the study of stationary processes and situations dealing with Toeplitz matrices such as the Szegö's orthogonal polynomials, trigonometric moment problems, geophysics, digital signal processing and filtering (see Landau, 1987; Pourahmadi, 2001), identification of ARMA models, the maximum likelihood estimation of their parameters (Jones, 1980) and simulating a random or “typical” ARMA model (Jones, 1987). The one-to-one correspondence between the stationary autocorrelation functions {ρ_k} and their PACF {π_k} (Barndorff-Nielsen and Schou, 1973; Ramsey, 1974) makes it possible to remove the positive-definiteness constraint on {ρ_k}, and work with {π_k} which are free to vary over the interval (−1,1) independently of each other.

We parameterize a (non-Toeplitz) correlation matrix R = (ρ_ij) in terms of the lag-1 correlations π_i,i+1 = ρ_i,i+1,i = 1, ⋯, p − 1 and the partial autocorrelations π_ij = ρ_{ij|i+1,⋯,j−1} for j − i ≥ 2, or the matrix Π = (π_ij). This allows swapping the constrained matrix R by the simpler matrix Π with ones on the diagonal and where for i ≠ j, the π_ij's can vary freely in the interval (−1,1). The key idea behind this reparameterization is the well-known recursion formula (Anderson, 2003, p.41), see also (1) below, for computing partial correlations in terms of the marginal correlations (ρ_ij). It also lies at the heart of Kurowicka and Cooke (2003, 2006) and Joe (2006) approaches to constructing Π; however the recursive Levinson-Durbin algorithm used by Dégerine and Lambert-Lacroix (2003) will be used in our presentation in Section 2.3.

Following Joe (2006), for j = 1, ⋯, p−k, k = 1, ⋯, p−1, let $r_1^{\prime }\left(j,k\right)=({\rho }_{j,j+1},\dots ,{\rho }_{j,j+k-1}),r_3^{\prime }\left(j,k\right)=({\rho }_{j+k,j+1},\dots ,{\rho }_{j+k,j+k-1})$, and R₂(j,k) be the correlation matrix corresponding to the components (j + 1,…, j + k−1). Then, the partial autocorrelations between Y_j and Y_j+k adjusted for the intervening variables, denoted by ${\pi }_{j,j+k}\equiv {\rho }_{j,j+k|j+1,\dots ,j+k-1}$, are computed using the expression

In what follows and in analogy with R, it is convenient to arrange these partial autocorrelations in a matrix Π where its (j,j + k)th entry π_j,j+k Note that the function g(·) in (1) that maps a correlation matrix R into the partial autocorrelation matrix Π, is indeed invertible, so that solving (1) for ρ_j,j+k, one obtains

where D_jk is the denominator of the expression in (1). Then, the formulae (1)-(2) clearly establish a one-to-one correspondence between the matrices R and Π. In the sequel, with a slight abuse of language and following the tradition in times series analysis we refer to the matrix Π or π_j,j+k viewed as a function of (j, k), as the partial autocorrelation function (PACF) of Y.

Evidently, when R is a stationary (Toeplitz) correlation matrix, then π_j,j+k depends only the lag k, see (1). Consequently, Π is a stationary (Toeplitz) matrix. Fortunately, the converse is also true and follows from (2). For ease of reference, we summarize these observation in Lemma 1 in Section 2.3. A correlation matrix R is stationary (Toeplitz) if and only if its associated PACF Π is a stationary (Toeplitz) matrix.

Moreover, for a stationary correlation matrix, R reduces precisely to the celebrated Levinson-Durbin formula (Pourahmadi, 2001, Theorem 7.3) for computing the PACF recursively.

2.2 An Alternative Reparameterization: Cholesky Decomposition

Next, we present an alternative reparameterization of a covariance matrix via its Cholesky decomposition or the idea of autoregression for the underlying random vector.

Consider a mean-zero random vector Y with the positive-definite covariance matrix Σ = (σ_st). For1 ≤ t ≤ p, let Ŷ_t be the linear least-squares predictor of Y_t based on its predecessors Y₁, ⋯, Y_t−1 and let ε_t = Y_t −Ŷ_t be its prediction error with variance ${\sigma }_t^2=\mathit{\text{Var}}\left({\epsilon }_t\right)$. Then, there are unique scalars φ_tj so that ${\widehat{Y}}_t=\text{∑}_{j=1}^{t-1}{\phi }_{tj}{Y}_j$ or

Let ε = (ε₁, ⋯, ε_p)′ be the vector of successive uncorrelated prediction errors with $\mathit{\text{Cov}}(\epsilon )=\mathit{\text{diag}}\left({\sigma }_1^2,\cdots ,{\sigma }_p^2\right)=D$. Then, (3) rewritten in matrix form becomes ε = TY, where T is a unit lower triangular matrix with 1's on the main diagonal and −φ_tj in the (t,j)th position for 2 ≤ t ≤ p, j = 1, ⋯, t −1. Note that ${\sigma }_t^2=\mathit{\text{Var}}\left({\epsilon }_t\right)$ is different from σ_tt = Var(Y_t) However, when the responses are independent, then φ_tj = 0 and ${\sigma }_t^2={\sigma }_{tt}$, so that the matrices T and D gauge the “dependence” and “heterogeneity” of Y, respectively.

Computing covariances using ε = TY, it follows that

The first factorization in (4), called the modified Cholesky decomposition of Σ, makes it possible to swap the p(p + 1)/2 constrained parameters of Σ with the unconstrained set of parameters φ_tj and log ${\sigma }_t^2$ of the same cardinality. In view of the similarity of (3) to a sequence of varying order autoregressions, we refer to the parameters φ_tj and ${\sigma }_t^2$ as the generalized autoregressive parameters (GARP) and innovation variances (IV) of Y or Σ (Pourahmadi, 1999). A major advantage of (4) is its ability to guarantee the positive-definiteness of the estimated covariance matrix given by T̂⁻¹D̂T̂′⁻¹ so long as the diagonal entries of D̂ are positive.

It should be noted that imposing structures on Σ will certainly lead to constraints on T and D in (4). For example, a correlation matrix R with 1's as its diagonal entries is structured with possibly p(p − 1)/2 distinct parameters. In this case, certain entries of T and D are either known, redundant or constrained. In fact, it is easy to see that the diagonal entries of the matrix D for a correlation matrix are monotone decreasing with ${\sigma }_1^2=1$ For this reason and others, it seems more prudent to rely on the ordered partial correlations when reparameterizing a correlation matrix R as in Section 2.1, than using its Cholesky decomposition.

2.3 A Multiplicative Determinantal Identity: Partial Autocorrelations

First, we study the role of partial autocorrelations in measuring the reduction in prediction error variance when a variable is added to the set of predictors in a regression model. Using this and the second identity in (4) we obtain a fundamental determinantal identity expressing |Σ| in terms of the partial autocorrelations and diagonal entries of Σ. Joe (2006) and Kurowicka and Cooke (2006) had obtained this identity using determinantal recursions and graph-theoretical methods based on (1), respectively. An earlier and a slightly more general determinantal identity for covariance matrices in the context of nonstationary processes was given by Dégerine and Lambert-Lacroix (2003, p.54), using an analogue of the Levinson-Durbin algorithm.

For u and v two distinct integers in {1, 2, ⋯, p}, let L be a subset of {1, 2, ⋯, p}\{u, v} and π_uv|L stand for the partial correlation between Y_u and Y_v adjusted for Y_ℓ, ℓεL. We denote the linear least squares predictor of Y_u based on Y_ℓ, ℓεL by Ŷ _u|L, and for v an integer Lv stands for the union of the set L and the singleton {v}.

Lemma 1. Let Y = (Y₁, ⋯, Y_p)′ be a mean-zero random vector with a positive-definite covariance matrix Σ. Then, with u,v and L as above, we have

1. $${\widehat{Y}}_{u|Lv}={\widehat{Y}}_{u|L}+{\alpha }_{uv}\left(Y_v-{\widehat{Y}}_{v|L}\right),{\alpha }_{uv}={\pi }_{uv|L}\sqrt{\frac{\mathit{\text{Var}}\left(Y_u-{\widehat{Y}}_{u|L}\right)}{\mathit{\text{Var}}\left(Y_v-{\widehat{Y}}_{v|L}\right)}}.$$

   (5)

2. $$\mathit{\text{Var}}\left(Y_u-{\widehat{Y}}_{u|Lv}\right)=\left(1-{\pi }_{uv|L}^2\right)\mathit{\text{Var}}\left(Y_u-{\widehat{Y}}_{u|L}\right).$$

   (6)

Proof. (a) Let sp{Y_u; uεL} stand for the linear subspace generated by the indicated random variables. Since Y_v − Ŷ_v|L is orthogonal to sp{Y_u; uεL} it follows that

and from the linearity of the orthogonal projection we have

where α_uv, the regression coefficient of Y_u on Y_v − Ŷ_v|L is given by

Since Ŷ_u|Lε sp{Y_i; iεL} is orthogonal to Y_v − Ŷ_v|L, the numerator of the expression above can be replaced by Cov(Y_u −Ŷ _u|L, Y_v −Ŷ_v|L), so that

(b) From the first identity in (5), it is immediate that

or

Computing variances of both sides of (7) and using the fact that Y_u −Ŷ_u|Lv is orthogonal to Y_v−Ŷ_v|L, it follows that

Now, substituting for α_uv from (a), the desired result follows.

Next, we use the recursion (6) to express the innovation variances or the diagonal entries of D in terms of the partial correlations or the entries of π. Similar expressions for φ_tj's, the entries of T, are not available. The next theorem sheds some light on this problem. Though the expressions are recursive and ideal for computation (Levinson-Durbin algorithm), they are not as explicit or revealing. The approach we use here is in the spirit of the Levinson-Durbin algorithm (Pouramadi, 2001, Corollary 7.4) as extended by Dégerine and Lambert-Lacroix (2003) to nonstationary processes.

Theorem 1. Let Y = (Y₁, ⋯, Y_p) be a mean-zero random vector with a positive-definite covariance matrix Σ which can be decomposed as in (4).

1. Then, for t = 2, ⋯, p; j = 1, ⋯, t − 1, we have

   $${\sigma }_t^2={\sigma }_{tt}\prod _{j=1}^{t-1}\left(1-{\pi }_{jt}^2\right),$$

2. $$\left|\sum \right|=\left(\prod _{t=1}^p{\sigma }_{tt}\right)\prod _{t=2}^p\prod _{j=1}^{t-1}\left(1-{\pi }_{jt}^2\right).$$

3. For t = 2, ⋯, p; L = {2, ⋯, t − 1}

   $${\phi }_{t1}={\pi }_{1t}\sqrt{\frac{\mathit{\text{Var}}\left(Y_t-{\widehat{Y}}_{t|L}\right)}{\mathit{\text{Var}}\left(Y_1-{\widehat{Y}}_{1|L}\right)}}\cdot$$

4. φ_tj = φ_tj|L − φ_t1φ_1,t−j|L, for j = 2, ⋯, t −1,

where φ_tj|L and φ_1,t−j|L are, respectively, the forward and backward predictor coefficients of Y_t and Y₁ based on {Y_k; kεL}, defined by

Proof (a) follows from the repeated use of (6) with u = t and v = j,j = 1, ⋯, t −1. (b) follows from (4) and (a).

Part (c) proved first in Pourahmadi (2001, p.102) shows that only the entries of the first column of T are multiples of the ordered partial correlations appearing in the first column of π. However, for j > 1, since φ_tj is a multiple of the partial correlation between Y_j and Y_t adjusted for {Y_i; i ε [1, t)\{j}}, (see Lemma 1(a)), it is not of the form of the entries of π. Note that these observations are true even when Y is stationary or Σ is Toeplitz (Pourahmadi, 2001, Lemma 7.8 and Theorem 7.3). In the search for connection with the entries of π, it is instructive to note that the nonredundant entries of T⁻¹ = (θ_tj) can be interpreted as the generalized moving average parameters (GMAP) or the regression coefficient of ε_j when Y_t is regressed on the innovations ε_t, ⋯, ε_j, ⋯, ε₁, see Pourahmadi (2001, p.103). An alternative interpretation of θ_tj as the coefficient of Y_j when Y_t is regressed on Y_j, Y_j−1, ⋯, Y₁ is presented in Wermuth et al. (2006 Sec. 2.2). Consequently, θ_tj is a multiple of the partial correlation between Y_t and Y_j adjusted for {Y₁, ⋯, Y_j−1}. We hope these connections and working with partial correlations will offer similar advantage to working with the GARP in terms of the autoregression interpretation (Pourahmadi, 2001, Sections 3.5.3 and 3.5.4).

2.4 An Attractive Property of the PACF parameterization

Parsimonious modeling of the GARP of the modified Cholesky decomposition often relies on exploring for structure as a function of lag; for example, fitting a polynomial to the GARP as a function of lag (Pourahmadi, 1999). For such models, the GARP, in a sense, have different interpretations within lag; i.e., the lag 1 coefficient from the regression of Y₃ on (Y₂, Y₁) is the Y₂ coefficient, when another variable, Y is also in the model; however, for the regression of Y₂ on Y₁, the lag 1 coefficient is the Y₁ coefficient with no other variable in the model. So, for a given lag k, the lag k coefficients all come from conditional regressions where the number of variables conditioned on are different. However, by construction, the lag k partial autocorrelations are all based on conditional regressions where the number of conditioning variables are the same, always conditioning on k − 1 intervening variables. This will facilitate building models for the partial autocorrelations as a function of lag. We discuss such model building in the next section.

2.5 Parsimonious Modeling of the PACF

In this section, we use the generalized partial correlogram, i.e. the plot of {π_j,j+k; j = 1, ⋯, p−k} versus k = 1, ⋯, p−1, as a graphical tool to formulate parsimonious models for the PACF in terms of the lags and other covariates. If necessary, we transform its range to the entire real line using z transform. Such modeling environment is much simpler and avoids working with the complex constraints on the correlation matrix R (Barnard et al., 2000; Daniels and Normand, 2006) or the matrix of full partial correlations constructed from Σ⁻¹ (Wong et al. 2003); note that the full partial correlations are defined as the correlation between two components conditional on all the other components.

Note that the partial autocorrelations π_j,j+k between successive variables Y_j and Y_j+k are grouped by their lags k = 1, ⋯, p−1, and heuristically, π_j,j+k gauges the conditional (in)dependence between variables k units apart conditional on the intervening variables, so one expects it to be smaller for larger k. In the Bayesian framework, this intuition suggests putting shrinkage priors on the partial autocorrelations that shrink the matrix Π toward certain simpler structures (Daniels and Kass, 2001).

2.6 Data illustration

To illustrate the capabilities of the generalized correlograms in revealing patterns, we use the cattle data (Kenward, 1987) which consists of p = 11 bi-weekly measurements of the weights of n = 30 cows. Table 1, displays the sample (partial) correlations for the cattle data in the lower (upper) triangular segment and the sample variances are along the main diagonal. It reveals several interesting features of the dependence in the data that the commonly used profile plot of the data cannot discern. For example, note that all the correlations are positive, they decrease monotonically within the columns (time-separation), they are not constant (nonstationary) within each subdiagonal. In fact, they tend to increase over time (learning effect). Furthermore, the partial autocorrelations of lags 2 or more are insignificant except for the entries 0.56 and 0.35.

Figure 1 presents the generalized correlograms corresponding to the sample correlation matrix of the data, the full partial correlations, the generalized partial correlogram and the Fisher's z transform of the PACF. Note that the first two correlograms suggest linear and quadratic patterns in the lag k, but in fitting such models one has to be mindful of the constraints on the coefficients so that the corresponding fitted correlation matrices are positive definite. Details of fitting such models and the ensuing numerical results can be found in Pourahmadi (2001). The generalized partial correlogram in (c) reveals a cubic polynomial in the lags, i.e. π_j,j+k = γ₀ + γ₁k + γ₂k² + γ₃k³; in fitting such models the only constraint to observe is that the entries of the matrix Π are required to be in (−1,1). However, the Fisher z transform of the entries of Π are unconstrained and Figure 1(d) suggests a pattern that can be approximated by an (exponential) function α + β exp(−k), k = 1, ⋯, p−1, with no constraints on (α, β) or another cubic polynomials in the lags. The least-squares fits of a cubic polynomial and an exponential function to the correlograms in Figure 1 (c)-(d) are summarized in Table 2. Note that fitting such models to the PACF amounts to replacing the entries of the kth subdiagonal of the matrix Π by a single number and hence R is approximated by a stationary (Toeplitz) matrix, see Lemma 1 and Theorem 2 below. In addition, this parameterization also allows the marginal variances to be similarly modelled parsimoniously (as a function of time) similar to the modelling of the prediction variances in Pourahmadi (1999). The maximum likelihood estimation of the parameters and their asymptotic properties will be pursued in a future work.

3.1 Independent priors on partial autocorrelations

Given that each partial autocorrelation is free to vary in the interval (−1,1), we may construct priors for R derived from independent priors on the PACF. For example, independent Unif[-1,1] priors (IU priors) on the partial autocorrelations (Jones, 1987; Joe, 2006) can be shown to induce the following prior on the correlation matrix R:

One can express the prior in (8) in terms of the marginal correlations, ρ_j,j+k by plugging in for π_j,j+k from (1). This prior induces a particular behavior on the marginal correlations. Specifically, the priors on the (marginal) correlations ρ_j,j+k, become gradually more peaked at zero as the lag k grows. As an illustration of this behavior, Figure 2 shows the histograms based on 10,000 simulations from the uniform prior on the partial autocorrelations for p = 5. As the dimension p of the correlation matrix grows, the priors become more peaked at zero for larger lags. This can also be seen by examining the prior probability of being in some interval, say [−.5, .5], as a function of the lag. For p = 15, the (averaged) probabilities, ordered from lag 1 to lag 14, are respectively, (.50, .59, .65, .70, .73, .76, .78, .80, .82, .83, .84, .85, .86, .87). This would appear to be a desirable behavior for longitudinal data which typically exhibits serial correlation decaying with increasing lags.

It is also evident from Figure 2 that the priors for ρ_j,j+k with k fixed, appear to be the same (see the subdiagonals). We state this observation more formally in the following theorem; see also Lemma 1.

Theorem 2. The partial autocorrelations, π_j,k have independent “stationary” priors, i.e

(or the priors are the same along the subdiagonals of Π), if and only if the marginal priors on the correlations ρ_jk are also “stationary”, i.e.

(or the priors are the same along the corresponding subdiagonal of R).

This implies that “stationary” priors on Π induce “stationary” priors on the marginal correlations and vice versa. Most priors we introduce here satisfy this property.

More generally independent linearly transformed Beta priors on the interval (-1,1) for partial correlations are a convenient and flexible way to specify a prior for a correlation matrix R. These priors, denoted by Beta(α,γ), have the density

Interestingly, the uniform prior on the correlation matrix (Barnard et al, 2000) corresponds to the following stationary Beta priors on the partial correlations:

where ${\alpha }_k=1+\frac{1}{2}\left(p-1-k\right)$; see Joe (2006). We will refer to this prior as Barnard Beta (BB). As noted in Barnard et al. (2000), such a prior on R results in the marginal priors for each of the marginal correlations being somewhat peaked around zero (same peakedness for all ρ_jk). Also note that the priors become more peaked as p grows.

In general, priors for the correlation matrix proportional to powers of the determinant of the correlation matrix,

are constructed by setting ${\alpha }_k={\alpha }_p+\frac{1}{2}\left(p-1-k\right)$ in (12) Priors so constructed are proper, so improper priors like Jeffreys' for a correlation matrix in a multivariate normal model, π (R) = |R|^−(p+1)/2, are not special cases.

3.2 Shrinkage behavior of the BB and IU priors

The IU priors on the PACF induce desirable behavior for longitudinal (ordered data) by ‘shrinking’ higher lag correlations toward zero. The Beta priors in (12), which induce a uniform prior for R (BB priors) place a uniform (-1,1) prior on the lag p − 1 partial autocorrelations and shrink the other partial autocorrelations toward zero with the amount of shrinkage being inversely proportional to lag. This induces the desired behavior on the marginal correlations, making their marginal priors equivalent, but it is counter-intuitive for ordered/longitudinal data with serial correlation; in addition, the shrinkage of the lag one partial autocorrelations for the BB prior increases with p (recall the form in (12)). In such data, we would expect lower lag correlations to be less likely to be zero and higher lag correlations to be more likely to be zero. Thus, the independent uniform priors are likely a good default choice for the partial autocorrelations in terms of inducing desirable behavior on the marginal correlations and not counter-intuitively shrinking the partial autocorrelations. We explore this shrinkage behavior further via some simulations in Section 3.5.

In addition, we expect many of the higher lag partial correlations to be close to zero for longitudinal data with serial correlation via conditional independence (see Table 1). To account for this, we could (aggressively) shrink the partial autocorrelations toward zero (with the shrinkage increasing with lag) using shrinkage priors similar to those proposed in Daniels and Pourahmadi (2002) and Daniels and Kass (2001) or by creating such priors based on the Beta distributions proposed here. We are currently exploring this.

3.3 Some other priors for a correlation matrix

Other priors for R have been proposed in the literature which cannot be derived based on independent priors on the partial autocorrelations. For example, the prior on R that induces marginal uniform (-1,1) priors on the ρ_ij's (Barnard et al., 2000) has the form:

where R[−i,−i] is the submatrix of R with the ith row and column removed and $|R\left[-i,-i\right]={\left[R\right]}_{ii}^{-1}\left|R\right|$. Such a prior might not be a preferred one for longitudinal (ordered) data where the same marginal priors on all correlations (irrespective of lag) may not be the best default choice.

Eaves and Chang (1992) derived some related reference priors for the set of partial correlations, π₁,j for j = 2,…, p; however, their priors are not natural for longitudinal data. Chib and Greenberg (1998) specified a truncated multivariate normal distribution on the marginal correlations. Liechty et al. (2004) placed normal distributions on the marginal correlations with the goals of grouping the marginal correlations into clusters. The latter two priors along with those in Wong et al. (2003) for the full partial correlations (p^ij) are highly constrained given that they model the marginal correlations directly.

3.4 Bayesian Computing

An additional issue with modeling the correlation matrix is computational. Our development here will focus on cases without covariates in the correlation matrix (this will be left for future work) under the class of independent priors on the PACF discussed in Section 3.1. The proposal here might be viewed as an alternative to the PX-RPMH algorithm in Liu and Daniels (2006) that explicitly exploits the fact that we are modeling the partial correlations themselves (a computational comparison will be left for future work). However, our approach will naturally allow structures in the partial correlations which cannot be done when using current versions of the PX-RPMH (or similar) algorithms; for example, if the partial autocorrelations are zero or constant within lag since then the correlation matrix is highly constrained.

In the following, we assume the data, {Y_i : i = 1,…, n}are independent, normally distributed p-vectors with mean X_iβ and with covariance matrix Σ = R (a correlation matrix). A natural way to sample the partial autocorrelations is via a Gibbs sampling algorithm in which we sample from the full conditional distributions of each of the partial autocorrelations. Given that the full conditional distributions of the partial autocorrelations are not available in closed form there are several options to sample them. We explore a simple one next.

We propose to use an auxiliary variable approach to sample each partial autocorrelation. Define the likelihood for the partial autocorrelation, π_jk as L(π_jk) and the prior as p(π_jk). As in Damien et al. (1999), introduce a positive latent variable U_jk such that

To sample π_jk we can proceed in two steps,

1. Sample U_jk ∼ Unif(0,L(π_jk)).

2. Sample p(π_jk) constrained to the set {π_jk : L(π_jk) > U_jk}.

Truncated versions of the priors proposed here, linearly transformed Beta distributions (of which the uniform is a special case), can easily be sampled using the approach in Damien and Walker (2001). The truncation region for step 2, given that the domain of π_jk is bounded, can typically be found quickly numerically.

The likelihood evaluations needed to find the truncation interval can be made simpler by taking advantage of the following facts.

Fact 1. If we factor R⁻¹ = CPC, where P is a correlation matrix and C is a diagonal matrix, the elements of P are the full partial correlations, p_ij|rest (Anderson, 1984, Chapter 15).

Fact 2. To isolate the likelihood contribution of π_j,j+k, we can factor the entire multivariate normal distribution into p(y_j,…, y_j+k)p(y_l : l < j or l > j + k | y_j,…, y_j+k). The (l,k) entry of inverse of the correlation matrix, R[j : j +k] for the first factor is related to the partial correlation of interest (recall Fact 1).

Fact 3. Using the determinantal identity in Theorem 1(b), the determinant of submatrices of R in terms of partial correlations can be written as a function of the partial correlations,

It can be shown that the likelihood in terms of one of the partial correlations of interest, say π_jk, can be written as

where S_y[j : k] is the submatrix of Σ_i(y_i−X_iβ)(y_i−X_iβ)′ and R[j : k] is the submatrix of R based on the first component of the factorized distribution in Fact 2. above.

Extensions of these computational procedures to modelling the correlation matrix when the matrix of interest is a covariance matrix is straightforward (see, e.g., Liu and Daniels, 2006).

3.5 Simulations

We now conduct some simulations in a longitudinal setting to

1. examine the mixing behavior of the auxiliary variable sampler here and

2. compare the risk of the IU prior to the BB prior, the standard default prior for a correlation matrix.

In terms of the mixing of the Markov chain, the auxiliary variable sampler on the partial autocorrelations works quite well, with the lag correlation in the chain dissipating quickly. For example, for p = 5, n = 20, the lag correlations for each partial autocorrelation was negligible by lag 10. Similar results were seen for other p/n combinations.

For our simulation, we consider three true matrices representing typical serial correlation, an AR(1) with lag 1 correlation of .8 and one with lag correlation .6. Both these matrices have all partial autocorrelation beyond lag 1 equal to 0. We also considered a matrix that had more non-zero partial autocorrelations with lags 1,2,3,4 equal to (.8, .4, .1,0) respectively. We consider two size matrices/sample size combinations, p = 5 with n = 10,25,50,100 and p = 10 with sample sizes, n = 20,50,100. We also consider several loss functions, log likelihood (LL) loss, tr(R̂R⁻¹) − log|R̂R⁻¹| −p, with Bayes estimator the inverse of the posterior expectation of R⁻¹ and squared error loss on Fisher's z-transform of the partial autocorrelations, π_jk (SEL-P) and the marginal correlations, ρ_jk (SEL-M), with Bayes estimator the posterior mean (of the z-transformed correlations).

For p = 5, the risk reductions from the IU prior are clear from Table 3, with percentage reductions as large as 30% for n = 10, 20% for n = 25 and 10% for n = 50. For p = 10, the risk reductions from the uniform prior are clear from Table 4, with percentage reductions as large as 40%. The largest risk reductions were for loss SEL-M (squared error loss on the Fisher's z-transform on the marginal correlations). The lower risk reductions for the first order autoregressive covariance matrices are related to only the lag 1 partial autocorrelations being non-zero; so the shrinkage of the BB prior for all the other partial autocorrelations is not unreasonable. Examination of squared error loss for the partial autocorrelations by lag indicates large reductions for the lag 1 partial autocorrelations and small increases for the other lag partial autocorrelations.

In addition, the estimates of the first order lag correlation themselves show large differences (not shown). For AR(.8) with p = 10 and n = 20, the means were .76 under the IU prior and about .71 under the BB prior with similar discrepancies for p = 5.

Some of the risk reductions from using the IU priors on the partial autocorrelations are small. However, they come at no computational cost (unlike some priors for covariance matrices proposed in the literature) and are consistent with prior beliefs about partial autocorrelations representing serial correlation. The BB prior is not a good default choice due to its dampening effects on the important lower order partial autocorrelations. Further risk improvement might be expected through the use of more targeted shrinkage (Daniels and Pourahmadi, 2002).