You can, in fact, do better as the answer to this Physics.SE question of yours indicates, assuming you're after an asymptotic expression. I haven't read the book you mention, but if you're familiar with the bra-ket notation, I recommend the discussion of variational methods in Sakurai's Modern Quantum Mechanics. The following will be based on Sakurai's approach, but adapted to address the questions I think you're asking.
Let $\{|\psi_n\rangle\}$ be an orthonormal basis, where $|\psi_0\rangle$ is the true ground state. We consider a trial wave function (ket) $|\phi\rangle=|\psi_0\rangle + |\delta \phi\rangle$, and express it in the energy eigenstate basis,
$$
|\phi\rangle=\sum_{n=0}^{N} c_n|\psi_n\rangle \tag{1}
$$
where $\sum_{n=0}^{N} |c_n|^2=1$, such that $\delta \phi= \sum_{n=1}^N c_n |\psi_n\rangle$. Here $N$ is the dimension of the Hilbert space. You may take $N\rightarrow \infty$ as in your question, but a finite dimension allows a clearer discussion of the difference between the wave functions.
The energy in the trial wave function is (since $|\phi\rangle$ is normalized)
$$\begin{align}\tag{2}
E &= \langle \phi | \hat{H} | \phi\rangle = \left( \langle \psi_0| + \langle \delta\phi| \right) \hat{H} \left( |\psi_0\rangle + | \delta \phi\rangle\right)\\\tag{3}
&= \langle \psi_0| \hat{H} |\psi_0\rangle + \langle \psi_0| \hat{H} | \delta \phi\rangle + \langle \delta\phi| \hat{H} |\psi_0\rangle + \langle \delta\phi| \hat{H} | \delta \phi\rangle\\
&= E_0 + \langle \delta\phi| \hat{H} | \delta \phi\rangle, \tag{4}
\end{align}
$$
where we used that $|\delta \phi\rangle$ is orthogonal to $|\psi_0\rangle$. Thus the difference between the approximate and exact ground state energies can be written
$$
\delta E = E -E_0 = \langle \delta\phi| \hat{H} | \delta \phi\rangle = \sum_{n=1}^N E_n |c_n|^2. \tag{5}
$$
Now, if we think about $|\phi\rangle$ as a (long?) vector, where each element corresponds to $c_n$ ($n=0,1,2,\dots$), the difference between $|\psi\rangle$ and the exact ground state can be written
$$
|\delta\phi\rangle = |\phi\rangle - |\psi_0\rangle = (0, c_1, c_2, c_3, \dots, c_N). \tag{6}
$$
You now have to decide how you want to define the meaning of "$\phi$ and $\psi_0$ differ by at most $O(\epsilon)$". (There are many distance metrics for Hilbert spaces.) One simple way is to directly say that the largest element ($c_{n_0}$, say) of the $|\delta \phi\rangle$ vector is $O(\epsilon)$. Alternatively, we can require
$$
\sqrt{\langle \delta \phi | \delta \phi \rangle} = \sqrt{\sum_{n=1}^N |c_n|^2} = \sqrt{1-|c_0|^2} = \mathcal{O} (\epsilon), \tag{7}
$$
which is achieved if $|c_n|\lesssim \mathcal{O}(\epsilon)\quad \forall n>0$. In both cases we find
$$
\delta E = \sum_{n=1}^N E_n |c_n|^2 \lesssim \sum_{n=1}^N E_n \mathcal{O}\left(\epsilon^2\right). \tag{8}
$$
If the spectrum of $\hat{H}$ is bounded (i.e. each $E_n$ is finite), $E_n$ can be treated as a constant that is unimportant for asymptotic properties, yielding
$$
\delta E \sim \sum_{n=1}^N \mathcal{O}\left(\epsilon^2\right) = \mathcal{O}\left(\epsilon^2\right). \tag{9}
$$
Thus, an error $\mathcal{O}(\epsilon)$ in the trial wave function leads to an energy error that scales as $\mathcal{O}(\epsilon^2)$. Note, however, that this expression is not itself useful to estimate the magnitude of the energy difference, which depends non-universally on $E_n$. A trivial way to see this is to simply multiply $\hat{H}$ by some large number, which will affect the energy difference in Eq. (5), but not the difference between states in Eq. (6). Neither will it change the error in Eq. (9). If you need the magnitude of $\delta E$, you'll want to know $E_0$ independently, e.g. from experiment or some other theoretical technique (but if you have another theoretical technique, using a variational one might be unnecessary). The real utility of the variational method comes from the fact that you can just guess trial wave functions until you're tired, pick the one with lowest energy $E$, and be sure that it's the best approximation of $|\psi_0\rangle$ that you've tried so far.
To avoid confusion, perhaps it's worth mentioning that the overall units of Eqs. (8,9) are the same. In (9) I've absorbed the units of energy into the $\mathcal{O}$, whereas $\mathcal{O}$ is dimensionless in (8). This is common practice, since $\mathcal{O}$ is used as a shorthand.