This paper is inspired by the work of Nagel and Stein in which the L p (1 < p < ∞) theory has been developed in the setting of the product Carnot-Carath ́eodory spaces M = M 1 ×· · ·×M n formed by vector fields satisfying Ḧormander's finite rank condition. The main purpose of this paper is to provide a unified approach to develop the multiparameter Hardy space theory on product spaces of homogeneous type. This theory includes the product Hardy space, its dual, the product BMO space, the boundedness of singular integral operators and Calder ́on-Zygmund decomposition and interpolation of operators. As a consequence, we obtain the endpoint estimates for those singular integral operators considered by Nagel and Stein (2004). In fact, we will develop most of our theory in the framework of product spaces of homogeneous type which only satisfy the doubling condition and some regularity assumption on the metric. All of our results are established by introducing certain Banach spaces of test functions and distributions, developing discrete Calder ́on identity and discrete Littlewood-Paley-Stein theory. Our methods do not rely on the Journ ́e-type covering lemma which was the main tool to prove the boundedness of singular integrals on the classical product Hardy spaces.