How to create an arbitrary covariance matrix For example, in R, the MASS::mvrnorm() function is useful for generating data to demonstrate various things in statistics.  It takes a mandatory Sigma argument which is a symmetric matrix specifying the covariance matrix of the variables. How would I create a symmetric $n\times n$ matrix with arbitrary entries? 
 A: I like to have control over the objects I create, even when they might be arbitrary.
Consider, then, that all possible $n\times n$ covariance matrices $\Sigma$ can be expressed in the form
$$\Sigma= P^\prime\ \text{Diagonal}(\sigma_1,\sigma_2,\ldots, \sigma_n)\ P$$
where $P$ is an orthogonal matrix and $\sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_n \ge 0$.
Geometrically this describes a covariance structure with a range of principal components of sizes $\sigma_i$.  These components point in the directions of the rows of $P$.  See the figures at Making sense of principal component analysis, eigenvectors & eigenvalues for examples with $n=3$.  Setting the $\sigma_i$ will set the magnitudes of the covariances and their relative sizes, thereby determining any desired ellipsoidal shape.  The rows of $P$ orient the axes of the shape as you prefer.
One algebraic and computing benefit of this approach is that when $\sigma_n \gt 0$, $\Sigma$ is readily inverted (which is a common operation on covariance matrices):
$$\Sigma^{-1} = P^\prime\ \text{Diagonal}(1/\sigma_1, 1/\sigma_2, \ldots, 1/\sigma_n)\ P.$$
Don't care about the directions, but only about the ranges of sizes of the the $\sigma_i$?  That's fine: you can easily generate a random orthogonal matrix.  Just wrap $n^2$ iid standard Normal values into a square matrix and then orthogonalize it.  It will almost surely work (provided $n$ isn't huge). The QR decomposition will do that, as in this code
n <- 5
p <- qr.Q(qr(matrix(rnorm(n^2), n)))

This works because the $n$-variate multinormal distribution so generated is "elliptical": it is invariant under all rotations and reflections (through the origin).  Thus, all orthogonal matrices are generated uniformly, as argued at How to generate uniformly distributed points on the surface of the 3-d unit sphere?.
A quick way to obtain $\Sigma$ from $P$ and the $\sigma_i$, once you have specified or created them, uses crossprod and exploits R's re-use of arrays in arithmetic operations, as in this example with $\sigma=(\sigma_1, \ldots, \sigma_5) = (5,4,3,2,1)$:
Sigma <- crossprod(p, p*(5:1))

As a check, the Singular Value decomposition should return both $\sigma$ and $P^\prime$.  You may inspect it with the command
svd(Sigma)

The inverse of Sigma of course is obtained merely by changing the multiplication by $\sigma$ into a division:
Tau <- crossprod(p, p/(5:1))

You may verify this by viewing zapsmall(Sigma %*% Tau), which should be the $n\times n$ identity matrix.  A generalized inverse (essential for regression calculations) is obtained by replacing any $\sigma_i \ne 0$ by $1/\sigma_i$, exactly as above, but keeping any zeros among the $\sigma_i$ as they were.
A: Create an $n\times n$  matrix $A$ with arbitrary values 
and then use $\Sigma = A^T A$ as your covariance matrix.  
For example 
n <- 4  
A <- matrix(runif(n^2)*2-1, ncol=n) 
Sigma <- t(A) %*% A

A: You can simulate random positive definite matrices from the Wishart distribution using the function "rWishart" from stats (included in base)
n <- 4
rWishart(1,n,diag(n))

From the documentation for rWishart:
Usage should be:
rWishart(n, df, Sigma)

where,

*

*n: the number of samples.

*df: the degrees of freedom, i.e. the number of dimensions of the matrix.

*Sigma: a positive definite scaling matrix.

