Finite Population Variance for a Changing Population How does the addition of one unit affect the population variance of a finite population if everything else remains unchanged? What are the conditions such that the new unit leaves the variance unchanged (increases/decreases it)?
I was able to find the following paper regarding sample variances for changing finite populations:
http://www.amstat.org/sections/srms/Proceedings/papers/1987_087.pdf.
But I am asking specifically about population variances. Any help is appreciated.
 A: I was unable to find the sample calculations that correspond to the specific problem here (as suggested by Glen_b), but I was able to confirm the following answer with numerical calculations in R at the bottom of this answer.
Let $N$ be the initial number of units in the population and $N + 1$ be the number of units in the population after the change. Denote the initial set of observations $X = \{x_1, \ldots, x_N\}$ (i.e., one observation corresponding to each population unit). Denote the set of observations after the change $Y = X \cup \{x_{N+1}\}$.
The mean of $X$ is
$\mu_X = \frac{\sum_{i=1}^N{x_i}}{N}$.
The mean of Y is
$\mu_Y = \frac{\sum_{i=1}^{N+1}{x_i}}{N+1}
= \mu_X \frac{N}{N+1} + \frac{x_{N+1}}{N+1}$
Define $x_{N+1}$ as the original mean, $\mu_X$, plus some $\varepsilon$. Then,
the mean of $Y$ is
$\mu_Y = \mu_X \frac{N}{N+1} + \frac{\mu_X
+ \varepsilon}{N+1} = \mu_X + \frac{\varepsilon}{N+1}$
The variance of $Y$ is 
$\sigma^2_Y = \frac{\sum_{i=1}^{N+1}
\left(x_i - \mu_Y \right)^2}{N+1} =
\frac{\sum_{i=1}^{N+1}
\left(x_i - \mu_X - \frac{\varepsilon}{N + 1} \right)^2}{N+1}$
$= \frac{\sum_{i=1}^{N} x_i^2 + \mu_X^2
+ \frac{\varepsilon^2}{\left(N+1\right)^2} - 2x_i\mu_X
- 2x_i\frac{\varepsilon}{N+1} + 2\mu_X\frac{\varepsilon}{N+1}}{N + 1}$
$\frac{\left(\mu_X + \varepsilon - \mu_X - \frac{\varepsilon}{N + 1}\right)}{N
+ 1} $
$ = \frac{N}{N+1}\sigma^2_X + \frac{N\varepsilon^2}{\left(N+1\right)^3}
    - \frac{2N\mu_X\varepsilon}{\left(N+1\right)^2}
    + \frac{2N\mu_X\varepsilon}{\left(N+1\right)^2}
    + \frac{N^2\varepsilon^2}{\left(N+1\right)^3}$
$ = \frac{N}{N+1}
\sigma^2_X + \frac{N}{\left(N+1\right)^2}\varepsilon^2$
When $x_{N+1}$ is equal to $\mu_X$, the variance of $Y$ is 
$\frac{N}{N+1}\sigma^2_X < \sigma^2_X $
Thus, when $\varepsilon$ is sufficiently small $\sigma^2_Y$ is less than
$\sigma^2_X$. To determine how large $\varepsilon$ should be so that the
variance of $Y$ is greater than the variance of $X$, I set the two variances
equal.
$ \frac{N}{N+1} \sigma^2_X
    + \frac{N}{\left(N+1\right)^2}\varepsilon^2 = \sigma^2_X$
$   \frac{N}{\left(N+1\right)^2}\varepsilon^2 = \frac{1}{N+1} \sigma^2_X$
$   \varepsilon^2 = \frac{N+1}{N} \sigma^2_X $
$ \varepsilon = \pm \sigma_X
    \sqrt{\frac{N+1}{N}}$
Thus, adding a unit who's observation is within $\sqrt{\frac{N+1}{N}}$ standard deviations
of the old mean will lead to a lower variance.

The following R script verifies the above conclusion:
N <- 10
X <- runif(N)
width <- sqrt((N+1)/N)
# on the boundary
var(c(X, mean(X) + width * sqrt(var(X)))) - var(X) == 0
# outside the boundary
var(c(X, mean(X) + width * sqrt(var(X)) + 1)) - var(X) > 0
# inside the boundary
var(c(X, mean(X))) - var(X) < 0

A: I have the feeling that you may be confusing a finite population with a sample from it. The fact that a population is finite, does not make it "equivalent" to a sample (which is always finite of course).  
When we examine populations that are comprised of "identically and independently distributed" random variables, we have the habit of talking about the population mean or variance: strictly speaking this is wrong language -what we mean is the mean/variance of the common marginal distribution that each member of the population follows.  
If this is your case, then, the expression $\frac{\sum_{i=1}^N{x_i}}{N}$ represents the sample mean of a specific sample from this finite population, namely, of a specific set of realizations of the random variables comprising this finite population. It is not the mean (expected value) of the population, i.e. it is not the common expected value of the i.i.d. variables comprising the population.  
And all your calculations are consistent with examining the sample mean and the sample variance, not their population counterparts.  
Viewed in this light, your calculations are correct and intuitive: if the additional observation is exactly equal to the sample mean of the previous observations included in $X$, then dispersion lessens and the sample variance of $Y$ will be smaller.  
I guess it is evident that adding an i.i.d. random variable to an i.i.d. population does not change "the" population mean or variance (i.e. the moments of the common marginal distribution).
A: Suppose we use the population variance expression that incorporates Bessel's correction:
$$S_N^2 = \frac{1}{N-1} \sum_{i=1}^N (X_i - \bar{X}_N)^2.$$
Using a simple variance decomposition result in O'Neill (2014) (Result 1, p. 283) you have:
$$S_{N+1}^2 = \frac{N-1}{N} \cdot S_N^2 + \frac{1}{N+1} \cdot (\bar{X}_N - X_{N+1})^2.$$
This gives you a simple recursive formula that allows you to update your population variance as you add an additional element $X_{N+1}$.  It is simple to create a function to update the variance by a single new value.  We will code this in R as the function update.var, which requires specification of the mean, variance and size of the existing sample, and the new value to add to the sample.
#Create function to update variance with one new value
update.var <- function(n, mean, var, value) {

  #Check input values
  if (!is.vector(n))          { stop('Error: n should be a single positive integer') }
  if (length(n) != 1)         { stop('Error: n should be a single positive integer') }
  if (!is.numeric(n))         { stop('Error: n should be a positive integer') }
  if (as.integer(n) != n)     { stop('Error: n should be a positive integer') }
  if (n < 1)                  { stop('Error: n should be a positive integer') }
  if (!is.vector(mean))       { stop('Error: mean should be a numeric value') }
  if (length(mean) != 1)      { stop('Error: mean should be a numeric value') }
  if (!is.numeric(mean))      { stop('Error: mean should be a numeric value') }
  if (!is.vector(var))        { stop('Error: var should be a numeric value') }
  if (length(var) != 1)       { stop('Error: var should be a numeric value') }
  if (!is.numeric(var))       { stop('Error: var should be a non-negative value') }
  if (var < 0)                { stop('Error: var should be a non-negative value') }
  if ((n == 1)&&(var != Inf)) { stop('Error: If n = 1 then var must be Inf') }
  if (!is.vector(value))      { stop('Error: value should be a numeric value') }
  if (length(value) != 1)     { stop('Error: value should be a numeric value') }
  if (!is.numeric(value))     { stop('Error: value should be a numeric value') }

  newvar <- ifelse(n == 1, ((mean - value)^2)/(n+1),
                           ((n-1)/n)*var + ((mean - value)^2)/(n+1));
  newvar; }

We can check that this method works by comparing the iterative and non-iterative computation.  The code below shows that the values from both methods are the same (to within a very small tolerance that occurs due to rounding error in the two different methods).
#Set initial seed and max pop size
set.seed(1);
M  <- 20;

#Generate random populations and compute the updating mean and variance values
VALUES <- rnorm(M);
MEANS  <- rep(0, M);
for (n in 1:M) { MEANS[n] <- mean(VALUES[1:n]); }
VARS1  <- rep(Inf, M);
for (n in 2:M) { VARS1[n] <- var(VALUES[1:n]); }

#Generate the variance values iteratively
VARS2  <- rep(Inf, M);
for (n in 2:M) { VARS2[n] <- update.var(n-1, MEANS[n-1], VARS2[n-1], value = VALUES[n]); }

#Check that results are the same (within rounding tolerance)
max(abs(VARS1[2:M]-VARS2[2:M]));
[1] 2.220446e-16

