Edit. For the sake of improving the quality of the post, I modified the proof to make it work for all $M>n$ after Prof. Tao’s comments (the previous version was admittedly way too loose).
In the end I just tried by hand, as I probably should have done from the beginning. The answer is yes. More precisely, if $M>n$, any function $f$ satisfying
$$ |f(x)|\leq\frac {C}{1+|x|^M} $$
belongs to the Hardy space $\mathcal H^1(\mathbb R^n)$ if and only if $\int_{\mathbb R^n} f(x)\,dx=0$, with $\mathcal H^1$ norm bounded by $C$. This is of course sharp in the constant $M$. It would be reasonable to expect that a similar proof works for weighted $L^p$ spaces or similar, as soon as there is enough regularity and a way of controlling the amount of mass of the function at infinity, but I will refrain from guessing.
Proof
We fix a function $\phi\in\mathscr S(\mathbb R^n)$ with non-zero mean. We know that $f$ belongs to $\mathcal H^1(\mathbb R^n)$ if and only if the maximal function defined as
$$ M_\phi f(x):=\sup_{t>0}|\phi_t*f(x)| $$
is in $L^1(\mathbb R^n)$, where $\phi_t(x):=\frac{1}{t^n}\phi(\frac{x}{t})$. In that case, we have
$$ \|f\|_{\mathcal H^1(\mathbb R^n)}\lesssim_{n,\phi} \|M_\phi f\|_{L^1(\mathbb R^n)}. $$
We then want to try to find an upper bound on $|\phi_t*f(x)|$ that does not depend on $t$. By linearity, we assume $C=1$.
First bound
Here we use the zero mean condition:
$$ |\phi_t*f(x)|=\left|\int_{\mathbb R^n} t^{-n}\phi(t^{-1}(x-y))f(y)\,dy\right|= $$
$$ =\left|\int_{\mathbb R^n} t^{-n}\left[\phi(t^{-1}(x-y))-\phi(t^{-1}x)\right]f(y)\,dy\right|, $$
since $f$ has zero mean. Then one can estimate the previous integral with
$$ \leq\left|\int_{\mathbb R^n} t^{-(n+1)}|y|\int_0^1|\nabla\phi(t^{-1}(x-sy))|\,ds\,f(y)\,dy\right|, $$
which, since $\phi$ has bounded gradient, can be estimated by
$$ \lesssim_{n,\phi} t^{-(n+1)}\|\,|\cdot|f\,\|_{L^1(\mathbb R^n)}\lesssim_{n,M,\phi} t^{-(n+1)}. $$
Second bound
We fix $N>0$. Now we look for a bound for all $t$ using the decay of $f$ and knowing that $\phi\lesssim \frac{1}{1+|x|^N}$.
$$ |\phi_t*f(x)|=\left|\int_{\mathbb R^n} t^{-n}\phi(t^{-1}y)f(x-y)\,dy\right|\lesssim $$
$$ \lesssim \int_{\mathbb R^n} \frac{t^{-n}}{1+|t^{-1}y|^N}\frac{1}{1+|x-y|^M}\,dy.$$
Now we split the integral into the regions $|y|\leq \frac{|x|}{2}$ and its complementary:
a)
$$ \left|\int_{|y|\leq \frac{|x|}{2}} \frac{t^{-n}}{1+|t^{-1}y|^N}\frac{1}{1+|x-y|^M}\,dy\right|\leq $$
$$ \leq 2^M \int_{|y|\leq \frac{|x|}{2}} \frac{t^{-n}}{1+|t^{-1}y|^N}\frac{1}{1+|x|^M}\,dy\lesssim$$
$$ \lesssim \frac{1}{1+|x|^M}. $$
b)
$$ \int_{|y|\geq \frac{|x|}{2}} \frac{t^{-n}}{1+|t^{-1}y|^N}\frac{1}{1+|x-y|^M}\,dy\leq $$
$$ \leq \int_{|y|\geq \frac{|x|}{2}} \frac{t^{-n}}{1+|t^{-1}y|^N} dy=$$
$$ = \int_{|z|\geq \frac{|x|}{2t}} \frac{1}{1+|z|^N}\lesssim $$
$$ \lesssim \frac{1}{1+|\frac{x}{t}|^{N-n}}. $$
Hence, the function is bounded up to a constant by
$$ \max\left\{\frac{1}{1+|x|^M},\frac{1}{1+|\frac{x}{t}|^{N-n}}\right\} $$
Final estimate
Now, we take the minimum of the two bounds we obtained and then we take the supremum over $t$:
$$ \sup_{t>0}\min\left\{t^{-(n+1)}, \max\left\{\frac{1}{1+|x|^M},\frac{1}{1+|\frac{x}{t}|^{N-n}}\right\}\right\} $$
Using the property $a\wedge(b\vee c)=(a\wedge b)\vee(a\wedge c)$ one can bring the inner $\max$ outside and split the bound into the maximum of two separate quantities, the first of which yelds trivially
$$ \frac{1}{1+|x|^M}, $$
and the second can be found explicitly simply by plugging $t$ such that $t^{-(n+1)}=\frac{1}{1+|x/t|^{N-n}}$, and can be estimated by
$$ \frac{1}{1+|x|^{-\alpha}}, $$
with
$$ \alpha=\frac{(N-n)(n+1)}{N+1}, $$
and we note that $\alpha>n$ (that is, the second bound is integrable in $x$) for $N$ large enough.
So, if $M>n$, then $M_\phi f$ belongs to $L^1(\mathbb R^n)$ and the claim is proved, namely $f\in\mathcal H^1(\mathbb R^n)$ with
$$ \|f\|_{\mathcal H^1(\mathbb R^n)}\lesssim_{n,M} C. $$
