First of all, there is a book by Lars Ahlfors available here which treats Moebius transformations in higher dimensional spaces in great detail.
One can define a general Moebius transformation $g: S^n\to S^n$ as follows. I will think of $S^n$ as the 1-point compactification of the Euclidean n-space $E^n$. First, you define the standard inversion in $S^n$
$$
I(x)= \frac{x}{|x|^2}.
$$
It sends the center of inversion (the origin) to infinity and fixes the unit sphere in $E^n$ pointwise. Composing this inversion with dilations and translations you get more inversions. In addition, reflections in hyperplanes in $E^n$ are also regarded as inversions.
Definition. The group $M(S^n)$ of Moebius transformations of $S^n$ consists of finite compositions of inversions defined as above.
For instance, the stabilizer of $\infty$ in this group is the group of Euclidean similarities, i.e. compositions of dilations, translations and orthogonal transformations of $R^n$.
There are various equivalent descriptions of this group. For instance, according to the Liouville's theorem, it equals the group of conformal transformations of $S^n$, i.e. the group of diffeomorphisms of the unit sphere in $E^{n+1}$ which preserve angles between tangent vectors.
Now, as for your question, as written it does not really make sense. To begin with, you have to specify what are $a, b, c, d,$ and $z$. Are they arbitrary quaternions? Then $z$ belongs to the $4$-dimensional real vector space, not $E^3$. Are they purely imaginary quaterions? But then $f(z)$ need not be purely imaginary, even if you define the division properly. Ok, you can ask if this describes the group $M(S^4)$.
You can identify the algebra of quaternions with $E^4$ by using the standard quaternionic norm to define distances. I think, you get an index 2 subgroup of $M(4)$ this way using the formula
$$
f(z)=(az+b)(cz+d)^{-1}.
$$