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$ – Chris Rackauckas Jul 13 '17 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 '17 at 6:58
  • $\begingroup$ Does $V$ depend on any of the $k_i$? $\endgroup$ – Feodoran Jul 13 '17 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 '17 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 '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}$

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

| cite | improve this answer | |
  • $\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 '17 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 '17 at 16:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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