The main purpose of this paper is to establish the theory of the multi-parameter Hardy spaces H p H^p ( 0 > p ≤ 1 0>p\leq 1 ) associated to a class of multi-parameter singular integrals extensively studied in the recent book of B. Street (2014), where the L p L^p ( 1 > p > ∞ ) (1>p>\infty ) estimates are proved for this class of singular integrals. This class of multi-parameter singular integrals are intrinsic to the underlying multi-parameter Carnot-Carathéodory geometry, where the quantitative Frobenius theorem was established by B. Street (2011), and are closely related to both the one-parameter and multi-parameter settings of singular Radon transforms considered by Stein and Street (2011, 2012a, 2012b, 2013). More precisely, Street (2014) studied the L p L^p ( 1 > p > ∞ ) (1>p>\infty ) boundedness, using elementary operators, of a type of generalized multi-parameter Calderón Zygmund operators on smooth and compact manifolds, which include a certain type of singular Radon transforms. In this work, we are interested in the endpoint estimates for the singular integral operators in both one and multi-parameter settings considered by Street (2014). Actually, using the discrete Littlewood-Paley-Stein analysis, we will introduce the Hardy space H p H^p ( 0 > p ≤ 1 0>p\leq 1 ) associated with the multi-parameter structures arising from the multi-parameter Carnot-Carathéodory metrics using the appropriate discrete Littlewood-Paley-Stein square functions, and then establish the Hardy space boundedness of singular integrals in both the single and multi-parameter settings. Our approach is much inspired by the work of Street (2014) where he introduced the notions of elementary operators so that the type of singular integrals under consideration can be decomposed into elementary operators.