Parallel MRI is a fast imaging technique that enables the acquisition of highly resolved images in space. It relies on $k$-space undersampling and multiple receiver coils with complementary sensitivity profiles in order to reconstruct a full Field-Of-View (FOV) image. The performance of parallel imaging mainly depends on the reconstruction algorithm, which can proceed either in the original $k$-space (GRAPPA, SMASH) or in the image domain (SENSE-like methods). To improve the performance of the widely used SENSE algorithm, 2D- or slice-specific regularization in the wavelet domain has been efficiently investigated. In this paper, we extend this approach using 3D-wavelet representations in order to handle all slices together and address reconstruction artifacts which propagate across adjacent slices. The extension also accounts for temporal correlations that exist between successive scans in functional MRI (fMRI). The proposed 4D reconstruction scheme is fully \emph{unsupervised} in the sense that all regularization parameters are estimated in the maximum likelihood sense on a reference scan. The gain induced by such extensions is first illustrated on EPI image reconstruction but also measured in terms of statistical sensitivity during a fast event-related fMRI protocol. The proposed 4D-UWR-SENSE algorithm outperforms the SENSE reconstruction at the subject and group-levels (15 subjects) for different contrasts of interest and using different parallel acceleration factors on $2\times2\times3$mm$^3$ EPI images.