Using the staggered-grid finite difference method, a numerical modeling algorithm for a 3D arbitrary anisotropic Earth is implemented based on magnetotelluric (MT) theory. After the validation of this algorithm and comparison with predecessors, it was applied to several qualitative and quantitative analyses containing electrical anisotropy and a simple 3D prism model. It was found that anisotropic parameters for ρ1, ρ2, and ρ3 play almost the same role in affecting 3D MT responses as in 1D and 2D without considering three Euler's angles αS, αD, and αL. Significant differences appear between the off-diagonal components of the apparent resistivity tensor and also between the diagonal components in their values and distributing features under the influence of 3D anisotropy, which in turn help to identify whether the MT data are generated from 3D anisotropic earth. Considering the deflecting effects arising from the inconsistency between the anisotropy axes and the measuring axes, some strategies are also provided to estimate the deflecting angles associated with anisotropy strike αS or dip αD, which may be used as initial values for the 3D anisotropy inversion.