You are misunderstanding the notation. Let us denote by $\sigma_i$, $i\in V\subset\mathbb{Z}^2$, the spins in your system.
The magnetization the takes the form
$$
M_V = \sum_{i\in V} \sigma_i.
$$
Therefore (note that the brackets denote expectation with respect to the Gibbs measure, not an average over spins as you seem to be doing),
$$
\langle M_V^2 \rangle
= \Bigl\langle\Bigl( \sum_{i\in V} \sigma_i \Bigr)^2\Bigr\rangle
= \Bigl\langle\sum_{i,j\in V} \sigma_i\sigma_j\Bigr\rangle
= \sum_{i,j\in V} \langle\sigma_i\sigma_j\rangle.
$$
Similarly,
$$
\langle M_V\rangle^2
= \Bigl( \Bigl\langle\sum_{i\in V} \sigma_i \Bigr\rangle\Bigr)^2
= \Bigl( \sum_{i\in V} \langle\sigma_i\rangle \Bigr)^2
= \sum_{i,j \in V} \langle\sigma_i\rangle\langle\sigma_j\rangle.
$$
Its variance is thus given by
$$
\operatorname{Var}(M_V) = \langle M_V^2 \rangle - \langle M_V \rangle^2
= \sum_{i,j\in V} \bigl( \langle\sigma_i\sigma_j\rangle - \langle\sigma_i\rangle\langle\sigma_j\rangle \bigr).
$$
In the thermodynamic limit, using translation invariance, the susceptibility can then be written as
$$
\chi = \lim_{V\to\mathbb{Z}^2} \frac{\beta}{|V|}\operatorname{Var}(M_V) = \beta \sum_{i\in\mathbb{Z}^2} \bigl( \langle\sigma_0\sigma_i\rangle - \langle\sigma_0\rangle\langle\sigma_i\rangle \bigr).
$$
Concerning your computations, note that, when $i\neq 0$, $\langle\sigma_0\sigma_i\rangle\neq 1$ except at $T=0$, while $\langle\sigma_0\rangle=\langle\sigma_i\rangle = 0$ if $T>T_{\rm c}$.
In fact, the truncated 2-point function $\langle\sigma_0\sigma_i\rangle - \langle\sigma_0\rangle\langle\sigma_i\rangle$ decays exponentially fast in $\|i\|$ at all $T\neq T_{\rm c}$ (in pure phases), but the decay becomes a non-summable power law at $T_{\rm c}$, which is the cause of the divergence of the susceptibility $\chi$.
