The problem was in how to write down the sum $\sum_{\{\vec{S}\}}$. Since this sum is over all the possible vectors $\vec{S}=(S_1,...,S_N)$, (where $S_i\in\{-1,1\}$), we can rewrite this sum like
$$\sum_{S_{1}\in\{-1,1\}}\cdots\sum_{S_{N}\in\{-1,1\}}$$
i.e.
$$\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)$$
Clearly this sum has $2^N$ elements of the type $\prod_{i=1}^{N}e^{\beta HS_{i}}$. Now, since
$$\prod_{i=1}^{N}e^{\beta HS_{i}}=e^{\beta HS_{1}}\cdots e^{\beta HS_{N}},$$
then (1) turns
$$\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)$$
Since each $S_i$ is independent of the others, we can "factorize" the $\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)$$
And finally:
$$\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.$$
Note. This result can be generalized for vectors $\vec{S}=(S_1,...,S_N)$ with components $S_i\in\{1,...,k\}$ for some integer $k$.