We review the results of Monte Carlo studies of chosen nonlinear optical effects in host-guest systems, using methods based on the bond-fluctuation model (BFM) for a polymer matrix. In particular, we simulate the inscription of various types of diffraction gratings in degenerate two wave mixing (DTWM) experiments (surface relief gratings (SRG), gratings in polymers doped with azo-dye molecules and gratings in biopolymers), poling effects (electric field poling of dipolar molecules and all-optical poling) and photomechanical effect. All these processes are characterized in terms of parameters measured in experiments, such as diffraction efficiency, nonlinear susceptibilities, density profiles or loading parameters. Local free volume in the BFM matrix, characterized by probabilistic distributions and correlation functions, displays a complex mosaic-like structure of scale-free clusters, which are thought to be responsible for heterogeneous dynamics of nonlinear optical processes. The photoinduced dynamics of single azopolymer chains, studied in two and three dimensions, displays complex sub-diffusive, diffusive and super-diffusive dynamical regimes. A directly related mathematical model of SRG inscription, based on the continuous time random walk (CTRW) formalism, is formulated and studied. Theoretical part of the review is devoted to the justification of the a priori assumptions made in the BFM modeling of photoinduced motion of the azo-polymer chains.