Suppose that $f:B(0,1)\setminus\left\{0\right\}\subset \mathbb{R}^n \to \mathbb{R}$ is a harmonic function. Clearly, the property that $\overline{f(B(0,\epsilon))}=\mathbb{R}$ for all $\epsilon>0$ is equivalent to $\limsup_{x \to 0 } f(x) = +\infty$ and $\liminf_{x \to 0 } f(x) = -\infty$ (through the intermediate value theorem).
Therefore, after considering $f$ on a smaller ball, adding a constant and replacing $f$ by $-f$ if necessary, my question is wether a harmonic function $f:B(0,1)\setminus\left\{0\right\}\subset \mathbb{R}^n \to [0,+\infty)$ is necessarily of the form \begin{equation} f(x)=h(x)+\frac{C}{|x|^{n-2}} \end{equation} where $h$ is a harmonic function on the entire ball $B(0,1)$ (this formula needs an appropriate modification when $n=2$ of course.).
NOTE: Using the Kelvin transform, we can restate this problem as saying that for a harmonic function $f:\mathbb{R}^n\setminus \overline{B(0,1)} \to [0,\infty)$ we have that $f$ converges to some $C$ at infinity and $|f(x)-C|=O(|x|^{...})$. Maybe the mean value property can do the trick from this point in some way?
