  In the framework of the Standard Model we present new theoretical predictions
for the branching ratios, double and single differential distributions and
forward -- backward leptonic asymmetries for the $\bar B_s \to \mu^+ \mu^- e^+
e^-$ decay. In our consideration we take into account the $\phi(1020)$ --
resonance contribution; the main contributions of four charmonium resonances :
$\psi(3770)$, $\psi(4040)$, $\psi(4160)$ and $\psi(4415)$; $u\bar{u}$ --
resonant contribution from $\rho(770)$ and $\omega(782)$; tails contributions
from $J/\psi$ and $\psi(2S)$ resonances; non -- resonant contribution of the
$b\,\bar{b}$ - pairs, bremsstrahlung and the contribution of the weak
annihilation. We provide the prediction for the branching ratio of $\bar B_s
\to \mu^+ \mu^- e^+ e^-$ decay with and without the $\phi(1020)$ -- resonance
contribution. We use the model of vector meson dominance (VMD) for calculation
of resonances contributions and take into account all substantive terms in
$\bar B_s \to \mu^+ \mu^- e^+ e^-$ amplitude that was not considered in the
previously papers.
