An interesting integral $\cos(x)\cos(x^2)\cos(x^3)...$
The upper bound is not horribly difficult to prove.
By the Weierstrass product for the cosine function we have
$$ -\log\cos x = \sum_{n\geq 1}\frac{(4^n-1)\zeta(2n)}{n \pi^{2n}}x^{2n} \tag{1}$$ hence $$ \sum_{m\geq 1} -\log\cos x^m = \sum_{n\geq 1}x^{2n}\sum_{d\mid n}\frac{(4^d-1)\zeta(2d)}{d\pi^{2d}}=\sum_{n\geq 1} c_n x^{2n}\tag{2}$$ while $$ -\log\sqrt{1-x^2} = \sum_{n\geq 1}\frac{1}{2n} x^{2n}\tag{3} $$ hence it is enough to show $c_n\geq \frac{1}{2n}$, which is trivial since $$ c_n\geq \frac{(4^1-1)\zeta(2\cdot 1)}{1\cdot \pi^{2\cdot 1}}=\frac{1}{2}.\tag{4}$$ Indeed, this inequality proves $$\forall x\in(-1,1),\qquad \prod_{m\geq 1}\cos(x^m)\approx\exp\left(\frac{x^2}{2(x^2-1)}\right) $$ such that the given integral turns out to be pretty close to
$$ 2\int_{0}^{1}\exp\left(\frac{x^2}{2(x^2-1)}\right)\,dx = \int_{0}^{+\infty}\frac{e^{-z/2}\,dz}{(z+1)\sqrt{z^2+z}}=\sqrt{\pi}\,U\left(\tfrac{1}{2},0,\tfrac{1}{2}\right)\approx \sqrt{2}$$ where $U$ is Tricomi's confluent hypergeometric function.
The infinite product under the integral converges extremely fast for $|x|<1$, since we have:
$$\cos x^n=1-\frac{x^{2n}}{2}+\frac{x^{4n}}{24}- \dots$$
So the following estimation for the infinite product can be made:
$$P(x)=\prod_{n=1}^\infty \cos x^n=\prod_{n=1}^N \cos x^n+\text{o} (x^{2N}) \tag{1}$$
To me it's quite interesting to prove that for $|x|>1$ the product approaches $0$. I think the probability argument can be made, since for large $n$ the factors start to behave like random variables evenly distributed on $[-1,1]$. Thus, their expected value would be $0$, and so $P(x)$ will vanish.
As the function $P(x)$ is continuous (and likely infinitely differentiable, though I'm not sure how to prove that), and we have $P(0)=1$ and $P(|x|>1)=0$, the argument can be made that only $|x|<1$ contribute to the value of the integral.
Thus, the Taylor series around $x=0$ should give a good approximation to $P(x)$ and integrating it by term we should obtain a good estimation of the integral.
As (1) shows we need to only consider $N$ factors in the product to obtain Taylor series up to $x^{2N}$.
For $N=6$ (with a little help of Wolfram Alpha) we obtain:
$$P(x)=1 - \frac{x^2}{2} - \frac{11 x^4}{24} - \frac{181 x^6}{720} - \frac{9239 x^8}{40320} - \frac{148681 x^{10}}{3628800} - \frac{49402979 x^{12}}{479001600}+\text{o} (x^{12})$$
The fact that all the terms except the first are negative may just be a coincidence (or not, see Jack D'Aurizio's answer).
Integrating the above expression by term gives us a lower bound on the integral:
$$2 \int_0^1 P(x) dx > \frac{832721237}{622702080}=1.3372706848835321 \dots$$
I believe that leaving our some of the terms actually gives us a better bound, but without careful consideration, and getting a good numerical value of the integral first, we might as well take only the first two terms:
$$2 \int_0^1 P(x) dx \approx 2-\frac{1}{3}=\frac{5}{3}=1.66666\dots$$
Which is poor estimation, and actually from above, not from below.
I'll try to check with other methods and maybe get a good numerical value with Mathematica later (if nobody else provides it).
Update
Even without Mathematica, I have been able to get a good numerical estimation, by doing quadratic spline interpolation on an equal grid with step size $0.1$. For the last interval I used a cubic spline because of the bend.
Here's the interpolation result (right side) compared to the truncated product from the OP (left side):
$$2 \int_0^1 P(x) dx \approx 1.4002420432 $$
This value agrees very well with the one provided by Jack D'Aurizio in the comments below.
I used WolframAlpha to compute the product for $9$ points and also the set the splines derivatives to $0$ at $x=0$ and $x=1$, as well as the usual conditions on the equality of splines and their derivatives at grid points.
Here are the splines explicitly.
-0.50460870178410*x^2+1
-0.56117523777048*(x-1/10)^2-0.10092174035682*(x-1/10)+0.99495391298216
-0.68453974237033*(x-2/10)^2-0.21315678791092*(x-2/10)+0.97924998656877
-0.90014478829508*(x-3/10)^2-0.35006473638498*(x-3/10)+0.95108891035398
-1.26006346202322*(x-4/10)^2-0.53009369404400*(x-4/10)+0.90708098883253
-1.86364442303217*(x-5/10)^2-0.78210638644864*(x-5/10)+0.84147098480790
-2.87493565303363*(x-6/10)^2-1.15483527105508*(x-6/10)+0.74462390193271
-4.30929593526087*(x-7/10)^2-1.72982240166180*(x-7/10)+0.60039101829687
-2.71129664596445*(x-8/10)^2-2.59168158871398*(x-8/10)+0.38431581877808
-117.324704896617*(x-9/10)^3+33.2684103240269*(x-9/10)^2-3.13394091790687*(x-9/10)+0.09803469344703
Using the numerical value, we find that the closest approximation is obtained by truncating the Taylor series at:
$$P(x) \approx 1 - \frac{x^2}{2} - \frac{11 x^4}{24} - \frac{181 x^6}{720}$$
Which gives for the integral:
$$2 \int_0^1 P(x) dx \approx \frac{3557}{2520}=1.4115\dots$$