It is true that no energy-momentum tensor for the gravitational field exists, however, it is easy to understand why not and then derive a perfectly correct formulation for conservation of energy and momenta of the gravitational field.
The invariance group of special relativity is the Poincare group. Energy and momentum combine in special relativity to form a 4-vector which belongs to a representation of the poincare group. The current of this four vector is the energy-momentum stress tensor whose divergence is zero.
When going from special to general relativity it is often assumed that any tensor quantity can be replaced by a similar one with ordinary derivatives replaced by covariant derivatives, this is not always the case. In general relativity the symmetry is the diffeomorphism group, not the Poincare group or Lorentz group. 4-vectors in GR only exist locally but energy and momentum are not local quantities so they cannot form a 4-vector. Instead they should form an object from a representation of the Diffeomorphism group.
If your spacetime is topologically equivalent (diffeomorphic) to $R^4$ then you can choose any global system of 4 co-ordinates and transform those co-ordinates using Poincare transforms. These are diffeomorphisms so this means you can embed the Poincare group in the diffeomorphism group by choosing such co-ordinates. For this reason it is possible to derive an energy-momentum pseudo-tensor for the gravitational field. It is co-ordinate dependent and not a tensor, but it works.
There is a better approach conceptually that is covariant and works for any topology. This is derived by applying Noether's first theorem directly to the Einstein-Hilbert action using the symmetry generators of the diffeomorphism group which are contravariant vector fields $k^\mu$. The result is a current with a linear dependence on the field $k^\mu$ which simplifies to the Komar Superpotential using the field equations
$J^{\mu} = (k^{\mu;\nu} - k^{\nu;\mu})_{;\nu}$
Using this formulation the energy and momenta belong to the dual of the adjoint representation of the diffeomorphism group.
Edit: I will add one more important point that is often misunderstood.
The matter and radiation part of the energy-momentum-stress tensor can be derived using the formula given in the question applied to the matter+radiation part of the Lagrangian
$T_{MR}^{\mu\nu} = -2 \frac{\delta L_{MR}}{\delta g_{\mu\nu}}$
If you use this expression on the full action as suggested in the reference it gives the gravitational equations of motion, which are dynamically zero. It is crucial to understand that this is not how to derive the Noether current which is correctly given by this expression (see Wikipedia for details)
$T_\mu{}^\nu =
\left( \frac{\partial L}{\partial \boldsymbol\phi_{,\nu}} \right) \cdot \boldsymbol\phi_{,\mu} - L\,\delta_\mu^\nu$
Some people cunfuse these two things and think that they give the same answer for the full Lagrangian, so that the Noether current must be zero under the field equations. This is certainly not the case. When the Noether current is derived correctly it gives the Komar Superpotential using the field equations and this is not zero. If you take a co-ordinate dependent approach you can alternatively use Noether's theorem to get pseudotensor expressions which again are not equal to zero.