In this paper, we propose a new method for automatic change detection in multi-temporal fully polarimetric synthetic aperture radar (PolSAR) images based on multivariate statistical wavelet subband modeling. The proposed method allows us to take into account the correlation structure between subbands by modeling the wavelet coefficients through multi-variate probability distributions. Three types of correlation are investigated: inter-scale, inter-orientation, and inter-polarization dependences. The multivariate generalized Gaussian distribution (MGGD) is used to model the interdependencies between wavelet coefficients at different orientations, scales, and polarizations. Kullback-Leibler similarity measures are computed and used to generate the change map. Simulated and real multilook PolSAR data are employed to assess the performance of the method and are compared to the multivariate Gaussian distribution (MGD) based method. We show that the information embedded in the correlation between subbands improves the accuracy of the change map, leading to better performance. Moreover, the MGGD represents better the correlations between wavelet coefficients and outperforms the MGD.