How accurate is $\prod_{p < y} \left(1-\frac{1}{p}\right)$ for the number of $y$-rough integers? Call a positive integer $y$-rough if it has no factors below $y$.  I am interested in accurate ways to estimate the number of positive integers less than or equal to $x$ that are $y$-rough. Assuming that $x > y^2$, one option is first to compute:
$$P=\prod_{p < y} \left(1-\frac{1}{p}\right)$$
where the $p$ are primes.
We can then give $xP$ as the estimate.
What can be said about how accurate this estimate is?

I don't know how much analytic number theory you know, so if this is way past your pay grade then let me know and perhaps there's an elementary method I can find. This is just the general method that almost always works with this sort of problem.
Suppose we want to count the amount of $y$-rough numbers. Let $i_n$ be the indicator function for $y$-rough numbers, i.e,
$$i_n=
\begin{cases}
1 & n \text{ is }y\text{-rough}\\
0 & \text{otherwise}.
\end{cases}$$
We construct the Dirichlet series
$$f(s)=\sum_{n=1}^{\infty}\frac{i_n}{n^s}.$$
It is clear that $nm$ is $y$-rough if and only if both $n$ and $m$ are $y$-rough, and hence $i_{nm}=i_{n}i_{m}$. Thus, we have an Euler product expansion
\begin{align*}
f(s)&=\prod_{p}\left(1+\frac{i_p}{p^s}+\frac{i_{p^2}}{p^{2s}}...\right)\\
&=\prod_{p\geq y}\left(1+\frac{1}{p^s}+\frac{1}{p^{2s}}...\right)\\
&=\prod_{p\geq y}\frac{1}{1-p^{-s}}\\
&=\zeta(s)\prod_{p<y}(1-p^{-s}).
\end{align*}
In particular, since $\lim_{s\to 1}(s-1)\zeta(s)=1$, we have
$$\lim_{s\to 1}(s-1)f(s)=\prod_{p<y}\left(1-\frac{1}{p}\right).$$
The Hardy-Littlewood Tauberian Theore gives us thus that
$$\lim_{x\to\infty}\frac{1}{x}\sum_{n<x}i_{n}=\prod_{p<y}\left(1-\frac{1}{p}\right).$$
This is exactly the statement that your approximation is correct, i.e, the proportion of $y$-rough numbers tends to $\prod_{p<y}\left(1-\frac{1}{p}\right)$. To get uniform error bounds, you could use stronger analytic number theory results. For instance, Perron's formula says that
$$\sum_{n<x}i_n=\frac{1}{2\pi i}\int_{c-i\infty}^{c+i\infty}f(z)\frac{x^z}{z}dz$$
for any $c>1$. Now, $f(s)$ has no poles apart from the one at $s=1$ and hence using Cauchy's integral formula we can push $c$ as close as we want to 0. Hence, we get that
$$\left(\text{proportion of }y\text{-rough numbers less than }x \right)=\prod_{p<y}\left(1-\frac{1}{p}\right)+O(1/(x^{1-\epsilon}))$$
for every $\epsilon>0$. The constant in the big-Oh depends on $\epsilon$ and $y$.
