El problema estaba en cómo escribir la suma $\sum_{\{\vec{S}\}}$. Dado que esta suma es sobre todos los posibles vectores $\vec{S}=(S_1,...,S_N)$, (donde $S_i\in\{-1,1\}$), podemos reescribir esta suma como
$$\sum_{S_{1}\in\{-1,1\}}\cdots\sum_{S_{N}\in\{-1,1\}}$$ es decir, $$\sum_{\{\vec{S}\}}\prod_{i=1}^{N}e^{\beta HS_{i}}=\sum_{S_{1}\in\{-1,1\}}\cdots\sum_{S_{N}\in\{-1,1\}}\prod_{i=1}^{N}e^{\beta HS_{i}} \qquad(1)$$
Claramente esta suma tiene $2^N$ elementos del tipo $\prod_{i=1}^{N}e^{\beta HS_{i}}. Ahora, dado que
$$\prod_{i=1}^{N}e^{\beta HS_{i}}=e^{\beta HS_{1}}\cdots e^{\beta HS_{N}},$$ entonces (1) se convierte en
$$\sum_{\{\vec{S}\}}\prod_{i=1}^{N}e^{\beta HS_{i}}=\sum_{S_{1}\in\{-1,1\}}\cdots\sum_{S_{N}\in\{-1,1\}}e^{\beta HS_{1}}\cdots e^{\beta HS_{N}} \qquad(2)$$ Dado que cada $S_i$ es independiente de los demás, podemos "factorizar" las $\Sigma$'s: $$\sum_{\{\vec{S}\}}\prod_{i=1}^{N}e^{\beta HS_{i}}=\left(\sum_{S_{1}\in\{-1,1\}}e^{\beta HS_{1}}\right)\cdots\left(\sum_{S_{N}\in\{-1,1\}}e^{\beta HS_{N}}\right)$$ Y finalmente: $$\sum_{\{\vec{S}\}}\prod_{i=1}^{N}e^{\beta HS_{i}}=\prod_{j=1}^{N}\sum_{S_{j}\in\{-1,1\}}e^{\beta HS_{j}} \qquad Q.E.D.$$
Nota. Este resultado puede ser generalizado para vectores $\vec{S}=(S_1,...,S_N)$ con componentes $S_i\in\{1,...,k\}$ para algún entero $k$.