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.

  • $\begingroup$ Write it in C/Cython/Fortran and use the FFI? $\endgroup$ Jul 13, 2017 at 6:52
  • $\begingroup$ @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 $\endgroup$
    – math
    Jul 13, 2017 at 6:58
  • $\begingroup$ Does $V$ depend on any of the $k_i$? $\endgroup$
    – Feodoran
    Jul 13, 2017 at 10:34
  • $\begingroup$ 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, ...]? $\endgroup$
    – Feodoran
    Jul 13, 2017 at 12:02
  • $\begingroup$ @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. $\endgroup$
    – math
    Jul 13, 2017 at 12:27

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

  • $\begingroup$ 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$. $\endgroup$
    – math
    Jul 14, 2017 at 16:40
  • $\begingroup$ 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. $\endgroup$
    – Feodoran
    Jul 16, 2017 at 16:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.