We consider the system of stochastic differential equation dXt = A(Xt−) dZt, X0 = x, driven by cylindrical α-stable process Zt in R d . We assume that A(x) = (aij(x)) is diagonal and aii(x) are bounded away from zero, from infinity and Hölder continuous. We construct transition density p A (t, x, y) of the process Xt and show sharp two-sided estimates of this density. We also prove Hölder and gradient estimates of x → p A (t, x, y). Our approach is based on the method developed by Chen and Zhang in [11].