Chase The Devilhttp://chasethedevil.blogspot.com/Technical blog for Fabienennoreply@blogger.com (Fabien)Sun, 16 Sep 2018 20:55:48 PDTBlogger http://www.blogger.com670125Controlling the SABR wings with Hagan PDE http://feedproxy.google.com/~r/chasethedevil/~3/FagxywMNtOI/controlling-sabr-wings-with-hagan-pde.htmlnoreply@blogger.com (Fabien)Thu, 17 Dec 2015 07:21:39 PSTtag:blogger.com,1999:blog-14415229.post-846622718474615961On the <a href="http://www.wilmott.com/messageview.cfm?catid=4&threadid=78001&FTVAR_MSGDBTABLE=&STARTPAGE=4">Wilmott forum</a>, Pat Hagan has recently suggested to cap the equivalent local volatility in order to control the wings and better match CMS prices. It also helps making the SABR approximation better behaved as the expansion is only valid when<br /><div>$$ 1 + 2\frac{\rho\nu}{\alpha}y(K)+\frac{\nu^2}{\alpha^2}y^2(K) $$</div><div>is close to 1. </div><div><br /><div>In the PDE approach (especially the non transformed one), it is very simple, one just needs to update the equivalent local vol as </div></div><div>$$\alpha K^\beta \min\left(M, \sqrt{1 + 2\frac{\rho\nu}{\alpha}y(K)+\frac{\nu^2}{\alpha^2}y^2(K)}\right)$$</div><div><br /></div><div>While it is straightforward to include in the PDE, it is more difficult to derive a good approximation. The zero-th order behaves as expected, but the first order formula has a unnatural kink, likely because of the non differentiability due to the min function. </div><div><br /></div><div>The following graphs presents the non capped PDE, the capped PDE with M=4*nu (PDEC4) and M=6*nu (PDEC6) as well as the approximation (Andersen Ratcliffe / Gatheral first order) where I have only taken care of the right wing. The SABR parameters are alpha = 0.0630, beta = 0.7, rho = -0.363, nu = 0.421, T = 10, f = 0.0439.</div><div><br /></div><div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-7N928DyGhHY/Vm_ibyWlVcI/AAAAAAAAIP0/YF7Mfcpm4w4/s1600/Screenshot%2Bfrom%2B2015-12-15%2B10-25-16.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="400" src="http://1.bp.blogspot.com/-7N928DyGhHY/Vm_ibyWlVcI/AAAAAAAAIP0/YF7Mfcpm4w4/s400/Screenshot%2Bfrom%2B2015-12-15%2B10-25-16.png" width="385" /></a></div><div><br /></div><div>We can see that the higher the cap is, the closer we are to the standard SABR PDE, and the lower the cap is, the flatter are the wings.<br /><br />The approximation matches well ATM (it is then equivalent to standard SABR PDE) but then has a discontinuous derivative for the K that reaches the threshold M. Far away, it matches very well again.</div><img src="http://feeds.feedburner.com/~r/chasethedevil/~4/FagxywMNtOI" height="1" width="1" alt=""/>0http://chasethedevil.blogspot.com/2015/12/controlling-sabr-wings-with-hagan-pde.htmlControlling the SABR wings with Hagan PDE http://feedproxy.google.com/~r/chasethedevil/~3/oaizABIeCks/wings-in-hagan-pde.htmlnoreply@blogger.com (Fabien)Tue, 15 Dec 2015 03:50:06 PSTtag:blogger.com,1999:blog-14415229.post-106851415660776087On the <a href="http://www.wilmott.com/messageview.cfm?catid=4&threadid=78001&FTVAR_MSGDBTABLE=&STARTPAGE=4">Wilmott forum</a>, Pat Hagan has recently suggested to cap the equivalent local volatility in order to control the wings and better match CMS prices. It also helps making the SABR approximation better behaved as the expansion is only valid when<br /><div>$$ 1 + 2\frac{\rho\nu}{\alpha}y(K)+\frac{\nu^2}{\alpha^2}y^2(K) $$</div><div>is close to 1. </div><div><br /><div>In the PDE approach (especially the non transformed one), it is very simple, one just needs to update the equivalent local vol as </div></div><div>$$\alpha K^\beta \min\left(M, \sqrt{1 + 2\frac{\rho\nu}{\alpha}y(K)+\frac{\nu^2}{\alpha^2}y^2(K)}\right)$$</div><div><br /></div><div>While it is straightforward to include in the PDE, it is more difficult to derive a good approximation. The zero-th order behaves as expected, but the first order formula has a unnatural kink, likely because of the non differentiability due to the min function. </div><div><br /></div><div>The following graphs presents the non capped PDE, the capped PDE with M=4*nu (PDEC4) and M=6*nu (PDEC6) as well as the approximation (Andersen Ratcliffe / Gatheral first order) where I have only taken care of the right wing. The SABR parameters are alpha = 0.0630, beta = 0.7, rho = -0.363, nu = 0.421, T = 10, f = 0.0439.</div><div><br /></div><div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-7N928DyGhHY/Vm_ibyWlVcI/AAAAAAAAIP0/YF7Mfcpm4w4/s1600/Screenshot%2Bfrom%2B2015-12-15%2B10-25-16.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="400" src="http://1.bp.blogspot.com/-7N928DyGhHY/Vm_ibyWlVcI/AAAAAAAAIP0/YF7Mfcpm4w4/s400/Screenshot%2Bfrom%2B2015-12-15%2B10-25-16.png" width="385" /></a></div><div><br /></div><div>We can see that the higher the cap is, the closer we are to the standard SABR PDE, and the lower the cap is, the flatter are the wings.<br /><br />The approximation matches well ATM (it is then equivalent to standard SABR PDE) but then has a discontinuous derivative for the K that reaches the threshold M. Far away, it matches very well again.</div><img src="http://feeds.feedburner.com/~r/chasethedevil/~4/oaizABIeCks" height="1" width="1" alt=""/>0http://chasethedevil.blogspot.com/2015/12/wings-in-hagan-pde.htmlBroken Internet?http://feedproxy.google.com/~r/chasethedevil/~3/jvJ8lZZjEIo/broken-internet_9.htmlnoreply@blogger.com (Fabien)Thu, 17 Dec 2015 07:21:39 PSTtag:blogger.com,1999:blog-14415229.post-1814672475461651820There is something funny going on with upcoming generic top level domains (gTLDs), they seem to be looked up in a strange manner (at least on latest Linux). For example:<br /><br /><span style="font-family: "Courier New",Courier,monospace;">ping chrome </span><br /><br />or<br /><br /><span style="font-family: "Courier New",Courier,monospace;">ping nexus </span><br /><br />returns 127.0.53.53.<br /><br />While existing <a href="https://www.name.com/new-gtld">official gTLD</a>s don't (<span style="font-family: "Courier New",Courier,monospace;">ping dental</span> returns "unknown host" as expected). I first thought it was a network misconfiguration, but as <a href="https://groups.google.com/forum/#!msg/public-dns-discuss/bzhTQnFqE6I/E9F46xhka98J">I am not the only one to notice this</a>, it's likely a genuine internet issue.<br /><br />Strange times.<img src="http://feeds.feedburner.com/~r/chasethedevil/~4/jvJ8lZZjEIo" height="1" width="1" alt=""/>0http://chasethedevil.blogspot.com/2015/11/broken-internet_9.htmlBroken Internet?http://feedproxy.google.com/~r/chasethedevil/~3/y_VIlpxNC3Q/broken-internet.htmlnoreply@blogger.com (Fabien)Mon, 09 Nov 2015 04:40:09 PSTtag:blogger.com,1999:blog-14415229.post-7221597406458069951There is something funny going on with upcoming generic top level domains (gTLDs), they seem to be looked up in a strange manner (at least on latest Linux). For example:<br /><br /><span style="font-family: "Courier New",Courier,monospace;">ping chrome </span><br /><br />or<br /><br /><span style="font-family: "Courier New",Courier,monospace;">ping nexus </span><br /><br />returns 127.0.53.53.<br /><br />While existing <a href="https://www.name.com/new-gtld">official gTLD</a>s don't (<span style="font-family: "Courier New",Courier,monospace;">ping dental</span> returns "unknown host" as expected). I first thought it was a network misconfiguration, but as <a href="https://groups.google.com/forum/#!msg/public-dns-discuss/bzhTQnFqE6I/E9F46xhka98J">I am not the only one to notice this</a>, it's likely a genuine internet issue.<br /><br />Strange times.<img src="http://feeds.feedburner.com/~r/chasethedevil/~4/y_VIlpxNC3Q" height="1" width="1" alt=""/>0http://chasethedevil.blogspot.com/2015/11/broken-internet.htmlHoliday's read - DFW - Everything and morehttp://feedproxy.google.com/~r/chasethedevil/~3/qnp6fEKnslg/holiday-read-dfw-everything-and-more.htmlnoreply@blogger.com (Fabien)Thu, 17 Dec 2015 07:21:39 PSTtag:blogger.com,1999:blog-14415229.post-1805129951152872170I am ambivalent towards David Foster Wallace. He can write the most creative sentences and make innocuous subjects very interesting. At the same time, i never finished his book <a href="https://en.wikipedia.org/wiki/Infinite_Jest">Infinite Jest</a>, partly because the characters names are too awkward for me so that i never exactly remember who is who, but also because the story itself is a bit too crazy.<br />I knew however that <a href="http://www.amazon.com/Everything-More-Compact-History-Infinity/dp/0393339289">a non fiction book on the subject of infinity</a> written by him would make for a very interesting read. And I have not been disappointed. It's in between maths and philosophy going back to the Greeks up to Gödel through a lot of Cantor following more or less the historical chronology.<br />Most of it is easy to read and follow, except the last part around sets and transfinite numbers. This last part is actually quite significant as it tries to explain why we still have no satisfying theory around the problems raised by infinity especially in the context of a Sets theory. I did not expect to learn much around the subject, I was disappointed. The book showed me how naive I was and how tricky the concept of infinity can be.<br />While I found the different explanations around <a href="https://en.wikipedia.org/wiki/Zeno%27s_paradoxes">Zeno's paradox of the arrow</a> very clever, there is one other view possible: the arrow really does not move at each instant (you could think of those as a snapshot) but an interval of time is just not a simple succession of instants. This is not so far of Aristotle attack, but the key here is around what is an interval really. DFW suggests slightly this interpretation as well p144 but it's not very explicit.<br />I had not heard about Kronecker's conception that only integers were mathematically real (against decimals, irrationals, infinite sets). I find it very appropriate in the frame of computer science. Everything ends up as finite integers (a binary representation) and we are always confronted to the process of transforming the continuous, that despite all its conceptual issues is often simpler to reason in to solve concrete problems, to the finite discrete.<img src="http://feeds.feedburner.com/~r/chasethedevil/~4/qnp6fEKnslg" height="1" width="1" alt=""/>0http://chasethedevil.blogspot.com/2015/11/holiday-read-dfw-everything-and-more.htmlHoliday's read - DFW - Everything and morehttp://feedproxy.google.com/~r/chasethedevil/~3/yhcQEPEBdlE/holidays-read-dfw-everything-and-more.htmlnoreply@blogger.com (Fabien)Sun, 01 Nov 2015 08:56:57 PSTtag:blogger.com,1999:blog-14415229.post-7474886190581667950I am ambivalent towards David Foster Wallace. He can write the most creative sentences and make innocuous subjects very interesting. At the same time, i never finished his book <a href="https://en.wikipedia.org/wiki/Infinite_Jest">Infinite Jest</a>, partly because the characters names are too awkward for me so that i never exactly remember who is who, but also because the story itself is a bit too crazy.<br />I knew however that <a href="http://www.amazon.com/Everything-More-Compact-History-Infinity/dp/0393339289">a non fiction book on the subject of infinity</a> written by him would make for a very interesting read. And I have not been disappointed. It's in between maths and philosophy going back to the Greeks up to Gödel through a lot of Cantor following more or less the historical chronology.<br />Most of it is easy to read and follow, except the last part around sets and transfinite numbers. This last part is actually quite significant as it tries to explain why we still have no satisfying theory around the problems raised by infinity especially in the context of a Sets theory. I did not expect to learn much around the subject, I was disappointed. The book showed me how naive I was and how tricky the concept of infinity can be.<br />While I found the different explanations around <a href="https://en.wikipedia.org/wiki/Zeno%27s_paradoxes">Zeno's paradox of the arrow</a> very clever, there is one other view possible: the arrow really does not move at each instant (you could think of those as a snapshot) but an interval of time is just not a simple succession of instants. This is not so far of Aristotle attack, but the key here is around what is an interval really. DFW suggests slightly this interpretation as well p144 but it's not very explicit.<br />I had not heard about Kronecker's conception that only integers were mathematically real (against decimals, irrationals, infinite sets). I find it very appropriate in the frame of computer science. Everything ends up as finite integers (a binary representation) and we are always confronted to the process of transforming the continuous, that despite all its conceptual issues is often simpler to reason in to solve concrete problems, to the finite discrete.<img src="http://feeds.feedburner.com/~r/chasethedevil/~4/yhcQEPEBdlE" height="1" width="1" alt=""/>0http://chasethedevil.blogspot.com/2015/11/holidays-read-dfw-everything-and-more.htmlCrank-Nicolson and Rannacher Issues with Touch optionshttp://feedproxy.google.com/~r/chasethedevil/~3/5nRgkVU4cWI/crank-nicolson-and-rannacher-issues_30.htmlnoreply@blogger.com (Fabien)Thu, 17 Dec 2015 07:21:39 PSTtag:blogger.com,1999:blog-14415229.post-8719439168906982425I just stumbled upon this particularly illustrative case where the Crank-Nicolson finite difference scheme behaves badly, and the Rannacher smoothing (2-steps backward Euler) is less than ideal: <a href="http://www.investopedia.com/terms/d/doubleonetouch.asp">double one touch</a> and <a href="http://www.investopedia.com/terms/d/doublenotouch.asp">double no touch</a> options.<br /><br />It is particularly evident when the option is sure to be hit, for example when the barriers are narrow, that is our delta should be around zero as well as our gamma. Let's consider a double one touch option with spot=100, upBarrier=101, downBarrier=99.9, vol=20%, T=1 month and a payout of 50K.<br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://3.bp.blogspot.com/-skVEtpSetds/VgvDcK5MycI/AAAAAAAAIIc/BPj70_3z4lo/s1600/Screenshot%2Bfrom%2B2015-09-30%2B13%253A11%253A13.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="326" src="http://3.bp.blogspot.com/-skVEtpSetds/VgvDcK5MycI/AAAAAAAAIIc/BPj70_3z4lo/s400/Screenshot%2Bfrom%2B2015-09-30%2B13%253A11%253A13.png" width="400" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Crank-Nicolson shows big spikes in the delta near the boundary</td><td class="tr-caption" style="text-align: center;"><br /></td></tr></tbody></table><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://3.bp.blogspot.com/-w0esoFRdaSA/VgvD6QgMDwI/AAAAAAAAIIk/-qdQ6BcTAmU/s1600/Screenshot%2Bfrom%2B2015-09-30%2B13%253A13%253A33.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="325" src="http://3.bp.blogspot.com/-w0esoFRdaSA/VgvD6QgMDwI/AAAAAAAAIIk/-qdQ6BcTAmU/s400/Screenshot%2Bfrom%2B2015-09-30%2B13%253A13%253A33.png" width="400" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Rannacher shows spikes in the delta as well</td></tr></tbody></table>Crank-Nicolson spikes are so high that the price is actually a off itself.<br /><br />The Rannacher smoothing reduces the spikes by 100x but it's still quite high, and would be higher had we placed the spot closer to the boundary. The gamma is worse. Note that we applied the smoothing only at maturity. In reality as the barrier is continuous, the smoothing should really be applied at each step, but then the scheme would be not so different from a simple Backward Euler.<br /><br />In contrast, with a proper second order finite difference scheme, there is no spike.<br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://2.bp.blogspot.com/-mj0mUfRCSJk/VgvGiUPP1nI/AAAAAAAAIIw/KKK9sXTrne4/s1600/Screenshot%2Bfrom%2B2015-09-30%2B13%253A24%253A27.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="325" src="http://2.bp.blogspot.com/-mj0mUfRCSJk/VgvGiUPP1nI/AAAAAAAAIIw/KKK9sXTrne4/s400/Screenshot%2Bfrom%2B2015-09-30%2B13%253A24%253A27.png" width="400" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Delta with the TR-BDF2 finite difference method - the scale goes from -0.00008 to 0.00008.</td></tr></tbody></table><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://1.bp.blogspot.com/-okMVRlfdJGw/VgvGsntsjbI/AAAAAAAAII4/xNVchYODHGU/s1600/Screenshot%2Bfrom%2B2015-09-30%2B13%253A24%253A42.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="325" src="http://1.bp.blogspot.com/-okMVRlfdJGw/VgvGsntsjbI/AAAAAAAAII4/xNVchYODHGU/s400/Screenshot%2Bfrom%2B2015-09-30%2B13%253A24%253A42.png" width="400" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Delta with the Lawson-Morris finite difference scheme - the scale goes from -0.00005 to 0.00005</td></tr></tbody></table>Both <a href="http://www.risk.net/journal-of-computational-finance/technical-paper/2330321/tr-bdf2-for-fast-stable-american-option-pricing">TR-BDF2</a> and Lawson-Morris (based on a local Richardson extrapolation of backward Euler) have a very low delta error, similarly, their gamma is very clean. This is reminiscent of the behavior on American options, but the effect is magnified here.<br /><br /><br /><img src="http://feeds.feedburner.com/~r/chasethedevil/~4/5nRgkVU4cWI" height="1" width="1" alt=""/>1http://chasethedevil.blogspot.com/2015/09/crank-nicolson-and-rannacher-issues_30.htmlCrank-Nicolson and Rannacher Issues with Touch optionshttp://feedproxy.google.com/~r/chasethedevil/~3/aSx7fmQgTS4/crank-nicolson-and-rannacher-issues.htmlnoreply@blogger.com (Fabien)Wed, 30 Sep 2015 04:34:09 PDTtag:blogger.com,1999:blog-14415229.post-601514439529362658I just stumbled upon this particularly illustrative case where the Crank-Nicolson finite difference scheme behaves badly, and the Rannacher smoothing (2-steps backward Euler) is less than ideal: <a href="http://www.investopedia.com/terms/d/doubleonetouch.asp">double one touch</a> and <a href="http://www.investopedia.com/terms/d/doublenotouch.asp">double no touch</a> options.<br /><br />It is particularly evident when the option is sure to be hit, for example when the barriers are narrow, that is our delta should be around zero as well as our gamma. Let's consider a double one touch option with spot=100, upBarrier=101, downBarrier=99.9, vol=20%, T=1 month and a payout of 50K.<br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://3.bp.blogspot.com/-skVEtpSetds/VgvDcK5MycI/AAAAAAAAIIc/BPj70_3z4lo/s1600/Screenshot%2Bfrom%2B2015-09-30%2B13%253A11%253A13.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="326" src="http://3.bp.blogspot.com/-skVEtpSetds/VgvDcK5MycI/AAAAAAAAIIc/BPj70_3z4lo/s400/Screenshot%2Bfrom%2B2015-09-30%2B13%253A11%253A13.png" width="400" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Crank-Nicolson shows big spikes in the delta near the boundary</td><td class="tr-caption" style="text-align: center;"><br /></td></tr></tbody></table><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://3.bp.blogspot.com/-w0esoFRdaSA/VgvD6QgMDwI/AAAAAAAAIIk/-qdQ6BcTAmU/s1600/Screenshot%2Bfrom%2B2015-09-30%2B13%253A13%253A33.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="325" src="http://3.bp.blogspot.com/-w0esoFRdaSA/VgvD6QgMDwI/AAAAAAAAIIk/-qdQ6BcTAmU/s400/Screenshot%2Bfrom%2B2015-09-30%2B13%253A13%253A33.png" width="400" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Rannacher shows spikes in the delta as well</td></tr></tbody></table>Crank-Nicolson spikes are so high that the price is actually a off itself.<br /><br />The Rannacher smoothing reduces the spikes by 100x but it's still quite high, and would be higher had we placed the spot closer to the boundary. The gamma is worse. Note that we applied the smoothing only at maturity. In reality as the barrier is continuous, the smoothing should really be applied at each step, but then the scheme would be not so different from a simple Backward Euler.<br /><br />In contrast, with a proper second order finite difference scheme, there is no spike.<br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://2.bp.blogspot.com/-mj0mUfRCSJk/VgvGiUPP1nI/AAAAAAAAIIw/KKK9sXTrne4/s1600/Screenshot%2Bfrom%2B2015-09-30%2B13%253A24%253A27.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="325" src="http://2.bp.blogspot.com/-mj0mUfRCSJk/VgvGiUPP1nI/AAAAAAAAIIw/KKK9sXTrne4/s400/Screenshot%2Bfrom%2B2015-09-30%2B13%253A24%253A27.png" width="400" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Delta with the TR-BDF2 finite difference method - the scale goes from -0.00008 to 0.00008.</td></tr></tbody></table><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://1.bp.blogspot.com/-okMVRlfdJGw/VgvGsntsjbI/AAAAAAAAII4/xNVchYODHGU/s1600/Screenshot%2Bfrom%2B2015-09-30%2B13%253A24%253A42.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="325" src="http://1.bp.blogspot.com/-okMVRlfdJGw/VgvGsntsjbI/AAAAAAAAII4/xNVchYODHGU/s400/Screenshot%2Bfrom%2B2015-09-30%2B13%253A24%253A42.png" width="400" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Delta with the Lawson-Morris finite difference scheme - the scale goes from -0.00005 to 0.00005</td></tr></tbody></table>Both <a href="http://www.risk.net/journal-of-computational-finance/technical-paper/2330321/tr-bdf2-for-fast-stable-american-option-pricing">TR-BDF2</a> and Lawson-Morris (based on a local Richardson extrapolation of backward Euler) have a very low delta error, similarly, their gamma is very clean. This is reminiscent of the behavior on American options, but the effect is magnified here.<br /><br /><br /><img src="http://feeds.feedburner.com/~r/chasethedevil/~4/aSx7fmQgTS4" height="1" width="1" alt=""/>1http://chasethedevil.blogspot.com/2015/09/crank-nicolson-and-rannacher-issues.htmlCloudshttp://feedproxy.google.com/~r/chasethedevil/~3/H-sWj0gn9Jg/clouds_2.htmlnoreply@blogger.com (Fabien)Thu, 17 Dec 2015 07:21:39 PSTtag:blogger.com,1999:blog-14415229.post-5799469505047569690I was wondering how to generate some nice cloudy like texture with a simple program. I first thought about using the Brownian motion, but of course if one uses it raw, with one pixel representing one movement in time, it's just going to look like a very noisy and grainy picture like this:<br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://4.bp.blogspot.com/-QfeLi5F03oQ/Vebs3C-XWyI/AAAAAAAAIG4/5q1kKip0PlA/s1600/normal_rng.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="200" src="http://4.bp.blogspot.com/-QfeLi5F03oQ/Vebs3C-XWyI/AAAAAAAAIG4/5q1kKip0PlA/s200/normal_rng.png" width="200" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Normal noise</td></tr></tbody></table><br />There is however a nice continuous representation of the Brownian motion : the Paley-Wiener representation<br /><div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-fkDLapfROaY/Vebtgtm2BCI/AAAAAAAAIHA/-EgKgEF_rEM/s1600/Screenshot%2Bfrom%2B2015-09-02%2B14-37-02.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="31" src="http://1.bp.blogspot.com/-fkDLapfROaY/Vebtgtm2BCI/AAAAAAAAIHA/-EgKgEF_rEM/s400/Screenshot%2Bfrom%2B2015-09-02%2B14-37-02.png" width="400" /></a></div><div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-JdA4jtT3J-A/VebtgjUKcfI/AAAAAAAAIHE/4HCx4BaAPDY/s1600/Screenshot%2Bfrom%2B2015-09-02%2B14-36-42.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="60" src="http://1.bp.blogspot.com/-JdA4jtT3J-A/VebtgjUKcfI/AAAAAAAAIHE/4HCx4BaAPDY/s400/Screenshot%2Bfrom%2B2015-09-02%2B14-36-42.png" width="400" /></a></div><br />This can produce an interesting smooth pattern, but it is just 1D. In the following picture, I apply it to each row (the column index being time), and then for each column (the row index being time). Of course this produces a symmetric picture, especially as I reused the same random numbers<br /><div class="separator" style="clear: both; text-align: center;"><a href="http://4.bp.blogspot.com/-CzAUnTU3V7w/VebuSjnBuLI/AAAAAAAAIHQ/t6a1FGCsIMA/s1600/constructive.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="200" src="http://4.bp.blogspot.com/-CzAUnTU3V7w/VebuSjnBuLI/AAAAAAAAIHQ/t6a1FGCsIMA/s200/constructive.png" width="200" /></a></div>If I use new random numbers for the columns, it is still symmetric, but destructive rather than constructive.<br /><div class="separator" style="clear: both; text-align: center;"><a href="http://3.bp.blogspot.com/-zb3pUtY1ZMo/VebvK0C1xdI/AAAAAAAAIHY/A0-vaxGMwdI/s1600/destructive.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="200" src="http://3.bp.blogspot.com/-zb3pUtY1ZMo/VebvK0C1xdI/AAAAAAAAIHY/A0-vaxGMwdI/s200/destructive.png" width="200" /></a></div><br />It turns out that spatial processes are something more complex than I first imagined. It is not a simple as using a N-dimensional Brownian motion, as it would produce a very similar picture as the 1-dimensional one. But <a href="http://www.maths.uq.edu.au/~kroese/ps/MCSpatial.pdf">this paper has a nice overview of spatial processes</a>. Interestingly they even suggest to generate a Gaussian process using a <a href="https://en.wikipedia.org/wiki/Precision_%28statistics%29">Precision matrix</a> (inverse of covariance matrix). I never thought about doing such a thing and I am not sure what is the advantage of such a scheme.<br /><br />There is a standard graphic technique to generate nice textures, originating from Ken Perlin for Disney, it is called simply <a href="https://en.wikipedia.org/wiki/Perlin_noise">Perlin Noise</a>. It turns out that several web pages in the top Google results <a href="https://en.wikipedia.org/wiki/Talk%3APerlin_noise">confuse</a> simple Perlin noise with fractal sum of noise that Ken Perlin also helped popularize (see his slides: <a href="http://www.noisemachine.com/talk1/20.html">standard Perlin noise</a>, <a href="http://www.noisemachine.com/talk1/21.html">fractal noise</a>). Those pages also believe that the later is simpler/faster. But there are two issues with fractal sum of noise: the first one is that it relies on an existing noise function - you need to first build one (it can be done with a random number generator and an interpolator), and the second one is that it ends up being more complex to program and likely to evaluate as well, see for example the code needed <a href="http://devmag.org.za/2009/04/25/perlin-noise/">here</a>. The fractal sum of noise is really a complementary technique.<br /><br />The insight of Perlin noise is to not generate random color values that would be assigned to shades of grey as in my examples, but to generate random gradients, and interpolate on those gradient in a smooth manner. In computer graphics they like the cosine function to give a little bit of non-linearity in the colors. A good approximation, usually used as a replacement in this context is <a href="http://codeplea.com/simple-interpolation">3x^2 - 2x^3</a>. It's not much more complicated than that, <a href="http://webstaff.itn.liu.se/~stegu/TNM022-2005/perlinnoiselinks/perlin-noise-math-faq.html">this web page</a> explains it in great details. It can be programmed in a few lines of code.<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="http://3.bp.blogspot.com/-Yo5dgKoD5P8/Veb6rAQ9jAI/AAAAAAAAIHo/uISQjpVzLfI/s1600/perlin_bw.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="200" src="http://3.bp.blogspot.com/-Yo5dgKoD5P8/Veb6rAQ9jAI/AAAAAAAAIHo/uISQjpVzLfI/s200/perlin_bw.png" width="200" /></a><a href="http://4.bp.blogspot.com/-MovHUFAGbnI/Veb6rsjhkBI/AAAAAAAAIHs/4raZRohDySs/s1600/perlin_color.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="200" src="http://4.bp.blogspot.com/-MovHUFAGbnI/Veb6rsjhkBI/AAAAAAAAIHs/4raZRohDySs/s200/perlin_color.png" width="200" /></a></div><br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://3.bp.blogspot.com/-__3z1MEB92o/Veb65IEgVZI/AAAAAAAAIH4/1B9npSW9XKo/s1600/Screenshot%2Bfrom%2B2015-09-02%2B15-33-11.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="512" src="http://3.bp.blogspot.com/-__3z1MEB92o/Veb65IEgVZI/AAAAAAAAIH4/1B9npSW9XKo/s640/Screenshot%2Bfrom%2B2015-09-02%2B15-33-11.png" width="640" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">very procedural and non-optimized Go code for Perlin noise</td></tr></tbody></table><br /><img src="http://feeds.feedburner.com/~r/chasethedevil/~4/H-sWj0gn9Jg" height="1" width="1" alt=""/>0http://chasethedevil.blogspot.com/2015/09/clouds_2.htmlCloudshttp://feedproxy.google.com/~r/chasethedevil/~3/9K5IstGMnWQ/clouds.htmlnoreply@blogger.com (Fabien)Wed, 02 Sep 2015 06:36:15 PDTtag:blogger.com,1999:blog-14415229.post-2668356434970943788I was wondering how to generate some nice cloudy like texture with a simple program. I first thought about using the Brownian motion, but of course if one uses it raw, with one pixel representing one movement in time, it's just going to look like a very noisy and grainy picture like this:<br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://4.bp.blogspot.com/-QfeLi5F03oQ/Vebs3C-XWyI/AAAAAAAAIG4/5q1kKip0PlA/s1600/normal_rng.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="200" src="http://4.bp.blogspot.com/-QfeLi5F03oQ/Vebs3C-XWyI/AAAAAAAAIG4/5q1kKip0PlA/s200/normal_rng.png" width="200" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Normal noise</td></tr></tbody></table><br />There is however a nice continuous representation of the Brownian motion : the Paley-Wiener representation<br /><div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-fkDLapfROaY/Vebtgtm2BCI/AAAAAAAAIHA/-EgKgEF_rEM/s1600/Screenshot%2Bfrom%2B2015-09-02%2B14-37-02.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="31" src="http://1.bp.blogspot.com/-fkDLapfROaY/Vebtgtm2BCI/AAAAAAAAIHA/-EgKgEF_rEM/s400/Screenshot%2Bfrom%2B2015-09-02%2B14-37-02.png" width="400" /></a></div><div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-JdA4jtT3J-A/VebtgjUKcfI/AAAAAAAAIHE/4HCx4BaAPDY/s1600/Screenshot%2Bfrom%2B2015-09-02%2B14-36-42.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="60" src="http://1.bp.blogspot.com/-JdA4jtT3J-A/VebtgjUKcfI/AAAAAAAAIHE/4HCx4BaAPDY/s400/Screenshot%2Bfrom%2B2015-09-02%2B14-36-42.png" width="400" /></a></div><br />This can produce an interesting smooth pattern, but it is just 1D. In the following picture, I apply it to each row (the column index being time), and then for each column (the row index being time). Of course this produces a symmetric picture, especially as I reused the same random numbers<br /><div class="separator" style="clear: both; text-align: center;"><a href="http://4.bp.blogspot.com/-CzAUnTU3V7w/VebuSjnBuLI/AAAAAAAAIHQ/t6a1FGCsIMA/s1600/constructive.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="200" src="http://4.bp.blogspot.com/-CzAUnTU3V7w/VebuSjnBuLI/AAAAAAAAIHQ/t6a1FGCsIMA/s200/constructive.png" width="200" /></a></div>If I use new random numbers for the columns, it is still symmetric, but destructive rather than constructive.<br /><div class="separator" style="clear: both; text-align: center;"><a href="http://3.bp.blogspot.com/-zb3pUtY1ZMo/VebvK0C1xdI/AAAAAAAAIHY/A0-vaxGMwdI/s1600/destructive.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="200" src="http://3.bp.blogspot.com/-zb3pUtY1ZMo/VebvK0C1xdI/AAAAAAAAIHY/A0-vaxGMwdI/s200/destructive.png" width="200" /></a></div><br />It turns out that spatial processes are something more complex than I first imagined. It is not a simple as using a N-dimensional Brownian motion, as it would produce a very similar picture as the 1-dimensional one. But <a href="http://www.maths.uq.edu.au/~kroese/ps/MCSpatial.pdf">this paper has a nice overview of spatial processes</a>. Interestingly they even suggest to generate a Gaussian process using a <a href="https://en.wikipedia.org/wiki/Precision_%28statistics%29">Precision matrix</a> (inverse of covariance matrix). I never thought about doing such a thing and I am not sure what is the advantage of such a scheme.<br /><br />There is a standard graphic technique to generate nice textures, originating from Ken Perlin for Disney, it is called simply <a href="https://en.wikipedia.org/wiki/Perlin_noise">Perlin Noise</a>. It turns out that several web pages in the top Google results <a href="https://en.wikipedia.org/wiki/Talk%3APerlin_noise">confuse</a> simple Perlin noise with fractal sum of noise that Ken Perlin also helped popularize (see his slides: <a href="http://www.noisemachine.com/talk1/20.html">standard Perlin noise</a>, <a href="http://www.noisemachine.com/talk1/21.html">fractal noise</a>). Those pages also believe that the later is simpler/faster. But there are two issues with fractal sum of noise: the first one is that it relies on an existing noise function - you need to first build one (it can be done with a random number generator and an interpolator), and the second one is that it ends up being more complex to program and likely to evaluate as well, see for example the code needed <a href="http://devmag.org.za/2009/04/25/perlin-noise/">here</a>. The fractal sum of noise is really a complementary technique.<br /><br />The insight of Perlin noise is to not generate random color values that would be assigned to shades of grey as in my examples, but to generate random gradients, and interpolate on those gradient in a smooth manner. In computer graphics they like the cosine function to give a little bit of non-linearity in the colors. A good approximation, usually used as a replacement in this context is <a href="http://codeplea.com/simple-interpolation">3x^2 - 2x^3</a>. It's not much more complicated than that, <a href="http://webstaff.itn.liu.se/~stegu/TNM022-2005/perlinnoiselinks/perlin-noise-math-faq.html">this web page</a> explains it in great details. It can be programmed in a few lines of code.<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="http://3.bp.blogspot.com/-Yo5dgKoD5P8/Veb6rAQ9jAI/AAAAAAAAIHo/uISQjpVzLfI/s1600/perlin_bw.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="200" src="http://3.bp.blogspot.com/-Yo5dgKoD5P8/Veb6rAQ9jAI/AAAAAAAAIHo/uISQjpVzLfI/s200/perlin_bw.png" width="200" /></a><a href="http://4.bp.blogspot.com/-MovHUFAGbnI/Veb6rsjhkBI/AAAAAAAAIHs/4raZRohDySs/s1600/perlin_color.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="200" src="http://4.bp.blogspot.com/-MovHUFAGbnI/Veb6rsjhkBI/AAAAAAAAIHs/4raZRohDySs/s200/perlin_color.png" width="200" /></a></div><br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://3.bp.blogspot.com/-__3z1MEB92o/Veb65IEgVZI/AAAAAAAAIH4/1B9npSW9XKo/s1600/Screenshot%2Bfrom%2B2015-09-02%2B15-33-11.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="512" src="http://3.bp.blogspot.com/-__3z1MEB92o/Veb65IEgVZI/AAAAAAAAIH4/1B9npSW9XKo/s640/Screenshot%2Bfrom%2B2015-09-02%2B15-33-11.png" width="640" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">very procedural and non-optimized Go code for Perlin noise</td></tr></tbody></table><br /><img src="http://feeds.feedburner.com/~r/chasethedevil/~4/9K5IstGMnWQ" height="1" width="1" alt=""/>0http://chasethedevil.blogspot.com/2015/09/clouds.htmlGo for Monte-Carlohttp://feedproxy.google.com/~r/chasethedevil/~3/vEFgFCwDEJ4/go-for-monte-carlo_22.htmlquantnoreply@blogger.com (Fabien)Thu, 17 Dec 2015 07:21:39 PSTtag:blogger.com,1999:blog-14415229.post-3522511126617972760I have <a href="http://chasethedevil.blogspot.fr/2015/04/modern-programming-language-for-monte.html">looked</a> a few months ago already at Julia, Dart, Rust and Scala programming languages to see how practical they could be for a simple Monte-Carlo option pricing.<br /><br />I forgot <a href="https://golang.org/">the Go language</a>. I had tried it 1 or 2 years ago, and at that time, did not enjoy it too much. Looking at Go 1.5 benchmarks on the<a href="http://benchmarksgame.alioth.debian.org/"> computer language shootout</a>, I was surprised that it seemed so close to Java performance now, while having a GC that guarantees pauses of less 10ms and consuming much less memory.<br /><br />I am in general a bit skeptical about those benchmarks, some can be rigged. A few years ago, I <a href="http://chasethedevil.blogspot.fr/2009/01/end-of-rings-around-plain-java-better.html">tried my hand at the thread ring</a> test, and found that it actually performed fastest on a single thread while it is supposed to measure the language threading performance. I looked yesterday at one Go source code (I think it was for pidigits) and saw that it just called a C library (gmp) to compute with big integers. It's no surprise then that Go would be faster than Java on this test.<br /><br />So what about my small Monte-Carlo test?<br />Well it turns out that Go is quite fast on it:<br /><span style="font-family: "Courier New",Courier,monospace;">Multipl. Rust Go<br />1 0.005 0.007<br />10 0.03 0.03<br />100 0.21 0.29</span><br /><span style="font-family: "Courier New",Courier,monospace;">1000 2.01 2.88</span><br /><br /><br />It is faster than Java/Scala and not too far off Rust, except if one uses FastMath in Scala, then the longest test is slighly faster with Java (not the other ones).<br /><br />There are some issues with the Go language: there is no operator overloading, which can make matrix/vector algebra more tedious and there is no generic/template. The later is somewhat mitigated by the automatic interface implementation. And fortunately for the former, complex numbers are a standard type. Still, automatic differentiation would be painful.<br /><br />Still it was extremely quick to grasp and write code, because it's so simple, especially when compared to Rust. But then, contrary to Rust, there is not as much safety provided by the language. Rust is quite impressive on this side (but unfortunately that implies less readable code). I'd say that Go could become a serious alternative to Java.<br /><br />I also found an interesting minor performance issue with the default Go <a href="https://golang.org/src/math/rand/rand.go">Rand.Float64</a>, the library convert an Int63 to a double precision number this way:<br /><br /><pre>func (r *Rand) Float64() float64 {</pre><pre> f := float64(r.Int63()) / (1 << 63)<br /> if f == 1 {<br /><span class="ln" id="L128"> </span>f = 0<br /><span class="ln" id="L129"> }</span><br /><span class="ln" id="L130"> </span>return f<br /> }</pre><br />I was interested in having a number in (0,1) and not [0,1), so I just used the conversion pattern from MersenneTwister 64 code:<br /><br /><pre>f := (float64(r.Int63() >> 11) + 0.5) * (1.0/4503599627370496.0)</pre><pre> </pre>The reasoning behind this later code is that the mantissa is 52 bits, and this is the most accuracy we can have between 0 and 1. There is no need to go further, this also avoids the issue around 1. It's also straightforward that is will preserve the uniform property, while it's not so clear to me that r.Int63()/2^63 is going to preserve uniformity as double accuracy is higher around 0 (as the exponent part can be used there) and lesser around 1: there is going to be much more multiple identical results near 1 than near 0.<br /><br />It turns out that the if check adds 5% performance penalty on this test, likely because of processor caching issues. I was surprised by that since there are many other ifs afterwards in the code, for the inverse cumulative function, and for the payoff.<img src="http://feeds.feedburner.com/~r/chasethedevil/~4/vEFgFCwDEJ4" height="1" width="1" alt=""/>0http://chasethedevil.blogspot.com/2015/08/go-for-monte-carlo_22.htmlGo for Monte-Carlohttp://feedproxy.google.com/~r/chasethedevil/~3/Yy95xIIOUks/go-for-monte-carlo.htmlquantnoreply@blogger.com (Fabien)Fri, 04 Sep 2015 12:45:48 PDTtag:blogger.com,1999:blog-14415229.post-7850200118684941881I have <a href="http://chasethedevil.blogspot.fr/2015/04/modern-programming-language-for-monte.html">looked</a> a few months ago already at Julia, Dart, Rust and Scala programming languages to see how practical they could be for a simple Monte-Carlo option pricing.<br /><br />I forgot <a href="https://golang.org/">the Go language</a>. I had tried it 1 or 2 years ago, and at that time, did not enjoy it too much. Looking at Go 1.5 benchmarks on the<a href="http://benchmarksgame.alioth.debian.org/"> computer language shootout</a>, I was surprised that it seemed so close to Java performance now, while having a GC that guarantees pauses of less 10ms and consuming much less memory.<br /><br />I am in general a bit skeptical about those benchmarks, some can be rigged. A few years ago, I <a href="http://chasethedevil.blogspot.fr/2009/01/end-of-rings-around-plain-java-better.html">tried my hand at the thread ring</a> test, and found that it actually performed fastest on a single thread while it is supposed to measure the language threading performance. I looked yesterday at one Go source code (I think it was for pidigits) and saw that it just called a C library (gmp) to compute with big integers. It's no surprise then that Go would be faster than Java on this test.<br /><br />So what about my small Monte-Carlo test?<br />Well it turns out that Go is quite fast on it:<br /><span style="font-family: "Courier New",Courier,monospace;">Multipl. Rust Go<br />1 0.005 0.007<br />10 0.03 0.03<br />100 0.21 0.29</span><br /><span style="font-family: "Courier New",Courier,monospace;">1000 2.01 2.88</span><br /><br /><br />It is faster than Java/Scala and not too far off Rust, except if one uses FastMath in Scala, then the longest test is slighly faster with Java (not the other ones).<br /><br />There are some issues with the Go language: there is no operator overloading, which can make matrix/vector algebra more tedious and there is no generic/template. The later is somewhat mitigated by the automatic interface implementation. And fortunately for the former, complex numbers are a standard type. Still, automatic differentiation would be painful.<br /><br />Still it was extremely quick to grasp and write code, because it's so simple, especially when compared to Rust. But then, contrary to Rust, there is not as much safety provided by the language. Rust is quite impressive on this side (but unfortunately that implies less readable code). I'd say that Go could become a serious alternative to Java.<br /><br />I also found an interesting minor performance issue with the default Go <a href="https://golang.org/src/math/rand/rand.go">Rand.Float64</a>, the library convert an Int63 to a double precision number this way:<br /><br /><pre>func (r *Rand) Float64() float64 {</pre><pre> f := float64(r.Int63()) / (1 << 63)<br /> if f == 1 {<br /><span class="ln" id="L128"> </span>f = 0<br /><span class="ln" id="L129"> }</span><br /><span class="ln" id="L130"> </span>return f<br /> }</pre><br />I was interested in having a number in (0,1) and not [0,1), so I just used the conversion pattern from MersenneTwister 64 code:<br /><br /><pre>f := (float64(r.Int63() >> 11) + 0.5) * (1.0/4503599627370496.0)</pre><pre> </pre>The reasoning behind this later code is that the mantissa is 52 bits, and this is the most accuracy we can have between 0 and 1. There is no need to go further, this also avoids the issue around 1. It's also straightforward that is will preserve the uniform property, while it's not so clear to me that r.Int63()/2^63 is going to preserve uniformity as double accuracy is higher around 0 (as the exponent part can be used there) and lesser around 1: there is going to be much more multiple identical results near 1 than near 0.<br /><br />It turns out that the if check adds 5% performance penalty on this test, likely because of processor caching issues. I was surprised by that since there are many other ifs afterwards in the code, for the inverse cumulative function, and for the payoff.<img src="http://feeds.feedburner.com/~r/chasethedevil/~4/Yy95xIIOUks" height="1" width="1" alt=""/>0http://chasethedevil.blogspot.com/2015/08/go-for-monte-carlo.htmlBumping Correlationshttp://feedproxy.google.com/~r/chasethedevil/~3/X0ChtnNrkZc/bumping-correlations_25.htmlquantnoreply@blogger.com (Fabien)Thu, 17 Dec 2015 07:21:39 PSTtag:blogger.com,1999:blog-14415229.post-5778247596999238613In his book "<i>Monte Carlo Methods in Finance</i>", P. Jäckel explains a simple way to clean up a correlation matrix. When a given correlation matrix is not positive semi-definite, the idea is to do a <a href="https://en.wikipedia.org/wiki/Singular_value_decomposition">singular value decomposition</a> (SVD), replace the negative eigenvalues by 0, and renormalize the corresponding eigenvector accordingly.<br /><br />One of the cited applications is "<i>stress testing and scenario analysis for market risk</i>" or "<i>comparative pricing in order to ascertain the extent of correlation exposure for multi-asset derivatives</i>", saying that "<i>In many of these cases we end up with a matrix that is no longer positive semi-definite</i>".<br /><br />It turns out that if one bumps an invalid correlation matrix (the input), that is then cleaned up automatically, the effect can be a very different bump. Depending on how familiar you are with SVD, this could be more or less obvious from the procedure,<br /><br />As a simple illustration I take the matrix representing 3 assets A, B, C with rho_ab = -0.6, rho_ac = rho_bc = -0.5.<br /><br /> 1.00000 -0.60000 -0.50000<br /> -0.60000 1.00000 -0.50000<br /> -0.50000 -0.50000 1.00000<br /><br />For those rho_ac and rho_bc, the correlation matrix is not positive definite unless rho_ab in in the range (-0.5, 1). One way to verify this is to use the fact that positive definiteness is equivalent to a positive determinant. The determinant will be 1 - 2*0.25 - rho_ab^2 + 2*0.25*rho_ab.<br /><br />After using P. Jaeckel procedure, we end up with: <br /><br /> 1.00000 -0.56299 -0.46745<br /> -0.56299 1.00000 -0.46745<br /> -0.46745 -0.46745 1.00000<br /><br />If we bump now rho_bc by 1% (absolute), we end up after cleanup with:<br /><br /> 1.00000 -0.56637 -0.47045<br /> -0.56637 1.00000 -0.46081<br /> -0.47045 -0.46081 1.00000<br /><br />It turns out that rho_bc has changed by only 0.66% and rho_ac by -0.30%, rho_ab by -0.34%. So our initial bump (0,0,1) has been translated to a bump (-0.34, -0.30, 0.66). In other words, it does not work to compute sensitivities.<br /><br />One can optimize to obtain the nearest correlation matrix in some norm. Jaeckel proposes a hypersphere decomposition based optimization, using as initial guess the SVD solution. <a href="https://nickhigham.wordpress.com/2013/02/13/the-nearest-correlation-matrix/">Higham proposed a specific algorithm</a> just for that purpose. It turns out that on this example, they will converge to the same solution (if we use the same norm). I tried out of curiosity to see if that would lead to some improvement. The first matrix becomes<br /><br /> 1.00000 -0.56435 -0.46672<br /> -0.56435 1.00000 -0.46672<br /> -0.46672 -0.46672 1.00000<br /><br />And the bumped one becomes<br /><br /> 1.00000 -0.56766 -0.46984<br /> -0.56766 1.00000 -0.46002<br /> -0.46984 -0.46002 1.00000<br /><br />We find back the same issue, rho_bc has changed by only 0.67%, rho_ac by -0.31% and rho_ab by -0.33%. We also see that the SVD correlation or the real near correlation matrix are quite close, as noticed by P. Jaeckel.<br /><br />Of course, one should apply the bump directly to the cleaned up matrix, in which case it will actually work as expected, unless our bump produces another non positive definite matrix, and then we would have correlation leaking a bit everywhere. It's not entirely clear what kind of meaning the risk figures would have then.<img src="http://feeds.feedburner.com/~r/chasethedevil/~4/X0ChtnNrkZc" height="1" width="1" alt=""/>0http://chasethedevil.blogspot.com/2015/07/bumping-correlations_25.htmlBumping Correlationshttp://feedproxy.google.com/~r/chasethedevil/~3/BX_HRNkwJrs/bumping-correlations.htmlquantnoreply@blogger.com (Fabien)Fri, 04 Sep 2015 12:45:48 PDTtag:blogger.com,1999:blog-14415229.post-6407946176461874924In his book "<i>Monte Carlo Methods in Finance</i>", P. Jäckel explains a simple way to clean up a correlation matrix. When a given correlation matrix is not positive semi-definite, the idea is to do a <a href="https://en.wikipedia.org/wiki/Singular_value_decomposition">singular value decomposition</a> (SVD), replace the negative eigenvalues by 0, and renormalize the corresponding eigenvector accordingly.<br /><br />One of the cited applications is "<i>stress testing and scenario analysis for market risk</i>" or "<i>comparative pricing in order to ascertain the extent of correlation exposure for multi-asset derivatives</i>", saying that "<i>In many of these cases we end up with a matrix that is no longer positive semi-definite</i>".<br /><br />It turns out that if one bumps an invalid correlation matrix (the input), that is then cleaned up automatically, the effect can be a very different bump. Depending on how familiar you are with SVD, this could be more or less obvious from the procedure,<br /><br />As a simple illustration I take the matrix representing 3 assets A, B, C with rho_ab = -0.6, rho_ac = rho_bc = -0.5.<br /><br /> 1.00000 -0.60000 -0.50000<br /> -0.60000 1.00000 -0.50000<br /> -0.50000 -0.50000 1.00000<br /><br />For those rho_ac and rho_bc, the correlation matrix is not positive definite unless rho_ab in in the range (-0.5, 1). One way to verify this is to use the fact that positive definiteness is equivalent to a positive determinant. The determinant will be 1 - 2*0.25 - rho_ab^2 + 2*0.25*rho_ab.<br /><br />After using P. Jaeckel procedure, we end up with: <br /><br /> 1.00000 -0.56299 -0.46745<br /> -0.56299 1.00000 -0.46745<br /> -0.46745 -0.46745 1.00000<br /><br />If we bump now rho_bc by 1% (absolute), we end up after cleanup with:<br /><br /> 1.00000 -0.56637 -0.47045<br /> -0.56637 1.00000 -0.46081<br /> -0.47045 -0.46081 1.00000<br /><br />It turns out that rho_bc has changed by only 0.66% and rho_ac by -0.30%, rho_ab by -0.34%. So our initial bump (0,0,1) has been translated to a bump (-0.34, -0.30, 0.66). In other words, it does not work to compute sensitivities.<br /><br />One can optimize to obtain the nearest correlation matrix in some norm. Jaeckel proposes a hypersphere decomposition based optimization, using as initial guess the SVD solution. <a href="https://nickhigham.wordpress.com/2013/02/13/the-nearest-correlation-matrix/">Higham proposed a specific algorithm</a> just for that purpose. It turns out that on this example, they will converge to the same solution (if we use the same norm). I tried out of curiosity to see if that would lead to some improvement. The first matrix becomes<br /><br /> 1.00000 -0.56435 -0.46672<br /> -0.56435 1.00000 -0.46672<br /> -0.46672 -0.46672 1.00000<br /><br />And the bumped one becomes<br /><br /> 1.00000 -0.56766 -0.46984<br /> -0.56766 1.00000 -0.46002<br /> -0.46984 -0.46002 1.00000<br /><br />We find back the same issue, rho_bc has changed by only 0.67%, rho_ac by -0.31% and rho_ab by -0.33%. We also see that the SVD correlation or the real near correlation matrix are quite close, as noticed by P. Jaeckel.<br /><br />Of course, one should apply the bump directly to the cleaned up matrix, in which case it will actually work as expected, unless our bump produces another non positive definite matrix, and then we would have correlation leaking a bit everywhere. It's not entirely clear what kind of meaning the risk figures would have then.<img src="http://feeds.feedburner.com/~r/chasethedevil/~4/BX_HRNkwJrs" height="1" width="1" alt=""/>0http://chasethedevil.blogspot.com/2015/07/bumping-correlations.htmlAndreasen Huge extrapolationhttp://feedproxy.google.com/~r/chasethedevil/~3/tVXYcK9g5Y4/andreasen-huge-extrapolation_13.htmlquantnoreply@blogger.com (Fabien)Thu, 17 Dec 2015 07:21:39 PSTtag:blogger.com,1999:blog-14415229.post-935558059905277933There are not many arbitrage free extrapolation schemes. Benaim et al. extrapolation is one of the few that claims it. However, despite the paper's title, it is not truely arbitrage free. The density might be positive, but the forward is not preserved by the implied density. It can also lead to wings that don't obey Lee's moments condition.<br /><br />On a Wilmott forum, <a href="http://www.wilmott.com/messageview.cfm?catid=4&threadid=95309">P. Caspers proposed</a> the following counter-example based on extrapolating SABR: \( \alpha=15\%, \beta=80\%, \nu=50\%, \rho=-48\%, f=3\%, T=20.0 \). He cut this smile at 2.5% and 6% and used the BDK extrapolation scheme with mu=nu=1.<br /><br />A truly arbitrage free extrapolation can be obtained through <a href="http://ssrn.com/abstract=1694972">Andreasen Huge volatility interpolation</a>, making sure the grid is wide enough to allow extrapolation. Their method is basically a one step finite difference implicit Euler scheme applied to a local volatility parameterization that has as many parameters than prices. The method is presented with piecewise constant local volatility, but actually used with piecewise linear local volatility in their example.<br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://3.bp.blogspot.com/-Cp1keeZvJjY/VaPaq-_ttVI/AAAAAAAAIFU/3M-0n5N4eqw/s1600/Screenshot-Untitled%2BWindow-5.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="283" src="http://3.bp.blogspot.com/-Cp1keeZvJjY/VaPaq-_ttVI/AAAAAAAAIFU/3M-0n5N4eqw/s400/Screenshot-Untitled%2BWindow-5.png" width="400" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Smile</td></tr></tbody></table><br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://2.bp.blogspot.com/-mG5xhhikk5Y/VaPWmDv875I/AAAAAAAAIFA/X8SU8NUmEAE/s1600/Screenshot-Untitled%2BWindow-4.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="283" src="http://2.bp.blogspot.com/-mG5xhhikk5Y/VaPWmDv875I/AAAAAAAAIFA/X8SU8NUmEAE/s400/Screenshot-Untitled%2BWindow-4.png" width="400" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Density with piecewise linear local volatility</td></tr></tbody></table>There is still a tiny oscillation that makes the density negative, but one understands why typical extrapolations fail on the example: the change in density must be very steep.<br />Note that moving the left extrapolation point even closer to the forward might fix BDK negative density, but we are already very close, and we can really wonder if going closer is really a good idea since we would effectively use a somewhat arbitrary extrapolation in most of the interpolation zone.<br /><br />It turns out that we can also use a cubic spline local volatility with linear extrapolation, and the density would look then: <br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://2.bp.blogspot.com/-UytrsK6ficA/VaPWr-BEHmI/AAAAAAAAIFI/BtzCxV0MCSI/s1600/Screenshot-Untitled%2BWindow-3.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="283" src="http://2.bp.blogspot.com/-UytrsK6ficA/VaPWr-BEHmI/AAAAAAAAIFI/BtzCxV0MCSI/s400/Screenshot-Untitled%2BWindow-3.png" width="400" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Density with cubic spline local volatility</td></tr></tbody></table>Interestingly, the right side of the density is much better captured.<br />The wiggle persists, although it is smaller. This is likely due to the fact that I am using a cubic spline on top of the finite difference prices (in order to have a C2 density). Using a better C2 convexity preserving interpolation would likely remove this artefact.<br /><br />Those figures also show why relying just on extrapolation to fix SABR is not necessarily a good idea: even a real arbitrage free extrapolation will make a not so meaningful density. The proper solution is to really use <a href="http://chasethedevil.blogspot.fr/2013/12/arbitrage-free-sabr-another-view-on.html">Hagan's arbitrage free SABR PDE</a>, which would be as nearly fast in this case.<img src="http://feeds.feedburner.com/~r/chasethedevil/~4/tVXYcK9g5Y4" height="1" width="1" alt=""/>0http://chasethedevil.blogspot.com/2015/07/andreasen-huge-extrapolation_13.htmlAndreasen Huge extrapolationhttp://feedproxy.google.com/~r/chasethedevil/~3/Jr8_JMwmaZg/andreasen-huge-extrapolation.htmlquantnoreply@blogger.com (Fabien)Fri, 04 Sep 2015 12:45:48 PDTtag:blogger.com,1999:blog-14415229.post-4036359162389843003There are not many arbitrage free extrapolation schemes. Benaim et al. extrapolation is one of the few that claims it. However, despite the paper's title, it is not truely arbitrage free. The density might be positive, but the forward is not preserved by the implied density. It can also lead to wings that don't obey Lee's moments condition.<br /><br />On a Wilmott forum, <a href="http://www.wilmott.com/messageview.cfm?catid=4&threadid=95309">P. Caspers proposed</a> the following counter-example based on extrapolating SABR: \( \alpha=15\%, \beta=80\%, \nu=50\%, \rho=-48\%, f=3\%, T=20.0 \). He cut this smile at 2.5% and 6% and used the BDK extrapolation scheme with mu=nu=1.<br /><br />A truly arbitrage free extrapolation can be obtained through <a href="http://ssrn.com/abstract=1694972">Andreasen Huge volatility interpolation</a>, making sure the grid is wide enough to allow extrapolation. Their method is basically a one step finite difference implicit Euler scheme applied to a local volatility parameterization that has as many parameters than prices. The method is presented with piecewise constant local volatility, but actually used with piecewise linear local volatility in their example.<br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://3.bp.blogspot.com/-Cp1keeZvJjY/VaPaq-_ttVI/AAAAAAAAIFU/3M-0n5N4eqw/s1600/Screenshot-Untitled%2BWindow-5.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="283" src="http://3.bp.blogspot.com/-Cp1keeZvJjY/VaPaq-_ttVI/AAAAAAAAIFU/3M-0n5N4eqw/s400/Screenshot-Untitled%2BWindow-5.png" width="400" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Smile</td></tr></tbody></table><br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://2.bp.blogspot.com/-mG5xhhikk5Y/VaPWmDv875I/AAAAAAAAIFA/X8SU8NUmEAE/s1600/Screenshot-Untitled%2BWindow-4.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="283" src="http://2.bp.blogspot.com/-mG5xhhikk5Y/VaPWmDv875I/AAAAAAAAIFA/X8SU8NUmEAE/s400/Screenshot-Untitled%2BWindow-4.png" width="400" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Density with piecewise linear local volatility</td></tr></tbody></table>There is still a tiny oscillation that makes the density negative, but one understands why typical extrapolations fail on the example: the change in density must be very steep.<br />Note that moving the left extrapolation point even closer to the forward might fix BDK negative density, but we are already very close, and we can really wonder if going closer is really a good idea since we would effectively use a somewhat arbitrary extrapolation in most of the interpolation zone.<br /><br />It turns out that we can also use a cubic spline local volatility with linear extrapolation, and the density would look then: <br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://2.bp.blogspot.com/-UytrsK6ficA/VaPWr-BEHmI/AAAAAAAAIFI/BtzCxV0MCSI/s1600/Screenshot-Untitled%2BWindow-3.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="283" src="http://2.bp.blogspot.com/-UytrsK6ficA/VaPWr-BEHmI/AAAAAAAAIFI/BtzCxV0MCSI/s400/Screenshot-Untitled%2BWindow-3.png" width="400" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Density with cubic spline local volatility</td></tr></tbody></table>Interestingly, the right side of the density is much better captured.<br />The wiggle persists, although it is smaller. This is likely due to the fact that I am using a cubic spline on top of the finite difference prices (in order to have a C2 density). Using a better C2 convexity preserving interpolation would likely remove this artefact.<br /><br />Those figures also show why relying just on extrapolation to fix SABR is not necessarily a good idea: even a real arbitrage free extrapolation will make a not so meaningful density. The proper solution is to really use <a href="http://chasethedevil.blogspot.fr/2013/12/arbitrage-free-sabr-another-view-on.html">Hagan's arbitrage free SABR PDE</a>, which would be as nearly fast in this case.<img src="http://feeds.feedburner.com/~r/chasethedevil/~4/Jr8_JMwmaZg" height="1" width="1" alt=""/>0http://chasethedevil.blogspot.com/2015/07/andreasen-huge-extrapolation.htmlUnintuitive behavior of the Black-Scholes formula - negative volatilities in displaced diffusion extrapolationhttp://feedproxy.google.com/~r/chasethedevil/~3/9WArbO4EffM/unintuitive-behavior-of-black-scholes_7.htmlquantnoreply@blogger.com (Fabien)Thu, 17 Dec 2015 07:21:39 PSTtag:blogger.com,1999:blog-14415229.post-6441762317525776208I am looking at various extrapolation schemes of the implied volatilities. An interesting one I stumbled upon is due to Kahale. Even if <a href="http://nkahale.free.fr/papers/Interpolation.pdf">his paper</a> is on interpolation, there is actually a small paragraph on using the same kind of function for extrapolation. His idea is to simply lookup the standard deviation \( \Sigma \) and the forward \(f\) corresponding to a given market volatility and slope:<br />$$ c_{f,\Sigma} = f N(d_1) - k N(d_2)$$<br />with<br />$$ d_1 = \frac{\log(f/k)+\Sigma^2/2}{\Sigma} $$<br />We have simply:<br />$$ c'(k) = - N(d_2)$$<br /><br />He also proves that we can always find those two parameters for any \( k_0 > 0, c_0 > 0, -1 < c_0' < 0 \)<br /><br />Then I had the silly idea of trying to match with a put instead of a call for the left wing (as those are out-of-the-money, and therefore easier to invert numerically). It turns out that it works in most cases in practice and produces relatively nice looking extrapolations, but it does not always work. This is because contrary to the call, the put value is bounded with \(f\).<br />$$ p_{f,\Sigma} = k N(-d_2) - f N(-d_1)$$ <br />Inverting \(p_0'\) is going to lead to a specific \(d_2\), and you are not guaranteed that you can push \(f\) high and have \(p_{f,\Sigma}\) large enough to match \(p_0\). As example we can just take \(p_0 \geq k N(-d_2)\) which will only be matched if \(f \leq 0\).<br /><br />This is slightly unintuitive as put-call parity would suggest some kind of equivalence. The problem here is that we would need to consider the function of \(k\) instead of \(f\) for it to work, so we can't really work with a put directly.<br /><br />Here are the two different extrapolations on Kahale own example:<br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://4.bp.blogspot.com/-UcomxGsx_r0/VZvSdzoBBVI/AAAAAAAAIEU/_V542xidbgc/s1600/Screenshot-Untitled%2BWindow.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="297" src="http://4.bp.blogspot.com/-UcomxGsx_r0/VZvSdzoBBVI/AAAAAAAAIEU/_V542xidbgc/s400/Screenshot-Untitled%2BWindow.png" width="400" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Extrapolation of the left wing with calls (blue doted line)</td></tr></tbody></table><br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://4.bp.blogspot.com/-sDl37fFAImE/VZvSd5eWCgI/AAAAAAAAIEQ/tOUG7nHrg_Q/s1600/Screenshot-Untitled%2BWindow-1.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="297" src="http://4.bp.blogspot.com/-sDl37fFAImE/VZvSd5eWCgI/AAAAAAAAIEQ/tOUG7nHrg_Q/s400/Screenshot-Untitled%2BWindow-1.png" width="400" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Extrapolation of the left wing with puts (blue doted line)</td></tr></tbody></table>Displaced diffusion extrapolation is sometimes advocated. It is not the same as Kahale extrapolation: In Kahale, only the forward variable is varying in the Black-Scholes formula, and there is no real underlying stochastic process. In a displaced diffusion setting, we would adjust both strike and forward, keeping put-call parity at the formula level. But unfortunately, it suffers from the same kind of problem: it can not always be solved for slope and price. When it can however, it will give a more consistent extrapolation.<br /><br />I find it interesting that some smiles can not be extrapolated by displaced diffusion in a C1 manner except if one allows negative volatilities in the formula (in which case we are not anymore in a pure displaced diffusion setting).<br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://4.bp.blogspot.com/-VSy7uTzu56U/VZwJUDT5iJI/AAAAAAAAIEo/lsi8EakZ-kA/s1600/Screenshot-Untitled%2BWindow-2.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="297" src="http://4.bp.blogspot.com/-VSy7uTzu56U/VZwJUDT5iJI/AAAAAAAAIEo/lsi8EakZ-kA/s400/Screenshot-Untitled%2BWindow-2.png" width="400" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Extrapolation of the left wing using negative displaced diffusion volatilities (blue dotted line)</td></tr></tbody></table><img src="http://feeds.feedburner.com/~r/chasethedevil/~4/9WArbO4EffM" height="1" width="1" alt=""/>0http://chasethedevil.blogspot.com/2015/07/unintuitive-behavior-of-black-scholes_7.htmlUnintuitive behavior of the Black-Scholes formula - negative volatilities in displaced diffusion extrapolationhttp://feedproxy.google.com/~r/chasethedevil/~3/KhZgkTfmBrs/unintuitive-behavior-of-black-scholes.htmlquantnoreply@blogger.com (Fabien)Fri, 04 Sep 2015 12:45:48 PDTtag:blogger.com,1999:blog-14415229.post-2350170856749440306I am looking at various extrapolation schemes of the implied volatilities. An interesting one I stumbled upon is due to Kahale. Even if <a href="http://nkahale.free.fr/papers/Interpolation.pdf">his paper</a> is on interpolation, there is actually a small paragraph on using the same kind of function for extrapolation. His idea is to simply lookup the standard deviation \( \Sigma \) and the forward \(f\) corresponding to a given market volatility and slope:<br />$$ c_{f,\Sigma} = f N(d_1) - k N(d_2)$$<br />with<br />$$ d_1 = \frac{\log(f/k)+\Sigma^2/2}{\Sigma} $$<br />We have simply:<br />$$ c'(k) = - N(d_2)$$<br /><br />He also proves that we can always find those two parameters for any \( k_0 > 0, c_0 > 0, -1 < c_0' < 0 \)<br /><br />Then I had the silly idea of trying to match with a put instead of a call for the left wing (as those are out-of-the-money, and therefore easier to invert numerically). It turns out that it works in most cases in practice and produces relatively nice looking extrapolations, but it does not always work. This is because contrary to the call, the put value is bounded with \(f\).<br />$$ p_{f,\Sigma} = k N(-d_2) - f N(-d_1)$$ <br />Inverting \(p_0'\) is going to lead to a specific \(d_2\), and you are not guaranteed that you can push \(f\) high and have \(p_{f,\Sigma}\) large enough to match \(p_0\). As example we can just take \(p_0 \geq k N(-d_2)\) which will only be matched if \(f \leq 0\).<br /><br />This is slightly unintuitive as put-call parity would suggest some kind of equivalence. The problem here is that we would need to consider the function of \(k\) instead of \(f\) for it to work, so we can't really work with a put directly.<br /><br />Here are the two different extrapolations on Kahale own example:<br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://4.bp.blogspot.com/-UcomxGsx_r0/VZvSdzoBBVI/AAAAAAAAIEU/_V542xidbgc/s1600/Screenshot-Untitled%2BWindow.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="297" src="http://4.bp.blogspot.com/-UcomxGsx_r0/VZvSdzoBBVI/AAAAAAAAIEU/_V542xidbgc/s400/Screenshot-Untitled%2BWindow.png" width="400" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Extrapolation of the left wing with calls (blue doted line)</td></tr></tbody></table><br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://4.bp.blogspot.com/-sDl37fFAImE/VZvSd5eWCgI/AAAAAAAAIEQ/tOUG7nHrg_Q/s1600/Screenshot-Untitled%2BWindow-1.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="297" src="http://4.bp.blogspot.com/-sDl37fFAImE/VZvSd5eWCgI/AAAAAAAAIEQ/tOUG7nHrg_Q/s400/Screenshot-Untitled%2BWindow-1.png" width="400" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Extrapolation of the left wing with puts (blue doted line)</td></tr></tbody></table>Displaced diffusion extrapolation is sometimes advocated. It is not the same as Kahale extrapolation: In Kahale, only the forward variable is varying in the Black-Scholes formula, and there is no real underlying stochastic process. In a displaced diffusion setting, we would adjust both strike and forward, keeping put-call parity at the formula level. But unfortunately, it suffers from the same kind of problem: it can not always be solved for slope and price. When it can however, it will give a more consistent extrapolation.<br /><br />I find it interesting that some smiles can not be extrapolated by displaced diffusion in a C1 manner except if one allows negative volatilities in the formula (in which case we are not anymore in a pure displaced diffusion setting).<br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://4.bp.blogspot.com/-VSy7uTzu56U/VZwJUDT5iJI/AAAAAAAAIEo/lsi8EakZ-kA/s1600/Screenshot-Untitled%2BWindow-2.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="297" src="http://4.bp.blogspot.com/-VSy7uTzu56U/VZwJUDT5iJI/AAAAAAAAIEo/lsi8EakZ-kA/s400/Screenshot-Untitled%2BWindow-2.png" width="400" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Extrapolation of the left wing using negative displaced diffusion volatilities (blue dotted line)</td></tr></tbody></table><img src="http://feeds.feedburner.com/~r/chasethedevil/~4/KhZgkTfmBrs" height="1" width="1" alt=""/>0http://chasethedevil.blogspot.com/2015/07/unintuitive-behavior-of-black-scholes.htmlLinux Desktops in 2015http://feedproxy.google.com/~r/chasethedevil/~3/jogGBSbBkWc/linux-desktops-in-2015_24.htmllinuxnoreply@blogger.com (Fabien)Thu, 17 Dec 2015 07:21:39 PSTtag:blogger.com,1999:blog-14415229.post-1087223011543300990I seem to <a href="http://chasethedevil.blogspot.fr/2014/05/kde-xfce-gnome-shell-in-2014.html">never be entirely happy with any of the linux desktops</a> these days. I have used XFCE on Ubuntu quite a bit in the past year, it mostly works, but I still had minor annoyances:<br />- sometimes (rarely) my laptop would not wake up from sleep.<br />- notifications sometimes keep popping up too much.<br />- on my desktop, experienced strong tearing issues with the Radeon graphic card, except with some very specific combination of video player settings and desktop settings (and then I had annoying redraw issue when pushing volume up/down in movies).<br /><br />I am satisfied with two different approaches since:<br />- OpenSuse 13.2 with KDE 4. I use that on my desktop, all issues are gone, and the integration of KDE in OpenSuse is clearly the best I have experienced. In contrast, KDE 5 on Ubuntu was a disaster for me. I also managed to fuck up the apt dependencies so much that I thought it would be simpler to reinstall a new distribution.<br />- Mate on Ubuntu 15.04. Very impressed so far. It's probably what Gnome should have been instead of going to 3.0. Even if there are nice aspects of the Gnome shell, Mate is fast, pretty, user friendly, much better than Cinnamon. There are even a few layouts to choose (most of them are good), here is "Eleven with Mate menu" (it installed and setup the Plank dock automatically for that layout, more traditional layouts without dock are available):<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-aCnzz5ujbA4/VYrfl3lhz1I/AAAAAAAAIDw/5qV5NrGn7a8/s1600/Mate-Eleven.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="358" src="http://1.bp.blogspot.com/-aCnzz5ujbA4/VYrfl3lhz1I/AAAAAAAAIDw/5qV5NrGn7a8/s640/Mate-Eleven.png" width="640" /></a></div><br /><img src="http://feeds.feedburner.com/~r/chasethedevil/~4/jogGBSbBkWc" height="1" width="1" alt=""/>0http://chasethedevil.blogspot.com/2015/06/linux-desktops-in-2015_24.htmlLinux Desktops in 2015http://feedproxy.google.com/~r/chasethedevil/~3/kaosf7GIs-M/linux-desktops-in-2015.htmllinuxnoreply@blogger.com (Fabien)Fri, 04 Sep 2015 12:46:04 PDTtag:blogger.com,1999:blog-14415229.post-911525716485680250I seem to <a href="http://chasethedevil.blogspot.fr/2014/05/kde-xfce-gnome-shell-in-2014.html">never be entirely happy with any of the linux desktops</a> these days. I have used XFCE on Ubuntu quite a bit in the past year, it mostly works, but I still had minor annoyances:<br />- sometimes (rarely) my laptop would not wake up from sleep.<br />- notifications sometimes keep popping up too much.<br />- on my desktop, experienced strong tearing issues with the Radeon graphic card, except with some very specific combination of video player settings and desktop settings (and then I had annoying redraw issue when pushing volume up/down in movies).<br /><br />I am satisfied with two different approaches since:<br />- OpenSuse 13.2 with KDE 4. I use that on my desktop, all issues are gone, and the integration of KDE in OpenSuse is clearly the best I have experienced. In contrast, KDE 5 on Ubuntu was a disaster for me. I also managed to fuck up the apt dependencies so much that I thought it would be simpler to reinstall a new distribution.<br />- Mate on Ubuntu 15.04. Very impressed so far. It's probably what Gnome should have been instead of going to 3.0. Even if there are nice aspects of the Gnome shell, Mate is fast, pretty, user friendly, much better than Cinnamon. There are even a few layouts to choose (most of them are good), here is "Eleven with Mate menu" (it installed and setup the Plank dock automatically for that layout, more traditional layouts without dock are available):<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-aCnzz5ujbA4/VYrfl3lhz1I/AAAAAAAAIDw/5qV5NrGn7a8/s1600/Mate-Eleven.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="358" src="http://1.bp.blogspot.com/-aCnzz5ujbA4/VYrfl3lhz1I/AAAAAAAAIDw/5qV5NrGn7a8/s640/Mate-Eleven.png" width="640" /></a></div><br /><img src="http://feeds.feedburner.com/~r/chasethedevil/~4/kaosf7GIs-M" height="1" width="1" alt=""/>0http://chasethedevil.blogspot.com/2015/06/linux-desktops-in-2015.htmlSquare Root Crank-Nicolsonhttp://feedproxy.google.com/~r/chasethedevil/~3/jjBn3fQCu3M/square-root-crank-nicolson_19.htmlquantnoreply@blogger.com (Fabien)Thu, 17 Dec 2015 07:21:39 PSTtag:blogger.com,1999:blog-14415229.post-1762372781056647052C. Reisinger kindly pointed out to me <a href="http://arxiv.org/abs/1210.5487">this paper around square root Crank-Nicolson</a>. The idea is to apply a square root of time transformation to the PDE, and discretize the resulting PDE with Crank-Nicolson. Two reasons come to mind to try this: <br /><ul><li>the square root transform will result in small steps initially, where the solution is potentially not so smooth, making Crank-Nicolson behave better.</li><li> it is the natural time of the Brownian motion.</li></ul>Interestingly, it has nicer properties than what those reasons may suggest. On the Fokker-Planck density PDE, it does not oscillate under some very mild conditions and <a href="http://papers.ssrn.com/sol3/papers.cfm?abstract_id=2605160">preserves density positivity at the peak</a>.<br /><br />Out of curiosity I tried it to price a one touch barrier option. Of course there is an analytical solution in my test case (Black-Scholes assumptions), but as soon as rates are assumed not constant or local volatility is used, there is no other solution than a numerical method. In the later case, finite difference methods are quite good in terms of performance vs accuracy.<br /><br />The classic Crank-Nicolson gives a reasonable price, but the strong oscillations near the barrier, at every time step are not very comforting.<br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://3.bp.blogspot.com/-XNr7vE77Dfo/VYQoBELx2KI/AAAAAAAAICw/FVVYrehW39Y/s1600/Screenshot%2Bfrom%2B2015-06-19%2B16%253A24%253A59.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="337" src="http://3.bp.blogspot.com/-XNr7vE77Dfo/VYQoBELx2KI/AAAAAAAAICw/FVVYrehW39Y/s640/Screenshot%2Bfrom%2B2015-06-19%2B16%253A24%253A59.png" width="640" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Crank-Nicolson Prices near the Barrier. Each line is a different time.</td></tr></tbody></table><br />Moving to square root of time removes nearly all oscillations on this problem, even with a relatively low number of time steps compared to the number of space steps.<br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://1.bp.blogspot.com/-JXVRuhLrMOQ/VYQodTFBDGI/AAAAAAAAIC4/hrsSdQbA5Wo/s1600/Screenshot%2Bfrom%2B2015-06-19%2B16%253A25%253A14.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="338" src="http://1.bp.blogspot.com/-JXVRuhLrMOQ/VYQodTFBDGI/AAAAAAAAIC4/hrsSdQbA5Wo/s640/Screenshot%2Bfrom%2B2015-06-19%2B16%253A25%253A14.png" width="640" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Square Root Crank-Nicolson Prices near the Barrier. Each line is a different time.</td></tr></tbody></table><br />We can see that the second step prices are a bit higher than the third step (the lines cross), which looks like a small numerical oscillation in time, even if there is no oscillation is space.<br /><br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://1.bp.blogspot.com/-d-xXHOvO1H0/VYQon2V5eiI/AAAAAAAAIDA/4g-YYby4R0A/s1600/Screenshot%2Bfrom%2B2015-06-19%2B16%253A25%253A06.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="338" src="http://1.bp.blogspot.com/-d-xXHOvO1H0/VYQon2V5eiI/AAAAAAAAIDA/4g-YYby4R0A/s640/Screenshot%2Bfrom%2B2015-06-19%2B16%253A25%253A06.png" width="640" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">TR-BDF2 Prices near the Barrier. Each line is a different time.</td></tr></tbody></table><br />As a comparison, the TR-BDF2 scheme does relatively well: oscillations are removed after the second step, even with the extreme ratio of time steps vs space steps used on this example so that illustrations are clearer - Crank-Nicolson would still oscillate a lot with 10 times less space steps but we would not see oscillation on the square root Crank-Nicolson and a very mild one on TR-BDF2.<br /><br />The LMG2 scheme (a local richardson extrapolation) does not oscillate at all on this problem but is the slowest:<br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://3.bp.blogspot.com/-HTxKtzO8au4/VYQo8t-wkdI/AAAAAAAAIDI/oPOjJWs0SI0/s1600/Screenshot%2Bfrom%2B2015-06-19%2B16%253A25%253A53.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="338" src="http://3.bp.blogspot.com/-HTxKtzO8au4/VYQo8t-wkdI/AAAAAAAAIDI/oPOjJWs0SI0/s640/Screenshot%2Bfrom%2B2015-06-19%2B16%253A25%253A53.png" width="640" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">LMG2 Prices near the Barrier. Each line is a different time.</td></tr></tbody></table><br />The square root Crank-Nicolson is quite elegant. It can however not be applied to that many problems in practice, as often some grid times are imposed by the payoff to evaluate, for example in a case of a discrete weekly barrier. But for continuous time problems (density PDE, Vanilla, American, continuous barriers) it's quite good.<br /><br />In reality, with a continuous barrier, the payoff is not discontinuous at every step, but it is only discontinuous at the first step. So Rannacher smoothing would work very well on that problem:<br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://4.bp.blogspot.com/-A-qqKczuefQ/VYQwf2ba_MI/AAAAAAAAIDY/2cYpi_3Y_pI/s1600/Screenshot%2Bfrom%2B2015-06-19%2B17%253A08%253A35.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="338" src="http://4.bp.blogspot.com/-A-qqKczuefQ/VYQwf2ba_MI/AAAAAAAAIDY/2cYpi_3Y_pI/s640/Screenshot%2Bfrom%2B2015-06-19%2B17%253A08%253A35.png" width="640" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Rannacher Prices near the Barrier. Each line is a different time.</td></tr></tbody></table>The somewhat interesting payoff left for the square root Crank-Nicolson is the American.<br /><img src="http://feeds.feedburner.com/~r/chasethedevil/~4/jjBn3fQCu3M" height="1" width="1" alt=""/>0http://chasethedevil.blogspot.com/2015/06/square-root-crank-nicolson_19.htmlSquare Root Crank-Nicolsonhttp://feedproxy.google.com/~r/chasethedevil/~3/WE4pRhPcZ50/square-root-crank-nicolson.htmlquantnoreply@blogger.com (Fabien)Fri, 04 Sep 2015 12:45:48 PDTtag:blogger.com,1999:blog-14415229.post-6232965918131989202C. Reisinger kindly pointed out to me <a href="http://arxiv.org/abs/1210.5487">this paper around square root Crank-Nicolson</a>. The idea is to apply a square root of time transformation to the PDE, and discretize the resulting PDE with Crank-Nicolson. Two reasons come to mind to try this: <br /><ul><li>the square root transform will result in small steps initially, where the solution is potentially not so smooth, making Crank-Nicolson behave better.</li><li> it is the natural time of the Brownian motion.</li></ul>Interestingly, it has nicer properties than what those reasons may suggest. On the Fokker-Planck density PDE, it does not oscillate under some very mild conditions and <a href="http://papers.ssrn.com/sol3/papers.cfm?abstract_id=2605160">preserves density positivity at the peak</a>.<br /><br />Out of curiosity I tried it to price a one touch barrier option. Of course there is an analytical solution in my test case (Black-Scholes assumptions), but as soon as rates are assumed not constant or local volatility is used, there is no other solution than a numerical method. In the later case, finite difference methods are quite good in terms of performance vs accuracy.<br /><br />The classic Crank-Nicolson gives a reasonable price, but the strong oscillations near the barrier, at every time step are not very comforting.<br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://3.bp.blogspot.com/-XNr7vE77Dfo/VYQoBELx2KI/AAAAAAAAICw/FVVYrehW39Y/s1600/Screenshot%2Bfrom%2B2015-06-19%2B16%253A24%253A59.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="337" src="http://3.bp.blogspot.com/-XNr7vE77Dfo/VYQoBELx2KI/AAAAAAAAICw/FVVYrehW39Y/s640/Screenshot%2Bfrom%2B2015-06-19%2B16%253A24%253A59.png" width="640" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Crank-Nicolson Prices near the Barrier. Each line is a different time.</td></tr></tbody></table><br />Moving to square root of time removes nearly all oscillations on this problem, even with a relatively low number of time steps compared to the number of space steps.<br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://1.bp.blogspot.com/-JXVRuhLrMOQ/VYQodTFBDGI/AAAAAAAAIC4/hrsSdQbA5Wo/s1600/Screenshot%2Bfrom%2B2015-06-19%2B16%253A25%253A14.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="338" src="http://1.bp.blogspot.com/-JXVRuhLrMOQ/VYQodTFBDGI/AAAAAAAAIC4/hrsSdQbA5Wo/s640/Screenshot%2Bfrom%2B2015-06-19%2B16%253A25%253A14.png" width="640" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Square Root Crank-Nicolson Prices near the Barrier. Each line is a different time.</td></tr></tbody></table><br />We can see that the second step prices are a bit higher than the third step (the lines cross), which looks like a small numerical oscillation in time, even if there is no oscillation is space.<br /><br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://1.bp.blogspot.com/-d-xXHOvO1H0/VYQon2V5eiI/AAAAAAAAIDA/4g-YYby4R0A/s1600/Screenshot%2Bfrom%2B2015-06-19%2B16%253A25%253A06.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="338" src="http://1.bp.blogspot.com/-d-xXHOvO1H0/VYQon2V5eiI/AAAAAAAAIDA/4g-YYby4R0A/s640/Screenshot%2Bfrom%2B2015-06-19%2B16%253A25%253A06.png" width="640" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">TR-BDF2 Prices near the Barrier. Each line is a different time.</td></tr></tbody></table><br />As a comparison, the TR-BDF2 scheme does relatively well: oscillations are removed after the second step, even with the extreme ratio of time steps vs space steps used on this example so that illustrations are clearer - Crank-Nicolson would still oscillate a lot with 10 times less space steps but we would not see oscillation on the square root Crank-Nicolson and a very mild one on TR-BDF2.<br /><br />The LMG2 scheme (a local richardson extrapolation) does not oscillate at all on this problem but is the slowest:<br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://3.bp.blogspot.com/-HTxKtzO8au4/VYQo8t-wkdI/AAAAAAAAIDI/oPOjJWs0SI0/s1600/Screenshot%2Bfrom%2B2015-06-19%2B16%253A25%253A53.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="338" src="http://3.bp.blogspot.com/-HTxKtzO8au4/VYQo8t-wkdI/AAAAAAAAIDI/oPOjJWs0SI0/s640/Screenshot%2Bfrom%2B2015-06-19%2B16%253A25%253A53.png" width="640" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">LMG2 Prices near the Barrier. Each line is a different time.</td></tr></tbody></table><br />The square root Crank-Nicolson is quite elegant. It can however not be applied to that many problems in practice, as often some grid times are imposed by the payoff to evaluate, for example in a case of a discrete weekly barrier. But for continuous time problems (density PDE, Vanilla, American, continuous barriers) it's quite good.<br /><br />In reality, with a continuous barrier, the payoff is not discontinuous at every step, but it is only discontinuous at the first step. So Rannacher smoothing would work very well on that problem:<br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://4.bp.blogspot.com/-A-qqKczuefQ/VYQwf2ba_MI/AAAAAAAAIDY/2cYpi_3Y_pI/s1600/Screenshot%2Bfrom%2B2015-06-19%2B17%253A08%253A35.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="338" src="http://4.bp.blogspot.com/-A-qqKczuefQ/VYQwf2ba_MI/AAAAAAAAIDY/2cYpi_3Y_pI/s640/Screenshot%2Bfrom%2B2015-06-19%2B17%253A08%253A35.png" width="640" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Rannacher Prices near the Barrier. Each line is a different time.</td></tr></tbody></table>The somewhat interesting payoff left for the square root Crank-Nicolson is the American.<br /><img src="http://feeds.feedburner.com/~r/chasethedevil/~4/WE4pRhPcZ50" height="1" width="1" alt=""/>0http://chasethedevil.blogspot.com/2015/06/square-root-crank-nicolson.htmlDecoding Hagan's arbitrage free SABR PDE derivationhttp://feedproxy.google.com/~r/chasethedevil/~3/D1NJAiZvLo8/decoding-hagan-arbitrage-free-sabr-pde.htmlquantnoreply@blogger.com (Fabien)Thu, 17 Dec 2015 07:21:39 PSTtag:blogger.com,1999:blog-14415229.post-7824797914935689209<br />Here are the main steps of Hagan derivation. Let's recall his notation for the SABR model where typically, \(C(F) = F^\beta\)<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-uq2IhJPDd7M/VUzC4Hh3xoI/AAAAAAAAH9k/sY034iAD38Y/s1600/Screenshot_2015-05-08_16-04-27.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="http://1.bp.blogspot.com/-uq2IhJPDd7M/VUzC4Hh3xoI/AAAAAAAAH9k/sY034iAD38Y/s1600/Screenshot_2015-05-08_16-04-27.png" /></a></div>First, he defines the moments of stochastic volatility:<br /><div class="separator" style="clear: both; text-align: center;"></div><div class="separator" style="clear: both; text-align: center;"></div><div class="separator" style="clear: both; text-align: center;"><a href="http://2.bp.blogspot.com/-RxzKVIrxbl8/VUzC4HGhMNI/AAAAAAAAH9c/7FmxbMuB4kw/s1600/Screenshot_2015-05-08_16-04-53.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="60" src="http://2.bp.blogspot.com/-RxzKVIrxbl8/VUzC4HGhMNI/AAAAAAAAH9c/7FmxbMuB4kw/s320/Screenshot_2015-05-08_16-04-53.png" width="320" /></a></div>Then he integrates the Fokker-Planck equation over all A, to obtain<br /><div class="separator" style="clear: both; text-align: center;"><a href="http://4.bp.blogspot.com/-mDsC-Zd6FMA/VUzC4L9dF_I/AAAAAAAAH9g/cAvRNw1VTkg/s1600/Screenshot_2015-05-08_16-05-36.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="http://4.bp.blogspot.com/-mDsC-Zd6FMA/VUzC4L9dF_I/AAAAAAAAH9g/cAvRNw1VTkg/s1600/Screenshot_2015-05-08_16-05-36.png" /></a></div>On the backward Komolgorov equation, he applies a Lamperti transform like change of variable:<br /><div class="separator" style="clear: both; text-align: center;"><a href="http://2.bp.blogspot.com/-QiYromr1SDU/VUzE5q_OfjI/AAAAAAAAH94/nfbr14tEnj0/s1600/Screenshot_2015-05-08_16-10-33.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="http://2.bp.blogspot.com/-QiYromr1SDU/VUzE5q_OfjI/AAAAAAAAH94/nfbr14tEnj0/s1600/Screenshot_2015-05-08_16-10-33.png" /></a></div>And then makes another change of variable so that the PDE has the same initial conditions for all moments: <br /><div class="separator" style="clear: both; text-align: center;"><a href="http://2.bp.blogspot.com/-1C8j58UD1lA/VUzE5hZJ95I/AAAAAAAAH-I/541DaJBFbAU/s1600/Screenshot_2015-05-08_16-12-34.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="http://2.bp.blogspot.com/-1C8j58UD1lA/VUzE5hZJ95I/AAAAAAAAH-I/541DaJBFbAU/s1600/Screenshot_2015-05-08_16-12-34.png" /></a></div> This leads to<br /><div class="separator" style="clear: both; text-align: center;"><a href="http://3.bp.blogspot.com/--wj3nYsi9g8/VUzE5glQ0VI/AAAAAAAAH98/NlVSXIc2NDI/s1600/Screenshot_2015-05-08_16-13-32.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="107" src="http://3.bp.blogspot.com/--wj3nYsi9g8/VUzE5glQ0VI/AAAAAAAAH98/NlVSXIc2NDI/s320/Screenshot_2015-05-08_16-13-32.png" width="320" /></a></div>It turns out that there is a magical symmetry for k=0 and k=2.<br /><div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-OStv8D0OSw0/VUzGJ6HKqDI/AAAAAAAAH-Y/jk85U57eHnY/s1600/Screenshot_2015-05-08_16-19-45.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="32" src="http://1.bp.blogspot.com/-OStv8D0OSw0/VUzGJ6HKqDI/AAAAAAAAH-Y/jk85U57eHnY/s320/Screenshot_2015-05-08_16-19-45.png" width="320" /></a></div><div class="separator" style="clear: both; text-align: center;"><a href="http://3.bp.blogspot.com/-67k-xTphi7Q/VUzGJwbwSWI/AAAAAAAAH-U/Faz235QNpKs/s1600/Screenshot_2015-05-08_16-20-09.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="32" src="http://3.bp.blogspot.com/-67k-xTphi7Q/VUzGJwbwSWI/AAAAAAAAH-U/Faz235QNpKs/s320/Screenshot_2015-05-08_16-20-09.png" width="320" /></a></div>Note that in the second equation, the second derivative applies to the whole.<br />Because of this, he can express \(Q^{(2)}\) in terms of \(Q^{(0)}\):<br /><div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-k-3DQd9MuIQ/VUzHl5TSuAI/AAAAAAAAH-s/CjTFtlX4upw/s1600/Screenshot_2015-05-08_16-25-36.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="29" src="http://1.bp.blogspot.com/-k-3DQd9MuIQ/VUzHl5TSuAI/AAAAAAAAH-s/CjTFtlX4upw/s320/Screenshot_2015-05-08_16-25-36.png" width="320" /></a></div>And he plugs that back to the integrated Fokker-Planck equation to obtain the arbitrage free SABR PDE:<br /><div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-aEwTOxulsWI/VUzHl3qfZAI/AAAAAAAAH-o/qQZPb8Vvz7o/s1600/Screenshot_2015-05-08_16-25-48.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="29" src="http://1.bp.blogspot.com/-aEwTOxulsWI/VUzHl3qfZAI/AAAAAAAAH-o/qQZPb8Vvz7o/s320/Screenshot_2015-05-08_16-25-48.png" width="320" /></a></div><br />There is a simple more common explanation in the world of local stochastic volatility for what's going on. For example, in the particle method paper from Guyon-Labordère, we have the following expression for the true local volatility.<br /><div class="separator" style="clear: both; text-align: center;"><a href="http://3.bp.blogspot.com/-uV792mNz4Xo/Ul6baEs2ktI/AAAAAAAAG1U/8Iv1i23oXok/s1600/Screenshot%2Bfrom%2B2013-10-16%2B15%3A51%3A46.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="127" src="http://3.bp.blogspot.com/-uV792mNz4Xo/Ul6baEs2ktI/AAAAAAAAG1U/8Iv1i23oXok/s320/Screenshot%2Bfrom%2B2013-10-16%2B15%3A51%3A46.png" width="320" /></a></div><br />In the first equation, the numerator is simply \(Q^{(2)}\) and the denominator \(Q^{(0)}\). Of course, the integrated Fokker-Planck equation can be rewritten as:<br /><br />$$Q^{(0)}_T = \frac{1}{2}\epsilon^2 \left[C^2(F) \frac{Q^{(2)}}{Q^{(0)}} Q^{(0)}\right]_{FF}$$<br /><br />Karlsmark uses that approach directly in his thesis, using the expansions of Doust for \(Q^{(k)}\). Looking a Doust expansions, the fraction reduces straightforwardly to the same expression as Hagan, and the symmetry in the equations appears a bit less coincidental. <br /><br /><br /><img src="http://feeds.feedburner.com/~r/chasethedevil/~4/D1NJAiZvLo8" height="1" width="1" alt=""/>0http://chasethedevil.blogspot.com/2015/05/decoding-hagan-arbitrage-free-sabr-pde.htmlDecoding Hagan's arbitrage free SABR PDE derivationhttp://feedproxy.google.com/~r/chasethedevil/~3/GmgtnJAUJ0Q/decoding-hagans-arbitrage-free-sabr-pde.htmlquantnoreply@blogger.com (Fabien)Fri, 04 Sep 2015 12:45:48 PDTtag:blogger.com,1999:blog-14415229.post-4603128956043803575<br />Here are the main steps of Hagan derivation. Let's recall his notation for the SABR model where typically, \(C(F) = F^\beta\)<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-uq2IhJPDd7M/VUzC4Hh3xoI/AAAAAAAAH9k/sY034iAD38Y/s1600/Screenshot_2015-05-08_16-04-27.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="http://1.bp.blogspot.com/-uq2IhJPDd7M/VUzC4Hh3xoI/AAAAAAAAH9k/sY034iAD38Y/s1600/Screenshot_2015-05-08_16-04-27.png" /></a></div>First, he defines the moments of stochastic volatility:<br /><div class="separator" style="clear: both; text-align: center;"></div><div class="separator" style="clear: both; text-align: center;"></div><div class="separator" style="clear: both; text-align: center;"><a href="http://2.bp.blogspot.com/-RxzKVIrxbl8/VUzC4HGhMNI/AAAAAAAAH9c/7FmxbMuB4kw/s1600/Screenshot_2015-05-08_16-04-53.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="60" src="http://2.bp.blogspot.com/-RxzKVIrxbl8/VUzC4HGhMNI/AAAAAAAAH9c/7FmxbMuB4kw/s320/Screenshot_2015-05-08_16-04-53.png" width="320" /></a></div>Then he integrates the Fokker-Planck equation over all A, to obtain<br /><div class="separator" style="clear: both; text-align: center;"><a href="http://4.bp.blogspot.com/-mDsC-Zd6FMA/VUzC4L9dF_I/AAAAAAAAH9g/cAvRNw1VTkg/s1600/Screenshot_2015-05-08_16-05-36.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="http://4.bp.blogspot.com/-mDsC-Zd6FMA/VUzC4L9dF_I/AAAAAAAAH9g/cAvRNw1VTkg/s1600/Screenshot_2015-05-08_16-05-36.png" /></a></div>On the backward Komolgorov equation, he applies a Lamperti transform like change of variable:<br /><div class="separator" style="clear: both; text-align: center;"><a href="http://2.bp.blogspot.com/-QiYromr1SDU/VUzE5q_OfjI/AAAAAAAAH94/nfbr14tEnj0/s1600/Screenshot_2015-05-08_16-10-33.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="http://2.bp.blogspot.com/-QiYromr1SDU/VUzE5q_OfjI/AAAAAAAAH94/nfbr14tEnj0/s1600/Screenshot_2015-05-08_16-10-33.png" /></a></div>And then makes another change of variable so that the PDE has the same initial conditions for all moments: <br /><div class="separator" style="clear: both; text-align: center;"><a href="http://2.bp.blogspot.com/-1C8j58UD1lA/VUzE5hZJ95I/AAAAAAAAH-I/541DaJBFbAU/s1600/Screenshot_2015-05-08_16-12-34.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="http://2.bp.blogspot.com/-1C8j58UD1lA/VUzE5hZJ95I/AAAAAAAAH-I/541DaJBFbAU/s1600/Screenshot_2015-05-08_16-12-34.png" /></a></div> This leads to<br /><div class="separator" style="clear: both; text-align: center;"><a href="http://3.bp.blogspot.com/--wj3nYsi9g8/VUzE5glQ0VI/AAAAAAAAH98/NlVSXIc2NDI/s1600/Screenshot_2015-05-08_16-13-32.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="107" src="http://3.bp.blogspot.com/--wj3nYsi9g8/VUzE5glQ0VI/AAAAAAAAH98/NlVSXIc2NDI/s320/Screenshot_2015-05-08_16-13-32.png" width="320" /></a></div>It turns out that there is a magical symmetry for k=0 and k=2.<br /><div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-OStv8D0OSw0/VUzGJ6HKqDI/AAAAAAAAH-Y/jk85U57eHnY/s1600/Screenshot_2015-05-08_16-19-45.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="32" src="http://1.bp.blogspot.com/-OStv8D0OSw0/VUzGJ6HKqDI/AAAAAAAAH-Y/jk85U57eHnY/s320/Screenshot_2015-05-08_16-19-45.png" width="320" /></a></div><div class="separator" style="clear: both; text-align: center;"><a href="http://3.bp.blogspot.com/-67k-xTphi7Q/VUzGJwbwSWI/AAAAAAAAH-U/Faz235QNpKs/s1600/Screenshot_2015-05-08_16-20-09.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="32" src="http://3.bp.blogspot.com/-67k-xTphi7Q/VUzGJwbwSWI/AAAAAAAAH-U/Faz235QNpKs/s320/Screenshot_2015-05-08_16-20-09.png" width="320" /></a></div>Note that in the second equation, the second derivative applies to the whole.<br />Because of this, he can express \(Q^{(2)}\) in terms of \(Q^{(0)}\):<br /><div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-k-3DQd9MuIQ/VUzHl5TSuAI/AAAAAAAAH-s/CjTFtlX4upw/s1600/Screenshot_2015-05-08_16-25-36.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="29" src="http://1.bp.blogspot.com/-k-3DQd9MuIQ/VUzHl5TSuAI/AAAAAAAAH-s/CjTFtlX4upw/s320/Screenshot_2015-05-08_16-25-36.png" width="320" /></a></div>And he plugs that back to the integrated Fokker-Planck equation to obtain the arbitrage free SABR PDE:<br /><div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-aEwTOxulsWI/VUzHl3qfZAI/AAAAAAAAH-o/qQZPb8Vvz7o/s1600/Screenshot_2015-05-08_16-25-48.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="29" src="http://1.bp.blogspot.com/-aEwTOxulsWI/VUzHl3qfZAI/AAAAAAAAH-o/qQZPb8Vvz7o/s320/Screenshot_2015-05-08_16-25-48.png" width="320" /></a></div><br />There is a simple more common explanation in the world of local stochastic volatility for what's going on. For example, in the particle method paper from Guyon-Labordère, we have the following expression for the true local volatility.<br /><div class="separator" style="clear: both; text-align: center;"><a href="http://3.bp.blogspot.com/-uV792mNz4Xo/Ul6baEs2ktI/AAAAAAAAG1U/8Iv1i23oXok/s1600/Screenshot%2Bfrom%2B2013-10-16%2B15%3A51%3A46.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="127" src="http://3.bp.blogspot.com/-uV792mNz4Xo/Ul6baEs2ktI/AAAAAAAAG1U/8Iv1i23oXok/s320/Screenshot%2Bfrom%2B2013-10-16%2B15%3A51%3A46.png" width="320" /></a></div><br />In the first equation, the numerator is simply \(Q^{(2)}\) and the denominator \(Q^{(0)}\). Of course, the integrated Fokker-Planck equation can be rewritten as:<br /><br />$$Q^{(0)}_T = \frac{1}{2}\epsilon^2 \left[C^2(F) \frac{Q^{(2)}}{Q^{(0)}} Q^{(0)}\right]_{FF}$$<br /><br />Karlsmark uses that approach directly in his thesis, using the expansions of Doust for \(Q^{(k)}\). Looking a Doust expansions, the fraction reduces straightforwardly to the same expression as Hagan, and the symmetry in the equations appears a bit less coincidental. <br /><br /><br /><img src="http://feeds.feedburner.com/~r/chasethedevil/~4/GmgtnJAUJ0Q" height="1" width="1" alt=""/>0http://chasethedevil.blogspot.com/2015/05/decoding-hagans-arbitrage-free-sabr-pde.htmlMatching Hagan PDE SABR with the one-step Andreasen-Huge SABRhttp://feedproxy.google.com/~r/chasethedevil/~3/5uxiTEU3oh4/matching-hagan-pde-sabr-with-one-step_30.htmlquantnoreply@blogger.com (Fabien)Thu, 17 Dec 2015 07:21:39 PSTtag:blogger.com,1999:blog-14415229.post-1951464426210018523I looked nearly two years ago already at <a href="http://chasethedevil.blogspot.co.uk/2013/05/sabr-with-new-hagan-pde-approach.html">the arbitrage free SABR of Andreasen-Huge in comparison to the arbitrage free PDE of Hagan</a> and showed how close the <a href="http://chasethedevil.blogspot.fr/2013/12/arbitrage-free-sabr-another-view-on.html">ideas were</a>: Andreasen-Huge relies on the normal Dupire forward PDE using a slightly simpler local vol (no time dependent exponential term) while Hagan works directly on the Fokker-Planck PDE (you can think of it as the Dupire Forward PDE for the density) and uses an expansion of same order as the original SABR formula (which leads to an additional exponential term in the local volatility).<br /><br />One clever idea from Andreasen-Huge is the use of a single step. It turns out that their idea is not completely new. Daniel Duffy sent me some old papers from Shishkin around fitted schemes (<a href="http://oai.cwi.nl/oai/asset/10209/10209A.pdf">here is one</a>). This is very much the same thing, except Shishkin concern is about a good handling of discontinuity in the initial condition, and therefore makes the association step function => cumulative density to fit the diffusion parameter. Andreasen-Huge work directly with the call prices as this is what they solve.<br /><br />One drawback of Andreasen-Huge one step method is the inability to match the standard SABR smile: the parameters don't have exactly the same meaning. It turns out that by just shifting proportionally the local volatility by a constant factor, Andreasen Huge matches Hagan PDE vols quite closely.<br /><br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://2.bp.blogspot.com/-BNmlkwTWu7Y/VUJGeLtiKII/AAAAAAAAH88/5ZeorrS4WRg/s1600/Screenshot_2015-04-30_17-12-43.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="http://2.bp.blogspot.com/-BNmlkwTWu7Y/VUJGeLtiKII/AAAAAAAAH88/5ZeorrS4WRg/s1600/Screenshot_2015-04-30_17-12-43.png" height="395" width="640" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Implied Black volatilities</td></tr></tbody></table><br /><br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="http://3.bp.blogspot.com/-rIUYFj3y7ro/VUJF6HOwyHI/AAAAAAAAH80/aOJrEeI1p9U/s1600/Screenshot_2015-04-30_17-09-53.png" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" src="http://3.bp.blogspot.com/-rIUYFj3y7ro/VUJF6HOwyHI/AAAAAAAAH80/aOJrEeI1p9U/s1600/Screenshot_2015-04-30_17-09-53.png" height="395" width="640" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">Probability density</td></tr></tbody></table><br />While this is interesting in itself, it's still not so simple to backup this factor without solving for it (and then the method looses appeal).<img src="http://feeds.feedburner.com/~r/chasethedevil/~4/5uxiTEU3oh4" height="1" width="1" alt=""/>0http://chasethedevil.blogspot.com/2015/04/matching-hagan-pde-sabr-with-one-step_30.html