A closed form expression is provided in the following paper.
SV Amari, RB Misra, Closed-form expressions for distribution of sum of exponential random variables, IEEE Transactions on Reliability, 46 (4), 519-522.
Note: My analysis assumes $\lambda=1$, but can be readily extended to the $\lambda \neq 1$ case.
We can write $\mathbb{P}(Y \leq y, N=k) = \mathbb{P}(Y \leq y | N=k)\mathbb{P}(N=k)$.
Let's use the fact that the PDF of the sum of two independent random variables is given by the convolution of the the respective PDFs. For $k=1$, Y|N=1 = $X_1 ~\sim exp(\lambda)$. For $k=2$, $Y|N=2=X_1+X_2$ and $f_{Y|N=2}(\cdot) = f_x(\cdot) \star f_x(\cdot)$. The convolution can be computed (by definition) as:
$f_{Y|N=2}(y) = \int_{0}^{y} e^{-x}e^{-(y-x)}dx = ye^{-y} \; \forall \; y \geq 0$.
Extending to $k=3$, the PDF of $Y|N=3 = X_1+X_2+X_3$ can be computed as:
$f_{Y|N=3}(y) = \int_0^y x e^{-x} e^{-(y-x)} dx = \frac{1}{2}y^2 e^{-y} \; \forall \; y \geq 0$.
This can be generalised to (and proved via mathematical induction) to: $f_{Y|N=k}(y) = \frac{1}{(k-1)!}y^{k-1}e^{-y} \; \forall \; y \geq 0$.
Now that we have the PDF, we can compute the CDF as:
$\mathbb{P}(Y \leq y | N=k) = \int_0^y f_{Y|N=k}(x)dx = \int_0^y \frac{1}{(k-1)!}x^{k-1}e^{-x} dx = \frac{1}{(k-1)!}M_k(y)$.
Let's use integration by parts to evaluate the above expression. We have the recursion: $M_k(y) = \int_0^y x^{k-1}e^{-x}dx = -x^{k-1}e^{-x} |_0^y + \int_0^y (k-1)x^{k-2}e^{-x} dx = -y^{k-1}e^{-y} + (k-1)M_{k-1}(y)$. With a bit of algebra one can show that $M_k(y) = (k-1)! - e^{-y} \sum_{j=0}^{k-1} \frac{(k-1)!}{j!} y^j $. Substituting in our previous expression, we get:
$\mathbb{P}(Y \leq y | N=k) = 1 - e^{-y}\sum_{j=0}^{k-1} \frac{y^j}{j!} $.
Finally, the desired joint distribution is given by:
$\mathbb{P}(Y \leq y, N=k) = \left[1 - e^{-y}\sum_{j=0}^{k-1} \frac{y^j}{j!} \right](1-p)^kp$.
You can compute $\mathbb{P}(Y=y) = \sum_{k=1}^\infty \mathbb{P}(Y=y, N=k)$ and plug in the expression computed above. I am not sure if there is a closed form analytical answer.
Best Answer
$$\begin{align} \int_{0}^{\infty}\int_{0}^{y} p(x,y) \,dxdy &=\int_{0}^{\infty}\frac{\gamma^ky^{k-1}e^{-\gamma y}}{(k-1)!}\int_{0}^{y} \lambda e^{-\lambda x} \,dxdy \\ &=\frac{\gamma^k}{(k-1)!}\int_{0}^{\infty}y^{k-1}e^{-\gamma y}(1-e^{-\lambda y}) \,dy\\ &=\frac{\gamma^k}{(k-1)!}\int_{0}^{\infty}y^{k-1}e^{-\gamma y}\,dy -\frac{\gamma^k}{(k-1)!}\int^\infty_0 y^{k-1}e^{-(\gamma+\lambda) y} \,dy\\ &=\frac{1}{(k-1)!}\int^\infty_0 t^{k-1}e^{-t}\,dt-\frac{\gamma^k}{(\gamma+\lambda)^k(k-1)!}\int^\infty_0t^{k-1}e^{-t}\,dt\\ &=\frac{\Gamma(k)}{(k-1)!}-\frac{\gamma^k\Gamma(k)}{(\gamma+\lambda)^k(k-1)!}=1-\frac{\gamma^k}{(\lambda+\gamma)^k} \end{align} $$