How can this multidimensional integral be efficiently implemented in python using Gauss-Hermite quadrature

Question

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 distrïbuted, 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.

Answer 1

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)$.