IIR Hilbert Transformer
This is achievable with two parallel all pass filters.
The two all pass filters synthesize an odd ordered low pass filter whose pass band extends from -90º to +90º in the z-domain. (I will discuss this below).
$$
G_{lowpass}(z) = \frac{A_0(z)+z^{-1}A_1(z)}{2}
$$
The low pass filter is then rotated by +90º so that its pass band extends from 0º to 180º, which approximates the Hilbert transform. Rotation mathematically is:
$$
H_{Hilbert}(z) =G_{lowpass}(-jz)
$$
As a consequence, one of the all pass filters becomes completely imaginary - this is the Hilbert transformed signal path. By necessity, since the filter is approximating this Hilbert transform, then the imaginary branch must be applying a 90º phase shift relative to the output of $A_0(z)$.
$$
H_{Hillbert}(z)=G_{lowpass}(-jz)=\frac{A_0(z)+jz^{-1}A_1(z)}{2}
$$

Example with Elliptical Filter
The following example was taken from Handbook for Digital Signal Processing pg. 920
You can also use this program that I interpreted to design a new half-band elliptic filter with a higher order.
Starting with a real, half-band elliptical low pass filter $G(z)$ whose frequency response satisfies
$$
1-\delta_1\leq|G(e^{j\omega})|\leq1\&gt;\&gt;\&gt;\&gt;\&gt;for\&gt;\&gt;0\leq\omega\leq\omega_p
$$
$$
\&gt;\&gt;\&gt;|G(e^{j\omega})|\leq\delta_2\&gt;\&gt;\&gt;\&gt;\&gt;for\&gt;\&gt;\omega_s\leq\omega\leq\pi
$$
To obtain half band we choose:
$$
\omega_s+\omega_p=\pi
$$
$$
(1-\delta_1)^{2}+\delta_2^{2}=1
$$
But I just used the examples given seventh order poles or an elliptical low pass half-band filter:
$$
z=0,\&gt;\&gt;\&gt;z=\pm j0.436688,\&gt;\&gt;\&gt;z=\pm j0.743707,\&gt;\&gt;\&gt;z=\pm j0.927758
$$
$$
G_{lowpass}(z) = \frac{A_0(z)+z^{-1}A_1(z)}{2}
$$
Distribute the poles between the two all pass filters like so:

$$
A_0(z)=\frac{z^{-2}+0.190696}{1+0.190696z^{-2}}\cdot\frac{z^{-2}+0.860735}{1+0.860735z^{-2}}
$$
$$
A_1(z)=\frac{z^{-2}+0.553100}{1+0.553100^{-2}}
$$
Rotate by 90º
$$
H_{Hilbert}(z)=G_{lowpass}(-jz)=\frac{A_0(-jz)+jz^{-1}A_1(-jz)}{2}
$$
$$
A_0(-jz)=\frac{0.190696-z^{-2}}{1-0.190696z^{-2}}\cdot\frac{0.860735-z^{-2}}{1-0.860735z^{-2}}
$$
$$
A_1(z)=\frac{0.553100-z^{-2}}{1-0.553100^{-2}}
$$
Phases of $A_0(z)$ and $A_1(z)$:

Magnitude Response of $H_{Hilbert}(z)$:

Phase Difference Between $A_0(z)$ and $A_1(z)$ Branches: