  In this paper, using the remarkable orthonormal wavelet basis constructed
recently by Auscher and Hyt\"onen, we establish the theory of product Hardy
spaces on spaces ${\widetilde X} = X_1\times X_2\times\cdot \cdot\cdot\times
X_n$, where each factor $X_i$ is a space of homogeneous type in the sense of
Coifman and Weiss. The main tool we develop is the Littlewood--Paley theory on
$\widetilde X$, which in turn is a consequence of a corresponding theory on
each factor space. We define the square function for this theory in terms of
the wavelet coefficients. The Hardy space theory developed in this paper
includes product~$H^p$, the dual $\cmo^p$ of $H^p$ with the special case $\bmo
= \cmo^1$, and the predual $\vmo$ of $H^1$. We also use the wavelet expansion
to establish the Calder\'on--Zygmund decomposition for product $H^p$, and
deduce an interpolation theorem. We make no additional assumptions on the
quasi-metric or the doubling measure for each factor space, and thus we extend
to the full generality of product spaces of homogeneous type the aspects of
both one-parameter and multiparameter theory involving the Littlewood--Paley
theory and function spaces. Moreover, our methods would be expected to be a
powerful tool for developing wavelet analysis on spaces of homogeneous type.
