Asymptotic value of a multivariate integral The following question is a simple case of a type of problem that occurs in combinatorial enumeration problems.
Define 
$$F(x_1,\ldots,x_n) = \frac{1}{(2\pi)^{n/2}}\exp\biggl( -\frac12\sum_{j=1}^n x_j^2 + i\,n^{-1/4} \sum_{j=1}^n x_j^3 \biggr), $$
where $i$ is the imaginary unit.
Let $\epsilon\gt 0$ be small enough. Let
$$I_0(n) = \int_{-n^\epsilon}^{n^\epsilon}\cdots \int_{-n^\epsilon}^{n^\epsilon}
  F(x_1,\ldots,x_n) ~ dx_1\cdots dx_n. $$
We can estimate $I_0(n)$ as $n\to\infty$ by factoring
$$F(x_1,\ldots,x_n) = \prod_{j=1}^n \frac{1}{\sqrt{2\pi}} 
  \exp\bigl(-x_j^2/2+in^{-1/4}x_j^3\bigr),$$
which separates the integral into a product of $n$ 1-dimensional integrals.  This very easily gives
$$I_0(n) = \exp\bigl( -n^{1/2}/2 + O(n^{-3})\bigr).$$
So far so good.  But now let $A=A(n)$ be a symmetric positive definite matrix. We can assume it is pretty nice, say all eigenvalues bounded between two positive constants independent of $n$.  Now define
$$G(x_1,\ldots,x_n) = \frac{1}{(2\pi)^{n/2}}\exp\biggl( -\frac12 \mathbf{x}^TA\mathbf{x} + i\,n^{-1/4} \sum_{j=1}^n x_j^3 \biggr). $$
How do we estimate
$$I_1(n) = \int_{-n^\epsilon}^{n^\epsilon}\cdots \int_{-n^\epsilon}^{n^\epsilon}
  G(x_1,\ldots,x_n) ~ dx_1\cdots dx_n? $$
The source of the difficulty is that $\int G$ is exponentially smaller than $\int |G|$, so as soon as you approximate the integrand the answer goes away.  What to do?
Note that in this case (and very commonly in practice) symmetry implies that $I_1(n)$ is real.  So we can discard the imaginary part of the integrand and integrate only the real part, which is
$$ \mathfrak{R}G(x_1,\ldots,x_n) = \frac{1}{(2\pi)^{n/2}}\exp\biggl( -\frac12 \mathbf{x}^TA\mathbf{x}\biggr)\cos\biggl(n^{-1/4} \sum_{j=1}^n x_j^3 \biggr). $$
Not sure that helps.
ADDED 2 Jan 2018.  Mikhail Isaev and I now have a method that can precisely estimate such integrals in many general cases.  It isn't completely written up yet.
 A: I'll consider the integrals extended from $-\infty$ to $\infty$ and show that the integral of $G$ is exponentially small compared to the 
integral of $|G|$ -- precisely, the oscillating integral is smaller by a factor of $C_0\exp(-C \sqrt{n})$ for some positive constants $C_0$ and $C$ 
that depend only on the eigenvalues of the quadratic form $A$.   One can get the same estimate for the truncated integrals.  It may be possible to refine this to get asymptotics but 
I have not thought about that.  
The problem asks for an estimate for 
$$ 
I = \frac{1}{(2\pi)^{n/2}} \int_{-\infty}^{\infty} \cdots \int_{-\infty}^{\infty} \exp\Big( -\frac{1}{2} A(x_1,\ldots,x_n) +\frac{i}{n^{\frac 14}} \sum_j x_j^3\Big) dx_1\ldots dx_n.
$$
where $A$ is a positive definite quadratic form.  We wish to show that $I$ is small compared to the same integral without the oscillating term; this is 
$$ 
\frac{1}{(2\pi)^{n/2} } \int_{-\infty}^{\infty}\cdots \int_{-\infty}^{\infty} \exp\Big(- \frac 12 A(x_1,\ldots,x_n )\Big) dx_1 \ldots dx_n = (\text{det} A)^{-\frac 12}.
$$ 
Think of the integrals in $I$ as contour integrals being integrated on the real axis.  The idea is to replace the integrals on the real axis by 
integrals along the line from $-\infty + i\alpha$ to $\infty +i\alpha$ for some real number $\alpha$ to be chosen carefully.   Thus 
$$ 
I =\frac{1}{(2\pi)^{n/2}} \int_{-\infty+i\alpha}^{\infty+i\alpha} \cdots \int_{-\infty+i\alpha}^{\infty+i\alpha}\exp\Big(-\frac12 A(z_1,\ldots,z_n) + \frac{i}{n^{\frac 14}} 
\sum z_j^3 \Big) dz_1 \ldots dz_n,
$$ 
and writing now $z_j =x_j +i\alpha$ this equals 
$$
\frac{1}{(2\pi)^{n/2}} \int_{-\infty}^{\infty}\cdots \int_{-\infty}^{\infty} \exp\Big( -\frac 12 A(x_1+i\alpha,\ldots, x_n+i\alpha) + \frac{i}{n^{\frac 14}} 
\sum_j (x_j+i\alpha)^3\Big) dx_1\ldots dx_n.
$$ 
We now estimate the integral above by just taking the absolute value of the integrand (and choosing $\alpha$ carefully).  The integrand is in modulus 
$$ 
\exp \Big( -\frac 12 A(x_1,\ldots, x_n) + \frac{\alpha^2 }{2} A(1,\ldots, 1) +\alpha^3 n^{\frac 34} -\frac{3\alpha}{n^{\frac 14} } \sum_j x_j^2\Big). 
$$ 
We will take $\alpha =\beta/n^{\frac 14}$ for a suitably small positive constant $\beta$.  Since the eigenvalues of $A$ are bounded between two positive constants independent of $n$, we have that $A(1,\ldots,1) \le 2C_1 n$ for some positive constant $C_1$, and that $3\sum_j  x_j^2 \ge \frac{C_2}{2} A(x_1,\ldots, x_n)$ for some positive constant $C_2$.
Thus the quantity above is 
$$ 
\le \exp\Big(C_1 \beta^2  \sqrt{n} + \beta^3 - \frac{1}{2} A(x_1,\ldots,x_n) \Big(1+\frac{C_2\beta}{\sqrt{n}}\Big) \Big).
$$ 
Integrating this, we find that 
$$ 
I \le \frac{\exp(C_1 \beta^2\sqrt{n} +\beta^3)}{(2\pi)^{n/2}} \int_{\infty}^{\infty}\cdots\int_{-\infty}^{\infty} \exp\Big(-\frac 12 A(x_1,\ldots,x_n) \Big(1+\frac{\beta C_2}{\sqrt{n}}\Big) \Big) dx_1 \ldots dx_n 
$$
which is readily seen to be 
$$ 
 \exp(C_1 \beta^2\sqrt{n} +\beta^3) (\text{det} A)^{-\frac 12} \Big( 1 +\frac{\beta C_2}{\sqrt{n}}\Big)^{-n/2}.
 $$ 
In other words we have shown that $I$ is smaller than the trivial bound by a factor of 
 $$ 
 \exp(C_1 \beta^2\sqrt{n} +\beta^3) \Big( 1 +\frac{\beta C_2}{\sqrt{n}}\Big)^{-n/2} \le 
 \exp(C_1 \beta^2 \sqrt{n} +\beta^3 - C_3 \beta\sqrt{n} \Big) 
 $$ 
 for a suitable positive constant $C_3$.  By choosing $\beta$ small enough, we find that this is $\le C_0 \exp(-C\sqrt{n})$ 
 for some positive constants $C_0$ and $C$.  This finishes the proof.
