Two operators $X$ and $Z$ in an infinite dimensional Hilbert space satisfying $X^2=Z^2=I$ and $\{X,Z\}= 0$

I assume that your operators are bounded in a Hilbert space $H$, otherwise I am not sure that the result I am proving is still true in the absence of other assumptions like the essential self-adjointness of the operator $\sum_{a=1}^3 X_a^2$.

Define $iY=ZX$ and next $X_1:=X$, $X_2:=Y$, $X_3:= Z$.

With this definition and your hypotheses, you easily see that $X,Y,Z$ are bounded self-adjoint operators such that $$\{X_a,X_b\}= 2\delta_{ab}I\tag{1}$$ $$[X_a,X_b] = 2 i\sum_c \epsilon_{abc} X_c\:.\tag{2}$$ (2) are the commutation relations of $su(2)$.

As the operators are bounded and self-adjoint, $\sum_{a=1}^3 X_a^2$ is bounded and self-adjoint as well, so in particular it is essentially self-adjoint. Nelson's theorem implies that the Hilbert spece supports a strongly-continuous unitary representation of $SU(2)$, whose Lie-algebra is represented by operators $-iX_a$.

At this point, since $SU(2)$ is compact, Peter-Weyl theorem says that $H$ decomposes into an orthogonal direct sum of finite dimensional irreducible subrepresentatons $H = \oplus_{j}H_j$. Above $j=0,1/2,1,3/2,2,\ldots$. The generators of the $j$-th subrepresentation are the restrictions of the $X_a$ to $H_j$.

Let us focus on these generators $X_{aj}$. As well known $H_j$ is the eigenspace of $X_j^2= \sum_{a=1}^3 X_{aj}^2$ with eigenvalue $4j(j+1)$. So $$X_j^2 = 4j(j+1)I_j$$ but the constraint (1) implies $$3 I_j= 4j(j+1)I_j\:,$$ thus $j=1/2$. The only possible representation appearing in the decomposition of $H$ is the one with $j=1/2$. It may appear infinitely times if $H$ is infinite dimensional. Thus $$H = H_{1/2}\otimes K$$ where $K$ is infinite dimensional if $H$ is. The representation of $X_a$ in $H_{1/2}$ are the ones given in terms of Pauli matrices. We end up with $$X_1 = \sigma_1 \otimes I\:, \quad X_2 = \sigma_2 \otimes I \:, \quad X_3 = \sigma_3 \otimes I\:,$$ where $I$ is the identity operator in $K$.