Find all polynomials $p: \mathbb{Z} \to \mathbb{Z}$ with $p(a)p(b) \in p(\mathbb{Z})$ for $a,b\in\mathbb Z$.

Edited. Now we present a complete solution to the original problem.
Edited again. Thanks for @lhf, who pointed out a flaw in the argument. The original answer missed some solutions by claiming that $a\mid b(b-1)\implies a\mid b$ or $a\mid (b-1)$, which is incorrect. Fortunately, that wasn't a serious problem, since we may just write the solutions with the restriction $a\mid b(b-1)$ instead of $a\mid b$ or $a\mid(b-1)$.

Answer. All possible $p\in\mathbb Q[x]$, with $p(\mathbb Z)\subseteq\mathbb Z$ and satisfying the given condition, can be expressed in one of the following form, with $a,b,n,k$ integers, $a\not=0,n>0$:
$$p(x)\equiv 0 \text{ or } p(x)\equiv 1$$ $$p(x)=(ax+b)^n,\quad \text{with }a\mid b(b-1)$$ $$p(x)=(ax+b)^{2n},\quad \text{with }a\mid b(b+1)$$

Let $n=\deg p$. We first prove that $p(x)$ must be of the form $(ax+b)^n$.

Claim. $p(k)$ is a perfect $n$-th power for any $k\in\mathbb Z$.

Proof. Suppose $p(k)\not=0$. For each $a$ there exists some $c(a)$ such that $p(k)=p(c(a))/(p(a)$. We may observe that for sufficiently large $a$, we have $|c(a)|\ge (1-\epsilon)|a|$ since $|p(k)|\ge 1$, hence $c(a)\to\infty$ as $a\to\infty$, and $(c(a)/a)^n\to p(c(a))/p(a)=p(k)$. Now suppose $\lambda=\sqrt[n]{p(k)}\not=0$ is irrational. First we assume $n$ is odd, so that we have $c/a\to\lambda$. We present some estimation $$\left|\left(\frac c a\right)^n-\frac{p(c)}{p(a)}\right|=\left|\left(\frac c a\right)^{n-1}\frac{a_{n-1}}{a_n}\left(\frac c a-1\right)\frac 1 a\right|+\mathcal O(a^{-2})$$ ($a_n,a_{n-1}$ is the coefficient of $x^n,x^{n-1}$ in $p(x)$, respectively) Denote by $\mu$ the ratio $a_{n-1}/a_n$, we have $$a\left|\left(\frac c a\right)^n-\frac{p(c)}{p(a)}\right|\to|\lambda|^{n-1}|\mu||\lambda-1|$$ Hence $$a\left|\left(\frac c a\right)^n-\lambda^n\right|=a|\lambda|^n\left|1-\left(1-\frac{\lambda-c/a}\lambda\right)^n\right|\to|\lambda|^{n-1}|\mu||\lambda-1|$$ Now note that $$\left|1-\left(1-\frac{\lambda-c/a}\lambda\right)^n\right|= \frac{n|\lambda-c/a|}{\lambda}+\mathcal O((\lambda-c/a)^2)$$ We conclude that $$|a\lambda-c|=|\mu||\lambda-1|/n+\mathcal O(a^{-1})$$ For the case where $n$ is even, the reasoning is not quite different: the negative part of $c/a$ tends to $-\lambda$ and the positive part tends to $\lambda$. The estimations remain almost the same except for some sign change. In other words, we have $$|a^-\lambda+c(a^-)|\to|\mu||\lambda+1|/n$$ $$|a^+\lambda-c(a^+)|\to|\mu||\lambda-1|/n$$ Here $\{a^+\}$ denotes the subset of $\mathbb Z$ where $c/a\ge0$, and the similar for $\{a^-\}$. The two sets form a partition of $\mathbb Z$. Obviously those limit properties are impossible for irrational $\lambda$ (by considering the fractional part of $a\lambda$), whereby we proved the rationality of $\lambda$. Since $p(k)\in\mathbb Z$, we may assert that $\lambda\in\mathbb{Z}$.

Now, as @u1571372 kindly pointed out in the comment, known results have shown that $p$ must be of the form $(ax+b)^n$(note that int the MO post it's assumed that $p(x)\in\mathbb Z[x]$, to which we may reduce the original problem by multiplying a factor $d^n$). In fact this can be done by more elementary method in our situation, where more restrictions are put on $p$, instead of the somehow overkill methods mentioned in the MO post. The main point is as following: since $|a^+\lambda-c(a^+)|$ is a integer sequence tending to a limit, for sufficiently large $a^+$ the expression is a constant. So for some $l$ the equality $p(\lambda x+l)=\lambda^n p(x)$ holds for infinitely many $x$. But it's a polynomial equation in $x$, so in fact it's an identity. Now it's not difficult to arrive at the conclusion by comparing coefficients(or more smartly, comparing derivatives of all orders).

At last we conclude the solution by considering all polynomials of the form $(ax+b)^n$. It's plain that $a,b\in\mathbb Z$ since $p(\mathbb Z)\subseteq\mathbb Z$. Note that $\sqrt[n]{p(u)p(v)}=a^2 uv+ab(u+v)+b^2$. When $n$ is odd, $(aw+b)^n=p(u)p(v)$ implies $aw+b=a^2 uv+ab(u+v)+b^2$, so the existence of $w$ is reduced to $b^2\equiv b\pmod{|a|}$, which is equivalent to $a\mid b(b-1)$. When $n$ is even, more possibility is allowed: we may have $aw+b=-(a^2 uv+ab(u+v)+b^2)$. It turns out that it's equivalent to $a\mid b(b-1)$ or $a\mid b(b+1)$ in this case. Now we're done.