• 1
• 2
• 3
• 4
• 5
• 6
• 8 ppauper
Posts: 70239
Joined: November 15th, 2001, 1:29 pm

### Re: Silly questions

so
$2x\,\mathrm{Heaviside}(x)+(1-2x)\,\mathrm{erfc}(-(x-1/2)/\epsilon)+2(x-1)\,\mathrm{Heaviside}(x-1)$

the area under that between 0 and 1 is I believe
$(-1/2+\epsilon^2)\,\mathrm{erf}(1/(2\epsilon))+1-\epsilon/\pi^{1/2}\exp(-1/(4\epsilon^2))$
What area when $\varepsilon = 0.01$? (exact area = 0.5)
I suppose a graph in Maple is easy?
$\varepsilon = 0.01$, area 0.5001000000
Attachments
area.pdf
area
curve.pdf
eps=0.01 Cuchulainn
Topic Author
Posts: 64440
Joined: July 16th, 2004, 7:38 am
Location: Drosophila melanogaster
Contact:

### Re: Silly questions

Very nice!
"Compatibility means deliberately repeating other people's mistakes."
David Wheeler

http://www.datasimfinancial.com
http://www.datasim.nl lovenatalya
Posts: 287
Joined: December 10th, 2013, 5:54 pm

### Re: Silly questions

I saw some confusions and misconceptions in several posts before Monday. I PM'd Cuchulainn, two days ago, the following clarified proposal along the same line as my solution #2. I did not post it here publicly because I would like to write out the detailed derivation and clarify the misconceptions I saw. But now it seems I do not have a stretch of time to do that quickly, so I am just going to post what I PM'ed Cuchulainn two days ago first.

The formula is simpler than all the previously suggested forms and is derived from the first principle thus possesses the property of the heat equation.

$f_l(t,x) = 2\Big(\frac12+\big(\frac12-x\big)\text{erf(u)} - \sqrt{\frac{2t}\pi}\,e^{-u^2} \Big)$

where $u:=\frac{x-\frac12}{\sqrt{2t}}$. To make sure it hits zero at $\{0,1\}$, you can take $f_l(t,x)-f_l(t,0)$. But I do not think the $f_l(t,0)$ makes any significant difference especially when $t$ is small. This is convolution over the wedge function over the whole real line as I had indicated before instead of the hat function Alan previously was integrating over. The difference is how one extends the hat function to the whole real line. The convergence rate is $\sim\sqrt t$. I will prove this also later.

Here is a plot of my function. Compare this to ppauper's function plot with the equivalent parameter. The former has a single concavity.  The value at $x=\frac12$ is lower than $1$. These are features a heat equation solution with a concave initial condition and zero boundary value should possess. We can not lose the exponential function since this accounts for the deviation of the target function or initial condition from a flat line. The latter curve switches between concavity and convexity, crisscrossing the original hat function and rounds above $(\frac12, 1)$, and below $(0,0)$ and $(1,0)$. These violate the properties of a heat equation solution with a concave initial condition and zero boundary value. Also, the former function needs only to evaluate 1 error function plus an exponential function, whereas the latter 3 error functions.

Of course, in the end, only actual computation can tell which is better. But I think the more faithful the approximation is to the features of the original equation the better. lovenatalya
Posts: 287
Joined: December 10th, 2013, 5:54 pm

### Re: Silly questions

This post is a bit long, because I present a full derivation of my solution mentioned in my last post and clear up the confusion and misconception in several previous posts. In the following $\Theta$ is the Heaviside step function.

The spatial domain is already compact. Saying a function is compact supported in the space dimension is a bit leading. Usually saying a function has compact support means the closure of the point of the domain where the function value is nonzero is a compact strict subset of the whole domain. If it is truly a strict compact subset and that function is analytic within the whole domain, that function can only be zero by the analytic continuation. The partial sum or truncated Fourier series

$f_{N}(x)=\Theta(|x-\frac12|)\big(\sum_{m=1}^{N}\frac{8}{m^2\pi^2}e^{-m^{2}\pi^{2}t}\sin\frac{m\pi}{2}\sin m\pi x\big)$

is only infinitely differentiable in $[0,1]$ and not in $[-\delta, 1+\delta]$ for some positive $\delta$. So saying $f_{N}(x)$ is compact-supported and infinitely differentiable is misleading.

According to Cuchulainn, the problem to be solved is a parabolic (heat) PDE $u$ define on $(t,x)\in[0,\infty)\times [0,1]$ with Dirichlet boundary condition
$u(t,x=0)=u(t,x=1)=0$. If the equation is $\Big(\frac{\partial }{\partial t}-\frac12\frac{\partial^2}{\partial x^2}\Big)u(t,x)=0$ We can solve it two ways. One is the Fourier series method as ppauper used to obtain

$u(t,x)=\sum_{m=1}^{\infty}\frac{8}{m^2\pi^2}e^{-m^{2}\pi^{2}t}\sin\frac{m\pi}{2}\sin m\pi x$

For $t=0$, the Fourier series can be shown to converge uniformly $\sim\frac1m$. The proof is a bit complicated. I may put it up later though. Thus the convergence is very slow for very small $t$ at least for the first partial sums but fast for large $t$. Besides, the partial sum is going to oscillate, which is not desirable to represent a heat equation solution with a concave initial condition and a vanishing boundary.

The other is the image (Green's function) method or the Poisson summation formula Alan and I used to obtain the following solution

$u(t,x)=\frac1{\sqrt{2\pi t}}\int_{-\infty}^\infty dy\, e^{-\frac{(x-y)^2}{2t}} \sum_{m=-\infty}^\infty (u(t=0,y+2m)-u(t=0,-y+2m)),$

where in this setting $u(t,y):=u(t,y)\Theta(|\frac12-y|)$ i.e., the extension of $u(t,\cdot)$ to the whole real axis with zero values assigned outside of $[0,1]$. This is extends the initial value function to an odd function of period $2$. This series converges very fast for small $t$ at the rate of $\sim e^{-\frac{m^2}{2t}}$ but slowly for large $t$. So we use this series.

The fast convergence allows us to just take the first term. We would like to enforce as much as we can without complicating the mathematics too much the vanishing boundary condition, so we take, say the quadratic Taylor expansion at the boundary and make an odd function extension with respect to that boundary point. In our particular case, it is even simpler. We just extend the wedge (not the hat) to infinity. We have

$u_{\text{l}}(t,x)=\frac1{\sqrt{2\pi t}}\int_{-\infty}^\infty dy\, e^{-\frac{(x-y)^2}{2t}}(1-2|y-\frac12|)=1-\frac2{\sqrt{2\pi t}}\int_{-\infty}^\infty dy\, e^{-\frac{y^2}{2t}}|x-\frac12-y|$

$=2\Big(\frac12-\big(\frac12-x\big)\text{erf}(v) - \sqrt{\frac {2t}\pi}\,e^{-v^2} \Big),$

where $v:=\frac{x-\frac12}{\sqrt{2t}}$. We see that the above $u_{\text{l}}(t,x)$ is the positively and analytically weighted $y$integral of concave functions $-|x-\frac12-y|$ which is then concave and analytic in $x$. Note the exponential function reflects the effect of $y$ term or the slope at each point, particularly the broken slope at the vertex, of the initial condition function. The concavity is provided by the exponential function. So it can not be dropped. The error function gives the base line shape and rounds off (making the resulting function analytic) the sharp corner of the vertex.

In fact, the error of the approximation is

$e\big(t,x+\frac12\big) = (1-2|x|)-u_{\text{l}}(t,x+\frac12)=\frac2{\sqrt{2\pi t}}\int_{-\infty}^\infty dy\, e^{-\frac{y^2}{2t}}(|x|-|x-y|).$

We can easily see from the graph of the integrand in the parenthesis that $e\big(t,x+\frac12\big)>0$ and its maximum occurs at $x=0$ i.e. the vertex with value $e\big(t,\frac12\big)=2\sqrt{\frac{2t}\pi}$. I have to emphasize that this expression particularly $\sqrt t$ is to give the functional dependency between $\max_x(e(t,x))$ and the parameterization of $u_{\text l}(t,x)$. One may just pick $\max_x(e(t,x))$ as the parameter and re-parameterize $u_{\text l}(\cdot,x)$. The approximation is always below the target function. This lowering is the effect of the exponential function. and drops the lowest at the vertex. So we know not only the convergence rate but the exact location of the error and its exact value. This is contributed solely by the exponential term. This gives yet another reason the exponential function cannot be dropped.

It is only natural that this approximation should possess all the properties of a heat equation because it IS a solution of a heat equation, while dropping one without the exponential function is NOT.

I have now derived my formula and its desired properties.
Last edited by lovenatalya on April 26th, 2018, 4:42 pm, edited 4 times in total. lovenatalya
Posts: 287
Joined: December 10th, 2013, 5:54 pm

### Re: Silly questions

if you want something with compact support, go with one of the truncated series
It depends on whether he wants the whole approximating function to be infinitely differentiable particularly at $\{0,1\}$.
For the Fourier series, we have
$f_{N}(x)=\sum_{m=1}^{N}\frac{8}{m^2\pi^2}e^{-c^{2}m^{2}\pi^{2}t}\sin\frac{m\pi}{2}\sin m\pi x$
which is infinitely differentiable for finite $N$ (as is the Chebyshev series)
Similarly, $g_{\epsilon}(x)=x\,\mathrm{erfc}(-x/\epsilon)+(1-2x)\,\mathrm{erfc}(-(x-1/2)/\epsilon)+(x-1)\,\mathrm{erfc}(-(x-1)/\epsilon)$ is infinitely differentiable for $\epsilon\ne 0$
Please see my last post for why this is ambiguous to say the least. We can not have analyticity (as these functions are) and of compact support at the same time.
Last edited by lovenatalya on April 27th, 2018, 4:20 am, edited 2 times in total. lovenatalya
Posts: 287
Joined: December 10th, 2013, 5:54 pm

### Re: Silly questions

Just to recap: the solution of the 1d heat equation in an infinite rod with initial condition $f$ whose support is the interval $(a,b)$ is given by the convolution integral

$u(x,t)= \dfrac {1} {2 \sqrt{\pi t}} \int _{a}^{b} f(v) \exp(- (x-v)^2/4t) dv$.
using the function cuch gave originally, this gives

$x\,\mathrm{erf}(x/(2t^{1/2}))+(1-2x)\mathrm{erf}((x-1/2)/(2t^{1/2})) +(x-1)\mathrm{erf}(x-1)/(2t^{1/2}))$
$+\frac{2t^{1/2}}{\pi^{1/2}}\left[e^{-x^{2}/(4t)}-2e^{-(x-1/2)^{2}/(4t)}+e^{-(x-1)^{2}/(4t)}\right]$

the solution I posted had erfc instead of erf and didn't have the exponential terms
Please see my earlier post on why these exponential functions cannot be dropped. Cuchulainn
Topic Author
Posts: 64440
Joined: July 16th, 2004, 7:38 am
Location: Drosophila melanogaster
Contact:

### Re: Silly questions

Thanks again, both.

The background is that most FDM schemes have issues at sharp corners (e.g. max(K-S,0)) and accuracy is affected. There are load of workarounds in discrete space to 'smooth' things but are kind of ad-hoc. The current discussion is to 'smooth' the PDE _before_ using FDM.
Scary thing is I have not seen this technique being used..

// Smoothing a discrete payoff by cubic splines gives overshoots. Is that serious in this case?
"Compatibility means deliberately repeating other people's mistakes."
David Wheeler

http://www.datasimfinancial.com
http://www.datasim.nl ppauper
Posts: 70239
Joined: November 15th, 2001, 1:29 pm

### Re: Silly questions

Thanks again, both.

The background is that most FDM schemes have issues at sharp corners (e.g. max(K-S,0)) and accuracy is affected. There are load of workarounds in discrete space to 'smooth' things but are kind of ad-hoc. The current discussion is to 'smooth' the PDE _before_ using FDM.
Scary thing is I have not seen this technique being used..

// Smoothing a discrete payoff by cubic splines gives overshoots. Is that serious in this case?
what is the PDE? that may affect the "best" approximation

if it's heat conduction, you've already got the solution on $0\le x\le 1$ and on $-\infty<x<\infty$ and use that as your initial condition

if it's Black-Scholes, if you give us the exact payoff it would only take a couple of minutes to work out the European solution and you could use that as your initial condition
I know you said
$f(x) = 2x, 0 \leq x \leq 1/2;$ $f(x) = 2(1-x), 1/2 \leq x \leq 1$.
but does it actually start from $x=0$ lovenatalya
Posts: 287
Joined: December 10th, 2013, 5:54 pm

### Re: Silly questions

There are indeed infinitely many ways to smooth things. The ultimate judge is trying them numerically directly.

The asymptotic error/speed of my previous method is $e^{-\frac{x^2}{2t}}$.

Here is my 3rd way to smooth your function: $1-2 t\ln\Big(2\cosh\big(\frac{x-\frac12}t\big)\Big)$ where $t>0$ is again some positive of your choosing to control the error $2t\ln2$ of the approximation. The asymptotics error/speed is $2t\,e^{-2\frac xt}$, slower than my previous method with respect to $x$, but still exponential. You can replace $t$ with $\sqrt{2t}$ for this expression to be in the same unit as my previous one for comparison. Cuchulainn
Topic Author
Posts: 64440
Joined: July 16th, 2004, 7:38 am
Location: Drosophila melanogaster
Contact:

### Re: Silly questions

Can the method of Separation of Variables be used solve the Black Scholes PDE?

1. If not, why not, where does it break down?
2. If yes, how do you compute a solution?

It's an interesting quizzie. There's a lot of maths in there.
"Compatibility means deliberately repeating other people's mistakes."
David Wheeler

http://www.datasimfinancial.com
http://www.datasim.nl ppauper
Posts: 70239
Joined: November 15th, 2001, 1:29 pm

### Re: Silly questions

you get solutions like $V_{n}(S,t)=\exp(a_{n}t)S^{b_{n}}$ and sum over them $V(S,t)=\sum_{n}C_{n}\exp(a_{n}t)S^{b_{n}}$
the issue is that there's a limited number of problems where that's useful Maybe barrier options or double barrier options.
If you go to the continuous spectrum, you might have a payoff $V(S,T)$ which you can write as an inverse Mellin transform
$V(S,T)=\frac{1}{2\pi i}\int_{\gamma-i\infty}^{\gamma+i\infty}S^{-k}v(k)dk$
then for $t<T$ you'd have an exponential in time in there
$V(S,t)=\frac{1}{2\pi i}\int_{\gamma-i\infty}^{\gamma+i\infty}\exp(a(k)(T-t))S^{-k}v(k)dk$
Last edited by ppauper on November 17th, 2018, 8:20 am, edited 3 times in total. Alan
Posts: 10612
Joined: December 19th, 2001, 4:01 am
Location: California
Contact:

### Re: Silly questions

Can the method of Separation of Variables be used solve the Black Scholes PDE?

1. If not, why not, where does it break down?
2. If yes, how do you compute a solution?

It's an interesting quizzie. There's a lot of maths in there.

2. Yes. Start with the easy driftless BS case, where $x = \log S$ and $\mu = r - \sigma^2/2 = 0$. Then, your question amounts to asking: can we use separation of variables to solve the heat equation for $u(t,x)$ on $x \in (-\infty,\infty)$ with $u(0,x) = f(x)$ given?

Ans: sure, tedious, but start by truncating the problem to $x \in (-L,L)$  with $u=0$ at boundaries. Use standard separation of variables. Following wikipedia, for example, the result has the form:

$u_L(t,x) = \int_{-L}^L G_L(t,x,y) f(y) \, dy$,

where $G_L$ is given by an infinite sum. Then, prove that, as $L \rightarrow \infty$,

$G_L \rightarrow \frac{e^{-(x-y)^2/(2 \sigma^2 t)}}{\sqrt{2 \pi \sigma^2 t}}$.

Surely doable and involves many of the same formulas already found in this thread.

Finally, generalize to non-zero drift, and you're done. ppauper
Posts: 70239
Joined: November 15th, 2001, 1:29 pm

### Re: Silly questions

to add to that, in Paul's "Derivatives" text there's a section "series solution" in the chapter on PDEs where he has an infinite sum
$\sum f_{n}(x)g_{n}(\tau)$
written in the dfffusion equation variables so $x$ is something like $x=\log(S/X)$ Cuchulainn Topic Author Posts: 64440 Joined: July 16th, 2004, 7:38 am Location: Drosophila melanogaster Contact: ### Re: Silly questions you get solutions like$V_{n}(S,t)=\exp(a_{n}t)S^{b_{n}}$and sum over them$V(S,t)=\sum_{n}C_{n}\exp(a_{n}t)S^{b_{n}}$the issue is that there's a limited number of problems where that's useful Maybe barrier options or double barrier options. If you go to the continuous spectrum, you might have a payoff V(S,T) which you can write as an inverse Mellin transform$V(S,T)=\frac{1}{2\pi i}\int_{\gamma-i\infty}^{\gamma+i\infty}S^{-k}v(k)dk$then for$t<T$you'd have an exponential in time in there$V(S,t)=\frac{1}{2\pi i}\int_{\gamma-i\infty}^{\gamma+i\infty}\exp(a(k)(T-t))S^{-k}v(k)dk$Mellin transform seems appropriate since the elliptic part of BSPDE is an Euler equation. Most methods' solution end up as an infinite series or integral when they run their course. At some stage they have to solved numerically (Maple, Mathematica?). The formula (22) in the paper looks doable but equation (27) not so because we have to sole for both the price and the free boundary S*. https://ijpam.eu/contents/2014-97-3/3/3.pdf Even if we had an extra condition in (27) it would be a big numerical challenge to compute P and S*? "Compatibility means deliberately repeating other people's mistakes." David Wheeler http://www.datasimfinancial.com http://www.datasim.nl Cuchulainn Topic Author Posts: 64440 Joined: July 16th, 2004, 7:38 am Location: Drosophila melanogaster Contact: ### Re: Silly questions Can the method of Separation of Variables be used solve the Black Scholes PDE? 1. If not, why not, where does it break down? 2. If yes, how do you compute a solution? It's an interesting quizzie. There's a lot of maths in there. 2. Yes. Start with the easy driftless BS case, where$x = \log S$and$\mu = r - \sigma^2/2 = 0$. Then, your question amounts to asking: can we use separation of variables to solve the heat equation for$u(t,x)$on$x \in (-\infty,\infty)$with$u(0,x) = f(x)$given? Ans: sure, tedious, but start by truncating the problem to$x \in (-L,L)$with$u=0$at boundaries. Use standard separation of variables. Following wikipedia, for example, the result has the form:$u_L(t,x) = \int_{-L}^L G_L(t,x,y) f(y) \, dy$, where$G_L$is given by an infinite sum. Then, prove that, as$L \rightarrow \infty$,$G_L \rightarrow \frac{e^{-(x-y)^2/(2 \sigma^2 t)}}{\sqrt{2 \pi \sigma^2 t}}$. Surely doable and involves many of the same formulas already found in this thread. Finally, generalize to non-zero drift, and you're done. I think this approach will also work in more general cases. Here are my ideas (BTW I have not implemented it): 1. Reduce BSPDE(S,t) to self-adjoint form (removes drift) 2. Define y = S/(1+S) to get BSPDE(y,t) on a unit interval [0,1]. 3. Separation of variables P = Y(y)T(t) to get two uncoupled ODEs (pde's coefficients independent of t). 4. Solve the Sturm-Liouville (SL) problem for the ODE for using FDM on a mesh of size N. Solve for the eigenvalues and eigenvectors. 5. Plug the stuff from 4 into the (truncated) series expansion for P = YT. remark a. For each eigenvalue we have to solve FDM for the corresponding eigenvector (N of them). b. What if the Dirichlet BC for SL are not homogeneous. c. Step a can be parallelised (?) d. We only discretise in space (S or y) but not in t (so, no Crank-Nicolson and such-like needed). e. Your domain truncation [-L, L] is probably better than domain transformation _provided_ we don't lose 'information'/loss of accuracy. Is this approach computationally feasible? // In the case of the heat equation the eigenvalues are$\lambda_n = -n^2\pi^2$and eigenvectors are$X_n(x) = sin(n\pi x)[$] but in general the are all different.
"Compatibility means deliberately repeating other people's mistakes."
David Wheeler

http://www.datasimfinancial.com
http://www.datasim.nl