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

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.

• Write it in C/Cython/Fortran and use the FFI? – Chris Rackauckas Jul 13 '17 at 6:52
• @ChrisRackauckas it's fine to have it python. My question is about the expression $\sum_{k_1=1}^m\sum_{k_2=1}^m\sum_{k_3=1}^m\dots\sum_{k_n=1}^m\dots$. How this sum of sum ... can be efficiently implemented – math Jul 13 '17 at 6:58
• Does $V$ depend on any of the $k_i$? – Feodoran Jul 13 '17 at 10:34
• So we can pull V before the sums and do $V\sum_{k_1\cdots k_n}w_{k_1}\cdots w_{k_n}$? Maybe you can use numpy.linalg.multi_dot(*w) where w=[w1, w2, w3, ...]? – Feodoran Jul 13 '17 at 12:02
• @Feodoran no you can't the nodes $q_{k_j}$ go in there. Sry, I did misread your last comment. I will delete my previous one. – math Jul 13 '17 at 12:27

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}$

$$\tilde L_{ij} = \sigma^2 L_{ij}q_{k_j}$$

### Step 2: precalculate $E_i$

$$E_i = \exp\left(\left(\mu_i-\frac{\sigma_i^2}{2}\right) +\sum_{j=1}^i\tilde L_{ij} \right)$$

### Step 3: calculate $X_n$

$$X_n = W_t\left(R_f(1-e^Tx_t)+\sum_i^n x_{ti}E_i\right)$$

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

• I think there is a problem in step 1. The $q_{k_j}$ vary because $k_j$ takes a possible range from $1\dots m$. – math Jul 14 '17 at 16:40
• I don't see the problem here, but I am not really a mathematician. Maybe you could build a simpler example equation and test/compare the issue with it. – Feodoran Jul 16 '17 at 16:21