We have applied ab initio based reflection principle to simulate photoelectron spectra of small water clusters, ranging from monomer to octamer. The role of quantum and thermal effects on the structure of the water photoelectron spectra is discussed within the ab initio path integral molecular dynamics (PIMD) framework. We have used the PIMD method with up to 40 beads to sample the ground state quantum distribution at temperature T = 180 K. We have thoroughly tested the performance of various density functionals (B3LYP, BHandHLYP, M06HF, BNL, LC-ωPBE, and CAM-B3LYP) for the ionization process description. The benchmarking based on a comparison of simulated photoelectron spectra to experimental data and high level equation-of-motion ionization potential coupled clusters with singles and doubles calculations has singled out the BHandHLYP and LC-ωPBE functionals as the most reliable ones for simulations of light induced processes in water. The good performance of the density functional theory functionals to model the water photoelectron spectra also reflects their ability to reliably describe open shell excited states. The width of the photoelectron spectrum converges quickly with the cluster size as it is controlled by specific interactions of local character. The peak position is, on the other hand, defined by long-range non-specific solvent effects; it therefore only slowly converges to the corresponding bulk value. We are able to reproduce the experimental valence photoelectron spectrum of liquid water within the combined model of the water octamer embedded in a polarizable dielectric continuum. We demonstrate that including the long-range polarization and the state-specific treatment of the solvent response are needed for a reliable liquid water ionization description.