I have a collection of $Z_1, ... Z_n$ which are all i.i.d standard normals. I have a collection of $X_1, ..., X_m$ such that each $X_i$ is a (potentially distinct) linear combination of $Z_1, ..., Z_n$, plus an extra constant term at the end.
i.e. $$X_i = a_{i1}Z_1 + ... + a_{in}Z_n + \mu_i$$
I am asked to find the joint m.g.f of $X_1, ..., X_m$. I applied the definition $$M_X(t_1, ..., t_m) = \mathbb{E}[exp(t_1X_1 + ... + t_mX_m)]$$
It seems that I can group the terms in the exponent by $Z_i$, and then find the joint m.g.f that way. However, I get a very nasty looking integral and am not sure how to proceed. With finding the distribution of $Z_1, ..., Z_n$, one would be able to do some refactoring so that one can obtain some $exp(f(t))$ times an integral of a normal distribution (which would sum to 1 and disappear, leaving a nice $exp(f(t))$). I'm not sure how that would be done here. Please let me know if perhaps another approach is required, thanks!
$$\begin{align*} \sum_{i=1}^m t_i X_i &= \sum_{i=1}^m t_i \left( \mu_i + \sum_{j=1}^n a_{ij} Z_j \right) \\ &= \sum_{i=1}^m t_i \mu_i + \sum_{i=1}^m \sum_{j=1}^n t_i a_{ij} Z_j \\ &= \sum_{i=1}^n t_i \mu_i + \sum_{j=1}^n Z_j \sum_{i=1}^m t_i a_{ij}. \end{align*}$$ Then exponentiation turns the outermost sum into a product over IID standard normals, the MGFs of which are easily computed. For example, $$M_{\boldsymbol X}(\boldsymbol t) \propto \prod_{i=1}^n M_{Z_i}\left( \sum_{i=1}^m t_i a_{ij} \right)$$ where the constant of proportionality is $\exp(\sum_{i=1}^n t_i \mu_i).$