What is the expected value of $\prod\limits_{i=1}^n (X_i +1)$?
Following is a closed form expression for the expected value in terms of hypergeometric functions. Unluckily, I am not able to deduce the asymptotic behavior of it for large $n$. I hope someone can followup and complete the task.
For any $m = (m_1, m_2,\ldots,m_n) \in \mathbb{N}^n$ and $t = (t_1,t_2,\ldots,t_n) \in \mathbb{R}^n$, let
$$ |m| = \sum_{i=1}^n m_i,\quad m! = \prod_{i=1}^n m_i!,\quad t \cdot X = \sum_{i=1}^n t_i X_i\quad\text{ and }\quad t \cdot m = \sum_{i=1}^n t_i m_i $$ It is easy to see the probability that integer $1$ appear $m_1$ times, integer $2$ appear $m_2$ times and so on is given by the formula: $$\mathbb{P}\big[X_i = m_i, \forall i \big] = \begin{cases} \frac{n!}{n^n m!}, & |m| = n\\ 0, & \text{ otherwise } \end{cases} $$ Using multinomial expansion, this leads to $$ \mathbb{E} \left[e^{t\cdot X}\right] = \sum_{m \in \mathbb{N}^n, |m| = n} \frac{n!}{n^n m!} e^{t\cdot m} = \sum_{m \in \mathbb{N}^n, |m| = n} \frac{n!}{m!} \prod_{i=1}^n \left(\frac{e^{t_i}}{n}\right)^{m_i} = \left(\frac1n\sum_{i=1}^n e^{t_i}\right)^n $$ For any $1 \le k \le n$, By taking partial derivatives of above expansion, we obtain $$\begin{align} \mathbb{E} \left[\prod_{i=1}^k X_i\right] &= \left.\frac{\partial^k}{\partial t_1\partial t_2\cdots \partial t_k} \mathbb{E}\left[ e^{t\cdot X}\right]\right|_{t=0}\\ &= \left. \prod_{i=1}^k \left((n+1-i)\frac{e^{t_i}}{n}\right)\times \left(\frac1n\sum_{i=1}^n e^{t_i}\right)^{n-k}\right|_{t=0}\\ &= \prod_{i=0}^{k-1} \left(1-\frac{i}{n}\right) = \frac{n!}{n^k(n-k)!}\end{align}$$ Using symmetry, we can express the expectation value at hand in terms of hypergeometric function: $$\mathbb{E}\left[\prod_{i=1}^n (X_i+1)\right] = 1 + \sum_{k=1}^n \binom{n}{k}\mathbb{E}\left[\prod_{i=1}^k X_i\right] = \sum_{k=0}^n \frac{(-n)^{\bar k}(-n)^{\bar k}}{k!}\left(\frac{1}{n}\right)^n\\ = {}_2F_0\left(-n,-n;\frac1n\right) = (-n)^{-n}U(-n,1,-n) $$
where
- $(\alpha)^{\bar k} = \prod_{i=0}^{k-1}(\alpha + i)$ is the rising factorial.
- ${}_2F_0(a,b;z)$ is one of the generalized hypergeometric function
- $U(a,b,z)$ is Tricomi's (confluent hypergeometry) function
Edited. I added a better estimate, see $\text{(3)}$ below.
Here is an asymptotic computation using achille hui's result. Writing $k = pn$ and $q = 1-p$, it follows that
$$ \lim_{n\to\infty} \frac{1}{n}\log \left( \binom{n}{k}\frac{n!}{n^k(n-k)!}\right) = -\left( p + p \log p + 2q \log q \right) =: I(p). $$
More precisely, from the quantitative version of the Stirling's approximation, we have the following uniform estimate
$$ c \frac{e^{I(p)n}}{\sqrt{pq^2 n}} \leq \binom{n}{k}\frac{n!}{n^k(n-k)!} \leq C \frac{e^{I(p)n}}{\sqrt{pq^2 n}} \tag{1}$$
for some constants $0 < c < C $, which holds for any $1 \leq k \leq n-1$.
By computation, we can check that $I(p)$ is maximized when $p^* = \frac{3-\sqrt{5}}{2}$. From this together with the estimate above, it is not hard to check that
$$ \lim_{n\to\infty} \frac{1}{n}\log\mathbb{E}\left[ \prod_{i=1}^{n}(1+X_i) \right] = I(p^*) \approx 0.58045763886910174320. \tag{2} $$
Added. In fact we have a better estimate
$$ a_n := \mathbb{E}\left[ \prod_{i=1}^{n}(1+X_i) \right] \sim c e^{I(p^*)n} \tag{3}$$
as $n\to\infty$, where $c = \frac{5^{1/4}+5^{-1/4}}{2} \approx 1.082044543\cdots$.
A heuristic argument is as follows: Plugging $k = p^*n + x_k\sqrt{n}$, it follows that $a_n e^{-I(p^*)n}$ is roughly the following Riemann sum
$$ a_n e^{-I(p^*)n} \approx \frac{1}{(1-p^*)\sqrt{2\pi p^*}} \sum_k e^{\frac{I''(p^*)}{2}x_k^2} \Delta x_k. $$
As $n\to\infty$, this converges to the integral
$$ \frac{1}{(1-p^*)\sqrt{2\pi p^*}} \int_{-\infty}^{\infty} e^{\frac{I''(p^*)}{2}x^2} \, dx = \frac{1}{(1-p^*)\sqrt{-p^* I''(p^*)}} = \frac{5^{1/4}+5^{-1/4}}{2}. $$
Turning this into a rigorous argument requires some efforts, but it is rather a matter of technicality.
Proof of $\text{(3)}$. Let $\delta > 0$ be sufficiently small so that $p^*\pm\delta \in (0, 1)$. Also for each $1 \leq k \leq n-1$ define $c_{n,k}$ as the ratio
$$ c_{n,k} = \frac{\binom{n}{k}\frac{n!}{n^k(n-k)!}}{e^{I(k/n)n}/\sqrt{\frac{k}{n}(1-\frac{k}{n})^2 n}}. $$
The estimate $\text{(1)}$ shows that $c_{n,k}$ is uniformly bounded away both from $0$ and from $\infty$. Moreover, Stirling's approximation tells that $c_{n,k} \to 1/\sqrt{2\pi} $ as long as $k/n$ stays uniformly away both from $0$ and from $1$.
Now let $\alpha = I(p^*) - \max\{ I(p^*-\delta), I(p^*+\delta)\} > 0$. Then truncating the sum using the window of size $n\delta$ around $p^*n$, we have
$$ a_n e^{-I(p^*)n} = \sum_{\left| \frac{k}{n} - p^* \right| \leq \delta} \frac{c_{n,k}}{\sqrt{\frac{k}{n}(1-\frac{k}{n})^2 n}} e^{\left( I(\frac{k}{n}) - I(p^*) \right)n} + \mathcal{O}(n e^{-\alpha n})$$
Now define $x_k$ by the relation $k = p^*n + x_k\sqrt{n}$. By the Taylor expansion, we know that
\begin{align*} I(\tfrac{k}{n}) - I(p^*) &= \frac{I''(p^*)}{2}(p^* - \tfrac{k}{n})^2 \left(1 + \mathcal{O}(p^* - \tfrac{k}{n}) \right) \\ &= \frac{I''(p^*) x_k^2}{2n} \left(1 + \mathcal{O}(p^* - \tfrac{k}{n}) \right) \end{align*}
Plugging this back and denoting by $C > 0$ the implicit bound of $\mathcal{O}(p^* - \tfrac{k}{n})$ in the expansion above, we obtain the following simple bound
\begin{align*} &\left( \min_{\left| p - p^* \right| \leq \delta} \frac{1}{\sqrt{p(1-p)^2}} \right) \sum_{\left| \frac{k}{n} - p^* \right| \leq \delta} \frac{c_{n,k}}{\sqrt{n}} e^{\frac{1}{2} (1 +C\delta ) I''(p^*) x_k^2} + \mathcal{O}(n e^{-\alpha n}) \\ &\hspace{3em} \leq a_n e^{-I(p^*)n} \\ &\hspace{6em} \leq \left( \max_{\left| p - p^* \right| \leq \delta} \frac{1}{\sqrt{p(1-p)^2}} \right) \sum_{\left| \frac{k}{n} - p^* \right| \leq \delta} \frac{c_{n,k}}{\sqrt{n}} e^{\frac{1}{2} (1 -C\delta )I''(p^*) x_k^2} + \mathcal{O}(n e^{-\alpha n}) \end{align*}
Letting $n\to\infty$, dominated convergence theorem yields
\begin{align*} &\left( \min_{\left| p - p^* \right| \leq \delta} \frac{1}{\sqrt{p(1-p)^2}} \right) \int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi}} e^{\frac{1}{2} (1 +C\delta ) I''(p^*) x^2} \, dx \\ &\hspace{3em} \leq \liminf_{n\to\infty} a_n e^{-I(p^*)n} \leq \limsup_{n\to\infty} a_n e^{-I(p^*)n} \\ &\hspace{6em} \leq \left( \max_{\left| p - p^* \right| \leq \delta} \frac{1}{\sqrt{p(1-p)^2}} \right) \int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi}} e^{\frac{1}{2} (1 -C\delta ) I''(p^*) x^2} \, dx \end{align*}
Letting $\delta \downarrow 0$ proves the desired result $\text{(3)}$.