Let $f(x; \sigma) = \frac{1}{\sigma\sqrt{2\pi}}\cdot e^{-\frac{x^2}{2\sigma^2}}$ be the probability density function of a normal distribution $\mathcal{N}(0, \sigma^2)$. We consider a *discrete* normal distribution over $\mathbb{Z}$, obtained by sampling from $\mathcal{N}(0, \sigma^2)$ and rounding the result to the nearest integer.
The cumulative distribution function of this discrete normal is then given by
\begin{equation*}
F(x;\sigma) = \int\limits_{x-\frac{1}{2}}^{x+\frac{1}{2}} f(\theta; \sigma)\ \mathrm{d}\theta.
\end{equation*}

From the symmetry of the normal around $0$, we can easily deduce that the expected value of our discrete normal is $0$ as in the continuous case.

My question however, is how one could estimate the variance of this discrete distribution, which is thus given by \begin{equation*} \sum\limits_{x=-\infty}^{\infty} F(x; \sigma)\cdot x^2 \;. \end{equation*} From some numerical simulations, it would seem that this value is very close to $\sigma^2$, the variance of the continuous distribution. Is there a way to either prove what the exact value is, or at least give some tight bound with respect to $\sigma^2$?