  In this paper, we investigate the optimal strong convergence rate of
numerical approximations for the Cox--Ingersoll--Ross model driven by
fractional Brownian motion with Hurst parameter $H\in(1/2,1)$. To deal with the
difficulties caused by the unbounded diffusion coefficient, we study an
auxiliary equation based on Lamperti transformation. By means of Malliavin
calculus, we prove that the backward Euler scheme applied to this auxiliary
equation ensures the positivity of the numerical solution, and is of strong
order one. Furthermore, a numerical approximation for the original model is
obtained and converges with the same order.
