Unfortunately, the [tex]\langle \hat{a}^\dagger \hat{a}^\dagger \hat{a} \hat{a} \rangle[/tex]- term does not factorize to [tex]\langle \hat{n} \hat{n} \rangle[/tex].
Starting with the equal time correlation, you have [tex]\hat{a}^\dagger \hat{a} =\hat{a}\hat{a}^\dagger -1[/tex], so that the above term factorizes to [tex]\langle \hat{n} (\hat{n}-1) \rangle[/tex].
This makes sense as the detection of a photon changes the light field by destroying that photon. However, if you are interested in the time dependence of the correlation function, the term [tex]\hat{a}^\dagger (t+\tau) \hat{a}(t)[/tex] can be anything from [tex]\hat{a} (t) \hat{a}^\dagger(t+\tau) -1[/tex] to [tex]\hat{a} (t) \hat{a}^\dagger(t+\tau)[/tex] depending on the magnitude of [tex]\tau[/tex] compared to the coherence time of the light.
Finding a solution for this problem is rather demanding and depends on the kind of light field you are interested in. I suppose you are interested in lasers. "Classical" atom lasers are for example discussed within a birth-death model in "Photon statistics of a cavity-QED laser: A comment on the laser–phase-transition analogy" by P.R. Rice and H.J. Carmichael, Phys. Rev. A 50, 4318–4329 (1994). Semiconductor lasers are treated using the cluster expansion method in "Semiconductor model for quantum-dot-based microcavity lasers" by C. Gies et al., Phys. Rev. A 75, 013803 (2007).
If you could tell me what kind of system or light field you have in mind, I might be able to come up with more suitable references for your case.