Group delay of the integrator
Your integrator is a perfect accumulator, which is an unstable (marginally stable) system with a pole at $z=1$:

$$H(z)=\frac{1}{1-z^{-1}}\tag{1}$$

This means that the value of the frequency response at $\omega=0$ (corresponding to the pole at $z=1$) is undefined. Consequently, the value of the group delay at $\omega=0$ is also undefined. Note that you can write down the frequency response of an ideal accumulator using a Dirac delta impulse (see this question and its answer), but this doesn't help in finding a value of the frequency response (or the group delay) at $\omega=0$.

For $\omega\neq 0$ (or, more generally, $\omega\neq 2\pi k$), we can easily determine the frequency response and group delay of the system:

$$H(e^{j\omega})=\frac{1}{1-e^{-j\omega}}=e^{j\omega/2}\frac{1}{e^{j\omega/2}-e^{-j\omega/2}}=e^{j\omega/2}\frac{1}{2j\sin(\omega/2)},\quad\omega\neq 2\pi k\tag{2}$$

From $(2)$, we see that the phase $\phi(\omega)$ equals $\omega/2$ plus a constant term which is irrelevant for the group delay. Consequently, the group delay is given by

$$\tau_g(\omega)=-\frac{d\phi(\omega)}{d\omega}=-\frac12,\quad\omega\neq 2\pi k\tag{3}$$

which is what you found using Matlab. However, Matlab doesn't tell you what's going on at $\omega=0$, where the value of the group delay is indeed undefined.