The answer is basically that there are formulas for these matrix elements, but they are too complicated to be useful. However, the expectation $\langle n|x^{2p}|n\rangle$ has a reasonably simple formula.
Generating Function
First of all we can calculate the exponential generating function referred to in the question and in Cosmas Zachos's answer.
The position operator is $ x = \sqrt{\hbar/2m\omega}(a^\dagger + a) $, so setting $ \lambda = t\sqrt{\hbar/2m\omega} $ gives $e^{tx} = e^{\lambda(a^\dagger + a)}$. We can use the Baker–Campbell–Hausdorff formula to normal order this operator,
\begin{align}
e^{\lambda a^\dagger}e^{\lambda a} = e^{\lambda(a^\dagger + a) + \frac{1}{2}[\lambda a^\dagger, \lambda a]} &\implies e^{\lambda(a^\dagger + a)} = e^{\frac{\lambda^2}{2}}e^{\lambda a^\dagger}e^{\lambda a}\\
&\implies \langle m|e^{tx}|n\rangle = e^{\frac{\lambda^2}{2}}\color{blue}{\langle m|e^{\lambda a^\dagger}}\color{red}{e^{\lambda a}|n\rangle}.
\end{align}
The annihilation operators satisfy
$$ a^k|n\rangle = \left(\sqrt{n-k+1}\ldots\sqrt{n}\right)\,|n-k\rangle = \left(\frac{n!}{(n-k)!}\right)^{1/2}\,|n-k\rangle, $$
so
\begin{align}
\color{blue}{\langle m|e^{\lambda a^\dagger}} &\color{blue}{= \sum_{l=0}^m \frac{\lambda^l}{l!}\left(\frac{m!}{(m-l)!}\right)^{1/2}\langle m-l|}, & \color{red}{e^{\lambda a}|n\rangle} &\color{red}{= \sum_{k=0}^n \frac{\lambda^k}{k!} \left(\frac{n!}{(n-k)!} \right)^{1/2} |n-k\rangle}.
\end{align}
Then
\begin{align}
\langle m|e^{tx}|n\rangle &= e^{\frac{\lambda^2}{2}}\color{blue}{\sum_{l=0}^m \frac{\lambda^l}{l!}\left(\frac{m!}{(m-l)!}\right)^{1/2}}\color{red}{\sum_{k=0}^n \frac{\lambda^k}{k!}\left(\frac{n!}{(n-k)!}\right)^{1/2}}\underbrace{\color{blue}{\langle m-l}|\color{red}{n-k\rangle}}_{=\delta_{l,k+m-n}}\\
&= e^{\frac{\lambda^2}{2}}\sum_{k=0}^n \frac{(m!\,n!)^{1/2}\;\lambda^{2k+m-n}}{(k+m-n)!\,(n-k)!\,k!}\tag{$\dagger$}\label{$\dagger$}\\
&= \left(\frac{n!}{m!}\right)^{1/2}\,\lambda^{m-n}e^{\frac{\lambda^2}{2}}\sum_{k=0}^n {m\choose n-k} \frac{\left(\lambda^2\right)^k}{k!}.
\end{align}
Using the definition of a generalized Laguerre polynomial,
$$ L_n^{(\alpha)}(x) = \sum_{k=0}^n (-1)^k { n+\alpha \choose n-k } \frac{x^k}{k!}, $$
we find that the generating function is
$$ \langle m|e^{tx}|n\rangle = \left(\frac{n!}{m!}\right)^{1/2}\,\lambda^{m-n}e^{\frac{\lambda^2}{2}}\,L_n^{(m-n)}\left(-\lambda^2\right).$$
If $n = m$, then this reduces to $ \langle n|e^{tx}|n\rangle = e^{\frac{\lambda^2}{2}}L_n\left(-\lambda^2\right)$.
Matrix Elements
The matrix elements $\langle m|x^p|n\rangle$ can be calculated by expanding $\langle m|e^{tx}|n\rangle$ and equating coefficients of $t^p$. Writing $ [f(t)]_{t^p}$ for the coefficient of $t^p$ in the expansion of $f(t)$, ($\ref{$\dagger$}$) gives
\begin{align}
\frac{1}{p!}\langle m|x^p|n\rangle &= \left[\color{orange}{e^{\frac{\lambda^2}{2}}}\sum_{k=0}^n \frac{(m!\,n!)^{1/2}\;\lambda^{2k+m-n}}{(k+m-n)!\,(n-k)!\,k!}\right]_{t^p}\\
&= \left( \frac{\hbar}{2m\omega} \right)^{p/2}\left[\color{orange}{\sum_{l=0}^\infty \frac{1}{l!}\left( \frac{\lambda^2}{2} \right)^l} \sum_{k=0}^n \frac{(m!\,n!)^{1/2}\;\lambda^{2k+m-n}}{(k+m-n)!\,(n-k)!\,k!}\right]_{\lambda^p}\\
&= \left( \frac{\hbar}{2m\omega} \right)^{p/2} \sum_{k=0}^N \frac{(m!\,n!)^{1/2}}{2^{(p+n-m)/2-k}\left(\frac{p+n-m}{2}-k\right)!\,(k+m-n)!\,(n-k)!\,k!}\tag{$\star$}\label{$\star$}
\end{align}
where $ N = \text{min}(n,\frac{p+n-m}{2})$. Mathematica can evaluate this sum for us to give the following expression involving the hypergeometric function:
$$ \langle m|x^p|n\rangle = \left( \frac{\hbar}{m\omega} \right)^{p/2}\cdot \frac{2^{(m-n)/2}p!}{2^p\left( \frac{p+n-m}{2} \right)!} \cdot {m \choose n} \cdot \left( \frac{n!}{m!} \right)^{1/2}{}_2F_1\left( -n,\tfrac{m-n-p}{2};1+m-n;2 \right), $$
(yuck!) which pretty much rules out the possibility that there is a nice closed-form expression.
For the most interesting case $n=m$, we get a fairly nice result (replacing $p\rightarrow 2p$ because expectations of odd moments of $x$ vanish):
$$ \boxed{\langle n| x^{2p} |n\rangle = \left(\frac{\hbar}{m\omega}\right)^p\cdot\frac{(2p)!}{4^p\,p!}\cdot{}_2F_1(-n,-p;1;2).} $$
The equivalent form
$$ \boxed{\langle n| x^{2p} |n\rangle = \left(\frac{\hbar}{m\omega}\right)^p\cdot\frac{(2p)!}{4^p\,p!}\cdot\sum_{k=0}^\text{min(n,p)} {n \choose k}{p \choose k} 2^k,}$$
which follows from ($\ref{$\star$}$), is probably more useful (e.g. for finding asymptotics). It's also possible to use this form to work out the expectation of $x^{2p}$ in the canonical ensemble, which has a neat answer:
$$ \langle x^{2p} \rangle_\beta = \frac{1}{Z}\sum_{n=0}^\infty \langle n|x^{2p}|n\rangle e^{-\beta E_n} = \left( \frac{\hbar}{m\omega} \right)^p\cdot \frac{(2p)!}{4^p\,p!}\cdot \coth^p \left( \frac{\beta \hbar\omega}{2} \right) $$
where $\beta = 1/(k_BT)$ and $Z = \sum e^{-\beta E_n}$ is the partition function.
