Why insist on dyadic decomposition of cubes when you don't need it?
Observe that the maximal function $\mathcal{M}f = \mathcal{M} |f|$ by definition. And observe that since $\phi$ is positive, $|\phi_t\ast f| \leq \phi_t \ast |f|$. Hence we can assume without loss of generality that $f$ is positive.
Let $\lambda_\phi(s)$ for $s \geq 0$ be the set $\{ \phi(x) \geq s\}$. We have that $\phi(x) = \int_0^P \chi_{\lambda_{\phi}(s)}(x) \mathrm{d}s$, where $P = \sup \phi$. Note that by assumption $\lambda_\phi(s)$ for a decreasing family of balls around the origin. Let $\lambda^*\phi(s)$ be the ball whose radius is 1 more than the radius of $\lambda_\phi(s)$. Now, let $|y-x| < 1$. We have that
$$ \int_{\lambda_\phi(s)} f(y-z) \mathrm{d}z \leq \int_{\lambda^*_\phi(s)} f(x - z) \mathrm{d}z $$
since $f$ is positive and $x + \lambda^*_\phi(s) \supset y + \lambda_\phi(s)$.
Let $P'$ be the smallest number such that $\lambda_\phi(P')$ has radius at most 1. By assumption (that $\phi$ is a bounded integrable function) we have that $\lambda_\phi(P')$ has positive radius $R'$.
$$\begin{align} \phi\ast f(y) &= \int_0^{P'} \int_{\lambda_\phi(s)} f(y-z) \mathrm{d}z \mathrm{d}s + \int_{P'}^P \int_{\lambda_\phi(s)} f(y-z) \mathrm{d}z \mathrm{d}s \\
&\leq \int_0^{P'} \int_{\lambda^*_\phi(s)} f(x-z) \mathrm{d}z \mathrm{d}s + \int_{P'}^P\int_{\lambda^*_\phi(P')} f(x-z) \mathrm{d}z \mathrm{d}s \\
& \leq \int_0^{P'} |\lambda_\phi^*(s)| \mathcal{M}f(x) \mathrm{d}s + |P' - P||\lambda_\phi^*(P')|\mathcal{M}f(x)
\end{align}$$
Now using that for $s \leq P'$ we have that $|\lambda_\phi^*(s)| \leq c^n |\lambda_\phi(s)|$ where $c = 1 + 1/R'$ we have
$$ \leq \mathcal{M}f(x)\left( c^n \int_0^P |\lambda_\phi(s)| \mathrm{d}s + |P-P'| |\lambda_\phi^*(P')|\right)\tag{*}$$
The first term inside the parenthesis gives $\int \phi(z) \mathrm{d}z$ from the distributional function characterisation of Lebesgue spaces. The second term is quite obviously a finite constant depending on $\phi$.
Now, if we replace $\phi$ by $\phi_t$ in the above argument, then $P \mapsto t^{-n} P$. We will need to consider $|y-x| < t$ and we let $\lambda^*_{\phi_t}(s)$ to have radius $t$ more than its counterpart without the star. We also let $P'$ be such that the corresponding ball has radius at least $t$: by the scaling property of $\phi_t$ we see that $P' \mapsto t^{-n} P'$, and $\lambda_\phi(P') \to t \lambda_\phi(P')$. So the above analysis goes to show that the constant derived above (the term inside the parentheses in (*)) does not depend on the scaling $t$. Hence we get the desired inequality.
Yes, this does say that $Hf$ is non-integrable, which is why for example $H$ is not bounded in $L^1$.
Now I'll try to explain the proof of the inequality.
Choose a radius $r_0$ such that $\int_{B(0,r_0)} |f(y)|\,dy = c > 0$ (this is possible so long as $f$ is not zero up to a set of measure zero). Then for $|x| > r_0$ we know that
$$
Hf(x) \ge \frac{1}{m(B(x,|x|+r_0)} \int_{B(x,|x|+r_0)} |f(y)|\,dy \gtrsim c (|x|+r_0)^{-n}
$$
because $B(0, r_0) \subset B(x, |x|+r_0)$. Furthermore, when $|x|$ is much larger than $r_0$, $(|x|+r_0)^{-n} \gtrsim |x|^{-n}$.
Best Answer
Wlog, take $f$ nonnegative. Since $f \neq 0$, find an $R > 0$ large enough so that $\int_{B(0;R)} f > 0$. Fix any $x \in \mathbb{R}^d$, and take $S$ large enough so that $B(0;R) \subset B(x;S)$ (say $S = R + |x|$). So $f^*(x) \ge \frac{c}{S^d} \int_{B(x;s)} f \ge \frac{c}{(R+|x|)^d} \int_{B(0;R)} f$, where $c^{-1}$ is the volume of the unit ball in $\mathbb{R}^d$. It should be easy to obtain the last part of the estimate from here.