Why do correlation functions diverge in field theory?

Note that $\langle\cdots\rangle$ is just the statistical mechanics notation for the expectation denoted by $\mathbb{E}\cdots$ in probability textbooks. So, in principle, $$ \langle\phi(x)\phi(y)\rangle=\mathbb{E}\phi(x)\phi(y) $$ which for centered random variables (i.e., assuming $\mathbb{E}\phi(x)=\mathbb{E}\phi(y)=0$) should be equal to the covariance ${\rm Cov}(\phi(x),\phi(y))$ rather than what statisticians would call the coefficient of correlation, i.e., ${\rm Cov}(\phi(x),\phi(y))/\sqrt{(\mathbb{E}\phi^2(x)\times \mathbb{E}\phi^2(y))}$.

Now the main problem related to your question is that $\phi(x)$ is not a random variable. This is because the field $(\phi(x))_{x\in\mathbb{R}^d}$ is not a function but a Schwartz distribution which in general cannot be evaluated at a point. What can be evaluated are local averages $$ \phi(f)=``\int_{\mathbb{R}^d} \phi(x)f(x)\ d^dx" $$ for smooth test functions $f$, e.g., $C^{\infty}$ with compact support or in Schwartz space $S(\mathbb{R}^d)$. Take for instance the Euclidean 2d Ising CFT. Then $\phi$ is a random distribution, i.e., a random variable $\omega\mapsto \phi[\omega]$ with values in $S'(\mathbb{R}^2)$ on some probability space $(\Omega,\mathcal{F},\mathbb{P})$. The maps $\Omega\rightarrow \mathbb{R}$, $\omega\mapsto \phi[\omega](f)$ indexed by fixed test functions $f$, are honest real-valued random variables with moments of all orders. In particular, the variance $$ {\rm Var}(\phi(f))=\mathbb{E} \phi(f)^2={\rm Cst}\times \int_{\mathbb{R}^2}\int_{\mathbb{R}^2} \frac{f(x)f(y)}{|x-y|^{\frac{1}{4}}}\ d^2x\ d^2y $$ is finite.

There is a Kolmogorov-Chentsov Theorem for distributions with negative Hölder regularity which generalizes the usual one for positive regularity as in standard textbooks on Brownian motion. See for instance this article by Furlan and Mourrat. The blow-up on the diagonal is a signature of the fact the field $\phi$ is a distribution rather than a function.

It's because $\langle \phi(x)\, \phi(y)\rangle$ is not perfectly analogous to $\operatorname{corr}(X,Y)$, but to $\operatorname{cov}(X,Y)$. Clearly, $\langle \phi(x)\,\phi(x)\rangle$ is analogous to the variance quanitity used to define $\operatorname{corr}(X,Y)$, so it would be more correct to say that $\operatorname{corr}(X,Y)$ is analogous to $$\frac{\langle \phi(x)\, \phi(y)\rangle}{\sqrt{\langle \phi^2(x) \rangle\langle \phi^2(y) \rangle}} \rightarrow 0\ \forall\ x\neq y.$$

As for why this happens, there are a number of theorems constraining the form that $\langle \phi(x)\, \phi(y)\rangle$ can take. First, $\langle \phi(x)\, \phi(y)\rangle$ is related to the Green's function (propagator) of the free part of the field theory, which ties it to the equations of motion/Lagrangian/Hamiltonian. The Hamiltonian, in particular, is required to be positive in order to ensure the theory is stable, requiring that the equations of motion can't have more than two derivatives. There a few mathematically equivalent theorems that have the same implication (Ostrogradsky's theorem, classically, and one outlined in Weinberg's Quantum Field Theory books from the quantum side). It's that limit on the number of derivatives, combined with the fact that space-time is four dimensional, that inevitably leads to the propagator containing divergences when space-time is continuous.

As observed, the product $\phi(x)\phi(y)$ is unambiguous and finite for $x \neq y$. However, as $y\to x$, we have that $\phi(x)^2$ possesses an ultraviolet divergence.

In order to make sense of this, we must renormalise this composite operator, call it $[\phi^2].$ In order to do so, we could consider evaluating $\langle T\phi(x)\phi(y)\frac12 \phi(z)^2\rangle$ and adding appropriate counterterms.

In doing so, you will encounter terms like, $m^2 - \frac16(p_1+p_2)^2$ which if interpreted as a Feynman rule would come from a term $(m^2+\frac16 \square)\phi$. Thus, we can replace the unrenormalised operator by a renormalised operator, which by necessity will contain lower dimensional operators to reproduce the result of renormalising the Green's function.

For completeness, the result is,

$$\frac12 [\phi^2] = \left[ 1+\frac{g^2}{64\pi^3 (d-6)}\right]\frac12 \phi^2 + \frac{g\mu^{d/2-3}}{64\pi^3(d-6)}(m^2 + \frac16 \square)\phi + \dots$$

Thus in field theory the divergence of the Green's function should not be shocking as $y\to x$, the root cause is the same as for the other quantities of the theory which diverge.