Spin-flip methods applied to excited-state approaches like the Bethe-Salpeter Equation allow access to the excitation energies of open-shell systems, such as molecules and defects in solids. The eigenstates of these solutions, however, are generally not eigenstates of the spin operator S
2. Even for simple cases where the excitation vector is expected to be, for example, a triplet state, the value of 〈S
2〉 may be found to differ from 2.00; this difference is called "spin contamination." The expectation values 〈S
2〉 must be computed for each excitation vector, to assist with the characterization of the particular excitation and to determine the amount of spin contamination of the state. Our aim is to provide for the first time in the spin-flip methods literature a comprehensive resource on the derivation of the formulas for 〈S
2〉 as well as its computational implementation. After a brief discussion of the theory of the Spin-Flip Bethe-Salpeter Equation and some examples further illustrating the need for calculating 〈S
2〉, we present the derivation for the general equation for computing 〈S
2〉 with the eigenvectors from an SF-BSE calculation, how it is implemented in a Python script, and timing information on how this calculation scales with the size of the SF-BSE Hamiltonian.