Let $X_1,X_2,\dots$ are independent r.v. such that $P(X_i=1)=p=1-P(X_i=\epsilon_i)$, $0<\epsilon_i<1$
$$Y=X_1+\frac{X_1}{X_2+\frac{X_2}{X_3+\frac{X_3}{X_4+\dots}}}$$
1.What is the distribution of $Y$?
2.What is the characteristic function of $Y$?
[If you have read up to this line you must go to "Added" part in the end of this post. Starting with $X_i$'s are i.i.d. would be less complicated]
I do not know how to find these type of problem, so I had to post it. Any idea or website link would be helpful.
When I thought about it it was very easy for $X_i$'s are same, i.e. $$Y^*=X_1+\frac{X_1}{X_1+\frac{X_1}{X_1+\frac{X_1}{X_1+\dots}}}$$
In this case $Y^*=\frac{X_1+\sqrt{X^2_1+4X_1}}{2}$, which is not very interesting.
In the end I think it will need condition on $\epsilon_1,\epsilon_2,\dots$ to $Y$ be a random variable. But these are very trivial observation.
What can we say about 1 and 2 and what can we say more than 1 and 2 about $Y$?
Added:
As @Henry said we can think of less complicated version:
Let $X_1,X_2,\dots$ are i.i.d. r.v. such that $P(X_1=1)=p=1-P(X_1=\epsilon)$, $0<\epsilon<1$
$$Y=X_1+\frac{X_1}{X_2+\frac{X_2}{X_3+\frac{X_3}{X_4+\dots}}}$$
Here we can write as follows: $Y=Y_1$. $Y_1=X_1(1+\frac{1}{Y_2})$, $Y_2=X_2(1+\frac{1}{Y_3}),\dots$
It can be easily seen that $Y_i$'s are identical but not independent. So first we should find question 1, 2 for this case.
Added:
Some reference I found here in Mathoverflow. This may help.
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.)