Complex systems are composed of mutually interacting components and the output values of these components are usually long-range cross-correlated. We propose a method to characterize the joint multifractal nature of such long-range cross correlations based on wavelet analysis, termed multifractal cross wavelet analysis (MFXWT). We assess the performance of the MFXWT method by performing extensive numerical experiments on the dual binomial measures with multifractal cross correlations and the bivariate fractional Brownian motions (bFBMs) with monofractal cross correlations. For binomial multifractal measures, the empirical joint multifractality of MFXWT is found to be in approximate agreement with the theoretical formula. For bFBMs, MFXWT may provide spurious multifractality because of the wide spanning range of the multifractal spectrum. We also apply the MFXWT method to stock market indexes and uncover intriguing joint multifractal nature in pairs of index returns and volatilities.