Let $\underset{\sim}{X}\sim N_k(\underset{\sim}{\mu},V)$ i.e. $X$ follows $k-$variate normal distribution with the variance-covariance matrix, $V_{k\times k}$. Then, pdf of $\underset{\sim}{X}$ is
$$f_\underset{\sim}{X}(\underset{\sim}{x})=\frac{e^{-\frac{1}{2}[(\underset{\sim}{x}-\underset{\sim}{\mu})^TV^{-1}(\underset{\sim}{x}-\underset{\sim}{\mu})]}}{(2\pi)^\frac{k}{2}\cdot|V|^\frac{1}{2}}$$
Then, the mgf of $\underset{\sim}{X}$ is
$$E(e^{\underset{\sim}{t}^T\underset{\sim}{X}})=\underset{x_k}{\int}…\underset{x_2}{\int}\underset{x_1}{\int}e^{\underset{\sim}{t}^T\underset{\sim}{x}}\cdot\frac{e^{-\frac{1}{2}[(\underset{\sim}{x}-\underset{\sim}{\mu})^TV^{-1}(\underset{\sim}{x}-\underset{\sim}{\mu})]}}{(2\pi)^\frac{k}{2}\cdot|V|^\frac{1}{2}}dx_1dx_2…dx_k$$
$$=\underset{x_k}{\int}…\underset{x_2}{\int}\underset{x_1}{\int}\frac{e^{-\frac{1}{2}[(\underset{\sim}{x}-\underset{\sim}{\mu})^TV^{-1}(\underset{\sim}{x}-\underset{\sim}{\mu})-2\underset{\sim}{t}^T\underset{\sim}{x}]}}{(2\pi)^\frac{k}{2}\cdot|V|^\frac{1}{2}}dx_1dx_2…dx_k$$
The trick to simplifying such integrals is making a square in the power of $e$. So, I consider
$$(\underset{\sim}{x}-\underset{\sim}{\mu})^TV^{-1}(\underset{\sim}{x}-\underset{\sim}{\mu})-2\underset{\sim}{t}^T\underset{\sim}{x}$$
$$=\underset{\sim}{x}^TV^{-1}\underset{\sim}{x}-\underset{\sim}{x}^TV^{-1}\underset{\sim}{\mu}-\underset{\sim}{\mu}^TV^{-1}\underset{\sim}{x}+\underset{\sim}{\mu}^TV^{-1}\underset{\sim}{\mu}-2\underset{\sim}{t}^T\underset{\sim}{x}$$
$\because\underset{\sim}{\mu}^TV^{-1}\underset{\sim}{x}$ is a scalar and $V$ is a symmetric matrix
$\therefore V^{-1}$ is also symmetric and $\underset{\sim}{\mu}^TV^{-1}\underset{\sim}{x}=(\underset{\sim}{\mu}^TV^{-1}\underset{\sim}{x})^T=\underset{\sim}{x}^T(V^{-1})^T\underset{\sim}{\mu}=\underset{\sim}{x}^TV^{-1}\underset{\sim}{\mu}$
Therefore,
$$(\underset{\sim}{x}-\underset{\sim}{\mu})^TV^{-1}(\underset{\sim}{x}-\underset{\sim}{\mu})-2\underset{\sim}{t}^T\underset{\sim}{x}$$
$$=\underset{\sim}{x}^TV^{-1}\underset{\sim}{x}-2\underset{\sim}{\mu}^TV^{-1}\underset{\sim}{x}+\underset{\sim}{\mu}^TV^{-1}\underset{\sim}{\mu}-2\underset{\sim}{t}^T\underset{\sim}{x}$$
$$=[\underset{\sim}{x}^TV^{-1}\underset{\sim}{x}-2\underset{\sim}{t}^T\underset{\sim}{x}]-2\underset{\sim}{\mu}^TV^{-1}\underset{\sim}{x}+\underset{\sim}{\mu}^TV^{-1}\underset{\sim}{\mu}$$
Let $\underset{\sim}{y}=V^{-\frac{1}{2}}\underset{\sim}{x}$. I consider
$(\underset{\sim}{y}-V^\frac{1}{2}\underset{\sim}{t})^T(\underset{\sim}{y}-V^\frac{1}{2}\underset{\sim}{t})$
$=\underset{\sim}{y}^T\underset{\sim}{y}-\underset{\sim}{y}^TV^\frac{1}{2}\underset{\sim}{t}-\underset{\sim}{t}^T(V^\frac{1}{2})^T\underset{\sim}{y}+\underset{\sim}{t}^T(V^\frac{1}{2})^TV^\frac{1}{2}\underset{\sim}{t}$
$=\underset{\sim}{y}^T\underset{\sim}{y}-2\underset{\sim}{t}^TV^\frac{1}{2}\underset{\sim}{y}+\underset{\sim}{t}^TV\underset{\sim}{t}\qquad\qquad $since, $\underset{\sim}{t}^TV^\frac{1}{2}\underset{\sim}{y}$ is a scalar and $V$ is symmetric
$=(V^{-\frac{1}{2}}\underset{\sim}{x})^TV^{-\frac{1}{2}}\underset{\sim}{x}-2\underset{\sim}{t}^TV^\frac{1}{2}V^{-\frac{1}{2}}\underset{\sim}{x}+\underset{\sim}{t}^TV\underset{\sim}{t}$
$=\underset{\sim}{x}^TV^{-1}\underset{\sim}{x}-2\underset{\sim}{t}^T\underset{\sim}{x}+\underset{\sim}{t}^TV\underset{\sim}{t}$
$\Rightarrow[\underset{\sim}{x}^TV^{-1}\underset{\sim}{x}-2\underset{\sim}{t}^T\underset{\sim}{x}]=(\underset{\sim}{y}-V^\frac{1}{2}\underset{\sim}{t})^T(\underset{\sim}{y}-V^\frac{1}{2}\underset{\sim}{t})-\underset{\sim}{t}^TV\underset{\sim}{t}$
Therefore,
$$(\underset{\sim}{x}-\underset{\sim}{\mu})^TV^{-1}(\underset{\sim}{x}-\underset{\sim}{\mu})-2\underset{\sim}{t}^T\underset{\sim}{x}$$
$$=(\underset{\sim}{y}-V^\frac{1}{2}\underset{\sim}{t})^T(\underset{\sim}{y}-V^\frac{1}{2}\underset{\sim}{t})-\underset{\sim}{t}^TV\underset{\sim}{t}-2\underset{\sim}{\mu}^TV^{-\frac{1}{2}}\underset{\sim}{y}+\underset{\sim}{\mu}^TV^{-1}\underset{\sim}{\mu}$$
$$=(\underset{\sim}{y}-V^\frac{1}{2}\underset{\sim}{t})^T(\underset{\sim}{y}-V^\frac{1}{2}\underset{\sim}{t})-2\underset{\sim}{\mu}^TV^{-\frac{1}{2}}(\underset{\sim}{y}-V^\frac{1}{2}\underset{\sim}{t})-2\underset{\sim}{\mu}^TV^{-\frac{1}{2}}V^\frac{1}{2}\underset{\sim}{t}-\underset{\sim}{t}^TV\underset{\sim}{t}+\underset{\sim}{\mu}^TV^{-1}\underset{\sim}{\mu}$$
$$=[(\underset{\sim}{y}-V^\frac{1}{2}\underset{\sim}{t})^T(\underset{\sim}{y}-V^\frac{1}{2}\underset{\sim}{t})-2\underset{\sim}{\mu}^TV^{-\frac{1}{2}}(\underset{\sim}{y}-V^\frac{1}{2}\underset{\sim}{t})]-2\underset{\sim}{\mu}^T\underset{\sim}{t}-\underset{\sim}{t}^TV\underset{\sim}{t}+\underset{\sim}{\mu}^TV^{-1}\underset{\sim}{\mu}$$
I am not sure how I can simplify this further to get a general form for the integrand- need help with the above.
Also, is there any easier way?
Note: I am not comfortable with multivariable calculus (in case it is absolutely required, please explain a little).
Best Answer
Suppose $X\sim N_k(\mu,\Sigma)$ where $\Sigma$ is positive definite.
Then by a definition of a multivariate normal distribution, any linear combination of $X$ has a univariate normal distribution. That is, $t^TX\sim N(t^T\mu,t^T\Sigma t)$ for any $t\in\mathbb R^k$.
Now moment generating function of some $Z\sim N(\mu,\sigma^2)$ is
$$M_Z(s)=E[e^{s Z}]=e^{\mu s+\sigma^2s^2/2}\quad,\,s\in\mathbb R$$
Using this fact, we have
$$M_X(t)=E[e^{t^TX}]=M_{t^TX}(1)=\exp\left(\mu^Tt+\frac12t^T\Sigma t\right)$$
Alternatively, for a direct proof you can decompose $\Sigma=BB^T$ for some nonsingular matrix $B$ since $\Sigma$ is positive definite. Transform $X\mapsto Y$ such that $Y=B^{-1}(X-\mu)$, i.e. $X=\mu+BY$. Then it follows that $Y=(Y_1,\ldots,Y_k)^T\sim N_k(0,I_k)$.
In other words $Y_1,\ldots,Y_k$ are independent standard normal. You can verify all of these using change of variables.
Therefore using MGF of standard normal distribution,
\begin{align} M_X(t)&=E[e^{t^T X}] \\&=E[e^{t^T(\mu+BY)}] \\&=e^{t^T\mu}E[e^{\ell^T Y}]\qquad\quad,\,\ell^T=t^TB \\&=e^{t^T\mu}E\left[e^{\sum_{i=1}^k \ell_i Y_i}\right]\quad,\,\ell=(\ell_1,\ldots,\ell_k) \\&=e^{t^T\mu}\prod_{i=1}^k E\left[e^{l_iY_i}\right] \\&=e^{t^T\mu}\prod_{i=1}^k e^{\ell_i^2/2} \\&=\exp\left(\mu^Tt+\frac12 \ell^T\ell\right) \\&=\exp\left(\mu^Tt+\frac12 t^T\Sigma t\right) \end{align}