When sampling from a finite population without replacement, the observations are negatively correlated with each other, and the sample variance $s^2 = \frac{1}{n-1} \sum_i \left( x_i - \bar{x} \right)^2$ is a slightly biased estimate of the population variance $\sigma^2$.
The derivation in this link from Robert Serfling provides a clear explanation of what's going on. The author first proves that if the observations in a sample have constant covariance (i.e. $\mathrm{Cov}\left(x_i, x_j \right) = \gamma$ for all $i\neq j$) that:
$$ E[s^2] = \sigma^2 - \gamma$$
For independent draws (hence $\gamma = 0$), you have $E[s^2] = \sigma^2$ and the sample variance is an unbiased estimate of the population variance. But the issue you have with sampling without replacement from a finite population is that your draws are negatively correlated with each other!
In the case of sampling without replacement from a population of size $N$:
$$ \text{For $i\neq j$ }\quad \mathrm{Cov}\left(x_i, x_j \right) = \frac{-\sigma^2}{N-1}$$
Hence:
$$ E\left[s^2\right] = \frac{N}{N-1}\sigma^2 $$
Yes, that statement is correct.
A convenient and rigorous way to analyze this uses the indicator function of the sample. That is, define the random variable $I_i$ to be $1$ when $x_i$ is in the sample and $0$ otherwise. Clearly the sampling procedure defines the $I_i$ and conversely.
Because the sample size is fixed at $n,$
$$\sum_{i=1}^N I_i = n$$
is constant. We can deduce everything we need from this simple fact.
We can easily compute the moments of these random variables, because they are interchangeable: in particular, the distributions of any $I_i$ and $I_j$ must be the same, $I_i^2 = I_i,$ and the joint distributions of any of the ordered pairs $(I_i,I_j),$ for $i\ne j,$ must all be the same. This justifies writing
$$E[I_i^2] = E[I_i] = p\text{ and } E[I_iI_j] = q$$
for some fixed numbers $p$ and $q$ (assuming $N\ge 2;$ the case $N=1$ is not in the scope of the question). Moreover, the sample size determines both these expectations, since
$$n = E[n] = E\left[\sum_{i=1}^N I_i\right] = \sum_{i=1}^N E[I_i] = Np$$
and
$$n^2 = E[n^2] = E\left[\left(\sum_{i=1}^N I_i\right)^2\right] = \sum_{i=1}^N E[I_i^2] + N(N-1)\sum_{i\ne j} E[I_iI_j] = Np + N(N-1)q$$
imply
$$p = \frac{n}{N}\quad\text{and}\quad q =\frac{n}{N}\frac{n-1}{N-1}.$$
The rest is algebra--essentially, re-arranging sums and products. Here are the details.
The sample variance is a random variable determined by the sample. It can be expressed in terms of the $I_i$ and population values $x_i$ as
$$\begin{aligned}
&\sum_{i=1}^N (x_i I_i )^2 - \frac{1}{n}\left(\sum_{i=1}^N x_iI_i\right)^2 \\
&= \sum_{i=1}^N x_i^2 I_i - \frac{1}{n}\sum_{i=1}^N x_i^2 I_i - \frac{1}{n}\sum_{i\ne j}^N x_i x_jI_iI_j .
\end{aligned}$$
all divided by $n-1.$
The expectation of this therefore equals
$$\begin{aligned}
\sum_{i=1}^N x_i^2 E[I_i] - \frac{1}{n}\sum_{i=1}^N x_i^2 E[I_i] - \frac{1}{n}\sum_{i\ne j}^N x_i x_jE[I_iI_j]\\
=\sum_{i=1}^N x_i^2 \frac{n}{N} - \frac{1}{n}\sum_{i=1}^N x_i^2 \frac{n}{N} - \frac{1}{N}\sum_{i\ne j}^N x_i x_j\frac{n}{N}\frac{n-1}{N-1}\\
=\frac{n-1}{N}\left[\sum_{i=1}^N x_i^2 - \frac{1}{N-1}\sum_{i\ne j}^N x_i x_j\right]\\
=\frac{n-1}{N-1}\left[\sum_{i=1}^N x_i^2 - \frac{1}{N}\left(\sum_{i=1}^N x_i\right)^2\right] \\
=(n-1) \frac{1}{N-1}\sum_{i=1}^N (x_i - \mu)^2.
\end{aligned}$$
where, as in the question, $N\mu$ is the sum of all the $x_i.$
Upon dividing by $n-1,$ we obtain the stated result.
Best Answer
First, if you're trying to estimate the population variance, σ², then the population variance is a perfect estimator of the population variance. I suspect that is not what you meant. Let's call Σ(x - x̄)²/n the maximum-likelihood estimator of σ².
Yes, when sampling without replacement from a finite population, s² tends to overestimate the true variance, σ², and using the maximum-likelihood estimator gives you a lower estimate than s². But doing so might over-correct or under-correct. There would still be a bias.
If you multiply s² by the factor (N-1)/N, where N is the size of the finite population, you will get an unbiased estimator of σ².