Perhaps the simplest approach is to consider the problem as the calculation of a functional derivative. Let's put $$ (a_1, \ldots, a_n, b_1, \ldots, b_n)\triangleq (\mathbf{a}, \mathbf{b})\in\Bbb C^{2n} $$ and, in an analogous fashion, $$ (\delta_1, \ldots, \delta_n, \gamma_1, \ldots, \gamma_n)\triangleq (\boldsymbol{\delta}, \boldsymbol{\gamma})\in\Bbb C^{2n} $$ Then, defining the complex valued functional $$ \mathfrak{F}\big[(\mathbf{a}, \mathbf{b})\big]= \frac{\prod_{j=1}^n a_j}{\prod_{j=1}^n b_j} $$ it is natural to write $$ \mathfrak{F}\big[(\mathbf{a}, \mathbf{b})+\epsilon(\boldsymbol{\delta}, \boldsymbol{\gamma})\big]= \frac{\prod_{j=1}^n (a_j+\epsilon\delta_j)}{\prod_{j=1}^n (b_j+\epsilon \gamma_j)}\quad \forall \epsilon >0 $$ and calculating the functional derivative as usual $$ {\delta \mathfrak{F}}\big[(\mathbf{a}, \mathbf{b})\big](\boldsymbol{\delta}, \boldsymbol{\gamma}) \triangleq \left.\frac{\mathrm{d}}{\mathrm{d}\epsilon}\mathfrak{F}[(\mathbf{a}, \mathbf{b})+\epsilon(\boldsymbol{\delta}, \boldsymbol{\gamma})]\right|_{\epsilon=0} $$$$ {\delta \mathfrak{F}}\big[(\mathbf{a}, \mathbf{b})\big](\boldsymbol{\delta}, \boldsymbol{\gamma}) \triangleq \left.\frac{\mathrm{d}}{\mathrm{d}\epsilon}\mathfrak{F}\big[(\mathbf{a}, \mathbf{b})+\epsilon(\boldsymbol{\delta}, \boldsymbol{\gamma})\big]\right|_{\epsilon=0} $$ This procedure occasionally shows up in problems of sensitivity analysis in signal processing. This functional derivative is linear respect to the variation $(\boldsymbol{\delta}, \boldsymbol{\gamma})$: if it is a Fréchet derivative, it is possibly relatively easy to obtain a nice estimate in the form $$ \left|\mathfrak{F}\big[(\mathbf{a}, \mathbf{b})+(\boldsymbol{\delta}, \boldsymbol{\gamma})\big]-\mathfrak{F}\big[(\mathbf{a}, \mathbf{b})\big]\right|\le c\left|{\delta \mathfrak{F}}\big[(\mathbf{a}, \mathbf{b})\big](\boldsymbol{\delta}, \boldsymbol{\gamma})\right| $$$$ \left|\mathfrak{F}\big[(\mathbf{a}, \mathbf{b})+(\boldsymbol{\delta}, \boldsymbol{\gamma})\big]-\mathfrak{F}\big[(\mathbf{a}, \mathbf{b})\big]\right|\le c\left|{\delta \mathfrak{F}}\big[(\mathbf{a}, \mathbf{b})\big](\boldsymbol{\delta}, \boldsymbol{\gamma})\right|\label{1}\tag{☆} $$ with a sharp constant $c$. However, even if it is a Gâteaux derivative (i.e. it depends on the structure of the variation $(\boldsymbol{\delta}, \boldsymbol{\gamma})$) you may able to obtain a sharp estimate.
Addendum. Without going in the details and assuming the (plausible) non degeneracy hypothesis $\left|{\delta \mathfrak{F}}\big[(\mathbf{a}, \mathbf{b})\big](\boldsymbol{\delta}, \boldsymbol{\gamma})\right| >0$ when $(\boldsymbol{\delta}, \boldsymbol{\gamma})\neq 0$, further fairly elementary steps for possibly arriving at \eqref{1} are the following ones: $$ \begin{split} \left|\mathfrak{F}\big[(\mathbf{a}, \mathbf{b})+(\boldsymbol{\delta}, \boldsymbol{\gamma})\big]-\mathfrak{F}\big[(\mathbf{a}, \mathbf{b})\big]\right| &= \int\limits_0^1\frac{\mathrm{d}}{\mathrm{d}\epsilon}\mathfrak{F}\big[(\mathbf{a}, \mathbf{b})+\epsilon(\boldsymbol{\delta}, \boldsymbol{\gamma})\big]\mathrm{d}\epsilon \\ & \le \int\limits_0^1\left| \frac{\mathrm{d}}{\mathrm{d}\epsilon}\mathfrak{F}\big[(\mathbf{a}, \mathbf{b})+\epsilon(\boldsymbol{\delta}, \boldsymbol{\gamma})\big] \right| \mathrm{d}\epsilon \\ & = \left| {\delta \mathfrak{F}}\big[(\mathbf{a}, \mathbf{b})\big](\boldsymbol{\delta}, \boldsymbol{\gamma})\right|\int\limits_0^1\left| \frac{\frac{\mathrm{d}}{\mathrm{d}\epsilon}\mathfrak{F}\big[(\mathbf{a}, \mathbf{b})+\epsilon(\boldsymbol{\delta}, \boldsymbol{\gamma})\big]}{\left.\frac{\mathrm{d}}{\mathrm{d}\epsilon}\mathfrak{F}\big[(\mathbf{a}, \mathbf{b})+\epsilon(\boldsymbol{\delta}, \boldsymbol{\gamma})\big]\right|_{\epsilon=0}} \right| \mathrm{d}\epsilon \\ & \le c\left|{\delta \mathfrak{F}}\big[(\mathbf{a}, \mathbf{b})\big](\boldsymbol{\delta}, \boldsymbol{\gamma})\right| \end{split} $$ and therefore the relation $$ \int\limits_0^1\left| \frac{\frac{\mathrm{d}}{\mathrm{d}\epsilon}\mathfrak{F}\big[(\mathbf{a}, \mathbf{b})+\epsilon(\boldsymbol{\delta}, \boldsymbol{\gamma})\big]}{\left.\frac{\mathrm{d}}{\mathrm{d}\epsilon}\mathfrak{F}\big[(\mathbf{a}, \mathbf{b})+\epsilon(\boldsymbol{\delta}, \boldsymbol{\gamma})\big]\right|_{\epsilon=0}} \right| \mathrm{d}\epsilon\le c \label{2}\tag{☆☆} $$ may possibly be a first guess for the sought for constant. Note that the integral on the right side of \eqref{2} can be difficult to estimate: however, despite requiring possibly tedious calculations, you can made explicit its form.