Positivity of correlation functions in the ferromagnetic Ising model

Yes, it is a simple consequence of the first Griffiths (or GKS) inequality, which states that $$ langle sigma_A rangle geq 0 $$ for any finite set of vertices $A$. Above, I have used the standard notation $sigma_A = prod_{iin A}sigma_i$.

Griffiths' first inequality holds at any temperature and for any nonnegative magnetic field. Actually, it also holds for (ferromagnetic) interactions of arbitrary (possibly infinite) range. The proof is very easy (one simply expands the Boltzmann weight in a Taylor series and sum the resulting expression over the spins) and can be found, for instance, in Section 3.8.1 of this book (the version given there actually covers a substantially more general situation than described here).

In fact, one can even show that $$ langle sigma_A rangle > 0 $$ at all finite temperatures, for any finite set $A$ containing an even number of vertices. This stronger version follows from the second Griffiths inequality stated below (see Exercise 3.25 in the book for the 2-point function and apply the second Griffiths inequality for general $A$ containing an even number of vertices).

---

Maybe one comment: it is crucial that the boundary condition is free, periodic or $+$ (or a mixture of $+$ and free, for instance). You cannot use a mixture of $+$ and $-$ boundary conditions (such as Dobrushin boundary conditions), as one can construct counterexamples to the claim in this case.

For pure $-$ boundary conditions, the result remains true when $h=0$ and $A$ contains an even number of vertices, since in that case the expectation coincides with the corresponding expectation under $+$ boundary condition by symmetry.

---

Just to be complete, the inequality is called the first Griffiths inequality, because there is a second one, which you might also find interesting: under the same assumptions, for any finite sets of vertices $A$ and $B$, $$ langle sigma_Asigma_B rangle geq langle sigma_A rangle langle sigma_B rangle , $$ which shows that the random variables $sigma_A$ and $sigma_B$ are positively correlated (their covariance is nonnegative).