Solving integral over the space of two lognormally distributed variables

I am interested in the closed form solution to the following problem:

$$\int_a^b\int_0^{cx+d}(x+e)f(x)f(y)dydx$$, where $$f(.)$$ is the pdf of the lognormal distribution with mean $$0$$ and variance $$\sigma^2$$. $$a,b,c,d,e$$ are constants. Say $$F(.)$$ is the corresponding cdf.

I got as far as $$\int_a^b\int_0^{cx+d}(x+e)f(x)f(y)dydx=\int_a^bxf(x)F(cx+d)dx+e\int_a^bf(x)F(cx+d)dx$$.

But I am not sure if generating a closed form solution to the above integral is possible. Given this, I am happy to approximate the lognormal cdf $$F(x)$$ to $$G(x)=\frac{x^\alpha}{1+x^\alpha}$$, with pdf $$g(x)=\frac{\alpha x^{\alpha-1}}{(1+x^\alpha)^2}$$, where $$\alpha>1$$ and $$\alpha$$ is smaller the larger the $$\sigma$$. The idea for this approximation came from https://www.jstor.org/stable/3621676?seq=1#metadata_info_tab_contents (Mihálykó, Csaba, and Tibor Blickle. “84.48 On the Approximation of the Lognormal Distribution.” The Mathematical Gazette, vol. 84, no. 500, 2000, pp. 308–311.).

Then my problem becomes one of finding closed form solutions to:

$$\int_a^bxg(x)G(cx+d)dx=\int_a^bx\frac{\alpha x^{\alpha-1}}{(1+x^\alpha)^2}\frac{(cx+d)^\alpha}{1+(cx+d)^\alpha}dx$$

and

$$\int_a^bg(x)G(cx+d)dx=\int_a^b\frac{\alpha x^{\alpha-1}}{(1+x^\alpha)^2}\frac{(cx+d)^\alpha}{1+(cx+d)^\alpha}dx$$.

But I am still lost as to how to exactly do this. Detailed steps would be appreciated, thanks!

• what you have is essentially an indefinite integral of a complicated algebraic function; the closed form answer is really complicated even for $\alpha=2$, I don't think this will be of any use for you; why not evaluate it numerically? – Carlo Beenakker Oct 12 '18 at 16:36
• @CarloBeenakker Are you aware of any "good enough" approximations to the lognormal distribution that would make finding closed forms easy? By "good enough", I mean similar in shape; I am approximating the distribution of earnings in the US: sc.cnbcfm.com/applications/cnbc.com/resources/files/2015/10/21/…) – Yay Oct 12 '18 at 16:46
• @CarloBeenakker I could evaluate the integral numerically, but my aim was to get at a neat closed form, if that was possible. – Yay Oct 12 '18 at 16:52

For $$c=0$$ there is a lengthy closed form, $$\int_a^b\int_0^{d}(x+e)f(x)f(y)dydx=\frac{1}{4} \left(\text{erf}\left(\frac{\ln d}{\sqrt{2} \sigma}\right)+1\right)$$ $$\qquad\qquad\times\left(e \left(\text{erf}\left(\frac{\ln b}{\sqrt{2} \sigma}\right)-\text{erf}\left(\frac{\ln a}{\sqrt{2} \sigma}\right)\right)+e^{\sigma^2/2} \left(\text{erf}\left(\frac{\sigma^2-\ln a}{\sqrt{2} \sigma}\right)-\text{erf}\left(\frac{\sigma^2-\ln b}{\sqrt{2} \sigma}\right)\right)\right)$$
For $$c \neq 0$$ it looks pretty hopeless...