short answer
The reason the phase is reversed is because you are computing $$\frac{\texttt{cross}(AB)}{\texttt{auto}(A)}$$ when you should be computing $$\frac{\texttt{cross}(BA)}{\texttt{auto}(A)}$$
long answer
The FRF defined by the OP is a often-used estimator of the true frequency response of a LTI system.
In a practical measurement scenario, there is always some noise (at least in the measured output), which is problematic if one is to use a straight-forward frequency-domain ratio such as $H = Y/X$, $Y$ being the noisy output and $X$ the perhaps-noisy input.
One way to mitigate the effect of noise in the output is to use transfer function estimators. A straight-forward one is:
$$\tilde{H}(\omega) = \frac{P_{yx}(\omega)}{P_{x}(\omega)}$$
where $P_{yx}(\omega)$ is the Cross Power Spectral Density of $x$ and $y$ and $P_{x}(\omega)$ is the Power Spectral Density of $x$. These can be computed in many different ways, such as the Welch method.
Notice the order $_{yx}$ in $P_{yx}$ for the CPSD (see short answer above).
So, why do we use this estimator?
In Frequency domain, we have: $$Y = X\cdot H$$ so you'd be tempted to
simply re-arrange to get the frequency response: $$H = \dfrac{Y}{X}$$
This is problematic if there's noise $N_0$ in the output: the resulting estimate is: $$\tilde{H} = \frac{Y+N_0}{X} = H +
\frac{N_0}{X} = H+\frac{1}{\texttt{SNR}}$$ The lower the $\texttt{SNR}$, the higher the estimator's bias.
Here comes the frequency averaging estimator: Let's look at the definition of the Cross Power Spectral Density: \begin{align} P_{yx}
&= \mathbb{E}\left[ \overline{X}\cdot(Y+N_0) \right]\\ &=
\mathbb{E}\left[ \overline{X}\cdot Y + \overline{X}\cdot N_0\right]\\
&= \mathbb{E}\left[ \overline{X}\cdot HX + \overline{X}\cdot
N_0\right]\\ &= H\cdot\mathbb{E}\left[\overline{X}X\right] +
\mathbb{E}\left[\overline{X}\cdot N_0\right]\\ &= HP_{x} +
\mathbb{E}\left[\overline{X}\cdot N_0\right] \tag{1} \end{align}
Assuming the noise in the output, $N_0$, is un-correlated with the input, we have $$\mathbb{E}\left[\overline{X}\cdot N_0\right] = 0$$ so (1) becomes:
$$P_{yx} = HP_{x}$$ and we get $$H =
\frac{P_{yx}}{P_{x}} = \tilde{H} $$ i.e., the estimator is perfect.
Of course this is just
theory, since the Expected Value Operator $\mathbb{E}$ is a
statistical operator based in infinite time. In a practical
measurement scenario (finite time, noisy), the cross term $\overline{X}\cdot N_0$ is never actually $0$, but by using frequency averaging (Welch's method) to compute $P_{yx}$, the cross term $\overline{X}\cdot N_0$ does go towards 0 since $X$ and $N_0$ are un-correlated, resulting in a better estimate of $H(\omega)$ than simply dividing the output spectrum by the input spectrum.
Note: When the noise is in the input, and un-correlated with the output, the estimator should be:
$$\tilde{H}(\omega) = \frac{P_y(\omega)}{P_{xy}(\omega)}$$