γ-γ correlation functions are mathematical expressions that describe the angular distribution of cascade γ-rays emitted from an atomic nucleus. Cascade transitions may occur in either a two-step deexcitation or through an excitation-deexcitation process of a particular energy level inside the nucleus. In both cases, the nucleus returns to its ground energy state. Spin and parity of the excited state can be determined experimentally using the asymmetry of the angular distribution of the emitted radiation. γ-γ correlation functions are only valid for point-like targets and detectors. In the real experiments, however, neither the target nor the detector is point-like. Thus, misassignment of the spin-parity of energy levels may easily take place if only the analytical equations are considered. Here, we develop a new Monte Carlo simulation method of the γ-γ correlation functions to account for the extended target and detector involved in spin-parity measurements using nuclear resonance fluorescence of nuclei. The proposed simulation tool can handle arbitrary geometries and spin sequences. Additionally, we provide numerical calculations of a parametric study on the influence of the detection geometry on the angular distribution of the emitted γ-rays. Finally, we benchmark our simulation by comparing the simulation-estimated asymmetry ratios with those measured experimentally. The present simulation can be employed as a kernel of an implementation that simulates the nuclear resonance fluorescence process.