Dependence of the multipole moments on the origin Multipole moments of a system are defined with an explicit refrence to the co-ordinate system, e.g.
$$\boldsymbol{d}=\int dV\, \rho\,\boldsymbol{r}\\
\boldsymbol{\mu}=\frac{1}{2c}\int dV\, [\boldsymbol{r}\times\boldsymbol{j}]\\
Q_{ij}=\int dV\, \rho(3r_ir_j-\delta_{ij}r^2)$$
Only the leading multipole moment is independent of this choice. Let's assume it's the dipole moment $\boldsymbol{d}$. Now consider a harmonic system with frequency $\omega$. The system will radiate energy, and the radiation intensity $I$ can be written in terms of the multipole moments (sorry if coefficients are not correct)
$$I=\frac{2}{3}\frac{\omega^4\boldsymbol{d}^2}{c^3}+\frac{2}{3}\frac{\omega^4\boldsymbol{\mu}^2}{c^3}+\frac{\omega^6Q_{ij}^2}{180c^5}+\dots$$
Now $\boldsymbol{d},\boldsymbol{\mu}$ and $Q_{ij}$ are the amplitudes of the corresponding oscillating moments. According to textbooks, if the dipole moment is present than the dipole term gives the largest contribution to the radiation intensity and other terms could be considered to be corrections. Suppose we want to increase the accuracy by including the magnetic and the quadrupole contributions. My problem is that these corrections do not seem to be well-defined, because the moments themselvs depend on our choice of the coorditate origin. Take for example the magnetic moment and shift the coordinate system by vector $\boldsymbol{a}$. The corresponding change of the magnetic moment is given by
$$\Delta\boldsymbol{\mu}=\frac{1}{2c}\int dV\, [\boldsymbol{a}\times\boldsymbol{j}]$$
It is reasonable to take $\boldsymbol{a}$ to be of the order of the whole system size $L$. Indeed, for a generic system with no special symmetry the uncertanty of where is "the center" is of the order of the the system size. But then, the $\Delta\boldsymbol{\mu}$ and $\boldsymbol{\mu}$ seems to be of the same order $\Delta\boldsymbol{\mu}=\boldsymbol{\mu}\propto L j/c$. So the qestion is whether one really can increase the accuracy of the multipole expansion by including only some of the higher order terms (the full sum should be invariant of course)? 

Edit: let's consider a specific example, always a good idea. Take a point electric dipole oscillating with frequency $\omega$ and amplitude $\boldsymbol{d}$, i.e. $\boldsymbol{d}(t)=e^{-i\omega t}\boldsymbol{d}$. The total radiation is given by the well-known expression 
$$I_0=\frac{2\omega^4\boldsymbol{d}(t)^2}{3c^3}$$ Of course, the same must be true for an electric dipole located not at the origin but at some position $\boldsymbol{a}$. The corresponding current and charge densities are
$$\rho = -(\boldsymbol{d}(t)\nabla)\delta(\boldsymbol{r}-\boldsymbol{a}),\qquad \boldsymbol{j}=-i\omega\boldsymbol{d}(t)\delta(\boldsymbol{r}-\boldsymbol{a})$$
The dipole moment is of course $\boldsymbol{d}(t)$. Yet, the magnetic moment is non-vanishing 
$$\boldsymbol{\mu}=\frac{1}{2c}\int\, dV [\boldsymbol{r}\times -i\omega\boldsymbol{d}(t)\delta(\boldsymbol{r}-\boldsymbol{a})]=-\frac{ik}{2}[\boldsymbol{a}\times \boldsymbol{d}(t)]$$
and for generic $\boldsymbol{a}$ will give a non-vanishing contribution to the radiation intensity according to the general formula. However, the total intensity should be given by plain dipole radiation $I_0$. This problem is even stronger than the one I've originally posed: not only arbitrariness is present in the subleading terms, but even inclusion of the higher-order terms won't save the mismatch.
 A: Consider first the expansion of the electrostatic potential $\Phi$ in a stationary situation
$$\Phi(\vec{r}) = \frac{Q}{4 \pi \epsilon_0 r} + \frac{\vec{d}\cdot \vec{r}}{4\pi \epsilon_0 r^3} + \mathcal{O}(r^{-3})$$
Now let us shift our origin by $\vec{a}, \vec{R} = \vec{r} - \vec{a}$, where $|\vec{a}| \ll r$. The monopole term then gets expanded as
$$\frac{Q}{4 \pi \epsilon_0 r} = \frac{Q}{4 \pi \epsilon_0 |\vec{R} + \vec{a}|} = \frac{Q}{4 \pi \epsilon_0 R} - \frac{Q(\vec{a}\cdot \vec{R})}{4 \pi \epsilon_0 R^3} + \mathcal{O}(R^{-3})$$
The dipole term is simply
$$\frac{\vec{d}\cdot \vec{r}}{4\pi \epsilon_0 r^3} = \frac{\vec{d}\cdot \vec{R}}{4\pi \epsilon_0 R^3} + \mathcal{O}(R^{-3})$$
Now let's take a look at the whole potential
$$\Phi(\vec{r}) = \Phi(\vec{R}) = \frac{Q}{4 \pi \epsilon_0 R} + \frac{(\vec{d}-Q\vec{a})\cdot \vec{R}}{4\pi \epsilon_0 R^3} + \mathcal{O}(R^{-3})$$
Note that this is (within the approximation) just a different description of the same potential as above. In other words, if you take a look at the same physical point, you will get the same value of $\Phi$ independent of whether you are in the $\vec{r}$ coordinates or the $\vec{R}$ coordinates. 
Now let us denote the dipole defined with respect to $\vec{R}$ as $\vec{D}$ and compute
$$\vec{D} = \int \rho \vec{R} d V = \int \rho \vec{r} d V - \int \rho \vec{a} d V = \vec{d} - Q\vec{a}$$
We then see that our shifted potential in $\vec{R}$ coordinates can be written as
$$\Phi(\vec{R}) = \frac{Q}{4 \pi \epsilon_0 R} + \frac{\vec{D}\cdot \vec{R}}{4\pi \epsilon_0 R^3} + \mathcal{O}(R^{-3})$$
which is exactly the multipole expansion that you would get if you started from the $\vec{R}$ coordinates in the first place.
I.e., multipole expansions are covariant with respect to coordinate shifts. It is possible to show that this applies to all orders with more and more terms trickling down from lower orders to higher ones, and you can even show covariance of the multipole expansion with respect to the whole Poincaré group (with small translations).

You are now probably asking what is happening with your radiation formula then. The trick is that the formulas you give will apply only in inertial frames. In particular, $\vec{a}$ will typically be a constant shift. However, your radiation formula applies to oscillation magnitudes to which a constant $\vec{a}$ will have subleading or outright vanishing contributions.
Consider the dipole moment. We have $\vec{D} = \vec{d} - Q\vec{a}$. Then you see that since $\dot{Q} = \dot{a} = 0$, we have $\dot{D} = \dot{d}$ and there is no extra term arising in the radiation formula from the shift.
As for the dipole magnetic moment, we have
$$\vec{\mu}' = \vec{\mu} + \frac{1}{c}\vec{a} \times \int \vec{j}\, dV$$
It is a slightly more complicated argument why the extra term will be subleading. 
Let us first rewrite as a double integral using the divergence theorem
$$\int j_i (\vec{R}) \, dV(\vec{R}) = \int_{-\infty}^\infty \left(\int_{x_i=const.} \! \!\!\! \vec{j} \cdot d \vec{S} \right) d x_i = \int_{-\infty}^\infty \left( \int_{-\infty}^{x_i =const.} \! \! \! \nabla \cdot \vec{j}(\vec{R}') d V(\vec{R}') \right) d x_i$$
I.e., we use the fact that the integral of $j_i$ over a surface of constant $x_i$ can be also written as a divergence in a volume bounded by $x_i = const.$ (assuming of course that currents vanish outside the body, so the contributions from the other boundaries are just zero). Now let us use the continuity equation $\nabla \cdot \vec{j} = - \partial \rho/\partial t$ to finally express
$$\vec{\mu}' = \frac{1}{c}\int \vec{R} \times \vec{j} dV = \frac{1}{c}\int \vec{r} \times \vec{j} dV - \frac{1}{c}\vec{a} \times \int_{-\infty}^\infty \left(\int_{-\infty}^{\vec{x}} \!  \frac{\partial \rho} {\partial t} d V \right) d \vec{x}$$ 
(If you are unsure about what the vector product means, just write out the expressions using the Levi-Civita symbol and components.) 
Now let us assume that we have $\rho = \rho_0 + \rho_{osc} e^{i\omega t}$ and $\vec{j} = \vec{j}_{0} + \vec{j}_{osc} e^{i\omega t}$, where $\rho_0, \rho_{osc},\vec{j}_0 ,\vec{j}_{osc}$ are functions of position only. Then we see that 
$$\vec{\mu}'_{osc} = \vec{\mu}_{osc} - \frac{\omega}{c}\vec{a}\times \vec{\Delta}$$
where $\vec{\mu}_{osc} = \int \vec{j}_{osc} \times \vec{r} dV$, and $\vec{\Delta} = \int (\int^\vec{x} \rho_{osc} dV) d\vec{x}$. Since $\vec{d}_{osc} = \int \rho_{osc} \vec{r} dV$, we will have $\vec{\Delta} \sim \vec{d}$ and we can finally write
$$\vec{\mu}'_{osc} = \vec{\mu}_{osc} + \mathcal{O}( \frac{\omega}{c} a \,d_{osc} )$$ 
I.e. the shift induces only a subleading correction.

The gist of the argument is that if you have a body from which no currents are leaving, then a nonzero $\int \vec{j} dV$ corresponds to changes of charge density somewhere inside the body ($\partial \rho/\partial t \neq 0$). In a stationary oscillation, however, this corresponds to a term of a higher $\omega$ power as compared to the $\mu$ oscillation. 
The reason is the following: current has the dimensions $[Charge \cdot distance / time]$, the charge scale is determined by the overall charge in the body, the distance by the size of the body, and the time by 1) the crossing time of the charged particle in the body, and 2) the oscillation time of the charges. The term $\int \vec{r} \times \vec{j} dV$ captures this "crossing current" whose magnitude is not dependent on $\omega$, but the term $\vec{a} \times \int \vec{j} dV$ captures only the $\omega$-proportional oscillatory current.
