Joint distribution of sums of exponential random variables

This is a problem from Ross’s Stochastic Processes [1]. Let $$\begin{array}{rl}{S}_1& =X_1\\ S_2& =X_1+X_2\\ S_3& =X_1+X_2+X_3\end{array}$$ where $X_1,X_2,X_3$ are i.i.d exponential random variables with rate $\lambda$. Find the joint distribution of $S_1,S_2,S_3$.

Let $f$ be the PDF of each $X_1,X_2,X_3$ (since identically distributed). Since $X_1,X_2,X_3$ are independent the joint PDF is $$f(x,y,z)=f(x)f(y)f(z)={\lambda }^3{e}^{-\lambda x}{e}^{-\lambda y}{e}^{-\lambda z}$$ Then we can find the joint CDF of $S_1,S_2,S_3$ $$\begin{array}{rl}P(S_1\le t_1,S_2\le t_2,S_3\le t_3)& ={\int }_0^{t_1}{\int }_0^{t_2-x}{\int }_0^{t_3-x-y}f(x,y,z)dzdydx\\ & ={\int }_0^{t_1}{\int }_0^{t_2-x}{\int }_0^{t_3-x-y}\lambda e^{-\lambda z}\lambda e^{-\lambda y}\lambda e^{-\lambda x}dzdydx\\ & ={\int }_0^{t_1}\lambda e^{-\lambda x}{\int }_0^{t_2-x}(1-e^{-\lambda (t_3-x-y)})\lambda e^{-\lambda y}dydx\\ & ={\int }_0^{t_1}\lambda e^{-\lambda x}[{\int }_0^{t_2-x}\lambda e^{-\lambda y}dy-{\int }_0^{t_2-x}{e}^{-\lambda (t_3-x)}dy]dx\\ & ={\int }_0^{t_1}\lambda e^{-\lambda x}[(1-e^{-\lambda (t_2-x)})-(t_2-x)e^{-\lambda (t_3-x)}]dx\\ & ={\int }_0^{t_1}\lambda e^{-\lambda x}dx-{\int }_0^{t_1}\lambda e^{-\lambda t_2}dx-{\int }_0^{t_1}(t_2-x)e^{-\lambda t_3}dx\\ & =1-e^{-\lambda t_1}-\lambda t_1{e}^{-\lambda t_2}-t_1{t}_2{e}^{-\lambda t_3}+\frac{1}{2}{t}_1^2{e}^{-\lambda t_3}\end{array}$$

Simulation

We can confirm these results with a simple simulation in R.

# Computed joint CDF
expect_joint <- function(t1, t2, t3, lambda = 1) {
  1 - exp(- lambda * t1) - lambda * t1 * exp(-lambda * t2) - 
    t1 * t2 * exp(-lambda * t3) + 1/2 * t1^2 *exp(- lambda * t3)
}

# Simulate 1000 realizations of S1, S2, S3
sim <- function(t1, t2, t3, rate = 1, n = 1000) {
  s1 <- rexp(n, rate)
  s2 <- s1 + rexp(n, rate)
  s3 <- s2 + rexp(n, rate)
  
  mean(s1 <= t1 & s2 <= t2 & s3 <= t3)
}

# P(S1 <= 1, S2 <= 2, S3 <= 3) with lambda = 1
t1 <- 1
t2 <- 2
t3 <- 3
rate <- 1
# Replicate the simulation 1000 times
sim_res <- replicate(1000, sim(t1, t2, t3, rate))

cat("Simulated mean:", round(mean(sim_res), 3))

Simulated mean: 0.422

cat("Expected joint distribution:", round(expect_joint(t1, t2, t3, rate), 3))

Expected joint distribution: 0.422

```

References

[1]

Ross, S.M. et al. 1996. Stochastic processes. Wiley New York.