# A Gaussian integral containing a logarithm with applications to finance.

This question is closely related to An integral involving a Gaussian and a logarithm. .

Let $$\sigma >0$$, $$A >0$$ and $$\lambda >0$$. Consider the following integral:

$$$${\mathcal I}(\lambda,A,\sigma) := \int\limits_{-\infty}^\infty \log\left[ 1+ \lambda \cdot (A \cdot e^{\sigma \xi} - 1)\right] \cdot \frac{\exp(-\xi^2/2)}{\sqrt{2\pi}} d\xi$$$$

In mathematical finance, this integral represents the expected value of the logarithm of wealth in a portfolio with a percentage $$\lambda$$ being invested in a risky asset with expected return $$\log(A)$$ and volatility $$\sigma$$ and a percentage $$(1-\lambda)$$ being invested in a risk-free asset that earns an interest rate equal to zero. Now, by rearranging the terms in the logarithm and then expanding the logarithm in a series , integrating term by term and then regularizing the result we have obtained a following closed form for the integral in question:

$$\begin{eqnarray} {\mathcal I}(\lambda,A,\sigma) = \log(1-\lambda + \lambda A) + \sum\limits_{p=1}^\infty \frac{(\sigma^2/2)^p}{p!} \cdot \left. \left(A \cdot \sum\limits_{k=1}^{2p-1} \{ \begin{array}{c} 2p-1 \\ k \end{array} \} \frac{(-1)^{k-1} k!}{(1+A)^{k+1}}\right)\right|_{A= A \cdot \frac{\lambda}{1-\lambda}} \tag{ii}\label{ii} \end{eqnarray}$$

Despite appearances the sum in the right hand side of \eqref{ii} converges quite fast as I have checked numerically below. There, I randomly sample some typical model parameters and then firstly compute the integral from the definition and secondly from the closed form by truncating the relevant sum at $$p=3$$ terms. As you can see the last term in the sum is at worst of the order of $$10^{-15}$$ which is negligible for practical purposes.

In[1692]:=
M = 3;
T2 = Table[
A Sum[ StirlingS2[2 p - 1, k] (-1)^(k - 1) k!/(1 + A)^(k + 1), {k,
1, 2 p - 1}], {p, 1, M}];
Clear[Terms];
Terms [A_, ll_, s_] =
Table[Simplify[( T2[[p]] /. A :> A ll/(1 - ll))] (s^2/2)^p/p!, {p, 2,
M}];

s = (RandomInteger[{0, 20}]/100)/Sqrt[256];
ll = RandomReal[{0, 1}, WorkingPrecision -> 50];
A = RandomReal[{0, 5}, WorkingPrecision -> 50];
N[NIntegrate[
Log[1 + ll (A Exp[s xi] - 1)] Exp[-xi^2/2]/
Sqrt[2 Pi], {xi, -Infinity, Infinity}, WorkingPrecision -> 14], 10]

N[Log[1 - ll] +
NIntegrate[
Log[1 + A ll/(1 - ll) Exp[s xi]] Exp[-xi^2/2]/
Sqrt[2 Pi], {xi, -Infinity, Infinity},
WorkingPrecision -> 14], 10]

(*This comes from expanding the Log[] in the integral above to second \
order in xi and integrating term by term.*)
Join[{N[Log[1 - ll + A ll] + ( (s^2) A (ll - ll^2) )/(
2 (1 - ll + A ll)^2), 10] }, Terms[A, ll, s]]

Out[1698]= 0.4001260648

Out[1699]= 0.4001260647

Out[1700]= {0.4001260652,
4.63239352341047007113417223619846949288159549561*10^-12, \
-1.63137128039387624683354129276554122034541505684*10^-15}


Now my question is twofold. Firstly can we prove that the series in \eqref{ii} do converge? Secondly can we use the result \eqref{ii} to find the optimal asset allocation $$\lambda$$ that maximizes the growth of log wealth in the long term?

In here we answer the second part of the question, namely how to find the optimal asset allocation that maximizes log-Wealth in the long term. First of all let us stress that this is a practical problem and as such all the parameters in question have a certain order of magnitude. Both $$\mu := \log(A)$$ and $$\sigma^2$$ are given in terms of percent per year (typically $$\mu = 0.04/256$$ and $$\sigma^2 = 0.03/256$$ where $$256$$ is the number of trading days per year). Therefore both $$\mu$$ and $$\sigma^2$$ are of the same order of magnitude and their ratio $$\theta := \mu/\sigma^2 = O(1)$$ is of the order of unity. We have:

$$\begin{eqnarray} {\mathcal I}(\lambda,A,\sigma) &=& \log\left[ 1- \lambda + \lambda A\right] + \frac{\sigma^2}{2} \cdot \frac{\lambda (1-\lambda) A}{(1-\lambda + \lambda A)^2} + O(\sigma^4) \\ &=& \frac{\sigma^2}{2} \left( \lambda\cdot (1+2 \theta)- \lambda^2 \right) + O(\sigma^2) \end{eqnarray}$$ where in the second line above we replaced $$A \rightarrow \exp(\theta \sigma^2)$$ and then expanded the whole expression in powers of $$\lambda$$ and neglected all terms containing powers of $$\sigma$$ higher than two. Clearly the expression above has a unique maximum at $$\lambda = \lambda^{(*)}$$ and therefore we write: $$\begin{eqnarray} \lambda^{(*)} = \left(\theta + \frac{1}{2}\right) + O(\sigma^2) \end{eqnarray}$$

What we do now we go to higher orders of the volatility $$\sigma$$. Since the calculations are very cumbersome rather than doing them on paper it is better to employ a Computer Algebra system like Mathematica. The algorithm is as follows.

1. Take a dummy variable and set it equal to two. $$ord = 2$$,
2. Define $${\mathcal S}(\theta,\sigma, {\mathfrak A}) := (\theta+1/2) + {\mathfrak A} \sigma^2$$. We tern that function a replacement function.
3. Take $${\mathcal J}(\lambda,\theta,\sigma) := d/(d \lambda ) {\mathcal I}(\lambda,\exp(\theta \sigma^2), \sigma)$$ and expand it in a series in the variable $$\lambda$$ to the order $$2 \cdot ord$$ inclusive.
4. In the expansion above replace $$\lambda$$ by $$\lambda \rightarrow {\mathcal S}(\theta,\sigma, {\mathfrak A})$$ and then expand the result and neglect all powers of the volatility that are strictly bigger than $$2 \cdot ord$$.
5. Simplify the resulting expression , set it to zero and solve it for $${\mathfrak A}^{(*)}(\theta)$$ -- it turns out that the expression in question is always linear in $${\mathfrak A}$$ and as such can be solved.
6. Update the replacement function as follows $${\mathcal S}(\theta,\sigma, {\mathfrak A}) := {\mathcal S}(\theta,\sigma, {\mathfrak A}^{(*)}(\theta)) + {\mathfrak A} \sigma^{2 ord}$$.
7. Increase $$ord$$ by unity ($$ord \rightarrow ord+1$$) and go back to 3.

The algorithm has been implemented in the code below:

th =.; s =.; rt =.; rt = th s^2; A =
1 + rt + rt^2/2; A1 =.; ll =.; s =.; th =.; M = 8;
Clear[myUtility];
myUtility[A_, ll_] :=
Log[1 - ll + A ll] +
Sum[(A ll (1 - ll) Sum[
StirlingS2[2 p - 1, k] ((-1)^(k - 1) k! (1 - ll)^(
k - 1))/(1 - ll + A ll)^(k + 1), {k, 1, 2 p - 1}]) (s^2/2)^p/
p!, {p, 1, M}];
mySol[th_, s_, A1_] := (th + 1/2) + A1 s^2;

For[ord = 2, ord <= 4, ord++,
eX = Expand[
Normal[Series[
D[myUtility[Sum[rt^j/j!, {j, 0, ord}], ll], ll], {ll, 0,
2 ord}]] /. ll :> mySol[th, s, A1]] /. s^n_ :> 0 /; n > 2 ord;
eX = Simplify[eX];
A11 = A1 /. Simplify[First@Solve[eX == 0, A1]];
Print[{ord, eX, A11}];

mySol[th_, s_, A1_] = mySol[th, s, A11] + A1 s^(2 ord);
];

{2,1/4 s^4 (-4 A1+th-4 th^3),1/4 (th-4 th^3)}

{3,-(1/48) s^6 (48 A1+5 th-32 th^3+48 th^5),-((5 th)/48)+(2 th^3)/3-th^5}

{4,-(1/32) s^8 (32 A1+th (1-4 th^2)^2 (-3+14 th^2)),-(1/32) th (1-4 th^2)^2 (-3+14 th^2)}


As such we can write that the optimal investment allocation reads:

$$\begin{eqnarray} &&\lambda^{(*)} = \\ &&\left(\theta+\frac{1}{2}\right) + \\ &&\frac{1}{4} (\sigma)^2 \left(\theta-4 \theta^3\right)+ \\ &&(\sigma)^4 \left(-\theta^5+\frac{2 \theta^3}{3}-\frac{5 \theta}{48}\right)+ \\ &&-\frac{1}{32} (\sigma)^6 \theta \left(1-4 \theta^2\right)^2 \left(14 \theta^2-3\right)+\\ &&O(\sigma^8) \end{eqnarray}$$

Now we have verified the formula above numerically in a following way. We generated five hundred random pairs $$(\mu,\sigma^2) = (i,j)/(25600)$$ where $$(i,j)$$ are iid integer random variables from one to twenty. Now for each $$(\mu,\sigma^2)$$ we firstly obtained the optimal allocation by maximizing the utility function $${\mathfrak I}(\lambda,\exp(\mu),\sigma)$$ numerically and secondly we used the the power expansion formula above. After that we computed the mean and the standard deviation of the set of relative differences between the two numbers all that over that whole sample. The resulting figures were always of the order of $$10^{-8}$$ and $$10^{-7}$$ respectively.