Define a unitary representation $U$ of $\mathbf R^{2n}$ on $L^2(\mathbf R^n)$ by $(U_{z,y}g)(w)=g(w+z-y)$. (Special case for $G=\mathbf R^n$ of the so-called two-sided regular representation of $G\times G$ on $L^2(G)$.) Then we have \begin{align} (g,U_{z,y}g)&=\int_{\mathbf R^n}\overline{g(w)}g(w+z-y)\,dw\\ &=\int_{\mathbf R^n}\overline{g(x-z)}g(x-y)\,dx\\ &=f(z,y). \end{align} ThereforeThus $f$ is a diagonal matrix coefficient of a unitary representation, and as such is well known to be positive-definite definite: indeed $a_{ij}= f(a_i-a_j)=(U_{a_j}g,U_{a_i}g)$$a_{ij}= f(a_j-a_i)=(U_{a_i}g,U_{a_j}g)$ gives $\sum_{i,j}a_{ij}\bar c_ic_j=\left\|\sum_ic_iU_{a_i}g\right\|^2\geqslant 0$.