Fix $n\in \mathbb N$. Let $F(n)$ be the space of discrete harmonic functions $f:\mathbb Z^n\to \mathbb R$. What is the minimal $d \in \mathbb N$ (if there is one) such that there exists $B\subset F(n)$ with $|B|=d$ that every $f\in F(n)$ can be written as a linear combination of functions from $B$ "subject to a translation", i.e. there is $m\in \mathbb N$ and $(a^1,f^1,x^1),\ldots,(a^m,f^m,x^m) \in \mathbb R \times B \times \mathbb Z^n$ such that $$f(x) = \sum_{j=1}^m a^j f^m(x-x^m) \quad \forall x\in \mathbb Z^n. \tag{*}$$
A function $f:\mathbb Z^n\to \mathbb R$ is harmonic iff $$ 2nf(x) = \sum_{i=1}^n [f(x-e_i)+f(x+e_i)], $$ for all $x\in \mathbb Z^n$, where $(e_1,\ldots,e_n)$ is the canonical basis of $\mathbb Z^n$. Note that $f$ being discrete harmonic can be interpreted as that $f$ is a martingale w.r.t. a symmetric random walk on $\mathbb Z^n$.
Case $n=1$: It is easy to see that $f(x)$ has to have a constant growth, thus $F(1)$ is the space of all the linear functions $f:\mathbb Z\to \mathbb R$.
Case $n=2$: Analogously, $F(2)$ also contains all the linear functions $f:\mathbb Z^2\to \mathbb R$. However, the function $f^2_{12}(x) = (x_1)^2 - (x_2)^2$ is also discrete harmonic. Does $f^2_{12}$ together with the basis of the space of linear functions $f:\mathbb Z^2\to \mathbb R$ span the whole $F(2)$? Clearly not, as @NeitherNor pointed out, there are other linearly independent functions in $F(2)$ and as @RyszardSzwarc brought up, their translations are also in $F(2)$, thus the complicated definition of "basis" $B$.
Case $n\geq 3$: Besides the linear functions also the functions $f^3_{ij}=(x_i)^2 - (x_j^2), i\neq j, \ i,j \in \{1,2,3\}$ are harmonic...
The answer is : the dimension $d$ is infinite for $n\ge 2.$
At the beginning I will put the problem into a more general framework of finitely generated groups, although it will not help solving it. Let $G=\mathbb{Z}^n.$ The group is generated by $e_1,e_2,\ldots, e_n.$ The generators lead to a natural distance on $G$ $$\|x\|=|x_1|+|x_2|+\ldots +|x_n|$$ Consider the operator acting on functions defined on $G$ by $$(Af)(x)={1\over 2n}\sum_{j=1}^n[f(x-e_j)+f(x+e_j)]={1\over 2n}\sum_{\|y\|=1} f(x+y)$$ Observe that $A$ is a convolution operator on the group $G$ with the function $$\varphi={1\over 2n}\sum_{j=1}^n[\delta_{e_j}+\delta_{-e_j}]={1\over 2n}\sum_{\|y\|=1}\delta_y$$ i.e. the normalized indicator function of elements with norm $1.$
Concerning the problem, we want to determine the dimension of the eigenspace of $A$ corresponding to the eigenvalue $1,$ i.e. $Af=f.$ We do not impose any norm conditions on $f.$ Some solutions are multiplicative (called usually the characters of the group, if $|f|=1$). Namely if $f(a+b)=f(a)f(b)$ then $$f(x)=(Af)(x)={1\over 2n}\sum_{j=1}^n[f(x)f(e_j)+f(x)f(e_j)^{-1}]\\ = {1\over 2n}\sum_{j=1}^n [f(e_j)+f(e_j)^{-1}]\,f(x)$$ This provides a solution (found by @reuns) $${1\over 2n}\sum_{j=1}^n[z_j+z_j^{-1}]=1,\quad z_j:=f(e_j)$$
We are going to show that the dimension of harmonic functions is infinite. Consider $n=2.$ For $z>1 $ let $f_z$ denote the function corresponding to $z_2=-z$ and $z_1>1,$ i.e. $$z_1=2+{1\over 2}(z+z^{-1})+\sqrt{[2+\textstyle{1\over 2}(z+z^{-1})]^2-1}\ge 1+(z+z^{-1})>z=|z_2|$$ Observe that for $w>z>1 $ we have $w_1>z_1>1.$ Let $z(k)>z(k-1)>\ldots >z(1)>1.$ We claim that the functions $f_{z(j)}$ are linearly independent. The corresponding parameters are $z_1(j)$ and $z_2(j).$ Assume $$a_1f_{z(1)}+a_2f_{z(2)}+\ldots +a_kf_{z(k)}=0$$ With no loss of generality we may assume $a_k\neq 0.$ Substituting $x=je_1$ gives $$a_1z_1(1)^j+a_2z_1(2)^j+\ldots +a_kz_1(k)^j=0$$ The left hand side is equal asymptotically to $a_kz_1(k)^j$ when $j\to \infty.$ The translation is irrelevant as the translation of $f_z$ is a multiple of $f_z.$ Therefore $a_k=0,$ a contradiction. The proof is made for $n=2,$ but it can be adapted to any $n\ge 2.$