Q:

Calculate PDF and CDF of a product of independent exponentially distributed random variables

Let x[1], x[2], ..., x[n] be independent random variables each of which is distributed exponentially, i.e. with the PDF 
f[x_,a_] := a Exp[-a x], a>0, x>=0

and let 
t = Product[x[i],{i,1,n}]

Calculate the PDF and the CDF of t for arbitrary n = 1, 2, 3, ...
Remark: As I have indicated in my answer to Probability: Calculating a multiple integral this is an interesting problem which I have solved in Mathematica. I propose it here for others to find possible different paths to the solution.

A:

The Problem 
Let $(X_1, \dots, X_n)$ denote independent and identically distributed variables, each with common Exponential pdf $f(x)$:
f = a Exp[-a x];       domain[f] = {x, 0, Infinity}  &&  {a > 0}; 

We seek the pdf of $\prod_{i=1}^n X_i$, for $n = 2, 3, \dots$
Solution
The pdf of the product of two Exponentials is simply:

(source: tri.org.au) 
where I am using the  TransformProduct function from the mathStatica package for Mathematica. 
The domain of support is:
 domain[g] = {x, 0, Infinity}  &&  {a > 0};

The product of 3, 4 and 5 Exponentials is obtained by iteratively applying the same function (here 3 times):

(source: tri.org.au) 
By induction, the pdf of the product of $n$ iid Exponentials is thus:
$$a^n \text{MeijerG}[\{ \{ \},  \{ \} \}, \{ \{0_1, \dots, 0_n \}, \{ \} \}, a^n x] \quad \quad \text{ for } x > 0 $$
Quick Monte Carlo check:  Here is a quick check comparing:

the theoretical pdf just obtained (when $n = 3$ and $a=2$):  RED DASHED curve,
to the empirical Monte Carlo pdf based:  squiggly BLUE curve

(source: tri.org.au) 
And all is well.
CDF
In the same way, we can find the cdf at any iteration $n$. For the current value $n=5$, the cdf $P_g(X<x)$ is:

(source: tri.org.au) 
... using mathStatica's Prob function. This takes the same form for each value of $n$. By induction, the general form for the cdf is:
MeijerG[{{1}, {}}, {lis, {0}}, a^n x]  where lis is an $n$-length vector of 1's.

A:

Although we have now the elegant solution by wolfies I shall post my solution which uses only standard functions of MMA version 8 and which might be of interest because of a strategy to overcome difficulties with integration.
We calculate the probability distribution functions PDF ($fp$) and CDF ($fc$) of a product of $n$ independent random variables $x_1, x_2, \text{...}, x_n$
distributed according to exponential distributions $a_i e^{-a_i x_i}$ for arbitrary integer $n > 0$.
First we note that it is sufficient to study random variables distributed with the same PDF and with the Parameter set to $a = 1$. 
In fact, the reader can easily show that for the CDF $fc$ with a set ${a_1, a_2, ..., a_n}$ of parameters we have 
$fc\left(n,t,\{a_1,a_2, ..., a_n\}\right)=fc\left(n,t*{\prod _{i=1}^n a_i},\{1_1,1_2, ..., 1_n\}\right)$
Now let
f[x_]:=Exp[-x]

The CDF for n = 2 is defined by
fci[2, s_] = 
 Integrate[f[x] f[y] UnitStep[s - x y], {x, 0, \[Infinity]}, {y, 0, \[Infinity]}, Assumptions -> s > 0]

(* 1 - 2 Sqrt[s] BesselK[1, 2 Sqrt[s]] *)

For n = 3 we have 
fci[3, s_] = 
 Integrate[f[x] f[y] f[z] UnitStep[s - x y z], {x, 0, \[Infinity]}, {y,0, \[Infinity]}, {z, 0, \[Infinity]}, Assumptions -> s > 0]

(* 
Integrate[
 E^-x (1 - 2 Sqrt[s/x] BesselK[1, 2 Sqrt[s/x]]), {x, 0, \[Infinity]}, 
 Assumptions -> s > 0] 
*)

Unfortunately, already for n = 3 Mathematica does not perform all the necessary integrations.
In order to make progress let's perform some "heuristic experimental mathematics" with Mathematica.
To this end we shall first study fci[3,s] for some small integer values of s and look for patterns.
Table[{s, fci[3, s]}, {s, 1, 5}]

$\left\{\left\{1,1-\text{MeijerG}\left[\{\{\},\{\}\},\left\{\left\{-\frac{1}{2},\frac{1}{2},\frac{1}{2}\right\},\{\}\right\},1\right]\right\},\left\{2,1-\sqrt{2} \text{MeijerG}\left[\{\{\},\{\}\},\left\{\left\{-\frac{1}{2},\frac{1}{2},\frac{1}{2}\right\},\{\}\right\},2\right]\right\},\left\{3,1-\sqrt{3} \text{MeijerG}\left[\{\{\},\{\}\},\left\{\left\{-\frac{1}{2},\frac{1}{2},\frac{1}{2}\right\},\{\}\right\},3\right]\right\},\left\{4,1-2 \text{MeijerG}\left[\{\{\},\{\}\},\left\{\left\{-\frac{1}{2},\frac{1}{2},\frac{1}{2}\right\},\{\}\right\},4\right]\right\},\left\{5,1-\sqrt{5} \text{MeijerG}\left[\{\{\},\{\}\},\left\{\left\{-\frac{1}{2},\frac{1}{2},\frac{1}{2}\right\},\{\}\right\},5\right]\right\}\right\}$
Here the special function MeijerG appears which is, fortunately, contained as a standard function in Mathematica.
The pattern is easily seen, and we "conclude" that for n = 3 we have
fcx[3, t_] :=  1 - Sqrt[t] MeijerG[{{}, {}}, {{-(1/2), 1/2, 1/2}, {}}, t]

Now let's increase n
fci[4, s_] := 
 Integrate[f[x] f[y] f[z] f[u] UnitStep[s - x y z u], {x, 0, \[Infinity]}, {y,
    0, \[Infinity]}, {z, 0, \[Infinity]}, {u, 0, \[Infinity]}, 
  Assumptions -> s > 0]

Again for some integer arguments
Table[fci[4, s], {s, 1, 5}]

$\left\{1-\text{MeijerG}\left[\{\{\},\{\}\},\left\{\left\{-\frac{1}{2},\frac{1}{2},\frac{1}{2},\frac{1}{2}\right\},\{\}\right\},1\right],1-\sqrt{2} \text{MeijerG}\left[\{\{\},\{\}\},\left\{\left\{-\frac{1}{2},\frac{1}{2},\frac{1}{2},\frac{1}{2}\right\},\{\}\right\},2\right],1-\sqrt{3} \text{MeijerG}\left[\{\{\},\{\}\},\left\{\left\{-\frac{1}{2},\frac{1}{2},\frac{1}{2},\frac{1}{2}\right\},\{\}\right\},3\right],1-2 \text{MeijerG}\left[\{\{\},\{\}\},\left\{\left\{-\frac{1}{2},\frac{1}{2},\frac{1}{2},\frac{1}{2}\right\},\{\}\right\},4\right],1-\sqrt{5} \text{MeijerG}\left[\{\{\},\{\}\},\left\{\left\{-\frac{1}{2},\frac{1}{2},\frac{1}{2},\frac{1}{2}\right\},\{\}\right\},5\right]\right\}$
We find very similar expressions to those for n = 3, only the list in the MeijerG-function is longer by one
fcx[4, s_] = 1 - Sqrt[s ] MeijerG[{{}, {}}, {{-(1/2), 1/2, 1/2, 1/2}, {}}, s]

We "conclude" that generally for integer n > 0 and real s > 0 the CDF is given by
fcxx[n_, s_] := 
 1 - Sqrt[s ] MeijerG[{{}, {}}, {Join[{-(1/2)}, Array[1/2 &, n - 1]], {}}, s]

Calling it leads to an automatic simplification:
fcxx[4, s]

(*
1 - MeijerG[{{}, {}}, {{0, 1, 1, 1}, {}}, s]
*)

Hence we have found good reasons to assume that the general CDF is given by
fc[n_, s_] := 1 - MeijerG[{{}, {}}, {Join[{0}, Array[1 &, n - 1]], {}}, s]

This is in fact the same as the results of wolfies
fc1[n_, s_] := MeijerG[{{1}, {}}, {Array[1 &, n], {0}}, s] (* wolfies, 07.02.15 *)

as can be shown by plotting both functions for some values of n.
The PDF is given by the derivative with respect to s
fpx[n_, s_] := D[fc[n, s], s]

The first few are
Table[fpx[n, s], {n, 1, 5}]

(*
{E^-s, BesselK[0, 2 Sqrt[s]] - BesselK[1, 2 Sqrt[s]]/Sqrt[s] + 
  BesselK[2, 2 Sqrt[s]], MeijerG[{{}, {}}, {{0, 0, 0}, {}}, s], 
 MeijerG[{{}, {}}, {{0, 0, 0, 0}, {}}, s], 
 MeijerG[{{}, {}}, {{0, 0, 0, 0, 0}, {}}, s]}
*)

We conclude for the PDF the expression
fp[n_, t_] := MeijerG[{{}, {}}, {Array[0 &, n], {}}, t]

This is exactly the same expression as that obtained by wolfies.
EDIT #1 Moments, Asymptotics, Complex Integral
Moments
Because the variables are independent we have for the k-th moment of t
$m\left(t^k\right) = m\left((x_1 x_2 ... x_n)^k\right) = m\left(x_1^k \right) m\left(x_2^k\right) ... m\left(x_n^k \right) = m\left(x^k \right)^n$
and because
$m\left(x^k \right) = \int_0^{\infty } \exp (-x) x^k \, dx = \Gamma (k+1) = k!$
the k-th moment of t is just $(k!)^n$.
This result is confirmed by direct calculation
Table[1/(k!)^2 Integrate[t^k fp[2, t], {t, 0, \[Infinity]}], {k, 1, 10}]

(*
{1, 1, 1, 1, 1, 1, 1, 1, 1, 1}
*)

Table[1/(k!)^3 Integrate[t^k fp[3, t], {t, 0, \[Infinity]}], {k, 1, 10}]

(*
{1, 1, 1, 1, 1, 1, 1, 1, 1, 1}
*)

and so on for 
n = 4, 5, ...

Asymptotics at infinity
For the first few n we have
Series[fp[1, t], {t, \[Infinity], 1}] // Normal

$e^{-t}$
Series[fp[2, t], {t, \[Infinity], 1}] // Normal

$e^{-2 \sqrt{t}} \left(\sqrt{\pi } \left(\frac{1}{t}\right)^{1/4}-\frac{1}{16} \sqrt{\pi } \left(\frac{1}{t}\right)^{3/4}\right)$
Series[fp[3, t], {t, \[Infinity], 1}] // Normal

$e^{-3 t^{1/3}} \left(\frac{4 \pi  \left(\frac{1}{t}\right)^{1/3}}{\sqrt{3}}-\frac{4 \pi  \left(\frac{1}{t}\right)^{2/3}}{9 \sqrt{3}}\right)$
Looking at the expressions up to, say n = 10, we can guess for the leading term the following expression
$\text{fp}(n,t\to \infty ) = \sqrt{\frac{2}{\pi  n t}}\left(\sqrt{2\pi  t}\right)^n \exp \left(-n t^{\frac{1}{n}}\right)$
We can see that for large n the exponential decay becomes very slow.
Asymptotics at zero
For the first few n we get
Series[fp[1, t], {t, 0, 1}] // Normal // TraditionalForm

$1-t$
Series[fp[2, t], {t, 0, 1}] // Normal // TraditionalForm

$t (-\log (t)-2 \gamma +2)-\log (t)-2 \gamma$
where $\gamma$ is the constant EulerGamma.
Series[fp[3, t], {t, 0, 1}] // Normal // TraditionalForm

$\frac{1}{4} t \left(-2 \log ^2(t)-12 \gamma  \log (t)+12 \log (t)-\pi ^2-18 \gamma ^2+36 \gamma -24\right)+\frac{1}{4} \left(2 \log ^2(t)+12 \gamma  \log (t)+\pi ^2+18 \gamma ^2\right)$
...
Hence the general leading term for n>1 is a power of $log(1/t)$:
$\text{fp}(n,t\to 0 ) = \frac{\log ^{n-1}\left(\frac{1}{t}\right)}{(n-1)!}$
MeijerG as a complex integral
Looking at the General deifition of this function in terms of a complex integral (from http://functions.wolfram.com/HypergeometricFunctions/MeijerG/02/0001/)
$\text{MeijerG}\left[\left\{\left\{a_1,a_n\right\},\left\{a_{1+n},a_p\right\}\right\},\left\{\left\{b_1,b_m\right\},\left\{b_{1+m},b_q\right\}\right\},z\right]==\frac{1}{2\pi  i } \text{ContourIntegrate}\left[\frac{z^{-s} \left(\prod _{k=1}^n \text{Gamma}\left[1-s-a_k\right]\right) \prod _{k=1}^m \text{Gamma}\left[s+b_k\right]}{\left(\prod _{k=1+n}^p \text{Gamma}\left[s+a_k\right]\right) \prod _{k=1+m}^q \text{Gamma}\left[1-s-b_k\right]},\{s,\mathcal{L}\}\right]$
We find that all $a_i$ are missing, ${b_1, ... , b_m} = {0_1, ... , 0_n}$ and all other b's are missing so that
$fp(n,t)$ has the representation
$fp(n,t) = \frac{1}{2 \pi  i}\int _{c - i \infty }^{c + i \infty }t^{-s} \Gamma (s)^nds$
where c is a real number >0. Changing the integration variable to s = 1 + I u we find
$fp(n,t) = \frac{1}{2 \pi }\int _{-\infty }^{+ \infty }t^{-(1+ i u)} \Gamma (1+ i u)^ndu$
The identity of this expression to the original one can be shown with examples using NIntegrate.
Observation: the expression for $fp(n,t)$ is nothing but the inverse Mellin transformation of the n-th power of the Mellin transform of the exponential Distribution, which is igiven by 
$M( f(t), s) = \int_0^{\infty } f(t) t^{s-1} \, dt$
specifically
$\int_0^{\infty } \text{Exp}(-t) t^{s-1} \, dt = \Gamma (s)$

