Full $NVT$ quantum statistics of the H$_3^+$ ion is simulated at low densities using the path integral Monte Carlo approach. For the first time, the molecular total energy, partition function, free energy, entropy and heat capacity are evaluated in temperatures relevant for planetary atmospheric physics. Temperature and density dependent dissociation recombination reaction balance of the molecule and its fragments above 4000 K is described, and also, the density dependence of thermal ionization above $10 000$ K is demonstrated. We introduce a new well-behaving analytical model for the molecular partition function of the H$_3^+$ ion for the temperature range below dissociation and fit the parameters to the energetics from our simulations. The approach presented here can be regarded as an extension of the traditional {\it ab initio} quantum chemistry beyond the Born--Oppenheimer approximation to description of nonadiabatic phenomena, and even further, account of nuclear quantum dynamics.