Sum of correlated exponential random variables

The pdf $f_Z(z)$ of the sum $Z=X+Y$ of any two jointly continuous random variables $X$ and $Y$ with joint pdf $f_{X,Y}(x,y)$ is as follows: $$text{For all } z, -infty < z < infty, ~~ f_Z(z) = int_{-infty}^infty f_{X,Y}(x,z-x) , mathrm dx.tag{1}$$

For the special case when $X$ and $Y$ are nonnegative random variables (including as a special case, exponential random variables) and so take on nonnegative values only, $f_{X,Y}(x,y)$ has value $0$ if at least one of $x$ and $y$ is smaller than $0$. Hence, in this case, the integrand $f_{X,Y}(x,z-x)$ in $(1)$ has value $0$ if $x < 0$ or if $z < x$. Consequently, if $z$ is a negative number, then the integrand in $(1)$ is always $0$ regardless of the value of $x$ and therefore so is the integral. All of which is just a long-winded way of saying that $f_Z(z)$ has value $0$ when $z<0$, that is, $Z$ takes on nonegative values only, which any idiot could have deduced from the fact that $Z=X+Y$ and both $X$ and $Y$ are nonnegative. But the approach is useful even for $z>0$ since now we have that the integrand in $(1)$ is zero when $x<0$ or when $x >z$ and so for nonnegative $X$ and $Y$, we can simplify $(1)$ to $$f_Z(z) = begin{cases}displaystyleint_0^z f_{X,Y}(x,z-x) , mathrm dx, & z geq 0,\quad\ 0, & z < 0end{cases} tag{2}$$

No further simplification of $(2)$ is possible in general.

For the special case when $X$ and $Y$ are independent random variables, $f_{X,Y}(x,y)$ factors into $f_X(x)f_Y(y)$ and so $(1)$ becomes the familiar convolution integral and $(2)$ the somewhat-less-familiar convolution integral for causal signals. But no such simplification is possible for nonindependent random variables $X$ and $Y$; we need the joint pdf to calculate $f_Z(z)$ and just knowing that $X$ and $Y$ are correlated random variables (whether exponential or Gaussian or whatever) is not enough.