Bounds on Cov(X, Y) given Var(X), Var(Y)? I'm generating random multivariate normal data using the rmultnorm() function in R, which allows users to specify a vector of $k$ population means and a $k \times k$ variance-covariance matrix. 
Given $\newcommand{\Var}{\mathrm{Var}}\Var(X)$, and $\Var(Y)$, what are the boundaries for the values of $\newcommand{\Cov}{\mathrm{Cov}}\Cov(X, Y)$ I can select for the variance-covariance matrix? 
Does it make sense to use the property: 
\begin{align}\Var(X) + \Var(Y) = \Var(X+Y) - 2\Cov(X,Y)\end{align} 
and the Pearson correlation coefficient equation:
\begin{align}-1 \leq \frac{\Cov(X,Y)}{\sqrt{\Var(X)\Var(Y)}} \leq 1 \end{align} 
to set upper/lower bounds on $\Cov(X,Y)$?
 A: For a single covariance you only need the bottom equation.  The bounds are that the covariance cannot be greater than the product of the standard deviations (and cannot be less than the negative of the same value).  However for a covariance matrix of more than 2 terms there is an additional limit, the matrix has to be positive semi-definite (or positive definite in some cases).  This eliminates cases like X and Y are strongly positively correlated and X and Z are strongly possitively correlated, but Y and Z are strongly negatively correlated (this does not work in the real world, but could still fit the pairwise bounds).  One way to check this is that if all the eigen values are positive then it is postive definite and if the eigen values are all non-negative then it is positive semi-definite.
A: In the multivariate case, you can use what is called the multivariate Cauchy-Schwarz inequality:
$$
   \newcommand{\Var}{Var} \newcommand{\Cov}{Cov}
      \Var(z) \ge \Cov(z,y) \Var(y)^{-1} \Cov(y,z)
$$
Where here the inequality sign must be interpreted in the sense of the partial order on the cone of positive-definite matrices: $A \le B$ means that $B-A$ is positive definite.
