How can population variance be estimated from a bivariate sample? Let's assume a bivariate population with a correlation $\rho$ and a common $\sigma$ so that $\Sigma = \sigma^2 \begin{pmatrix}1 & \rho \\ \rho & 1\end{pmatrix}$.
I would like to know the different ways that could be used to estimate $\sigma^2$ from bivariate samples of size $n$ taken from this population? and what is their distribution? @whuber here said he could see at least three ways. I assume that the parameter $\rho$ is unknown although the sample correlation $r$ can be used; also, the result should be an unbiased estimator of $\sigma^2$.
I can see two ways only:

*

*Use the pooled variance $S_p^2$, that is, the mean of the two variances, $S_1^2$ for the first elements of the pairs, $S_2^2$ for the second element of the pairs. Ben here demonstrated that $S_p^2/ \sigma^2 \sim \chi^2_v / v$ where $v = \frac{2(n-1)}{1+\rho^2}$.


*A second way would be to flatten all the $n$ pairs into a single sample of size $2 n$ and compute the variance from this list. What is the distribution of this estimator, I do not know... It mixes uncorrelated items (from distinct pairs) with correlated items.


*I have no idea how else $\sigma^2$ could be estimated.

Assuming that $\rho$ is known (as opposed to being a nuisance parameter), I would think that the first method you mention would give quite reasonable results.  However, I would caution against assuming normality of the population, since inferences for the variance are heavily affected by skewness and kurtosis.  In particular, the variability of the sample variance is determined by the kurtosis of the underlying distribution, and so your inference would be highly sensitive to your normality assumption.  To deal with this, you can adopt a method set out in O'Neill (2014) to adjust a confidence interval for the variance to take account of the (known or estimated) kurtosis in the underlying population.  This method gives the standard results for a normal population, but generalises to allow for non-mesokurtic populations.
Using this proposed method, suppose we let $\kappa<\infty$ denote the kurtosis of the underlying distribution (here assumed to be finite), which then gives the asymptotic approximation:
$$\frac{S_\text{pooled}^2}{\sigma^2} \overset{\text{Approx}}{\sim} \frac{\chi_{DF_n}^2}{DF_n} 
\quad \quad \quad \quad \quad DF_n \equiv \frac{2}{1+\rho^2} \cdot \frac{2n}{\kappa - (n-3)/(n-1)}.$$
This asymptotic approximation should be reasonable for any population with a finite kurtosis, so long as $n$ is not too small.  We can re-write this asymptotic approximation in the more useful form:
$$\frac{\sigma^2}{DF_n S_\text{pooled}^2} \sim \text{InvGamma} \bigg( \text{Shape} = \frac{DF_n}{2}, \text{Scale} = 2 \bigg).$$
At any stipulated confidence level $1-\alpha$ we can obtain optimal values $0 \leqslant L_{\alpha} \leqslant U_{\alpha}$ for the critical points of the inverse-gamma distribution that minimise the interval width (implicitly depending on $DF_n$).  We have:
$$\begin{align}
1-\alpha
&= \mathbb{P} \bigg( L_{\alpha} \leqslant \frac{\sigma^2}{DF_n S_\text{pooled}^2} \leqslant U_{\alpha} \bigg) \\[6pt]
&= \mathbb{P} \bigg( L_{\alpha} DF_n \cdot S_\text{pooled}^2 \leqslant \sigma^2 \leqslant U_{\alpha} DF_n \cdot S_\text{pooled}^2 \bigg). \\[6pt]
\end{align}$$
Substitution of the observed value of the pooled sample variance then gives an optimal confidence interval:
$$\text{CI}_n(1-\alpha) \equiv \Bigg[ L_{\alpha} DF_n \cdot s_\text{pooled}^2 , \ U_{\alpha} DF_n \cdot s_\text{pooled}^2 \Bigg].$$
If the kurtosis is known then you use this as a substitute value in $DF_n$ and if not you can substitute an appropriate estimator.  (Adding an estimated kurtosis makes the interval slightly more variable, so you can make some crude adjustments to the degrees-of-freedom in this case if you want.)

Programming the CI: You can program this confidence interval in a manner similar to the function CONF.var in the stat.extend package.  Here is some basic code programmed in a similar style, but without the full generality and input checks of these other confidence intervals.  This function takes in two data vectors and stipulated values for the correlation and population kurtosis and it returns the confidence interval at a specified confidence level.
CONF.var.correlated <- function(conf.level, x0, x1, corr = 0, kurt = 3) {
  
  #Check input data
  n  <- length(x0); nn <- length(x1)
  if (n != nn)       stop('Error: Inputs x0 and x1 must have the same length')
  if (n < 3)         stop('Error: Method requires at least three data points')
  
  #Calculated degrees-of-freedom and optimal critical points
  DF <- (4*n)/((kurt - (n-3)/(n-1))*(1+corr^2))
  HH <- stat.extend::HDR.invgamma(cover.prob = conf.level, shape = DF/2, scale = 2)
  LL <- min(HH); UU <- max(HH)
  
  #Compute the pooled variance
  pooled.variance <- var(c(x0, x1))
  
  #Generate the confidence interval
  CONF     <- pooled.variance*DF*sets::interval(l = LL, r = UU, bounds = 'closed')
  METHOD   <- attributes(HH)$method
  DATADESC <- paste0('Interval uses ', n, 
                 ' pairs of data points with pooled sample variance = ',
                 sprintf(pooled.variance, fmt = '%#.4f'), ', assumed correlation = ',
                 sprintf(corr, fmt = '%#.4f'), ' and assumed kurtosis = ',
                 sprintf(kurt, fmt = '%#.4f'))
  
  #Add class and attributes and return confidence interval
  class(CONF) <- c('ci', 'interval');
  attr(CONF, 'confidence') <- conf.level
  attr(CONF, 'data')       <- DATADESC
  attr(CONF, 'method')     <- METHOD
  attr(CONF, 'parameter')  <- 'variance parameter for infinite population'
  attr(CONF, 'pooled.var') <- pooled.variance
  attr(CONF, 'DF')         <- DF
  CONF }

Here is an example where we use this confidence interval function on some mock data.    We generate $n=100$ pairs of data points with population variance $\sigma^2=9$ and correlation $\rho = 0.2$.  In the present case we can see that the 95% confidence interval includes the true value of the population variance.
#Load libraries
library(invgamma)
library(stat.extend)

#Set parameters
n    <- 40
VAR  <- 9
CORR <- 0.2

#Generate mock data
set.seed(1)
VAR.MATRIX <- matrix(VAR*c(1, CORR, CORR, 1), nrow = 2, ncol = 2)
XX <- mvtnorm::rmvnorm(n, sigma = VAR.MATRIX)
x0 <- XX[,1]
x1 <- XX[,2]

#Generate 95% confidence interval for variance
CONF <- CONF.var.correlated(conf.level = 0.95, x0, x1, corr = 0.2, kurt = 3)


    Confidence Interval (CI) 
 
95.00% CI for variance parameter for infinite population 
Interval uses 40 pairs of data points with pooled sample variance = 7.3948, 
assumed correlation = 0.2000 and assumed kurtosis = 3.0000 
Computed using nlm optimisation with 6 iterations (code = 1) 

[5.29264696608237, 10.153814031653]

