Variance-covariance matrix of a single variable I want to calculate the variance-covariance matrix of a single variable:
$$
\begin{align*}
Var({\bf y}) & = E({\bf yy'}) - E({\bf y})E({\bf y'}) \\
  & = \left( \begin{array}{cccc}
\sigma_{y_{1}}^{2} & \sigma_{y_{1}y_{2}} & \cdots & \sigma_{y_{1}y_{n}} \\
\sigma_{y_{1}y_{2}} & \sigma_{y_{2}}^{2} & \cdots & \sigma_{y_{2}y_{n}} \\
\vdots             & \vdots              & \ddots & \vdots \\
\sigma_{y_{1}y_{n}} & \sigma_{y_{2}y_{n}} & \cdots & \sigma_{y_{n}}^{2}
\end{array} \right)\\
\end{align*}$$
for a single variable ${\bf y}$:
$$
\begin{align*}
{\bf y} = \left( \begin{matrix}
y_1 \\ y_2 \\ \vdots \\ y_n \end{matrix} \right)
\end{align*}
$$
but I do not understand how to calculate $\sigma_{y_{1}y_{2}}$ since $y_{1}$ and $y_{2}$ are just two numbers. 
Any help?
Update:
Take for example the residuals of an OLS:
$$ \begin{align*}
{\bf e} &= ({\bf I} - {\bf X} ({\bf X}^\prime {\bf X})^{-1} {\bf X}^\prime) {\bf y}\\
&= ({\bf I} - {\bf H}) {\bf y}
\end{align*}
$$
For 10 observations, ${\bf e}$ is a $10 \times 1$ matrix. And ${\bf e} = (e_1 \quad e_2 \quad \ldots \quad  e_{10})^\prime$. 
$e_1$ is just a number, like $e_2$. How do I calculate the second term in first row of the variance-covariance matrix of the residuals, $Cov(e_1, e_2)$? 
The residual variance-covariance matrix is $({\bf I} - {\bf H}){\bf \Omega}({\bf I} - {\bf H})^{\prime}$, but I would like to learn a more "manual" way if possible.
 A: You're confusing the index of a variable with an index of an observation. In your first equation, there are no indices of observations: $y_1$ refers to the first variable of the vector of variables $y_i$. Hence, you have $\sigma_{y_1y_2}$ refers to the covariance of the variables $y_1,y_2$.
In your update, the variable $e$ is just one variable, and $e_1,e_2,\dots$ refer to the observations of this variable. So, $e_1$ is the first observation of the single variable $e$, it's not the first variable in the vector of variables. Hence, the your first equation does not apply to this case. In fact, the only thing that's related to this case is the variance $\sigma^2_{y_1}$ - the variance of the variable $y_1$.
To compute the variance of a variable $e$ you don't use the first equation, you use other relevant equations. For instance the standard estmator of the sample variance $s^2=\frac{1}{n-1}\sum_{i=1}^n(e_i-\bar e)^2$
A: They are not numbers, they are random variables. For example, they might represent different draws from the population $Y$, so that you get multiple random variables even if you just consider a single variable. (This seems to be one of the instances in which it might have been useful to distinguish between capital letters for the random variables and lowercase letters for their realizations.)
Suppose $n=2$, $\mathbf{Y}=(Y_1,Y_2)'$ is multivariate normal with means of zero, variances of 1 and a correlation of $\rho$.
Then, $E(Y_1)=E(Y_2)=0$, $Var(Y_1)=Var(Y_2)=1$ and $Corr(Y_1,Y_2)=\rho$. Since the variances are 1, the correlation is also the covariance $\sigma_{Y_1,Y_2}$ in this example.
