Molecular dynamics (MD) simulation and subsequent analysis of the macroscopic polarization developed in response to "a posteriori'' applied electric fields or of spontaneous fluctuations in the instantaneous polarization under zero applied field is used to assess the nonlinear optical properties of a polarizable liquid. Three strategies are proposed for the electrostatic analysis, all using as input static "gas phase'' (hyper)polarizabilities, obtained from ab initio calculations. All three strategies are shown to accurately reproduce the experimentally measured refractive index and second nonlinear susceptibility of liquid benzene. The simulation also predicts the distribution of orientations and magnitudes of the local electric fields experienced by the molecules in the liquid, and the nonlinear contributions to the local fields. This approach gives an 8% higher estimate of the second nonlinear susceptibility of liquid benzene than the Lorentz local field factor approach, in better agreement with experimental values.