Harmonic functions on $\mathbf{Z}^2$ 
Problem 1: Find all functions $f:\mathbf{Z}^2 \to \mathbf{R}$ which are harmonic in the sense that $$f(x,y) = \frac{f(x+1,y) + f(x-1,y) + f(x,y+1) + f(x,y-1)}{4}$$for all $(x,y)\in\mathbf{Z}^2$, and which are also Lipschitz in the sense that the gradients $$f(x+1,y)-f(x,y)\\f(x-1,y)-f(x,y)\\f(x,y+1)-f(x,y)\\f(x,y-1)-f(x,y)$$are all globally bounded.

Obvious examples: $f = 1$, $f=x$, $f=y$, and linear combinations of these. Is this all?

Problem 2: What further examples do we get if we weaken the Lipschitz condition, say by allowing the gradients to grow at most linearly (with respect to distance from $(0,0)$)?
Problem 3: How much do the character of our examples change if we replace the generating set $S = \lbrace (1,0),(-1,0),(0,1),(0,-1)\rbrace$ by another (symmetric) set $S$ which generates $\mathbf{Z}^2$? For example, what if we require that $$f(x) = \frac{1}{|S|}\sum_{s\in S}f(x+s),$$for, say, $S = \lbrace s: \|s\|_2\leq 100\rbrace$? Does the dimension of the space of Lipschitz harmonic functions change?

[Background: I'm trying to understand Kleiner's theorem, which states that if a finitely generated group $G$ has polynomial growth then the space of Lipschitz harmonic functions on $G$ has finite dimension. The simplest example $G=\mathbf{Z}$ is pretty simple, but the second simplest example $G=\mathbf{Z}^2$ already seems nonobvious to me.]
 A: Lemma. Let $f$ be a bounded harmonic function on $\mathbb Z^n$ with respect to a  random walk that is not restricted to any proper subset of the grid. Then $f$ is constant.
Proof (from Principles of Random Walk by F. Spitzer). If $f$ is nonconstant, then the function $g(x)=f(x+a)-f(x)$ has positive supremum $M$ for some $a$. Let $x_n$ be a sequence such that $g(x_n)\to M$. Let $g_n(x)=g(x+x_n)=f(x+a+x_n)-f(x+x_n)$. Using the Cantor diagonal argument, pick a subsequence of $g_n$ that converges at every point of the grid. Let $h$ be the limit of this subsequence. Since $h$ is harmonic and attains its maximum at $0$, we have $h\equiv M$. Due to pointwise convergence, for any positive integer $N$ there exists $n$ such that $g_n>M/2$ at the points $0,a, 2a, \dots, Na$. It follows that $f(x_n+Na)-f(x_n)=\sum_{k=0}^{N-1} g_n(ka)>MN/2$ which  contradicts the boundedness of $f$. $\Box$
Another, more probabilistic proof is here, but it uses the recurrence of random walk and therefore works only in two dimensions (not for general dimension as claimed there). 
Answer 1: If a harmonic function on $\mathbb Z^2$ is Lipschitz, then it's of the form $f(x,y)=ax+by+c$. Indeed, $g(x,y)=f(x+1,y)-f(x,y)$ is bounded and harmonic, therefore constant by the Lemma. Similarly, $f(x,y+1)-f(x,y)$ is constant and thus $f$ is linear. 
Partial Answer 2: Linear growth on derivatives allows for 2nd degree harmonic polynomials $f(x,y)=xy$ and $f(x,y)=x^2-y^2$. I think these are all (i.e., the space is 5-dimensional) but don't have a proof. Higher-order polynomial bounds will allow for higher degree polynomials, which are similar, but not identical to the harmonic polynomials on $\mathbb R^2$: see the expository article Discrete analytic functions by Lovász. 
Answer 3: The space will not change, because the proof from Answer 1 applies here as well.
A: Answer to Problem 2:
@user31373 already answered Problem 1, and the same argument can be applied to Problem 2.
If the gradient $g(x,y)=f(x+1,y)-f(x,y)$ can grow at most linearly then the second-order gradient $g^{(2)}(x,y)=g(x+1,y)-g(x,y)$ is bounded. Since $g^{(2)}$ is harmonic, thus by the @user31373's lemma (discussed in this article) $g^{(2)}(x,y)$ is constant. We conclude that $g(x,y)$ whose all second-order gradients are bounded is a polynomial of at most second order.
We can inductively show that if $f(x,y)$ is harmonic, then all of its (mixed) gradients are harmonic as well, and we conclude that:

If all the $n$-th gradients of a discrete harmonic function $f(x,y)$ are bounded, then $f(x,y)$ is a polynomial of at most $n$-th order.

