Continued Fraction and Random Variable

This isn't a complete solution, but hopefully it will give you some useful information on the problem. I'll work in the case where the $X_i$ are iid; i.e., $\epsilon_i=\epsilon$ for all $i$.

First of all, setting aside the probability $p$ for the moment, note that this question is really about compositions of two maps $M_1: t\mapsto1+\frac1t$ and $M_\epsilon: t\mapsto\epsilon(1+\frac1t)$. Since $\frac1t$ is decreasing on $t\gt 0$, so are both of the maps $M_1, M_\epsilon$; this means that the minimum and maximum possible values satisfy $t_\min = M_\epsilon(t_\max)$ and $t_\max=M_1(t_\min)$. But these equations together yield quadratic equations for $t_\min$ and $t_\max$, with solutions $t_\min=\frac12(2\epsilon-1+\sqrt{4\epsilon^2+1})$ and $t_\max = \frac1{2\epsilon}(1+\sqrt{4\epsilon^2+1})$; thus, the maximal domain of the distribution is $[\frac12(2\epsilon-1+\sqrt{4\epsilon^2+1}),\frac1{2\epsilon}(1+\sqrt{4\epsilon^2+1})]$. For instance, when $\epsilon=\frac12$, the maximum interval of support is $[\frac{\sqrt{2}}2,1+\sqrt{2}]$.

Given this, we can see that the range of the map $M_1$ on this domain is $[1+\frac{1}{2\epsilon}(\sqrt{4\epsilon^2+1}-1),\frac1{2\epsilon}(1+\sqrt{4\epsilon^2+1})]$, and the range of the map $M_\epsilon$ is exactly $\epsilon$ times this, $[\frac12(2\epsilon-1+\sqrt{4\epsilon^2+1}),\frac12(1+\sqrt{4\epsilon^2+1})]$; e.g., when $\epsilon=\frac12$, the range of $M_1$ is $[\sqrt{2},1+\sqrt{2}]$ and the range of $M_\epsilon$ is $[\frac{\sqrt2}{2},\frac12(1+\sqrt{2})]$.

Now, notice that in this example the ranges of $M_1$ and $M_\epsilon$ are disjoint: $\frac12(1+\sqrt{2})\lt \sqrt{2}$, so there's no overlap. In fact, this holds true in general: we have $\frac12(1+\sqrt{4\epsilon^2+1})\lt1+\frac{1}{2\epsilon}(\sqrt{4\epsilon^2+1}-1)$ for all $0\lt\epsilon\lt 1$, so the maps $M_1$ and $M_\epsilon$ are disjoint on their mutual domain. This implies that the support of the distribution is actually a Cantor set: totally disconnected, perfect, nowhere dense, etc. (With a lot of algebra you could probably find the Hausdorff dimension, but I don’t trust my algebra or my transcription skills to and from Alpha enough to get that one right.)

Since the structure of the domain is so complicated, an explicit description of any probability measure on the domain is going to be pretty messy. In your case, there's an additional complication: because the maps are 'inversive' (decreasing), for any given fixed $p$ the probability is neither monotonically increasing nor monotonically decreasing on the (non-trivial) overlap of any interval with the domain. For instance, letting $D$ be the overall domain, we can see that the order(s) of iterates of the two maps $M_i$ applied to this domain are $\langle M_\epsilon(D)\lt M_1(D)\rangle$; $\langle M_\epsilon(M_1(D))\lt M_\epsilon(M_\epsilon(D)\lt M_1(M_1(D)\lt M_1(M_\epsilon(D)\rangle$; $\langle M_\epsilon(M_1(M_\epsilon(D)))$ $\lt M_\epsilon(M_1(M_1(D)))$ $\lt M_\epsilon(M_\epsilon(M_\epsilon(D)))$ $\lt M_\epsilon(M_\epsilon(M_1(D)))$ $\lt M_1(M_1(M_\epsilon(D)))$ $\lt M_1(M_1(M_1(D)))$ $\lt M_1(M_\epsilon(M_\epsilon(D)))$ $\lt M_1(M_\epsilon(M_1(D)))\rangle$; etc. Thus, the probabilities that $Y$ lies in each of the third-order intervals (listed by increasing $Y$) are $\langle p^2(1-p), p(1-p)^2,p^3, p^2(1-p),p(1-p)^2,(1-p)^3, p^2(1-p), p(1-p)^2\rangle$; for instance, letting $p=\frac13$, these are $\frac1{27}\langle2,4,1,2,4,8,2,4\rangle$. These sequences have a close relation with dragon curves, and the problem as a whole is very closely coupled to (weighted) iterated function systems; I'd encourage looking up more information on those.

(Also, it's possibly worth noting that in the 'limit' case where $\epsilon =0$ — i.e., the $X_i$ are either $0$ or $1$ — the support consists of exactly the convergents to $\phi=\frac12(1+\sqrt{5})$, and the probability of attaining the $n$th convergent is geometric, $p^n(1-p)$. This is because any values of $X_i$ after the first zero are irrelevant.)