Why shouldn't the denominator of the covariance estimator be n-2 rather than n-1? The denominator of the (unbiased) variance estimator is $n-1$ as there are $n$ observations and only one parameter is being estimated.
$$
\mathbb{V}\left(X\right)=\frac{\sum_{i=1}^{n}\left(X_{i}-\overline{X}\right)^{2}}{n-1}
$$
By the same token I wonder why shouldn't the denominator of covariance be $n-2$ when two parameters are being estimated?
$$
\mathbb{Cov}\left(X, Y\right)=\frac{\sum_{i=1}^{n}\left(X_{i}-\overline{X}\right)\left(Y_{i}-\overline{Y}\right)}{n-1}
$$
 A: A quick and dirty answer... Let’s consider first $\text{var}(X)$; if you had $n$ observations with known expected value $E(X) = 0$ you would use ${1\over n}\sum_{i=1}^n X_i^2$ to estimate the variance. 
The expected value being unknown, you can transform your $n$ observations into $n-1$ observations with known expected value by taking $A_i = X_i - X_1$ for $i = 2, \dots,n$. You will get a formula with a $n-1$ in the denominator — however the $A_i$ are not independent and you’d have to take this into account; at the end you’d find the usual formula.
Now for the covariance you can use the same idea: if the expected value of $(X,Y)$ was $(0,0)$, you’d had a ${1\over n}$ in the formula. By subtracting $(X_1,Y_1)$ to all other observed values, you get $n-1$ observations with known expected value... and a ${1\over n-1}$ in the formula — once again, this introduces some dependence to take into account.
P.S. The clean way to do that is to choose an orthonormal basis of $\big\langle (1, \dots, 1)' \big\rangle^{\perp}$, that is $n-1$ vectors $c_1, \dots, c_{n-1} \in \mathbb R^n$ such that 


*

*$\sum_j c_{ij}^2 = 1$ for all $i$,

*$\sum_j c_{ij} = 0$ for all $i$, 

*$\sum_j c_{i_1j} c_{i_2j} = 0$ for all $i_1 \ne i_2$.


You can then define $n-1$ variables $A_i = \sum_j c_{ij} X_j$ and $B_i = \sum_j c_{ij} Y_j$. The $(A_i,B_i)$ are independent, have expected value $(0,0)$, and have same variance/covariance than the original variables.
All the point is that if you want to get rid of the unknown expectation, you drop one (and only one) observation. This works the same for both cases.
A: Here is a proof that the p-variate sample covariance estimator with denominator $\frac{1}{n-1}$ is an unbiased estimator of the covariance matrix: 
$ x' = (x_1,...,x_p) $.
$\Sigma= E((x-\mu)(x-\mu)') $
$S = \frac{1}{n} \sum (x_i - \bar{x})(x_i - \bar{x})'$
To show: $E(S) = \frac{n-1}{n}\Sigma$
Proof: $S= \frac{1}{n}\sum x_ix_i' - \bar{x}\bar{x}'$
Next:
(1) $ E(x_ix_i') = \Sigma + \mu\mu'$
(2) $E(\bar{x}\bar{x}') = \frac{1}{n} \Sigma+ \mu\mu' $
Therefore: $E(S) = \Sigma + \mu\mu' - (\frac{1}{n} \Sigma+ \mu\mu') = \frac{n-1}{n} \Sigma $
And so $S_u = \frac{n}{n-1}S $, with the final denominator $\frac{1}{n-1}$, is unbiased. The off-diagonal elements of $S_u$ are your individual sample covariances. 
Additional remarks:


*

*The n draws are independent. This is used in (2) to calculate the covariance of the sample mean.

*Step (1) and (2) use the fact that $Cov(x)= E[xx']-\mu\mu'$

*Step (2) uses the fact that $Cov(\bar{x})= \frac{1}{n}\Sigma$
A: I guess one way to build intuition behind using 'n-1' and not 'n-2' is - that for calculating co-variance we do not need to de-mean both X and Y, but either of the two, i.e.

A: Covariances are variances. 
Since by the polarization identity
$$\newcommand{\c}{\text{Cov}}\newcommand{\v}{\text{Var}}
\c(X,Y) = \v\left(\frac{X+Y}{2}\right) - \v\left(\frac{X-Y}{2}\right),$$
the denominators must be the same.
A: 1) Start $df=2n$. 
2) Sample covariance is proportional to $\Sigma_{i=1}^n(X_i-\bar{X})(Y_i-\bar{Y})$. Lose two $df$; one from $\bar{X}$, one from $\bar{Y}$ resulting in $df=2(n-1)$. 
3) However, $\Sigma_{i=1}^n(X_i-\bar{X})(Y_i-\bar{Y})$ only contains $n$ separate terms, one from each product. When two numbers are multiplied together the independent information from each separate number disappears. 
As a trite example, consider that 
$24=1*24=2*12=3*8=4*6=6*4=8*3=12*2=24*1$, 
and that does not include irrationals and fractions, e.g. $24=2\sqrt{6}*2\sqrt{6}$, so that when we multiply two number series together and examine their product, all we see are the $df=n-1$ from one number series, as we have lost half of the original information, that is, what those two numbers were before the pair-wise grouping into one number (i.e., multiplication) was performed. 
In other words, without loss of generality we can write 
$(X_i-\bar{X})(Y_i-\bar{Y})=z_i-\bar{z}$ for some $z_i$ and $\bar{z}$, 
i.e., $z_i=X_iY_i-\bar{X}Y_i-X_i\bar{Y}$, and, $\bar{z}=\bar{X}\bar{Y}$. From the $z$'s, which then clearly have $df=n-1$, the covariance formula becomes 
$\Sigma_{i=1}^n\frac{z_i-\bar{z}}{n-1}=$
$\Sigma_{i=1}^n\frac{[(X_i-\bar{X})(Y_i-\bar{Y})]}{n-1}=$
$\frac{1}{n-1}\Sigma_{i=1}^n(X_i-\bar{X})(Y_i-\bar{Y})$. 
Thus, the answer to the question is that the $df$ are halved by grouping.
A: A special case ought to give you an intuition; think about the following:
$$\hat{\mathbb{Cov}}\left(X, X\right)=
\hat{\mathbb{V}}\left(X\right)$$
You are happy that the latter is $\frac{\sum_{i=1}^{n}\left(X_{i}-\overline{X}\right)^{2}}{n-1}$ due to the Bessel correction.
But replacing $Y$ by $X$ in $\hat{\mathbb{Cov}}\left(X, Y\right)$ for the the former gives $\frac{\sum_{i=1}^{n}\left(X_{i}-\overline{X}\right)\left(X_{i}-\overline{X}\right)}{\text{mystery denominator}}$, so what do you now think might best fill in the blank?
A: My simple way of thinking about the intuition behind the 'n-1' degree of freedom is as follow. An observation is a pair $\{x_i, y_i\}$. Therefore, we can always infer the last observation $\{x_n, y_n\}$ if we know 'n-1' previous observations and the sample mean $\{\overline{x}, \overline{y}\}$. The last observation can be considered redundant.
$df = n-2$ should be wrong if we follow the same logic as above.
