<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/using-rcpparmadillo-to-price-european-put-options/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://gallery.rcpp.org//articles/black-scholes-three-ways/"> Rcpp Gallery</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<h3>Introduction</h3>
<p>In the quest for ever faster code, one generally begins exploring ways to integrate C++<br>
with R using <a rel="nofollow" href="http://www.rcpp.org/" target="_blank">Rcpp</a>. This post provides an example of multiple<br>
implementations of a European Put Option pricer. The implementations are done in pure R,<br>
pure <a rel="nofollow" href="http://www.rcpp.org/" target="_blank">Rcpp</a> using some <a rel="nofollow" href="http://www.rcpp.org/" target="_blank">Rcpp</a> sugar functions,<br>
and then in <a rel="nofollow" href="http://www.rcpp.org/" target="_blank">Rcpp</a> using<br>
<a rel="nofollow" href="http://dirk.eddelbuettel.com/code/rcpp.armadillo.html" target="_blank">RcppArmadillo</a>, which exposes the<br>
incredibly powerful linear algebra library, <a rel="nofollow" href="http://arma.sourceforge.net/" target="_blank">Armadillo</a>.</p>
<p>Under the <a rel="nofollow" href="https://en.wikipedia.org/wiki/Black%E2%80%93Scholes_model" target="_blank">Black-Scholes model</a> The value of a European put option has the closed form solution:</p>
<p></p>
<p>where</p>
<p></p>
<p>and</p>
<p></p>
<p>Armed with the formulas, we can create the pricer using just R.</p>
<figure>
<pre><code data-lang="r"><span>put_option_pricer</span><span> </span><span>&lt;-</span><span> </span><span>function</span><span>(</span><span>s</span><span>,</span><span> </span><span>k</span><span>,</span><span> </span><span>r</span><span>,</span><span> </span><span>y</span><span>,</span><span> </span><span>t</span><span>,</span><span> </span><span>sigma</span><span>)</span><span> </span><span>{</span><span>
</span><span>d</span><span>1</span><span> </span><span>&lt;-</span><span> </span><span>(</span><span>log</span><span>(</span><span>s</span><span> </span><span>/</span><span> </span><span>k</span><span>)</span><span> </span><span>+</span><span> </span><span>(</span><span>r</span><span> </span><span>-</span><span> </span><span>y</span><span> </span><span>+</span><span> </span><span>sigma</span><span>^</span><span>2</span><span> </span><span>/</span><span> </span><span>2</span><span>)</span><span> </span><span>*</span><span> </span><span>t</span><span>)</span><span> </span><span>/</span><span> </span><span>(</span><span>sigma</span><span> </span><span>*</span><span> </span><span>sqrt</span><span>(</span><span>t</span><span>))</span><span>
</span><span>d</span><span>2</span><span> </span><span>&lt;-</span><span> </span><span>d</span><span>1</span><span> </span><span>-</span><span> </span><span>sigma</span><span> </span><span>*</span><span> </span><span>sqrt</span><span>(</span><span>t</span><span>)</span><span>
</span><span>V</span><span> </span><span>&lt;-</span><span> </span><span>pnorm</span><span>(</span><span>-</span><span>d</span><span>2</span><span>)</span><span> </span><span>*</span><span> </span><span>k</span><span> </span><span>*</span><span> </span><span>exp</span><span>(</span><span>-</span><span>r</span><span> </span><span>*</span><span> </span><span>t</span><span>)</span><span> </span><span>-</span><span> </span><span>s</span><span> </span><span>*</span><span> </span><span>exp</span><span>(</span><span>-</span><span>y</span><span> </span><span>*</span><span> </span><span>t</span><span>)</span><span> </span><span>*</span><span> </span><span>pnorm</span><span>(</span><span>-</span><span>d</span><span>1</span><span>)</span><span>
</span><span>V</span><span>
</span><span>}</span><span>
</span><span># Valuation with 1 stock price</span><span>
</span><span>put_option_pricer</span><span>(</span><span>s</span><span> </span><span>=</span><span> </span><span>55</span><span>,</span><span> </span><span>60</span><span>,</span><span> </span><span>.01</span><span>,</span><span> </span><span>.02</span><span>,</span><span> </span><span>1</span><span>,</span><span> </span><span>.05</span><span>)</span></code></pre>
</figure>
<pre>
[1] 5.52021
</pre>
<figure>
<pre><code data-lang="r"><span># Valuation across multiple prices</span><span>
</span><span>put_option_pricer</span><span>(</span><span>s</span><span> </span><span>=</span><span> </span><span>55</span><span>:</span><span>60</span><span>,</span><span> </span><span>60</span><span>,</span><span> </span><span>.01</span><span>,</span><span> </span><span>.02</span><span>,</span><span> </span><span>1</span><span>,</span><span> </span><span>.05</span><span>)</span></code></pre>
</figure>
<pre>
[1] 5.52021 4.58142 3.68485 2.85517 2.11883 1.49793
</pre>
<p>Let’s see what we can do with <a rel="nofollow" href="http://www.rcpp.org/" target="_blank">Rcpp</a>. Besides explicitely stating the<br>
types of the variables, not much has to change. We can even use the sugar function,<br>
<code>Rcpp::pnorm()</code>, to keep the syntax as close to R as possible. Note how we are being<br>
explicit about the symbols we import from the <code>Rcpp</code> namespace: the basic vector type, and<br>
well the (vectorized) ‘Rcpp Sugar’ calls <code>log()</code> and <code>pnorm()</code> calls. Similarly, we use<br>
<code>sqrt()</code> and <code>exp()</code> for the calls on an atomic <code>double</code> variables from the C++ Standard<br>
Library. With these declarations the code itself is essentially identical to the R code<br>
(apart of course from requiring both static types and trailing semicolons).</p>
<figure>
<pre><code data-lang="cpp"><span>#include &lt;Rcpp.h>
</span>                                        
<span>using</span> <span>Rcpp</span><span>::</span><span>NumericVector</span><span>;</span>
<span>using</span> <span>Rcpp</span><span>::</span><span>log</span><span>;</span>
<span>using</span> <span>Rcpp</span><span>::</span><span>pnorm</span><span>;</span>
<span>using</span> <span>std</span><span>::</span><span>sqrt</span><span>;</span>
<span>using</span> <span>std</span><span>::</span><span>log</span><span>;</span>

<span>// [[Rcpp::export]]
</span><span>NumericVector</span> <span>put_option_pricer_rcpp</span><span>(</span><span>NumericVector</span> <span>s</span><span>,</span> <span>double</span> <span>k</span><span>,</span> <span>double</span> <span>r</span><span>,</span> <span>double</span> <span>y</span><span>,</span> <span>double</span> <span>t</span><span>,</span> <span>double</span> <span>sigma</span><span>)</span> <span>{</span>

    <span>NumericVector</span> <span>d1</span> <span>=</span> <span>(</span><span>log</span><span>(</span><span>s</span> <span>/</span> <span>k</span><span>)</span> <span>+</span> <span>(</span><span>r</span> <span>-</span> <span>y</span> <span>+</span> <span>sigma</span> <span>*</span> <span>sigma</span> <span>/</span> <span>2.0</span><span>)</span> <span>*</span> <span>t</span><span>)</span> <span>/</span> <span>(</span><span>sigma</span> <span>*</span> <span>sqrt</span><span>(</span><span>t</span><span>));</span>
    <span>NumericVector</span> <span>d2</span> <span>=</span> <span>d1</span> <span>-</span> <span>sigma</span> <span>*</span> <span>sqrt</span><span>(</span><span>t</span><span>);</span>
    
    <span>NumericVector</span> <span>V</span> <span>=</span> <span>pnorm</span><span>(</span><span>-</span><span>d2</span><span>)</span> <span>*</span> <span>k</span> <span>*</span> <span>exp</span><span>(</span><span>-</span><span>r</span> <span>*</span> <span>t</span><span>)</span> <span>-</span> <span>s</span> <span>*</span> <span>exp</span><span>(</span><span>-</span><span>y</span> <span>*</span> <span>t</span><span>)</span> <span>*</span> <span>pnorm</span><span>(</span><span>-</span><span>d1</span><span>);</span>
    <span>return</span> <span>V</span><span>;</span>
<span>}</span></code></pre>
</figure>
<p>We can call this from R as well:</p>
<figure>
<pre><code data-lang="r"><span># Valuation with 1 stock price</span><span>
</span><span>put_option_pricer_rcpp</span><span>(</span><span>s</span><span> </span><span>=</span><span> </span><span>55</span><span>,</span><span> </span><span>60</span><span>,</span><span> </span><span>.01</span><span>,</span><span> </span><span>.02</span><span>,</span><span> </span><span>1</span><span>,</span><span> </span><span>.05</span><span>)</span></code></pre>
</figure>
<pre>
[1] 5.52021
</pre>
<figure>
<pre><code data-lang="r"><span># Valuation across multiple prices</span><span>
</span><span>put_option_pricer_rcpp</span><span>(</span><span>s</span><span> </span><span>=</span><span> </span><span>55</span><span>:</span><span>60</span><span>,</span><span> </span><span>60</span><span>,</span><span> </span><span>.01</span><span>,</span><span> </span><span>.02</span><span>,</span><span> </span><span>1</span><span>,</span><span> </span><span>.05</span><span>)</span></code></pre>
</figure>
<pre>
[1] 5.52021 4.58142 3.68485 2.85517 2.11883 1.49793
</pre>
<p>Finally, let’s look at<br>
<a rel="nofollow" href="http://dirk.eddelbuettel.com/code/rcpp.armadillo.html" target="_blank">RcppArmadillo</a>. Armadillo has a<br>
number of object types, including <code>mat</code>, <code>colvec</code>, and <code>rowvec</code>. Here, we just use<br>
<code>colvec</code> to represent a column vector of prices. By default in Armadillo, <code>*</code> represents<br>
matrix multiplication, and <code>%</code> is used for element wise multiplication. We need to make<br>
this change to element wise multiplication in 1 place, but otherwise the changes are just<br>
switching out the types and the sugar functions for Armadillo specific functions.</p>
<p>Note that the <code>arma::normcdf()</code> function is in the upcoming release of<br>
<a rel="nofollow" href="http://dirk.eddelbuettel.com/code/rcpp.armadillo.html" target="_blank">RcppArmadillo</a>, which is<br>
<code>0.8.400.0.0</code> at the time of writing and still in CRAN’s incoming. It also requires the<br>
<code>C++11</code> plugin.</p>
<figure>
<pre><code data-lang="cpp"><span>#include &lt;RcppArmadillo.h>
</span>
<span>// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::plugins(cpp11)]]
</span>
<span>using</span> <span>arma</span><span>::</span><span>colvec</span><span>;</span>
<span>using</span> <span>arma</span><span>::</span><span>log</span><span>;</span>
<span>using</span> <span>arma</span><span>::</span><span>normcdf</span><span>;</span>
<span>using</span> <span>std</span><span>::</span><span>sqrt</span><span>;</span>
<span>using</span> <span>std</span><span>::</span><span>log</span><span>;</span>


<span>// [[Rcpp::export]]
</span><span>colvec</span> <span>put_option_pricer_arma</span><span>(</span><span>colvec</span> <span>s</span><span>,</span> <span>double</span> <span>k</span><span>,</span> <span>double</span> <span>r</span><span>,</span> <span>double</span> <span>y</span><span>,</span> <span>double</span> <span>t</span><span>,</span> <span>double</span> <span>sigma</span><span>)</span> <span>{</span>
  
    <span>colvec</span> <span>d1</span> <span>=</span> <span>(</span><span>log</span><span>(</span><span>s</span> <span>/</span> <span>k</span><span>)</span> <span>+</span> <span>(</span><span>r</span> <span>-</span> <span>y</span> <span>+</span> <span>sigma</span> <span>*</span> <span>sigma</span> <span>/</span> <span>2.0</span><span>)</span> <span>*</span> <span>t</span><span>)</span> <span>/</span> <span>(</span><span>sigma</span> <span>*</span> <span>sqrt</span><span>(</span><span>t</span><span>));</span>
    <span>colvec</span> <span>d2</span> <span>=</span> <span>d1</span> <span>-</span> <span>sigma</span> <span>*</span> <span>sqrt</span><span>(</span><span>t</span><span>);</span>
    
    <span>// Notice the use of % to represent element wise multiplication
</span>    <span>colvec</span> <span>V</span> <span>=</span> <span>normcdf</span><span>(</span><span>-</span><span>d2</span><span>)</span> <span>*</span> <span>k</span> <span>*</span> <span>exp</span><span>(</span><span>-</span><span>r</span> <span>*</span> <span>t</span><span>)</span> <span>-</span> <span>s</span> <span>*</span> <span>exp</span><span>(</span><span>-</span><span>y</span> <span>*</span> <span>t</span><span>)</span> <span>%</span> <span>normcdf</span><span>(</span><span>-</span><span>d1</span><span>);</span> 

    <span>return</span> <span>V</span><span>;</span>
<span>}</span></code></pre>
</figure>
<p>Use from R:</p>
<figure>
<pre><code data-lang="r"><span># Valuation with 1 stock price</span><span>
</span><span>put_option_pricer_arma</span><span>(</span><span>s</span><span> </span><span>=</span><span> </span><span>55</span><span>,</span><span> </span><span>60</span><span>,</span><span> </span><span>.01</span><span>,</span><span> </span><span>.02</span><span>,</span><span> </span><span>1</span><span>,</span><span> </span><span>.05</span><span>)</span></code></pre>
</figure>
<pre>
        [,1]
[1,] 5.52021
</pre>
<figure>
<pre><code data-lang="r"><span># Valuation across multiple prices</span><span>
</span><span>put_option_pricer_arma</span><span>(</span><span>s</span><span> </span><span>=</span><span> </span><span>55</span><span>:</span><span>60</span><span>,</span><span> </span><span>60</span><span>,</span><span> </span><span>.01</span><span>,</span><span> </span><span>.02</span><span>,</span><span> </span><span>1</span><span>,</span><span> </span><span>.05</span><span>)</span></code></pre>
</figure>
<pre>
        [,1]
[1,] 5.52021
[2,] 4.58142
[3,] 3.68485
[4,] 2.85517
[5,] 2.11883
[6,] 1.49793
</pre>
<p>Finally, we can run a speed test to see which comes out on top.</p>
<figure>
<pre><code data-lang="r"><span>s</span><span> </span><span>&lt;-</span><span> </span><span>matrix</span><span>(</span><span>seq</span><span>(</span><span>0</span><span>,</span><span> </span><span>100</span><span>,</span><span> </span><span>by</span><span> </span><span>=</span><span> </span><span>.0001</span><span>),</span><span> </span><span>ncol</span><span> </span><span>=</span><span> </span><span>1</span><span>)</span><span>
</span><span>rbenchmark</span><span>::</span><span>benchmark</span><span>(</span><span>R</span><span> </span><span>=</span><span> </span><span>put_option_pricer</span><span>(</span><span>s</span><span>,</span><span> </span><span>60</span><span>,</span><span> </span><span>.01</span><span>,</span><span> </span><span>.02</span><span>,</span><span> </span><span>1</span><span>,</span><span> </span><span>.05</span><span>),</span><span>
</span><span>Arma</span><span> </span><span>=</span><span> </span><span>put_option_pricer_arma</span><span>(</span><span>s</span><span>,</span><span> </span><span>60</span><span>,</span><span> </span><span>.01</span><span>,</span><span> </span><span>.02</span><span>,</span><span> </span><span>1</span><span>,</span><span> </span><span>.05</span><span>),</span><span>
</span><span>Rcpp</span><span> </span><span>=</span><span> </span><span>put_option_pricer_rcpp</span><span>(</span><span>s</span><span>,</span><span> </span><span>60</span><span>,</span><span> </span><span>.01</span><span>,</span><span> </span><span>.02</span><span>,</span><span> </span><span>1</span><span>,</span><span> </span><span>.05</span><span>),</span><span> </span><span>order</span><span> </span><span>=</span><span> </span><span>&quot;relative&quot;</span><span>,</span><span> </span><span>replications</span><span> </span><span>=</span><span> </span><span>100</span><span>)[,</span><span>1</span><span>:</span><span>4</span><span>]</span></code></pre>
</figure>
<pre>
  test replications elapsed relative
2 Arma          100   6.409    1.000
3 Rcpp          100   7.917    1.235
1    R          100   9.091    1.418
</pre>
<p>Interestingly, <a rel="nofollow" href="http://arma.sf.net/" target="_blank">Armadillo</a> comes out on top here on this (multi-core)<br>
machine (as Armadillo uses OpenMP where available in newer versions). But the difference<br>
is slender, and there is certainly variation in repeated runs. And the nicest thing about<br>
all of this is that it shows off the “embarassment of riches” that we have in the R and<br>
C++ ecosystem for multiple ways of solving the same problem.</p>
<p><img width="1" alt="" src="http://feeds.feedburner.com/~r/rcpp-gallery/~4/E-r3V2pUjg8" height="1"></p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://gallery.rcpp.org//articles/black-scholes-three-ways/"> Rcpp Gallery</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/KufAubJLdYY" height="1"></a></div>