The only method I know for computing covariant derivatives (which works, at least in principal, all the time), is the Koszul formula. This formula states that \begin{align*} g(\nabla_X(Y),Z) &= \frac{1}{2}\bigg(\partial_X(g(Y,Z)) + \partial_Y(g(X,Z)) - \partial_Z(g(X,Y)) \\ &+ g([X,Y],Z) - g([X,Z],Y) -g([Y,Z],X)\bigg). \end{align*} (Here, $X$, $Y$, and $Z$ can be any three vectors, not just the three in your question.)
Using the fact that your $X$,$Y$, and $Z$ form a basis, if we can compute $g(\nabla_X Y, W)$ where $W$ ranges over each element of the basis, and this completely determines $\nabla_X Y$.
Actually carrying this out is a bit tedious, but here goes. I'll compute $\nabla_X(Y)$. To simplify matters, note that $[X,Y] = 2Z$, $[Y,Z] = 2X$, $[Z,X]=2Y$ so that simplifies things a bit. Further, since $X$, $Y$, and $Z$ actually form an orthonormal basis, all possible inner products of basis elements are constant (either $0$ or $1$), so all their derivatives are $0$. This means the first three terms in the Koszul formula are all $0$ no matter what.
We have $g(\nabla_X(Y), X) = \frac{1}{2}\bigg( g([X,Y],X) - g([X,X],Y) - g([Y,X],X) \bigg) =0$ since each Lie bracket is proportional to $Z$.
An almost identical argument shows $g(\nabla_X(Y),Y) = 0$.
The last calculation is \begin{align*} g(\nabla_X(Y),Z) &= \frac{1}{2}\bigg(g([X,Y],Z) - g([X,Z],Y) -g([Y,Z],X)\bigg) \\ &= \frac{1}{2}\bigg(g(2Z,Z) - g(-2Y,Y) -g(2X,X)\bigg) \\ &= 1. \end{align*}
Putting this all together, we have computed that $\nabla_X(Y) = Z$.
To answer your question in the comments about Killing fields, probably one of the easiest ways to show something is a Killing field is to find a family of isometries whose derivative is that something. In this case, to show $X$ is Killing, one possible choice is to use the matrix $\begin{bmatrix} \cos t & \sin t & 0 & 0 \\\ -\sin t & \cos t & 0 & 0 \\\ 0 & 0 & \cos t & \sin t \\\ 0 & 0 &-\sin t & \cos t\end{bmatrix}$ which is, for each $t$, an isometry of $\mathbb{R}^4$, and thus, also of $S^3$. Applying this to the point $p=(x,y,z,w)^t$ in $S^3$ and taking the derivative at $t=0$, one should get the vector field $X$ at the point $p$. (I confess I didn't work out the details, but something like this will definitely work.)
An alternative characterization is that $X$ is Killing iff $g(\nabla_Y(X),Z) + g(\nabla_Z(X), Y) = 0$ for any choice of $Y$ and $Z$ (which can be checked on a basis). You could do all the computations as above and then just verify that this works out ok.