We consider the system of stochastic differential equation $dX_t = A(X_{t-}) \, dZ_t$, $ X_0 = x$, driven by cylindrical $α$-stable process $Z_t$ in $\mathbb{R}^d$. We assume that $A(x) = (a_{ij}(x))$ is diagonal and $a_{ii}(x)$ are bounded away from zero, from infinity and Hölder continuous. We construct transition density $p^A(t,x,y)$ of the process $X_t$ and show sharp two-sided estimates of this density. We also prove Hölder and gradient estimates of $x \to p^A(t,x,y)$. Our approach is based on the method developed by Chen and Zhang.