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?
 A: 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)



