Negative eigenvalues in principle component analysis in the presence of missing data For the purpose of dimension reduction I have performed an eigen analysis (using Jacobi-iteration) on a correlation matrix R of 163 variables (based on 1500 cases). The scree plot is attached. 
The sum of all eigenvalues is (to 6 places after the decimal) equal to the number of variables, as expected. However, the last 54 eigenvalues are negative, and the cumulative explained variance goes to 125% before droping to 100%. What is going on here, and does it still make sense to use the first 15 or so eigenvectors to calculate scores and communalities?
There are missing data in the set, therefore a robust calculation of correlation coefficients was used:
set sums and n to 0
for i := 1 to Cases do
  if ((IsNaN(x[i]) or (IsNaN(y[i]))
    then do nothing
    else add $x, y, xy, x^2, y^2$ to their respective sums and inc(n)
$r = \frac{n * SumXY - SumX * SumY}{\sqrt{[n * SumX^2 - (SumX)^2] [n* SumY^2 - (SumY)^2]}}$
The idea is to use the available information completely without inventing anything by imputation (at this state, for the calculation of scores by multiplication of data and eigenvectors I see no way other than imputation of $\bar{x}$).
 A: If you calculate pairwise correlation coefficients in presence of missing values, you correlation matrix may end up being non positive definite. In fact it's a very common phenomenon in quant finance. One way to deal with this issue is Ledoit Wolf procedure, see here. They developed a method for a different issue, but it's used for missing value issue too. An author has MATLAB code here.
Suppose you have three variables x,y and z. In observation 1 value of x is missing, but y and z are present. One way to calculate the correlation matrix is to skip the observation 1. 
Another, seemingly better way is to skip observation 1 only when calculating pairwise correlations xy and xz, and use it for YZ correlation. Values of Y and Z are available in observation 1, why not use them? If you go this way then obtained correlation matrix may not be a good estimate of the true correlation matrix, surprisingly. Particularly, your matrix may not be positive definite. Again, this is a common situation in many finance applications such as portfolio optimization and PCA.
I would skip the observations with missing values IF the data size allows it. This is not always possible, e.g. sometimes we have hundreds variables and about as many observations. If we skip an observation when at least one variable value is missing, easily a half of the observations may be flagged. In this case, it's worth the trouble to do pairwise correlations using "all available" data, then shrink the matrix using Ledoit Wolf procedure. Otherwise, if it was just a couple of rows dropping out, then I wouldn't bother and skip them.
A: After the discussion here, I can provide at least a partial answer. 
Apparently, pairwise calculation of correlation coefficients, especially if the data matrix has missing data, leads to a correlation matrix that is only positive semi-definite, not positive definite. Eigenvalues that are too extreme. According to http://epublications.bond.edu.au/cgi/viewcontent.cgi?article=1099&context=ejsie it is possible to counteract this effect by "shrinking" the matrix. To this purpose, a weighted average of the correlation matrix $R$ with a form matrix $F$ is calculated: $\hat{r}_{ij} = (1-\omega) r_{ij} + \omega f_{ij}$ with $\omega$ a weight factor from [0..1] (in the literature, it is often called $\lambda$, but as this symbol is used for eigenvalues already I have re-christend it). According to https://cssanalytics.files.wordpress.com/2013/10/shrinkage-simpler-is-better.pdf it doesn't really matter which form matrix one uses, I have tried both the identity matrix $I$ and the average matrix $\bar{R}$ (for each variable $i$ calculate the average correlation with all other variables $\bar{r}_i$, then $\bar{r}_{ij} = (\bar{r}_i + \bar{r}_j)/2$). The results are indeed very similar.
What I still don't understand is why $R$ is so affected by missing data. I have done simulations with 10,000 data from $y = 1 + 2x + \epsilon$, with $\epsilon$ a Gaussian random number so that $r = 0.825$. Then I have set random elements of this matrix to NaN and recalculated $r$. For up to 20% missing data, the maximal deviation of $r$ was 0.01, and most of the data fell within ± 0.005. Even for 50% missing data the maximum deviation was 0.015, for 70% 0.02. Surely this is not such a big effect?
