What is the expected value of $\prod\limits_{i=1}^n (X_i +1)$? 
Let $S_n$ be a  sample of $n$ integers taken uniformly with replacement from $\{1,\dots, n\}$.  Let $X_i$ by the number of times the integer $i$ occurs in $S_n$. What is
  $$\mathbb{E} \left(\prod_{i=1}^n (X_i +1)\right)\;?$$

From numerical experiments, it seems that $$\lim_{n \to \infty}\frac1n\log_2\mathbb{E} \left(\prod_{i=1}^n (X_i +1)\right) \approx 0.83$$
 A: 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
A: 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)}$.
