Question
Let $n\in N$, and let $S$ be a finite set in $Z^n$. What is the dimension $d(n,S)$ of the space $F(n,S)$ of bounded functions $f:\mathbb Z^n\to \mathbb R$ that are discrete harmonic on $Z^n\setminus S$?
Definition
A function $f:\mathbb Z^n\to \mathbb R$ is discrete harmonic iff $$ 2nf(x) = \sum_{i=1}^n [f(x-e_i)+f(x+e_i)], \tag{*}\label{H} $$ for all $x\in \mathbb Z^n$, where $(e_1,\ldots,e_n)$ is the canonical basis of $\mathbb Z^n$.
Note that the condition \eqref{H} can be interpreted as that $f$ is a martingale w.r.t. a symmetric random walk on $\mathbb Z^n$.
The case $n=3, S=\{0\}$
This question was motivated by the question about existence of a bounded non-constant function $f:\mathbb Z^3 \to \mathbb R$ that is discrete harmonic except at the origin $(0,0,0)$. An example of $f\in F(3,\{0\})$ is $f_1(x)$ equal to the probability that a symmetric 3D random walk starting from $x$ hits the origin. It can then be asked if all functions in $F(3,\{0\})$ can be obtained as linear combinations of $f_1(x)$ and the constant function $f_0(x)=1$.
The case $S=\{0\}$
For $n=1$ and $n=2$ the question is closely related to the related question with $S=\emptyset$. A random walk in 1D and 2D is guaranteed to return to the origin. This suggest that $d(1,\{0\})=d(2,\{0\})=1$.
For $n\geq 3$ it is no more guaranteed that the random walk will be absorbed in at the origin and so there is at least one more dimension associated the "value of never returning". Can we conclude that $d(n,\{0\}) = 2$ for $n\geq 3$?
The case $n=1$
Consider $S=\{a,b\}$. If $b=a+1$, then $f(a)=f(a-1)=\ldots$ and $f(b)=f(b+1)=\ldots$, but $f(a),f(b)$ can be arbitrary. What is more, $f(x)$ has to be linear between the points $a$ and $b$. For the general case of $S\subset \mathbb Z$ with $|S|\geq 2$, $f$ is constant below $\min(S)$ and above $\max(S)$, and is linear in between the points of $S$.
I try to answer to this interesting question.
Consider the discrete PDE $$ \left\{\begin{matrix} \Delta f(x) = \rho & \text{ on } G \\ f(x) = F(x)& \text{ on } \partial G \end{matrix}\right. $$ where $G$ is a graph. Then the problem is linear, i.e. if we have solutions for$f_1$ and $f_2$ for data $(\rho_1,F_1)$ and $(\rho_2,F_2)$ then $\alpha f_1 + \beta f_2 $ is the unique solution for $(\alpha \rho_1 + \beta \rho_2, \alpha F_1+ \beta F_2)$. We are interested in solutions of the data $(0,F)$. Which says that $f$ is harmonic on $G$. Let $y \in \partial G$ $$ \left\{\begin{matrix} \Delta f(x) = 0 & \text{ on } G \\ f(x) = \delta_{x=y}& \text{ on } \partial G \end{matrix}\right. $$ the discrete harmonic measure $H(x,y)$ is the unique solutions to this problem. Consider $S \subset \partial G$ and the PDE $$ \left\{\begin{matrix} \Delta f(x) = 0 & \text{ on } G \\ f(x) = \delta_{x \in S}& \text{ on } \partial G \end{matrix}\right. $$ we have that $H(x,S)= \sum_{ y \in S} H(x,y)$ is the unique solutions by linearity. In particular every solution $f$ to the general equation with data $(0,F)$ is of the form $$ f(x) = \sum_{y \in \partial G} H(x,y) F(y).$$
If we want to give probabilistic interpretation to this, if we consider $X:=(X_n)_{n \geq 0}$ to be a simple random walk on $\mathbb{Z}^d$ which start at $x $ and let $\tau $ to be the first time that $X$ visits the boundary then $H(x,S)= \mathbb{P}_x(X_{\tau} \in S)$, i.e. is the probability that the simple random walk starting at $x$ exits in $S$. And for the general solution is $$f(x) = \mathbb{E}[ F(X_{\tau}) ].$$
Consider the case $G= \mathbb{Z}^n \setminus S$ and $\partial G = S$, where $S$ is finite. For $n=1,2$ we have that simple random walks are recurrent, and so $f$ is fully determined by its values on $S$. Hence $d(1,S)=d(2,S) = |S|$. In contrast, for $n\geq 3$ there is a positive probability that a given random walk will never hit any of the points in $S$, and thus we can also choose a constant associate with the random walk never hitting $S$, and so $d(n,S)=|S|+1$.
Conclusion: $$d(n,S) = \begin{cases} |S| & \text{if } \ n\leq 2 \\ |S|+1 & \text{otherwise}. \end{cases}$$