Abstract: Cox-Ingersoll-Ross (CIR) process is an important tool to study stochastic interest rate and stochastic volatility in financial market. The statistical behavior of fractional CIR process is mainly simulated and discussed in this paper. Since there is no analytical expression of the CIR process, two different functions wfbm and fbmld are used to simulate the fractional Brownian motion, and the Euler-Maruyama (EM) method is used to investigate the expectation and variance of the fractional CIR process. Because the distribution of fractional CIR process can not be expressed by the solution of Fokker-Planck equation, the empirical distribution of fractional CIR process is simulated, and the change of empirical distribution with time is obtained. In order to further verify the algorithm and compare the advantages of the two different algorithms, a backward Euler type scheme of the CIR model and the fractional Ornstein-Uhlenbeck (OU) model with analytical solution is carried out. By comparing figure and table, it is found that simulation by the function fbmld have a very high fitting precision with the theoretical analytical solution with expectation and variance.