Let us clarify these two issues first. The expansion of the exponential is:
$$e^x = 1 + x + \mathcal{O}(x^2)\tag{1}\label{eq:exponential}$$
Second let us remind ourselves about integration by parts:
$$\int \mathrm{d}^4\, x\; g(x)\square f(x)  = \int \mathrm{d}^4\,x\; g(x)\partial_\mu\partial^\mu f = \sum_\mu g(x)\partial^\mu f(x)\bigg|_{\partial\mathbb{R}^4}-\int \mathrm{d}^4x\; \partial_\mu g(x)\partial^\mu f(x)\tag{2}\label{eq:intParts}$$
As it is common in this topic, fields are assume to vanish at infinity fast enough that boundary terms don't contribute. Now we can use \eqref{eq:intParts} two times, to place all the differential operators to act on $\phi$ instead of $\epsilon$, so that:
$$-\frac{1}{2}\int\mathrm{d}^4 x\;\phi(x)\square\epsilon(x) = \frac{1}{2}\int\mathrm{d}^4x\;\partial_\mu\epsilon(x)\partial^\mu\phi(x) =-\frac{1}{2}\int\mathrm{d}^4x\;\epsilon(x)\square\phi(x).\tag{3}\label{eq:epsParts}$$

Starting with the third equation of the OP's question can use first equation \eqref{eq:exponential} to deal with the second factor:
$$\langle \hat{\phi}(x)\rangle = \frac{1}{Z[0]}\int\mathcal{D}\phi e^{iS[\phi]}\left(1-\frac{i}{2}\int\mathrm{d}^4z\;(\phi(z)\square\epsilon(z)+\epsilon(z)\square\phi(z)) +\mathcal{O}(\epsilon^2)\right)\left[\phi(x) + \epsilon(x) \right] $$
where just the exponential having the $\epsilon$ field was expanded and $S$ stands for the kinetic term action of $\phi$. Now we can use the integration by parts as shown in \eqref{eq:epsParts}:
$$\langle \hat{\phi}(x)\rangle = \frac{1}{Z[0]}\int\mathcal{D}\phi e^{iS[\phi]}\left(1-i\int\mathrm{d}^4z\;\epsilon(z)\square\phi(z)+\mathcal{O}(\epsilon^2)\right)\left[\phi(x) + \epsilon(x) \right]$$
Keeping only terms up to linear order in $\epsilon$ one obtains Schwartz's result:
$$\langle \hat{\phi}(x)\rangle = \frac{1}{Z[0]}\int\mathcal{D}\phi e^{iS[\phi]}\left(\phi(x) + \epsilon(x) -i\phi(x)\int\mathrm{d}^4z\;\epsilon(z)\square\phi(z)\right)$$