Why is $\int\int f(x)f(y) |x-y|dxdy$ negative? The Setup
Let $f:\mathbb{R} \to\mathbb{R}$ be a smooth function with support in the 
interval $[-R,R]$ and satisfying $\int f = 0$.  By manipulating some integrals, I found the surprising inequality
$$
\int\int f(x)f(y) |x-y| \,dxdy \leq 0.
$$
My questions are 


*

*Is this inequality true? Is my derivation below correct?

*Is there a different reason why this is true?

*Are there other nontrivial functions $g(x,y)$ for which 
$$
\int \int f(x)f(y)g(x,y)\,dxdy \leq 0?
$$
What properties should I expect of these functions $g(x,y)$? 


The Derivation
Let $H:\mathbb{R}\to\mathbb{R}$ be the Heaviside function, so $H(x) = 1$
for $x>0$, and $H(x) = 0$ for $x\leq 0$.  I was interested in the convolution $H\ast f(x)$ defined by 
$$
H\ast f(x) = \int f(y) H(x-y)\,dy = \int_{-\infty}^x f(y)\,dy.
$$
Notice that since $f$ has support in $[-R,R]$ and $\int f = 0$, 
$H\ast f$ also has support in $[-R,R]$.
Now we're prepared to start the derivation of the inequality, beginning
with the simple observation
$$
0 \leq \int |H\ast f(x)|^2\,dx = \int_{-R}^R |H\ast f(x)|^2\,dx.
$$
First expand the square and the convolution, and then
rearrange the order of the integrals:
\begin{align*}
\int_{-R}^R |H\ast f(x)|^2\,dx 
&= \int_{-R}^R \left(\int_{-R}^R f(y)H(x-y)\,dy\right) \left(\int_{-R}^R f(z)H(x-z)\,dz\right) \,dx \\
&= \int_{-R}^R\int_{-R}^R f(y)f(z) \left(\int_{-R}^R H(x-y)H(x-z)\,dx\right)\,dydz
\end{align*}
The integrand $H(x-y)H(x-z)$ is $1$ when both $x>y$ and $x>z$, and $0$ otherwise.  Thus the integral comes out to $\min\{R-y,R-z\}$. We plug this
back into the integral, and use again the fact that $\int f = 0$:
\begin{align*}
\int_{-R}^R\int_{-R}^R f(y)f(z) \min\{R-y,R-z\} \,dydz
&= \int_{-R}^R \int_{-R}^R f(y)f(z)(\min\{-y,-z\} - R)\,dydz
\\&= \int_{-R}^R\int_{-R}^R f(y)f(z)\min\{-y,-z\}\,dydz.
\end{align*}
Now split up the domain according to which of $-y$ or $-z$ is smaller:
\begin{align*}
\int_{-R}^R\int_{-R}^R f(y)f(z)\min\{-y,-z\}\,dydz
&= -\int_{y=-R}^R f(y)\left(y\int_{z=-R}^y f(z)\,dz + \int_{z=y}^R zf(z)\,dz\right)\,dy. 
\end{align*}
Since $\int_{-R}^R f(z)\,dz = 0$, $\int_{-R}^yf(z)\,dz = -\int_y^R f(z)\,dz$, so we can combine the integrals to conclude that 
\begin{align*}
\int |H\ast f(x)|^2\,dx &= \int_{-R}^R f(y) \int_y^R f(z) (y-z)\,dz\,dy\\
&=  - \int_{-R}^R f(y) \int_y^R f(z) |y-z|\,dz\,dy 
\end{align*}
To get to the integral above, swap the names of the dummy variables and
then swap the order of the integrals:
\begin{align*}
\int |H\ast f(x)|^2\,dx &= -\int_{-R}^R f(z) \int_z^R f(y) |z-y|\,dy\,dz\\
&= -\int_{-R}^R f(y) \int_{-R}^y f(z) |z-y|\,dy\,dz.
\end{align*}
In conclusion, adding both of these formulas together, we obtain
$$
0\leq 2\int |H\ast f(x)|^2\,dx = -\int\int f(y)f(z)|y-z|\,dy\,dz.
$$
The derivation is quite long so there's a good chance I've made a mistake.
 A: You don't need smoothness: indeed, integrability is probably enough. Perhaps the "most illuminating" proof is to use Fourier series: let $f \in L^2[-R,R]$, then
$$ \tilde{f}(k) = \frac{1}{2R}\int_{-R}^R e^{i\pi kx/R}f(x) \, dx $$
all exist, and Parseval's identity says that
$$ C\sum_{k=-\infty}^{\infty} \overline{\tilde{f}(k)} \tilde{g}(k) = \int_{-R}^R \overline{f(x)}g(x) \, dx, $$
where $C$ is an uninteresting normalisation constant (most likely $2R$).
Also, we have the convolution theorem, so taking $g=f \star h $, where $h(x)=-\lvert x \rvert$, we find that
$$ \tilde{g}(k) = \tilde{f}(k) \tilde{h}(k), $$
and a calculation shows that
$$ \tilde{h}(k) = \begin{cases} -R/2 & k=0 \\ \frac{(1-(-1)^k) R}{k^2 \pi^2} & k \neq 0 \end{cases}, $$
and in particular, all of these but $\tilde{h}(0)$ are nonnegative. Also, $\tilde{f}(0)=\int f = 0$. So therefore we have
$$ -\iint f(x) f(y) \lvert x -y \rvert \, dx \, dy = 2R \sum_{k=-\infty}^{\infty} \lvert \tilde{f}(k) \rvert^2 \tilde{h}(k) = 2R\sum_{k \neq 0} \frac{(1-(-1)^k)R}{k^2\pi^2} \lvert \tilde{f}(k) \rvert^2 \geq 0 $$

More generally, it is apparent that if $g(x,y)=h(x-y)$, you require that all the Fourier coefficients of $h$ be nonnegative save for $k=0$. Beyond that, it's considerably more difficult (and in fact I am writing a paper that includes some such functions at the moment).
You can also generalise this to infinite intervals, if you have $\int \lvert x\rvert f(x) \, dx < \infty$, by a similar argument using the Fourier transform.
A: Setting $F(x)=\int_{-R}^xf(y)\mathrm dy=H\star f(x)$ and $G(x)=\int_{-R}^xyf(y)\mathrm dy$, we have $F'(x)=f(x)$, $G'(x)=xf(x)$,
$F(-R)=F(R)=0$. Proceed with integration 
$$\begin{split}A&=\iint f(x)f(y)|x-y|\,\mathrm dy\\
&=\int_{-R}^Rf(x)\left[x\int_{-R}^xf(y)\mathrm dy-\int_{-R}^xyf(y)\mathrm dy-x\int_x^Rf(y)\mathrm dy+\int_x^Ryf(y)\mathrm dy\right]\,\mathrm dx\\
&=\int_{-R}^Rf(x)\left(xF(x)-G(x)+xF(x)+G(R)-G(x)\right)\mathrm dx\\
&=\big[F(x)\left(2xF(x)-2G(x)+G(R)\right)\big]_{-R}^R-\int_{-R}^RF(x)\left(2xf(x)+2F(x)-2xf(x)\right)\mathrm dx\\
&=-2\int_{-R}^RF(x)^2\mathrm dx\le0.
\end{split}$$
We have use integration by parts and $\int_x^Rf(y)\mathrm dy=F(R)-\int_{-R}^xf(y)\mathrm dy=-F(x)$. 
