# Validity of confidence interval for $\rho$ when $X\sim N_3(0,\Sigma)$ with $\Sigma_{ij}=\rho^{|i-j|}$

Suppose $$X\sim N_3(0,\Sigma)$$, where $$\Sigma=\begin{pmatrix}1&\rho&\rho^2\\\rho&1&\rho\\\rho^2&\rho&1\end{pmatrix}$$.

On the basis of one observation $$x=(x_1,x_2,x_3)'$$, I have to obtain a confidence interval for $$\rho$$ with confidence coefficient $$1-\alpha$$.

We know that $$X'\Sigma^{-1}X\sim \chi^2_3$$.

So expanding the quadratic form, I get

$$x'\Sigma^{-1}x=\frac{1}{1-\rho^2}\left[x_1^2+(1+\rho^2)x_2^2+x_3^2-2\rho(x_1x_2+x_2x_3)\right]$$

To use this as a pivot for a two-sided C.I with confidence level $$1-\alpha$$, I set $$\chi^2_{1-\alpha/2,3}\le x'\Sigma^{-1}x\le \chi^2_{\alpha/2,3}$$

I get two inequalities of the form $$g_1(\rho)\le 0$$ and $$g_2(\rho)\ge 0$$, where

$$g_1(\rho)=(x_2^2+\chi^2_{\alpha/2,3})\rho^2-2(x_1x_2+x_2x_3)\rho+x_1^2+x_2^2+x_3^2-\chi^2_{\alpha/2,3}$$

and $$g_2(\rho)=(x_2^2+\chi^2_{1-\alpha/2,3})\rho^2-2(x_1x_2+x_2x_3)\rho+x_1^2+x_2^2+x_3^2-\chi^2_{1-\alpha/2,3}$$

Am I right in considering a both-sided C.I.? After solving the quadratics in $$\rho$$, I am guessing that the resulting C.I would be quite complicated.

Another suitable pivot is $$\frac{\mathbf1' x}{\sqrt{\mathbf1'\Sigma \mathbf 1}}\sim N(0,1)\quad\,,\,\mathbf1=(1,1,1)'$$

With $$\bar x=\frac{1}{3}\sum x_i$$, this is same as saying $$\frac{3\bar x}{\sqrt{3+4\rho+2\rho^2}}\sim N(0,1)$$

Using this, I start with $$\left|\frac{3\bar x}{\sqrt{3+4\rho+2\rho^2}}\right|\le z_{\alpha/2}$$

Therefore, $$\frac{9\bar x^2}{3+4\rho+2\rho^2}\le z^2_{\alpha/2}\implies 2(\rho+1)^2+1\ge \frac{9\bar x^2}{z^2_{\alpha/2}}$$

That is, $$\rho\ge \sqrt{\frac{9\bar x^2}{2z^2_{\alpha/2}}-\frac{1}{2}}-1$$

Since the question asks for any confidence interval, there are a number of options available here. I could have also squared the standard normal pivot to get a similar answer in terms of $$\chi^2_1$$ fractiles. I am quite sure that both methods I used are valid but I am not certain whether the resulting C.I. is a valid one. I am also interested in other ways to find a confidence interval here.

@JimB had commented that the intervals might not give real numbers within $$[-1,1]$$. Is there a way I can actually ensure that it does not happen?

• I think you want to use $x' \Sigma^{-1} x \leq \chi^2_{1-\alpha,3}$ and then solve the quadratic for the lower and upper bounds for $\rho$. Try it with actual samples from the trivariate normal and you'll see that the quadratic solution when using the bound on the left-hand side of your equation typically results in imaginary numbers and occasionally real numbers but outside of the range -1 to 1.
– JimB
Commented Mar 31, 2019 at 21:04
• @JimB I had not considered that. Will check again. Commented Mar 31, 2019 at 21:10
• You should also check the other confidence interval approach with values for $\bar{x}$. With values close enough to zero you're going to get imaginary numbers and large values of $\bar{x}$ will get you bounds way outside -1 and +1. Estimating $\rho$ with just a single sample of 3 correlated values is not going to have very "tight" coverage.
– JimB
Commented Mar 31, 2019 at 21:22
• Even using the $x'\Sigma^{-1}x \leq \chi^2_{1-\alpha,3}$ approach results in between 1% and 3% of the estimates being imaginary numbers. (At least my klutzy implementation ends up with those figures.)
– JimB
Commented Mar 31, 2019 at 21:45
• By the way, in general terms, each pivot gives rise to a CI right (apart from these issues). How do we know than what is the "best" CI ? I take the occasion to make this question here :) Commented Jan 12, 2022 at 8:37

Normally, inverting the likelihood ratio test statistic or Rao's score test statistic is a nice technique to keep the parameter within its parameter space. However, since we have a sample size of 1 and the exact distribution of these test statistics is unknown in this case, these methods will fail miserably.

It's times like these when the frequentist approach should give way to the Bayesian approach. The Bayesian approach will develop a credible interval, not a confidence interval, but it is guaranteed to be within its parameter space. Since we have a sample size of 1, I recommend using a uniform prior for $$\rho$$, namely $$\rho \sim Uni(-1,1)$$. Hence the posterior density of $$\rho$$ will be equivalent to the likelihood, up to a multiplicative constant. That is, $$\begin{eqnarray*} p(\rho | \boldsymbol{x} ) = \frac{1}{C} \exp\left(-\frac{1}{2} \left[\log |\boldsymbol{\Sigma}| + \boldsymbol{x}^{\prime} \boldsymbol{\Sigma}^{-1} \boldsymbol{x}\right]\right), \end{eqnarray*}$$ where the constant $$C$$ is given by $$C = \int_{-1}^1 p(\rho | \boldsymbol{x} ) \mbox{d} \rho$$. This may be found by numerical integration.

In order to construct a two-sided $$100(1-\alpha)\%$$ credible interval, we can use a univariate root solver to find the values of $$\rho_l$$ and $$\rho_u$$ such that $$\alpha/2 =\int_{-1}^{\rho_l} p(\rho | \boldsymbol{x} ) \mbox{d} \rho$$ and $$\alpha/2 =\int_{\rho_u}^{1} p(\rho | \boldsymbol{x} ) \mbox{d} \rho$$. The $$100(1-\alpha)\%$$ credible interval is then $$(\rho_l, \rho_u)$$. Here is some R code to accomplish this.

library(MASS)

p = .5
Sigma = matrix(c(1,p,p^2,p,1,p,p^2,p,1),3,3)
x = mvrnorm(1,rep(0,3),Sigma)

prop.post = function(p,x){
sigma = matrix(c(1,p,p^2,p,1,p,p^2,p,1),3,3)
h=exp(-.5*(as.numeric(determinant(sigma, logarithm = TRUE)$$mod) + as.numeric(t(x)%*%solve(sigma)%*%x))) return(h) } v.prop.post = Vectorize(prop.post,vectorize.args="p") con = integrate(prop.post,-1,1,x=x)$$value

true.post = function(p,x,con){
sigma = matrix(c(1,p,p^2,p,1,p,p^2,p,1),3,3)
h=exp(-.5*(as.numeric(determinant(sigma, logarithm = TRUE)$mod) + as.numeric(t(x)%*%solve(sigma)%*%x))) return(h/con) } v.true.post = Vectorize(true.post,vectorize.args="p") a = .05 ci.l = function(x2,x,con,a){ return(integrate(v.true.post,-1,x2,x=x,con=con)$$value - a/2) } ci.u = function(x2,x,con,a){ return(integrate(v.true.post,x2,1,x=x,con=con)value - a/2) } v1=uniroot(ci.l,c(-.9999,.9999),x=x,con=con,a=a)$$root v2=uniroot(ci.u,c(-.9999,.9999),x=x,con=con,a=a)$root

s = seq(-.99,.99,length=1000)
y = v.true.post(s,x,con)
plot(y~s,type="l")
abline(v=v1,lty=2)
abline(v=v2,lty=2)