Parallel MRI is a fast imaging technique that allows reconstruction of full Field-of-View (FoV) images based on under-sampled k-space data acquired using multiple receiver coils with complementary sensitivity profiles. It enables the acquisition of highly resolved images either in space or in time, which is of particular interest in applications like functional neuroimaging. These improvements are counterbalanced by a degraded SNR and the presence of artifacts that depend on the reconstruction algorithm. To improve the performance of the widely used SENSE algorithm, 2D regularization in the wavelet domain has recently been efficiently investigated. In this paper, we extend this work to 3D-wavelet decompositions in order to manipulate all slices together. We illustrate the gain induced by such extension in terms of statistical impact on functional MRI (fMRI) data analysis using a fast-event related protocol. Our results show that our 3D reconstruction algorithm outperforms its 2D counterpart and the SENSE algorithm in several statistical respects at the group-level: peak localization, local maxima, cluster extent, robustness to high acceleration factors.