Is every covariance matrix positive definite? I guess the answer should be yes, but I still feel something is not right. There should be some general results in the literature, could anyone help me?
 A: No. 
Consider three variables, $X$, $Y$ and $Z = X+Y$. Their covariance matrix, $M$, is not positive definite, since there's a vector $z$ ($= (1, 1, -1)'$) for which $z'Mz$ is not positive.
Population covariance matrices are positive semi-definite.
(See property 2 here.)
The same should generally apply to covariance matrices of complete samples (no missing values), since they can also be seen as a form of discrete population covariance.
However due to inexactness of floating point numerical computations, even algebraically positive definite cases might occasionally be computed to not be even positive semi-definite; good choice of algorithms can help with this.
More generally, sample covariance matrices - depending on how they deal with missing values in some variables - may or may not be positive semi-definite, even in theory. If pairwise deletion is used, for example, then there's no guarantee of positive semi-definiteness. Further, accumulated numerical error can cause sample covariance matrices that should be notionally positive semi-definite to fail to be.
Like so:
 x <- rnorm(30)
 y <- rnorm(30) - x/10 # it doesn't matter for this if x and y are correlated or not
 z <- x+y
 M <- cov(data.frame(x=x,y=y,z=z))
 z <- rbind(1,1,-1)
 t(z)%*%M%*%z
              [,1]
[1,] -1.110223e-16

This happened on the first example I tried (I probably should supply a seed but it's not so rare that you should have to try a lot of examples before you get one).
The result came out negative, even though it should be algebraically zero. A different set of numbers might yield a positive number or an "exact" zero.
--
Example of moderate missingness leading to loss of positive semidefiniteness via pairwise deletion:
z <- x + y + rnorm(30)/50  # same x and y as before.
xyz1 <- data.frame(x=x,y=y,z=z) # high correlation but definitely of full rank 

xyz1$x[sample(1:30,5)] <- NA   # make 5 x's missing  

xyz1$y[sample(1:30,5)] <- NA   # make 5 y's missing  

xyz1$z[sample(1:30,5)] <- NA   # make 5 z's missing  

cov(xyz1,use="pairwise")     # the individual pairwise covars are fine ...

           x          y        z
x  1.2107760 -0.2552947 1.255868
y -0.2552947  1.2728156 1.037446
z  1.2558683  1.0374456 2.367978

 chol(cov(xyz1,use="pairwise"))  # ... but leave the matrix not positive semi-definite

Error in chol.default(cov(xyz1, use = "pairwise")) : 
  the leading minor of order 3 is not positive definite

 chol(cov(xyz1,use="complete")) # but deleting even more rows leaves it PSD

          x          y          z
x 0.8760209 -0.2253484 0.64303448
y 0.0000000  1.1088741 1.11270078
z 0.0000000  0.0000000 0.01345364

A: As the other answer note, the covariance matrix is positive semi-definite (which I prefer to call non-negative definite), but not necessarily positive definite.  We can show that the covariance matrix is positive semi-definite from first principles using its definition.  To do this, suppose we consider a random vector $\mathbf{X}$ with mean vector $\boldsymbol{\mu}$ and covariance matrix $\mathbf{\Sigma}_\mathbf{X}$.  For any conformable vector $\mathbf{z}$ we can define the corresponding vector:
$$\mathbf{Y} = (\mathbf{X} - \boldsymbol{\mu}_\mathbf{X})^\text{T} \mathbf{z}.$$
Since $||\mathbf{Y}|| \geqslant 0$ we then have:
$$\begin{aligned}
\mathbf{z}^\text{T} \mathbf{\Sigma}_\mathbf{X} \mathbf{z}
&= \mathbf{z}^\text{T} \mathbb{E}((\mathbf{X} - \boldsymbol{\mu}_\mathbf{X}) (\mathbf{X} - \boldsymbol{\mu}_\mathbf{X})^\text{T}) \mathbf{z} \\[6pt]
&= \mathbb{E}(\mathbf{z}^\text{T} (\mathbf{X} - \boldsymbol{\mu}_\mathbf{X}) (\mathbf{X} - \boldsymbol{\mu}_\mathbf{X})^\text{T} \mathbf{z}) \\[6pt]
&= \mathbb{E}(\mathbf{Y}^\text{T} \mathbf{Y}) \\[6pt]
&= \mathbb{E}(||\mathbf{Y}||^2) \geqslant 0. \\[6pt]
\end{aligned}$$
This establishes that the covariance matrix $\mathbf{\Sigma}_\mathbf{X}$ is positive semi-definite.  Moreover, we can see that $\mathbf{z}^\text{T} \mathbf{\Sigma}_\mathbf{X} \mathbf{z} = 0$ if and only if $\mathbf{Y}=(\mathbf{X} - \boldsymbol{\mu}_\mathbf{X})^\text{T} \mathbf{z}=\mathbf{0}$ almost surely.
A: Well, to understand why the covariance matrix of a population is always positive semi-definite, notice that:
$$
\sum_{i,j =1}^{n} y_i \cdot y_j \cdot Cov(X_i, X_j) = Var(\sum_{i=1}^n y_iX_i) \geq 0
$$
where $y_i$ are some real numbers, and $X_i$ are some real valued random variables.
This also explains why in the example given by Glen_b the covariance matrix was not positive definite . We had $y_1 =1 , y_2 = 1, y_3 = -1$, and $X_1 = X, X_2 = Y, X_3 = Z = X+Y$, so $\sum_{i=1}^{3} y_iX_i = 0$, and the variance of a random variable which is constant is $0$.
A: As the other answers already make clear, a covariance matrix is not necessarily positive definite, but only positive semi-definite.
However, a covariance matrix is generally positive definite unless the space spanned by the variables is actually a linear subspace of lower dimension. This is exactly why in the example with X, Y and Z=X+Y the result is only positive semi-definite, but not positive definite. Although the variables span a three-dimensional space, they actually describe only a two-dimensional linear subspace (because they are not linearly independent).
A: $$\begin{array}{l}theory:\left\{ {{{\bf{\Sigma }}_{\bf{X}}}{\rm{ is positive semi - definite}}} \right.\\proof::\\set:\left\{ {{\bf{a}} = {\rm{vector }}\left( {p \times 1} \right){\rm{ }}\left( {{\mathop{\rm const}\nolimits} } \right) \ne \vec 0} \right.\\{{\bf{a}}^T}\Sigma {\bf{a}} = {\left[ {\begin{array}{*{20}{c}}{{a_1}}\\{{a_2}}\\ \vdots \\{{a_p}}\end{array}} \right]^T}\left[ {\begin{array}{*{20}{c}}{{\sigma _{11}}}&{{\sigma _{12}}}& \cdots &{{\sigma _{1p}}}\\{{\sigma _{21}}}&{{\sigma _{22}}}& \cdots &{{\sigma _{2p}}}\\ \vdots & \vdots & \ddots & \vdots \\{{\sigma _{p1}}}&{{\sigma _{p2}}}& \cdots &{{\sigma _{pp}}}\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{{a_1}}\\{{a_2}}\\ \vdots \\{{a_p}}\end{array}} \right]\\ = \left[ {\begin{array}{*{20}{c}}{{a_1}{\sigma _{11}} + {a_2}{\sigma _{21}} +  \cdots  + {a_p}{\sigma _{p1}}}& \cdots & \cdots &{{a_1}{\sigma _{1p}} + {a_2}{\sigma _{2p}} +  \cdots  + {a_p}{\sigma _{pp}}}\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{{a_1}}\\{{a_2}}\\ \vdots \\{{a_p}}\end{array}} \right]\\ = \left[ {\begin{array}{*{20}{c}}{\sum\limits_{i = 1}^p {{a_i}{\sigma _{i1}}} }& \cdots & \cdots &{\sum\limits_{i = 1}^p {{a_i}{\sigma _{ip}}} }\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{{a_1}}\\{{a_2}}\\ \vdots \\{{a_p}}\end{array}} \right]\\ = \left[ {\begin{array}{*{20}{c}}{\sum\limits_{i = 1}^p {{a_i}{\mathop{\rm Cov}\nolimits} \left( {{X_i},{X_1}} \right)} }& \cdots & \cdots &{\sum\limits_{i = 1}^p {{a_i}{\mathop{\rm Cov}\nolimits} \left( {{X_i},{X_p}} \right)} }\end{array}} \right]\left[ {\begin{array}{*{20}{c}}{{a_1}}\\{{a_2}}\\ \vdots \\{{a_p}}\end{array}} \right]\\ = \sum\limits_{j = 1}^p {{a_j}\sum\limits_{i = 1}^p {{a_i}{\mathop{\rm Cov}\nolimits} \left( {{X_i},{X_j}} \right)} }  = \sum\limits_{j = 1}^p {\sum\limits_{i = 1}^p {{a_i}{a_j}{\mathop{\rm Cov}\nolimits} \left( {{X_i},{X_j}} \right)} } \\\left[ {rule:\left\{ {{\mathop{\rm var}} \left( {\sum\limits_{i = 1}^p {{a_i}{X_i}} } \right) = \sum\limits_{j = 1}^p {\sum\limits_{i = 1}^p {{a_i}{a_j}{\mathop{\rm cov}} \left( {{X_i},{X_j}} \right)} } } \right.{\rm{ }}\left[ {see{\rm{ }}below} \right]} \right.\\ = {\mathop{\rm Var}\nolimits} \left( {\sum\limits_{i = 1}^p {{a_i}{X_i}} } \right) = {\mathop{\rm Var}\nolimits} \left( {{{\bf{a}}^T}{\bf{X}}} \right) \ge 0\end{array}
$$

Following is followed by the answer from sjm.majewski
$$\begin{array}{l}rule:\left\{ {{\mathop{\rm var}} \left( {\sum\limits_{i = 1}^p {{a_i}{X_i}} } \right) = \sum\limits_{j = 1}^p {\sum\limits_{i = 1}^p {{a_i}{a_j}{\mathop{\rm cov}} \left( {{X_i},{X_j}} \right)} } } \right.\\proof,eg:\\\left\{ \begin{array}{l}{\mathop{\rm var}} \left( {{a_1}{X_1} + {a_2}{X_2} + {a_3}{X_3}} \right) = {\mathop{\rm var}} \left( {\sum\limits_{i = 1}^3 {{a_i}{X_i}} } \right)\\ = {\mathop{\rm var}} \left( {{a_1}{X_1} + {a_2}{X_2}} \right) + 2{\mathop{\rm cov}} \left( {{a_1}{X_1} + {a_2}{X_2},{a_3}{X_3}} \right) + {\mathop{\rm var}} \left( {{a_3}{X_3}} \right)\\ = \left( {{\mathop{\rm var}} \left( {{a_1}{X_1}} \right) + 2{\mathop{\rm cov}} \left( {{a_1}{X_1},{a_2}{X_2}} \right) + {\mathop{\rm var}} \left( {{a_2}{X_2}} \right)} \right) + 2\left( {{\mathop{\rm cov}} \left( {{a_1}{X_1},{a_3}{X_3}} \right) + {\mathop{\rm cov}} \left( {{a_2}{X_2},{a_3}{X_3}} \right)} \right) + {\mathop{\rm var}} \left( {{a_3}{X_3}} \right)\\ = \left( {{a_1}^2{\mathop{\rm var}} \left( {{X_1}} \right) + 2{a_1}{a_2}{\mathop{\rm cov}} \left( {{X_1},{X_2}} \right) + {a_2}^2{\mathop{\rm var}} \left( {{X_2}} \right)} \right) + 2\left( {{a_1}{a_3}{\mathop{\rm cov}} \left( {{X_1},{X_3}} \right) + {a_2}{a_3}{\mathop{\rm cov}} \left( {{X_2},{X_3}} \right)} \right) + {a_3}^2{\mathop{\rm var}} \left( {{X_3}} \right)\\ = {a_1}^2{\mathop{\rm var}} \left( {{X_1}} \right) + 2{a_1}{a_2}{\mathop{\rm cov}} \left( {{X_1},{X_2}} \right) + {a_2}^2{\mathop{\rm var}} \left( {{X_2}} \right) + 2{a_1}{a_3}{\mathop{\rm cov}} \left( {{X_1},{X_3}} \right) + 2{a_2}{a_3}{\mathop{\rm cov}} \left( {{X_2},{X_3}} \right) + {a_3}^2{\mathop{\rm var}} \left( {{X_3}} \right)\\ = \sum\limits_{j = 1}^3 {\sum\limits_{i = 1}^3 {{a_i}{a_j}{\mathop{\rm Cov}\nolimits} \left( {{X_i},{X_j}} \right)} } \end{array} \right.\end{array}%
$$
