First of all, that formula for the phase has a sign error, at least if $\alpha$ is defined as in the diagram of the filter structure. Second, your phase plot shows a jump from $-\pi$ to $\pi$, which is no jump at all because a phase of $\pi$ is the same as a phase of $-\pi$. The jump occurs only because the arctangent function returns the principal value of the phase in the interval $[-\pi,\pi]$. You can get rid of the jump by using the phase unwrapping function unwrap
in Matlab/Octave. But note that this is just a cosmetic thing, there is no phase jump in the physical sense of the word.
Now for the derivation of the formula for the phase. From the given filter structure we get the following equation for the $\mathcal{Z}$-transforms of the input and output signals:
$$Y(z)=X(z)z^{-2}+\alpha\left(X(z)-Y(z)z^{-2}\right)\tag{1}$$
From this the transfer function can be computed as
$$H(z)=\frac{Y(z)}{X(z)}=\frac{\alpha+z^{-2}}{1+\alpha z^{-2}}=\frac{\alpha z+z^{-1}}{z+\alpha z^{-1}}=\frac{B(z)}{B\left(\frac{1}{z}\right)}\tag{2}$$
with $B(z)=\alpha z+z^{-1}$.
$H(z)$ as given by $(2)$ is indeed an allpass transfer function with $|H(z)|=1$ for $|z|=1$. Note that the system described by $(2)$ is only stable for $-1<\alpha <1$.
On the unit circle $z=e^{j\omega}$, $H(z)$ can be written as
$$H(e^{j\omega})=\frac{B(e^{j\omega})}{B^*(e^{j\omega})}\tag{3}$$
where $^*$ denotes complex conjugation. Consequently, the phase of $H(e^{j\omega})$ is given by
$$\phi_H(\omega)=2\phi_B(\omega)\tag{4}$$
where $\phi_B(\omega)$ is the phase of the numerator $B(e^{j\omega})$:
$$\begin{align}\phi_B(\omega)&=\arg\big\{\alpha e^{j\omega}+e^{-j\omega}\big\}\\&=\arg\big\{(1+\alpha)\cos(\omega)-j(1-\alpha)\sin(\omega)\big\}\\&=-\arctan\left(\frac{1-\alpha}{1+\alpha}\tan(\omega)\right)\end{align}\tag{5}$$
The phase of $H(e^{j\omega})$ is thus given by
$$\phi_H(\omega)=-2\arctan\left(\frac{1-\alpha}{1+\alpha}\tan(\omega)\right)\tag{6}$$
which differs from the given formula in the sign of $\alpha$.