4
$\begingroup$

UPDATED: I am constructing a correlation matrix for an MA(1) process, which would look something like...

$$ C = \left( \begin{array}{cccccccccccccccccc} 1 & \rho & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ \rho & 1 & \rho & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & \rho & 1 & \rho & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & \rho & 1 & \rho & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & \rho & 1 & \rho & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & \rho & 1 & \rho & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & \rho & 1 & \rho & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & \rho & 1 & \rho & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \rho & 1 & \rho & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \rho & 1 & \rho & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \rho & 1 & \rho & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \rho & 1 & \rho & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \rho & 1 & \rho & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \rho & 1 & \rho & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \rho & 1 & \rho & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \rho & 1 & \rho & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \rho & 1 & \rho\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \rho & 1\end{array} \right)$$

However, I'm using JAGS to fit the model and need to invert this matrix to sample from a multivariate normal distribution (after multiplying it by $\sigma^2$ of course). $C$ must be positive definite for it to be valid and I recognize that $C$ is a tridiagonal Toeplitz matrix whose eigenvalues have a closed form solution:

$$ \lambda_k = 1 + 2\rho\cos\left(\frac{k\pi}{n+1}\right), \mbox{ for }k=1,\cdots,n $$

where $n$ is the number of rows or columns (as the matrix is symmetric). I can use the prior distribution on $\rho$ to ensure that $C$ is positive definite and the fact that I can calculate the eigenvalues for $C$ to determine the bounds of that prior, but I am stuck and need help. Some values between $[-1,1]$ result in a matrix that is not positive definite. I would also like to extend the method to a pentadiagonal matrix for an MA(2) process.

$\endgroup$
9
  • 1
    $\begingroup$ Your formula for $\lambda_k$ cannot be right, because multiplying $C$ by $\sigma^2$ must multiply all eigenvalues by $\sigma^2$ rather than adding $\sigma^2$ to them. Apart from correcting that error, where exactly do you need help? $\endgroup$
    – whuber
    Nov 18, 2014 at 18:27
  • $\begingroup$ You misapplied the formula: when $C$ is multiplied by $\sigma^2$, the entries along the bands become $\sigma^2$ (along the diagonal) and $\sigma^2\rho$ (just off the diagonal). You seem to have forgotten to multiply the off-diagonal entries by $\sigma^2$. But if that's what you intended, then you need to redefine $C$ in your question. And if you are trying to apply any of this to an AR-1 process, then $C$ needs to be a different matrix altogether. $\endgroup$
    – whuber
    Nov 18, 2014 at 19:13
  • 1
    $\begingroup$ @user777 I basically need to determine what are the lowest and highest values of $\rho$ that lead to a positive definite matrix $C$, as $\sigma^2$ will always be positive. $\endgroup$
    – user13317
    Nov 18, 2014 at 21:11
  • 1
    $\begingroup$ Your formula for $\lambda_k$ is wrong, as @whuber has already pointed out. One cannot agree to disagree on that! Wikipedia says that the eigenvalues are $\lambda_k = a+2\sqrt{bc}\cos(k\pi/(n+1))$. In your case, $a=\sigma^2$ and $b=c=\rho\sigma^2$, meaning that eigenvalues are $\lambda_k=\sigma^2+2\sigma^2\rho\cos(k\pi/(n+1))$. Compare this with what you wrote. $\endgroup$
    – amoeba
    Nov 19, 2014 at 15:14
  • 1
    $\begingroup$ @amoeba you're right and so is @whuber! My equation only pertains to the correlation matrix and not the covariance matrix, my mistake. I edited the post and deleted my comments to reflect this. Thanks! $\endgroup$
    – user13317
    Nov 21, 2014 at 14:32

2 Answers 2

4
$\begingroup$

Regardless of whether the given covariance matrix correctly models the covariance matrix of an AR(1) process (or an MA(1) process) or not, the sum of all the entries in a covariance matrix is the variance of the sum of the $n$ random variables. Since this variance must be nonnegative, we get that in order for your matrix to be a valid covariance matrix, it must be that $$n + 2(n-1)\rho \geq 0 ~ \Rightarrow \rho \geq -\frac{n}{2(n-1)} \approx -\frac 12.$$ So it is certainly true that some choices of $\rho \in [-1,1]$ will not result in valid covariance matrices.


In hindsight, the OP's problem has an even simpler solution. Suppose that $Y_0, Y_1, Y_2, \cdots, Y_n$ are iid random variables with variance $\sigma^2$, and define $$X_i = aY_{i-1} + bY_i, ~ i = 1, 2, \cdots, n.$$ where $a^2+b^2 = 1$. It follows that $\operatorname{var}(X_i) = \sigma^2$ for $1 \leq i \leq n$, and more generally that $$\operatorname{cov}(X_i,X_{i+k}) = \operatorname{cov}(aY_{i-1}+bY_i, aY_{i+k-1}+bY_{i+k}) = \begin{cases}\sigma^2,& \text{if}~ k=0,\\ ab\sigma^2,& \text{if}~ k = \pm 1,\\ 0,&\text{if}~ |k| > 1, \end{cases}$$ that is, the covariance matrix of the $X_i$'s is what the OP wants.

So, subject to the constraint that $a^2+b^2 = 1$, what can we say about $ab = \rho_{X_i,X_{i+1}}$, the correlation between adjacent $X_i$? It is easy to deduce that $ab \in \left[-\frac 12, \frac 12\right]$ which is what the OP determined by experiment.

$\endgroup$
4
  • $\begingroup$ Just to check my understanding of your answer, does this imply that the constraint on $\rho$ should be $\rho\ge-\frac{n}{2(n-1)}$? $\endgroup$
    – Sycorax
    Nov 18, 2014 at 23:52
  • $\begingroup$ Thanks Dilip, that was exactly what I was looking for! I had determined the interval needed to be $[-0.5,0.5]$ by trial and error, it's nice to see the mathematical reasoning behind it though! $\endgroup$
    – user13317
    Nov 19, 2014 at 14:09
  • $\begingroup$ @DilipSarwate, could you please clarify why $ab\notin [−1,1]$? Since $\mbox{cov}(X_i,X_{i\pm 1})=ab\sigma^2$ and $\mbox{std}(X_i)=\sigma$, shouldn't $−1\leq \rho \leq 1$, where $\rho = \mbox{cov}(X_i,X_{i\pm 1})/(\mbox{std}(X_i)\mbox{std}(X_i\pm1))=ab\sigma^2/(\sigma^2)=ab$? $\endgroup$ Sep 19, 2019 at 12:54
  • $\begingroup$ @VivekSubramanian You are forgetting the constraint that $a^2+b^2=1$. So what we really have is $\rho = ab = a\sqrt{1-a^2}$ and we are seeking the possible values of $a\sqrt{1-a^2}$ subject to the constraint that $|a|\leq 1$. $\endgroup$ Sep 19, 2019 at 14:11
2
$\begingroup$

There are two problems here:

  1. The given matrix is not the covariance matrix of an AR(1) process. If it were, there would be no problem since a covariance matrix is always positive definite. The matrix you have given resembles the inverse of the AR(1) covariance matrix, but even this is not an exact match.
  2. You don't need to form the covariance matrix in order to specify an AR(1) process in JAGS.

In JAGS, you can specify an AR(1) model directly as:

for(i in 2:n) {
  x[i] ~ dnorm(theta*x[i-1],precision) 
}
$\endgroup$
3
  • $\begingroup$ I understand that a covariance matrix will always be positive definite. The question more so relates on how to ensure that during estimation through bounds via a prior distribution. $\endgroup$
    – user13317
    Nov 18, 2014 at 19:45
  • $\begingroup$ I think @Tom just gave you what you needed with the JAGS line. You will need an extra few lines for priors on theta and or the precision. $\endgroup$ Nov 18, 2014 at 20:50
  • $\begingroup$ @Tom_Minka OP has clarified that this is a question about MA(1) processes. Just FYI. $\endgroup$
    – Sycorax
    Nov 18, 2014 at 20:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.