[guided]The variance calculation has two purposes: use independence to gain the factor $1/n$, and use the $h^{-d}$ scaling of the kernel to determine the size of one summand. For $x,y\in\mathcal X$, write
\begin{align*}
k_{h,x}(y)
=
\sum_{m\in\mathbb Z^d}
h^{-d}K\left(\frac{x-y+m}{h}\right).
\end{align*}
This is the function of the observation variable $y$ that is evaluated at $X_j$. Hence
\begin{align*}
\hat f_h(x)=\frac{1}{n}\sum_{j=1}^n k_{h,x}(X_j).
\end{align*}
Because the random variables $X_1,\dots,X_n$ are independent and identically distributed, the variance of the average is
\begin{align*}
\operatorname{Var}_f(\hat f_h(x))
=
\frac{1}{n}\operatorname{Var}_f(k_{h,x}(X_1)).
\end{align*}
We then bound the variance by the second moment:
\begin{align*}
\operatorname{Var}_f(k_{h,x}(X_1))
\le
\mathbb E_f[|k_{h,x}(X_1)|^2].
\end{align*}
Since $X_1$ has density $f$ with respect to $\mathcal L^d$ on $\mathcal X$,
\begin{align*}
\mathbb E_f[|k_{h,x}(X_1)|^2]
=
\int_{\mathcal X}|k_{h,x}(y)|^2 f(y)\,d\mathcal L^d(y).
\end{align*}
Now integrate over $x$ and use the uniform density bound $\|f\|_\infty\le B$:
\begin{align*}
\int_{\mathcal X}\operatorname{Var}_f(\hat f_h(x))\,d\mathcal L^d(x)
\le
\frac{1}{n}
\int_{\mathcal X}\int_{\mathcal X}
|k_{h,x}(y)|^2 f(y)
\,d\mathcal L^d(y)\,d\mathcal L^d(x).
\end{align*}
Using $f(y)\le B$ for $\mathcal L^d$-almost every $y\in\mathcal X$, we get
\begin{align*}
\int_{\mathcal X}\operatorname{Var}_f(\hat f_h(x))\,d\mathcal L^d(x)
\le
\frac{B}{n}
\int_{\mathcal X}\int_{\mathcal X}
|k_{h,x}(y)|^2
\,d\mathcal L^d(y)\,d\mathcal L^d(x).
\end{align*}
The remaining integral is the only place where periodization needs care. Squaring the periodized kernel creates cross terms, so we do not identify its squared mass with the squared mass of $K$. Instead choose $R_K<\infty$ such that $\operatorname{supp}K\subset B(0,R_K)$ and define
\begin{align*}
N_K
=
\sup_{u\in\mathbb R^d}\#\left\{m\in\mathbb Z^d:u+m\in B(0,R_K)\right\}.
\end{align*}
This number is finite because $B(0,R_K)$ is bounded and the lattice $\mathbb Z^d$ is locally finite. For $0<h\le 1$, at most $N_K$ terms in $k_{h,x}(y)$ can be non-zero for any fixed $(x,y)$. The finite-sum Cauchy-Schwarz inequality therefore gives
\begin{align*}
|k_{h,x}(y)|^2
\le
N_K\sum_{m\in\mathbb Z^d}h^{-2d}\left|K\left(\frac{x-y+m}{h}\right)\right|^2.
\end{align*}
The right-hand side is non-negative, so Tonelli's theorem allows us to interchange the sum and the integral over $x$. For a fixed $m\in\mathbb Z^d$, use the substitution $z=(x-y+m)/h$; the measure transforms as
\begin{align*}
d\mathcal L^d(x)=h^d\,d\mathcal L^d(z).
\end{align*}
The image of $x\in\mathcal X$ is $h^{-1}(\mathcal X-y+m)$. The translated cubes $\mathcal X-y+m$, indexed by $m\in\mathbb Z^d$, tile $\mathbb R^d$ up to boundary sets of $\mathcal L^d$-measure zero; scaling by $h^{-1}$ preserves this tiling property up to null sets. Hence the sum of the changed-variable integrals is exactly the integral over $\mathbb R^d$:
\begin{align*}
\sum_{m\in\mathbb Z^d}h^{-2d}
\int_{\mathcal X}\left|K\left(\frac{x-y+m}{h}\right)\right|^2
\,d\mathcal L^d(x)
=
h^{-d}\int_{\mathbb R^d}|K(z)|^2\,d\mathcal L^d(z).
\end{align*}
Combining this with the finite-overlap inequality gives
\begin{align*}
\int_{\mathcal X}|k_{h,x}(y)|^2\,d\mathcal L^d(x)
\le
N_Kh^{-d}\|K\|_{L^2(\mathbb R^d)}^2.
\end{align*}
Integrating in $y$ over $\mathcal X$ contributes the factor $\mathcal L^d(\mathcal X)=1$, so
\begin{align*}
\int_{\mathcal X}\operatorname{Var}_f(\hat f_h(x))\,d\mathcal L^d(x)
\le
\frac{BN_K\|K\|_{L^2(\mathbb R^d)}^2}{nh^d}.
\end{align*}
For a boundary-corrected estimator, the same argument is encoded in the assumed uniform $L^2$ kernel-size bound: the corrected kernel sections have squared $L^2$ mass at most a constant times $h^{-d}$ uniformly in the boundary location. Thus the variance remains of order $(nh^d)^{-1}$.[/guided]