Empirical studies suggest that asset price fluctuations exhibit “long memory”, “volatility smile”, “volatility clustering” and asset prices present “jump”. To fit the above empirical characteristics of the market, this paper proposes a fractional stochastic volatility jump-diffusion model by combining two fractional stochastic volatilities with mixed-exponential jumps. The characteristic function of the log-return is expressed in terms of the solution of two-dimensional fractional Riccati equations of which closed-form solution does not exist. To obtain the explicit characteristic function, we approximate the pricing model by a semimartingale and convert fractional Riccati equations into a classic PDE. By the multi-dimensional Feynman-Kac theorem and the affine structure of the approximate model, we obtain the solution of the PDE with which the explicit characteristic function and its cumulants are derived. Based on the derived characteristic function and Fourier cosine series expansion, we obtain approximate European options prices. By differential evolution algorithm, we calibrate our approximate model and its two nested models to S&P 500 index options and obtain optimal parameter estimates of these models. Numerical results demonstrate the pricing method is fast and accurate. Empirical results demonstrate our approximate model fits the market best among the three models.