Since it is not in my habit of just regurgitating stuff I see on the internet, here is my added value. How can one estimate the quadrupole gravitational radiation of a binary? How close does the binary have to be for it to coalesce within the age of the universe?

Let's begin. Suppose we have two black holes of mass M rotating around each other at a distance of R. Can we estimate the gravitational radiation power be? This is a classical back of the envelope problem, which we can do using dimensional analysis.

We are looking for the power P of the radiation. We know that the result can depend on the gravitational constant G, the speed of light c, the mass of the BHs M, the separation R and the frequency with which they rotate around each other ω.

From Buckingham's Pi theorem, we know that from the 6 variables above (P, G, c, M, R, ω) and 3 units (cm, gr, sec), we can construct 6-3 = 3 dimensionless numbers. To do better (that is, find a unique ratio between P and the other variables), we need some extra information.

First, the solution will only include the quadrupole moment Q. To see this, we need to look at the leading multipole order of two masses rotating.

The masses have a monopole, but it does not vary and therefore does not radiate.

The masses don't have a dipole moment

The masses do have a quadrupole moment. One way of seeing this is to describe the two masses as a sum of a cross with 1/2 M on the edges, and a second cross with 1/2 M on two opposing edges, and -1/2 M on the other (perpendicular) edges of the cross. At large distances, the first cross looks like a monopole (single mass), which the second cross is a classic quadrupole.

Thus, instead of any combination of M and R in the expression for P, we know that only the expression including a quadrupole should appear. The mass M is a monopole. Mass times distance is the dipole. Mass times area is a quadrupole. So, we'll have in the solution $Q = M R^2$.
Next, we also know that in the weak field limit, the gravitational field is linear in the mass, and thus also to the quadrupole moment. We also know that the energy in the gravitational waves is proportional to the field squared (as is the case with all linear waves). This means that the power radiated will be proportional to the quadrupole squared, that is: $P\propto Q^2$.

We now need to construct a dimensionless ratio between the following variables:

variable

units

$P/Q^2$

$1/(sec^3gr~cm^2)$

$c$

$cm/sec$

$\omega$

$1/sec$

$G$

$cm^3/(sec^2gr)$

To get rid of the gr in $P/Q^2$, we can take the combination $P/Q^2 G$ which has the units of $1/cm^5 sec$. To get rid of the cm, we can take the combination $P c^5/ Q^2 G$, which has the units of $1/sec^6$. Thus, the dimensionless ratio we are looking for is:
$$ \Pi = {P c^5 \over \omega^6 Q^2 G}. $$
In other words, the radiated power is going to be:
$$ P \sim {Q^2 G \omega^6 \over c^5}. $$
Furthermore, in circular orbits we have that that the centrifugal acceleration $\omega^2 r$ is equal to the gravitational acceleration $GM/r^2$. Thus, we have that
$$ \omega^2 = {GM \over R^3}.$$
The emitted power is therefore:
$$ P \sim {Q^2 G^4 M^3 \over c^5 R^3} \sim {G^4 M^5 \over c^5 R^5}. $$
The actual pre-factor is $16/\pi$, something dimensional analysis cannot pick up.

Next, if we want to estimate how long will it take a binary at radius R to in-spiral, we have to compare the power to the binding energy of the binary:
$$ \tau \sim {U_{grav} \over P} \sim {G M^2/ R \over G^4 M^5 / c^5 R^5} \sim { c^5 R^4 \over G^3 M^3 } .$$
If we care about R, we can invert the relation and get:
$$ R \sim \left({G^3 M^3 \tau \over c^5}\right)^{1/4}$$
For a Hubble time (the typical age of the universe) and solar type BHs, the radius has to be:
$$ R \sim 10^{11} cm \sim R_\odot, $$
that is, similar to the solar radius. If we're talking about two massive black holes, with say $10^6 M_\odot$ each, we're talking about a radius which is 3 104 times larger, or about 100 AU. Namely, for the massive blackholes to merge, there first must be a process which brings them close to each other, to within the size of the solar system. This can occur, for example, through three body collisions (i.e., giving energy and angular momentum to nearby stars, and blowing them to kingdom come).

Incidentally, these events are very interesting because they mostly likely occur, and they could in principle be observed with terrestrial gravitational wave detectors. The stellar types in LIGO and the massive BH's in LISA (which to be precise is not terrestrial...).