Error on Wikipedia: Nelson's proof of Liouville's theorem works only for bounded modulus?

On Wikipedia, it is stated:

If $f$ is a harmonic function defined on all of $\mathbb{R}^n$ which is bounded above or bounded below, then $f$ is constant...Edward Nelson gave a particularly short proof of this theorem, using the mean value property mentioned above:

Given two points, choose two balls with the given points as centers and of equal radius. If the radius is large enough, the two balls will coincide except for an arbitrarily small proportion of their volume. Since f is bounded, the averages of it over the two balls are arbitrarily close, and so f assumes the same value at any two points.

But this is an error, right? Because Nelson's proof works if the harmonic function is bounded in modulus, not just bounded above or below. What is a simple way to exetend the proof to the case that the function is bounded above or below?

If such doesn't exist, where can I find a different proof of this fact?

I guess one way (for $n=2$) might be to invoke complex analysis, create $f$ entire such that $\text{Re}(f)$ is our harmonic function, and then reason about $e^{f}$.

I think you are correct that the Wikipedia article's explanation only really works when $f$ is bounded in modulus; however, a slight change adapts the proof to the case when $f$ is bounded below or above only. The proof below can be found on page 45 of Harmonic Function Theory by Axler et al., for instance. I give another solution at the bottom that works when the dimension is equal to $2$.

Suppose $f$ is harmonic on $\mathbb R^n$ and bounded above or below. By replacing $f$ with $-f$ if necessary, we can assume that $f$ is bounded below, and then by adding a suitable constant we can assume that $f\geq 0$. It is clear that $|f(x) - f(y)|$ is bounded by the integral of $f$ over the symmetric difference of the balls $B(x,r)$ and $B(y,r)$ divided by the volume of a ball of radius $r$. But because $f$ is nonnegative, the integral of $f$ over the symmetric difference of $B(x,r)$ and $B(y,r)$ is bounded by the integral of $f$ over the annulus $B(x,r+|x-y|)\setminus B(x,r-|x-y|)$ (this is easiest to see if you draw a picture). Thus \begin{align*} |f(x) - f(y)| & \leq {1\over|B(0,r)|}\int_{B(x,r+|x-y|)\setminus B(y,r-|x-y|)}f\\ &= {1\over |B(0,r)|}\int_{B(x,r+|x-y|)}f - {1\over |B(0,r)|}\int_{B(x,r-|x-y|)}f \\[0.2em] & = {(r+|x-y|)^n - (r-|x-y|)^n\over r^n}f(x). \end{align*} The expression above is $O(r^{-1})$, so making $r\to\infty$ gives $f(x) = f(y)$.

Here's another proof for dimension $n=2$. Again suppose $f\geq0$ is harmonic. For $\epsilon>0$, put $f_\epsilon(z) = f(z) + \epsilon\log{|z|}$. Let $m = \inf{f}$. Now $f_\epsilon$ is harmonic on $|z|> 1$ and continuous on the boundary $|z| = 1$. Moreover, since $f\geq 0$, we know that $f_\epsilon(z)\geq \epsilon\log{|z|} \to\infty$ with $|z|$. Thus for all $\epsilon>0$, the maximum principle tells us that $\inf{f_\epsilon}$ is attained on the unit circle $|z|=1$. Since $\log{|z|} = 0$ on the circle, it follows that this infimum is independent of $\epsilon$; call it $k$, and note that $k>m$ because $f$ attains the value $k$ somewhere on the unit circle, while $f$ cannot attain $m$ anywhere on the plane by the maximum principle. Then for each $z$ with $|z|\geq 1$, we have $f_\epsilon(z)\geq k >m$ for all $\epsilon$. It follows by making $\epsilon\to0$ that $f(z) \geq k$. But if $f(z)\geq k$ whenever $|z|\geq1$, then $f(z)\geq k$ for $|z| < 1$ as well by the maximum principle. But now we have found that $f(z)\geq k>m$ for all $z$, contradicting the definition of $m$ as the infimum of $f$.