Asymptotic distribution of multinomial I'm looking for the limiting distribution of multinomial distribution over d outcomes. IE, the distribution of the following
$$\lim_{n\to \infty} n^{-\frac{1}{2}} \mathbf{X_n}$$
Where $\mathbf{X_n}$ is a vector value random variable with density $f_n(\mathbf{x})$ for $\mathbf{x}$ such that $\sum_i x_i=n$, $x_i\in \mathbb{Z}, x_i\ge 0$ and 0 for all other $\mathbf{x}$, where
$$f_{n}(\mathbf{x})=n!\prod_{i=1}^d\frac{p_i^{x_i}}{x_i!}$$
I found one form in Larry Wasserman's "All of Statistics" Theorem 14.6, page 237 but for limiting distribution it gives Normal with a singular covariance matrix, so I'm not sure how to normalize that. You could project the random vector into (d-1)-dimensional space to make covariance matrix full-rank, but what projection to use?
Update 11/5
Ray Koopman has a nice summary of the problem of singular Gaussian. Basically, singular covariance matrix represents perfect correlation between variables, which is not possible to represent with a Gaussian. However, one could get a Gaussian distribution for the conditional density, conditioned on the fact that the value of random vector is valid (components add up to $n$ in the case above). 
The difference for the conditional Gaussian, is that inverse is replaced with pseudo-inverse, and normalization factor uses "product of non-zero eigenvalues" instead of "product of all eigenvalues". Ian Frisce gives link with some details.
There's also a way to express normalization factor of conditional Gaussian without referring to eigenvalues, 
here's a derivation
 A: The covariance is still non-negative definite (so is a valid multivariate normal distribution), but not positive definite: what this means is that (at least) one element of the random vector is a linear combination of the others.
As a result, any draw from this distribution will always lie on a subspace of $R^d$. As a consequence, this means it is not possible to define a density function (as the distribution is concentrated on the subspace:  think of the way a univariate normal will concentrate at the mean if the variance is zero).
However, as suggested by Robby McKilliam, in this case you can drop the last element of the random vector. The covariance matrix of this reduced vector will be the original matrix, with the last column and row dropped, which will now be positive definite, and will have a density (this trick will work in other cases, but you have to be careful which element you drop, and you may need to drop more than one). 
A: It looks to me like Wasserman's covariance matrix is singular, to see, multiply it by a vector of $d$ ones, i.e. $[1,1,1,\dots,1]^\prime$ of length $d$.
Wikipedia gives the same covariance matrix anyway. If we restrict ourselves to just a binomial distribution then the standard central limit theorem tells us that the binomial distribution (after appropriate scaling) converges to the normal as $n$ gets big (see wikipedia again). Applying similar ideas you should be able to show that an appropriately scaled mulinomial is going to converge in distribution to the multivariate normal, i.e. each marginal distribution is just a binomial and converges to the normal distribution, and the variance between them is known.
So, I am very confident you will find that the distribution of
$$\frac{X_n - np}{\sqrt{n}}$$
converges to the multivariate normal with zero mean and covariance 
$$\frac{C}{n}$$
where $C$ is the covariance matrix of the multinomial in question and $p$ is the vector of probabilities $[p_1,\dots,p_d]$.
A: There no inherent problem with the singular covariance here.  Your asymptotic distribution is the singular normal.  See http://fedc.wiwi.hu-berlin.de/xplore/tutorials/mvahtmlnode34.html which gives the density of the singular normal.
A: Is it not the case that $|S_{-i}|=|S_{-j}|$ for all $i,j$ where $S_{-i}$ is the Multinomial covariance matrix with the $i$-th row and column removed?  Since this is the case, I don't understand what you mean by "freedom of choice" as any "choice" is equivalent.
A: Another way to think of this is your multinomial distribution parameters are embedded in the simplex, so your asymptotic distribution should also be embedded in the simplex. You can get to this by trying to instead characterize the distribution of $y$ given by the transform
$$
y = ilr(p)=V^T\log p
$$
where $ilr(p)$ is the isometric log-ratio transform and $V$ is an orthonormal basis of dimension $D-1\times D$.  It doesn't actually matter which basis you use, but you could use the top D-1 eigenvectors from the eigendecomposition of $\log p$ for starters.
It can then be shown that the asympotic distribution can be given by
$$
\hat{y} \sim \mathcal{N}_{D-1}\big(y, \frac{1}{n} V^TD^{-1}V\big)
$$
where $D=diag(p)$
See Ortego et al 2015 for more details, the paper can be found below
https://www.researchgate.net/publication/327646436_The_asymptotic_distribution_of_log-ratio_transformed_proportions_of_multinomial_count_data
