I'm playing around with dynamic programming and need to calculate a multidimensional integral $E[V(W)]$ where we assume $W$ has a log normal distribution. I was looking at the following example in this pdf, section 9.4.3 on page 83.To give some background (I summarize from the section): The example is from economics and is about asset allocation. Assume $R$ is a return vector of dimension $n$ and is log-normally distrbuted, i.e. $\log(R) = (\log(R_1),\dots, \log(R_n))$ has a multivariate normal distribution with given mean and covariance matrix, i.e. $\log(R)\sim\mathcal{N}((\mu-\frac{\sigma^2}{2})\Delta t,(\Lambda \Sigma\Lambda)\Delta t)$. The exact structure is not that important. $\Sigma$ is the correlation matrix and $\Lambda$ is simple a diagonal matrix with the standard deviation $\sigma_1,\dots,\sigma_n$ on its diagonal. Using Choleski decomposition we can write $\Sigma = LL^T$. Then one has$$\log(R_1) = (\mu_1-\frac{\sigma^2_1}{2})\Delta t + (L_{11}z_1)\sigma_1\sqrt{\Delta t}$$$$\log(R_2) = (\mu_2-\frac{\sigma^2_2}{2})\Delta t + (L_{21}z_1+ L_{22}z_2)\sigma_2\sqrt{\Delta t}$$and so on, where $z_i$ are independent normal distribution. So that we have$$R_i=\exp{((\mu_i-\frac{\sigma^2_i}{2})\Delta t+\sigma_i\sqrt{\Delta t}\sum_{j=1}^iL_{ij}z_j)}$$Let for simplicity no $\Delta t=1$ then we are interested in the quantity$$ W_{t+1}= W_t(R_f(1-e^Tx_t) + \sum_{i=1}^n\exp{((\mu_i-\frac{\sigma_i^2}{2})+\sigma^2\sum_{j=1}^iL_{ij}z_j)}x_{ti}) $$where $e$ is a $n$ dimensional vector of $1$ and $x_t$ is some $n$ dimensional vector. For a given function $V$ the conditional expectatino of $V(W_{t+1})$ given $W_t, x_t$ can be calculated using Gauss Hermite quadrature $$\sum_{k_1,\dots,k_n=1}^m w_{k_1}\cdot\cdot\cdot w_{k_n} V\left(W_t(R_f(1-e^Tx_t) + \sum_{i=1}^n\exp{((\mu_i-\frac{\sigma_i^2}{2})+\sigma^2\sum_{j=1}^iL_{ij}q_{k_j})}x_{ti})\right)$$ where $w_{k_i}$ are the Gauss-Hermite weights and $q_i$ the corresponding nodes. My question is how can the above sum of sums be efficiently implemented in python? The exponential part can be precalculated using a cummulative sum if I'm not wrong.
How can this multidimensional integral be efficiently implemented in python using Gauss-Hermite quadrature
python;numerical;integration
The summation over $q_{k_i}$ in $V$ is independent of the sums over the products $w_{k_1}w_{k_2}\cdots w_{k_n}$. So we can calculate $V$ first. To avoid double evaluations of same operations we can further decompose the equation:Step 1: precalculate $\tilde L_{ij}$\begin{equation}  \tilde L_{ij} = \sigma^2 L_{ij}q_{k_j}\end{equation}Step 2: precalculate $E_i$\begin{equation}  E_i = \exp\left(\left(\mu_i-\frac{\sigma_i^2}{2}\right) +\sum_{j=1}^i\tilde L_{ij} \right)\end{equation}Step 3: calculate $X_n$\begin{equation}  X_n = W_t\left(R_f(1-e^Tx_t)+\sum_i^n x_{ti}E_i\right)\end{equation}For the sums over the products $w_{k_1}w_{k_2}\cdots w_{k_n}$ we can use numpy.linalg.multi_dot(*w), where w=[w1,w2,...] is the list of all vectors $w_{k_i}$. In a final step we multiply this result by $V(X_n)$.