I understand that when sampling from a finite population and our sample size is more than 5% of the population, we need to make a correction on the sample's mean and standard error using this formula:
$\hspace{10mm} FPC=\sqrt{\frac{N-n}{N-1}}$
Where $N$ is the population size and $n$ is the sample size.
I have 3 questions about this formula:
The threshold is chosen such that it ensures convergence of the hypergeometric distribution ($\sqrt{\frac{N-n}{N-1}}$ is its SD), instead of a binomial distribution (for sampling with replacement), to a normal distribution (this is the Central Limit Theorem, see e.g., The Normal Curve, the Central Limit Theorem, and Markov's and Chebychev's Inequalities for Random Variables). In other words, when $n/N\leq 0.05$ (i.e., $n$ is not 'too large' compared to $N$), the FPC can safely be ignored; it is easy to see how the correction factor evolves with varying $n$ for a fixed $N$: with $N=10,000$, we have $\text{FPC}=.9995$ when $n=10$ while $\text{FPC}=.3162$ when $n=9,000$. When $N\to\infty$, the FPC approaches 1 and we are close to the situation of sampling with replacement (i.e., like with an infinite population).
To understand this results, a good starting point is to read some online tutorials on sampling theory where sampling is done without replacement (simple random sampling). This online tutorial on Nonparametric statistics has an illustration on computing the expectation and variance for a total.
You will notice that some authors use $N$ instead of $N-1$ in the denominator of the FPC; in fact, it depends on whether you work with the sample or population statistic: for the variance, it will be $N$ instead of $N-1$ if you are interested in $S^2$ rather than $\sigma^2$.
As for online references, I can suggest you
As for the derivation, the setup goes as follows. Suppose we have a population of size $N$, with mean $\mu$ and variance $\sigma^2$, where each element can assume values $v_k$ for $k = 1, 2, \dots, m$. Let $n_k$ be the number of times that the value $v_k$ occurs in the population, such that the probability that we draw the value $v_k$ at random from the population is: $$P(X=v_k) = \frac{n_k}{N}$$
We pick a sample of size $n$, without replacement, and we estimate the global mean $\mu$ with the estimator $\overline{X}=\frac{\sum_{i=1}^n X_i}{n}$.
We are going to find the formula of the finite population correction factor by looking at the variance of the estimator:
$$ \tag{1} \label{variance} \mathrm{Var}(\overline{X}) = \mathrm{Var} \left( \frac{\sum_{i=1}^n X_i}{n} \right) = \frac{1}{n^2} \sum_{i=1}^n\sum_{j=1}^n \mathrm{Cov}(X_i, X_j) $$
Notice that if we were doing sampling with replacement, the variables $X_i$ would be completely independent of each other, meaning that there would be no covariance between them: $$ \mathrm{Cov}(X_i, X_j) = 0, \quad i \ne j $$ This would imply that we could discard all terms where $i \ne j$. Also, when $i$ and $j$ are equal the covariance is just the variance: $$\mathrm{Cov}(X_i, X_i) = \mathrm{Var}(X_i) = \sigma^2$$ Which would mean we could work the variance like so: \begin{align*} \frac{1}{n^2} \sum_{i=1}^n\sum_{j=1}^n \mathrm{Cov}(X_i, X_j) &= \frac{1}{n^2} \sum_{i=1}^n \mathrm{Var}(X_i) \\ &= \frac{1}{n^2} \sum_{i=1}^n \sigma^2 \\ &= \frac{n \sigma^2}{n^2} \\ \mathrm{Var}(\overline{X}) &= \frac{\sigma^2}{n} \end{align*}
This is the variance for sampling with replacement (or with an infinite population).
However, since we are doing sampling without replacement, the random variables $X_i$ aren't independent (considering we can't get a given element more than once, the probability that we get a certain value for a given $X_i$ depends on the values of the remaining ones). We treat the summation above by splitting the indices where $i=j$ and where $i\ne j$, similarly like we did the covariance for sampling with replacement: \begin{align*} \frac{1}{n^2} \sum_{i=1}^n\sum_{j=1}^n \mathrm{Cov}(X_i, X_j) &= \frac{1}{n^2} \left( \sum_{i=1}^n\sum_{j=i} \mathrm{Cov}(X_i, X_j) + \sum_{i=1}^n\sum_{j \ne i} \mathrm{Cov}(X_i, X_j) \right) \\ &= \frac{1}{n^2} \left( \sum_{i=1}^n \mathrm{Var}(X_i) + \sum_{i=1}^n\sum_{j \ne i} \mathrm{Cov}(X_i, X_j) \right) \\ &= \frac{1}{n^2} \left( \sum_{i=1}^n \sigma^2 + \sum_{i=1}^n\sum_{j \ne i} \mathrm{Cov}(X_i, X_j) \right) \\ &= \frac{1}{n^2} \left( n \sigma^2 + \sum_{i=1}^n\sum_{j \ne i} \mathrm{Cov}(X_i, X_j) \right) \\ &= \frac{\sigma^2}{n} + \frac{1}{n^2} \sum_{i=1}^n\sum_{j \ne i} \mathrm{Cov}(X_i, X_j) \tag{2} \label{covariance} \end{align*}
Straightaway we must find the covariance between $X_i$ and $X_j$ when $i \ne j$. Recall the definition of covariance: $$ \mathrm{Cov} (X_i, X_j) = \mathrm{E}[X_i X_j] - \mathrm{E}[X_i]\mathrm{E}[X_j] $$ Since $\mathrm{E}[X_i] = \mathrm{E}[X_j] = \mu$, this yields: $$\tag{3} \label{covariance expectation} \mathrm{Cov} (X_i, X_j) = \mathrm{E}[X_i X_j] - \mu^2 $$ Immediately we proceed to calculate $\mathrm{E}[X_i X_j]$, which is defined as:
$$ \tag{4} \label{covariance expectation summation} \mathrm{E}[X_i X_j] = \sum_{k=1}^m \sum_{l=1}^m v_k v_l \ P(X_i=v_k , X_j=v_l) $$
The tricky part is calculating the probability above, because this probability changes depending on whether $k=l$ or not. The whole thing becomes clearer using Bayes' theorem: \begin{align*} P(X_i=v_k, X_j=v_l) &= P(X_i = v_k)P(X_j=v_l | X_i=v_k)\end{align*}
Consider the case where $k=l$: this is equivalent to drawing the same value $v_k$ twice. The probability of drawing $v_k$ is $P(X_i=v_k)=\frac{n_k}{N}$, and doing so again (given that we already drew $v_k$) is: $$P(X_j=v_k | X_i=v_k) = \frac{n_k-1}{N-1}$$
On the other hand, considering the case where $k \ne l$, if we draw $v_l$ given that we already drew $v_k$, we find that the number of occurences of $v_l$ in the population is unchanged ($n_l$). However, the total size of our population is now $N-1$. Hence: $$P(X_j=v_l | X_i=v_k) = \frac{n_l}{N-1}, \quad k \ne l$$ Therefore, our probability is: $$P(X_i=v_k, X_j=v_l) = \begin{cases} \dfrac{n_k (n_k - 1)}{N(N - 1)}, & \quad k=l\\ \dfrac{n_k n_l}{N(N-1)}, & \quad k \ne l \end{cases} $$
Because of this, we must split the summation at $\eqref{covariance expectation summation}$ on the indices where $k=l$ and $k \ne l$, as such: \begin{align*} \mathrm{E}[X_i X_j] &= \sum_{k=1}^m \sum_{l=1}^m v_k v_l \ P(X_i=v_k, X_j=v_l) \\ &= \sum_{k=1}^m \sum_{k=l} v_k^2 P(X_i=v_k, X_j=v_k) + \sum_{k=1}^m \sum_{k \ne l} v_k v_l P(X_i=v_k, X_j=v_l) \\ &= \sum_{k=1}^m v_k^2 \frac{n_k (n_k - 1)}{N(N-1)} + \sum_{k=1}^m \sum_{k \ne l} v_k v_l \frac{n_k n_l}{N(N-1)} \end{align*}
Now we can pull the $N(N-1)$ factor out and do some manipulation on these sums: \begin{align*} \mathrm{E}[X_iX_j] &= \frac{1}{N(N-1)} \left( \sum_{k=1}^m v_k^2 n_k(n_k -1) + \sum_{k=1}^m \sum_{k \ne l} v_k n_k v_l n_l \right) \\ &= \frac{1}{N(N-1)} \left( \sum_{k=1}^m v_k^2 n_k^2 - \sum_{k=1}^m v_k^2 n_k + \sum_{k=1}^m \sum_{k \ne l} v_k n_k v_l n_l \right) \\ &= \frac{1}{N(N-1)} \left( \sum_{k=1}^m v_k^2 n_k^2 + \sum_{k=1}^m \sum_{k \ne l} v_k n_k v_l n_l - \sum_{k=1}^m v_k^2 n_k \right) \tag{5} \label{expanded summation} \end{align*} We must realize that there is a way to simplify this expression, by recalling that: \begin{align*} \left( \sum_i a_i \right)^2 &= \sum_i \sum_j a_i a_j \\ &= \sum_i a_i^2 + \sum_i \sum_{j \ne i} a_i a_j \end{align*}
That is, if we square a sum, we can write the result splitting its indices. This means that the apparently intractable sum above is just: $$ \sum_{k=1}^m v_k^2 n_k^2 + \sum_{k=1}^m \sum_{k \ne l} v_k n_k v_l n_l = \left( \sum_{k=1}^m v_k n_k \right)^2 $$ So we simplify $\eqref{expanded summation}$ to: $$ \tag{6} \label{simplified expectation} \mathrm{E}[X_iX_j] = \frac{1}{N(N-1)} \left( \left( \sum_{k=1}^m v_k n_k \right)^2 - \sum_{k=1}^m v_k^2 n_k \right) $$
We're almost done. Our task now is to represent the sums above in terms of known constants. Let us remember that, in the case where we have repeating values in our domain, the expected value $E[\cdot]$ can be written as: $$ \mathrm{E}[X_i] = \frac{1}{N} \sum_{k=1}^m v_k n_k = \mu $$ The $n_k$ term accounts for the fact that we have more than one ocurrence of the value $v_k$. From this, it follows that: \begin{gather*} \sum_{k=1}^m v_k n_k = N \mu \\ \left( \sum_{k=1}^m v_k n_k \right)^2 = N^2 \mu^2 \tag{6.1} \label{square of sum} \end{gather*} Likewise, the expected value of the square of the variable can be written as: $$ \mathrm{E}[X_i^2] = \frac{1}{N} \sum_{k=1}^m v_k^2 n_k $$ From the definition of variance this simplifies to another expression: \begin{align*} \mathrm{Var}(X_i) &= \mathrm{E}[X_i^2] - \mathrm{E}^2[X_i] \\ \sigma^2 &= \mathrm{E}[X_i^2] - \mu^2 \\ \mathrm{E}[X_i^2] &= \mu^2 + \sigma^2 \end{align*} And it follows immediately that: \begin{gather*} \sum_{k=1}^m v_k^2 n_k = N \ \mathrm{E}[X_i^2] \\ \sum_{k=1}^m v_k^2 n_k = N(\mu^2 + \sigma^2) \tag{6.2} \label{sum of squares} \end{gather*}
Substituting $\eqref{square of sum}$ and $\eqref{sum of squares}$ back into $\eqref{simplified expectation}$, we get: \begin{align*} \mathrm{E}[X_iX_j] &= \frac{1}{N(N-1)} \left( \left( \sum_{k=1}^m v_k n_k \right)^2 - \sum_{k=1}^m v_k^2 n_k \right) \\ &= \frac{1}{N(N-1)} \left( N^2\mu^2 - N(\mu^2 + \sigma^2) \right) \\ &= \frac{N^2\mu^2 - N\mu^2 - N\sigma^2}{N(N-1)} \\ &= \frac{\mu^2N(N-1) - N\sigma^2}{N(N-1)} \\ &= \mu^2 - \frac{\sigma^2}{N-1} \end{align*}
We substitute back yet again into $\eqref{covariance expectation}$ to find our covariance: \begin{align*} \mathrm{Cov}(X_i, X_j) &= \mathrm{E}[X_i X_j] - \mu^2 \\ &= \left( \mu^2 - \frac{\sigma^2}{N-1} \right) - \mu^2 \\ &= - \frac{\sigma^2}{N-1} \end{align*}
At last: $$ \tag{7} \label{covariance for i not j} \boxed{\mathrm{Cov}(X_i, X_j) = - \dfrac{\sigma^2}{N-1}}$$ Lastly, we substitute $\eqref{covariance for i not j}$ into $\eqref{covariance}$ to find the variance of the estimator $\overline{X}$: \begin{align*} \mathrm{Var}(\overline{X}) &= \frac{1}{n^2} \left( n\sigma^2 + \sum_{i=1}^n \sum_{j \ne i} \mathrm{Cov}(X_i, X_j) \right) \\ &= \frac{1}{n^2} \left( n\sigma^2 - \sum_{i=1}^n \sum_{j \ne i} \frac{\sigma^2}{N-1} \right) \\ &= \frac{1}{n^2} \left( n\sigma^2 - \frac{n(n-1)\sigma^2}{N-1} \right) \\ &= \frac{\sigma^2}{n} - \frac{(n-1)\sigma^2}{(N-1)n} \end{align*}
Wraping everything up, if we pull out the common $\frac{\sigma^2}{n}$ term, we find our desired correction factor for the variance: $$ \boxed{\mathrm{Var}(\overline{X}) = \frac{\sigma^2}{n} \left( 1 - \frac{n-1}{N-1} \right)}$$
$$ \boxed{ \mathrm{FCF} = 1 - \frac{n-1}{N-1} } $$
So there it is. Also, $1 - \frac{n-1}{N-1} = \frac{N-n}{N-1}$, just in case anyone missed it.
Here is an alternative setup within the framework of the superpopulation model of sampling theory. It differs in notation and conception to classical sampling theory, but I think it is quite simple and intuitive.
Let $X_1,X_2,X_3,...$ be an exchangeable "superpopulation" of values. Take the first $N$ values to be the finite population of interest and the first $n \leqslant N$ values as a sample from this population. (The exchangeability of the superpopulation means that the sample is a simple-random-sample from the population.) Now, consider the mean-difference $\bar{X}_n - \bar{X}_N$ measuring the difference between the sample mean and population mean. This quantity can be written in the form:
$$\begin{align} \bar{X}_n - \bar{X}_N &= \frac{1}{n} \sum_{i=1}^n X_i - \frac{1}{N} \sum_{i=1}^N X_i \\[6pt] &= \Big( \frac{1}{n} - \frac{1}{N} \Big) \sum_{i=1}^n X_i - \frac{1}{N} \sum_{i=n+1}^N X_i \\[6pt] &= \frac{N-n}{nN} \sum_{i=1}^n X_i - \frac{1}{N} \sum_{i=n+1}^N X_i \\[6pt] &= \frac{1}{n} \Bigg[ \frac{N-n}{N} \sum_{i=1}^n X_i - \frac{n}{N} \sum_{i=n+1}^N X_i \Bigg]. \\[6pt] \end{align}$$
We clearly have $\mathbb{E}(\bar{X}_n - \bar{X}_N) = 0$, so we can use the sample mean as an unbiased estimator for the population mean. If we denote the variance of the superpopulation by $\sigma^2$ then our quantity has variance:
$$\begin{align} \mathbb{V}(\bar{X}_n - \bar{X}_N) &= \frac{1}{n^2} \Bigg[ \Big( \frac{N-n}{N} \Big)^2 \sum_{i=1}^n \mathbb{V}(X_i) + \Big( \frac{n}{N} \Big)^2 \sum_{i=n+1}^N \mathbb{V}(X_i) \Bigg] \\[6pt] &= \frac{1}{n^2} \Bigg[ \Big( \frac{N-n}{N} \Big)^2 n \sigma^2 + \Big( \frac{n}{N} \Big)^2 (N-n) \sigma^2 \Bigg] \\[6pt] &= \frac{1}{n^2 N^2} \Bigg[ (N-n)^2 n \sigma^2 + n^2 (N-n) \sigma^2 \Bigg] \\[6pt] &= \frac{1}{n^2 N^2} \cdot (N-n) N n \sigma^2 \\[6pt] &= \frac{N-n}{N} \cdot \frac{\sigma^2}{n}. \\[6pt] \end{align}$$
Suppose we let $S_N^2$ and $S_{N*}^2$ denote the variance values for the population, where the first uses Bessel's correction and the second does not (so we have $S_N^2 = \frac{N}{N-1} S_{N*}^2$). In classical sampling theory the latter quantity is considered to be "the variance" of the population. (Formally it is the variance of the empirical distribution of the population.) However, the first of these quantities is an unbiased estimator for the superpopulation variance, so we can estimate the variance of our mean-difference quantity by:
$$\begin{align} \hat{\mathbb{V}}(\bar{X}_n - \bar{X}_N) &= \frac{N-n}{N} \cdot \frac{S_N^2}{n}. \\[6pt] \end{align}$$
Consequently, using the central limit theorem we can establish the following confidence interval for the population mean $\bar{X}_N$:
$$\text{CI}_N(1-\alpha) = \Bigg[ \bar{X}_n \pm \sqrt{\frac{N-n}{N}} \cdot \frac{t_{n-1, \alpha/2}}{\sqrt{n}} \cdot S_{N} \Bigg].$$
This is the form of the confidence interval that I find the most natural. However, with this form, you will notice that we use a finite population correction term that is different to your expression. The expression you are using occurs when we convert to the variance estimator $S_{N*}^2$ that does not use Bessel's correction (purportedly "the variance" of the population). In this case we have the equivalent expression:
$$\text{CI}_N(1-\alpha) = \Bigg[ \bar{X}_n \pm \sqrt{\frac{N-n}{N-1}} \cdot \frac{t_{n-1, \alpha/2}}{\sqrt{n}} \cdot S_{N*} \Bigg].$$
As you can see, framed in this latter form, the finite population correction term is the one in your question. You can see that the finite population correction term appears in the confidence interval formula in order to "correct" for the finite population. Taking $N \rightarrow \infty$ (so that the population of interest is the superpopulation) we get $FPC \rightarrow 1$, yielding the standard confidence interval for the mean parameter of a "large" population.
Now, as to the "5% rule", that is an arbitrary rule, and I don't recommend it. In my view it is best to always include the FPC when you have a finite population. If the sample proportion is small then the FPC is close to one, so it does not change the interval much, but I find it silly to remove it. Practitioners who offer these rules-of-thumb evidently think that with an FPC close to one they should remove the term, but I see no sense in that; it is an approximation for approximation's sake.
I hope this alternative presentation of the matter elucidates the correction term within the broader framework of the superpopulation model. I have always preferred this model of sampling theory, since it makes it simpler to distinguish between the finite population case and the infinite population case. As you can see, within this framework the correction term pops out fairly simply in the course of attempting to estimate the mean of the finite population.
This deserves a simple derivation. The idea is to separate analysis of the population from analysis of the purely combinatorial process of sampling without replacement.
Let the population values be $(x_1,x_2,\ldots, x_N)$ with mean $\mu$ and variance $\sigma^2.$ We are interested in the standard error of the sample mean.
Let $I_j$ be the indicator of whether $x_j$ is in a random sample of size $K.$ These random variables are identically distributed (but not independent!) Bernoulli random variables and since they sum to $K,$ $\Pr(I_j=1) = K/N$ whence for $j=1,2,\ldots, N,$
$$\operatorname{Var}(I_j) = \frac{K}{N}\left(1 - \frac{K}{N}\right) = \frac{K(N-K)}{N^2}.$$
Moreover, since the $I_j$ sum to the constant random variable equal to $K,$
$$\begin{aligned} 0 &= \operatorname{Var}(K) = \operatorname{Var}\left(\sum_{j=1}^N I_j\right) = \sum_{j=1}^N \operatorname{Var}(I_j) + \sum_{i\ne j}^N \operatorname{Cov}(I_i,I_j). \end{aligned}$$
From the identical distributions of the $I_j$ all $N$ variances are equal to each other and all $N(N-1)$ covariances are equal to each other. Solving this equation yields
$$\operatorname{Cov}(I_i,I_j) = -\frac{1}{N-1}\operatorname{Var}(I_i) = -\frac{K(N-K)}{N^2(N-1)}.$$
That's the crux of the matter: sampling without replacement induces computable negative correlations among the indicators.
The rest is straightforward algebra. The sample mean is
$$\bar X = \frac{1}{K} \sum_{j=1}^N x_jI_j,$$
whence it has expectation $E[\bar X]=\mu$ (easily) and variance
$$\begin{aligned} \operatorname{Var}(\bar X) &= \frac{1}{K^2}\sum_{i=1}^N\sum_{j=1}^N x_ix_j\operatorname{Cov}(I_i,I_j) \\ &= \frac{1}{K^2}\sum_{j=1}^N x_j^2 \frac{K(N-K)}{N^2} - \frac{1}{K^2}\sum_{i\ne j}^N x_ix_j \frac{K(N-K)}{N^2(N-1)}\\ &= \frac{N-K}{K(N-1)}\,\frac{1}{N}\left(\sum_{j=1}^N x_j^2 - \frac{1}{N}\left(\sum_{j=1}^N x_i\right)^2\right)\\ &= \frac{N-K}{K(N-1)}\sigma^2. \end{aligned}$$
That's all there is to it. (Please notice that $\sigma^2$ is the actual population variance, not an estimator. It differs by using a fraction $1/N$ rather than $1/(N-1)$ in its definition.)
The usual "infinite" population formula for the standard error of $\bar X$ is $\sqrt{\sigma^2/K}.$ The foregoing result multiplies this by $\sqrt{(N-K)/(N-1)},$ the so-called "finite population correction factor" for random samples without replacement. A first-order approximation is
$$\sqrt{\frac{N-K}{N-1}} \ \approx\ 1 - \frac{K-1}{2(N-1)}\ \approx \ 1 - \frac{K}{2N}.$$
Thus, the relative sample size $K/N$ governs the accuracy of the approximation. The 5% threshold merely assumes you want to estimate the standard error of the mean to within 5% accuracy.
