Gerber Statistic Implementation in Rcpp and OpenMP

Posted on May 25, 2016 by Next Level Analytics in R bloggers | 0 Comments

[This article was first published on Next Level Analytics' Blog, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

By Yi, Nico, and Chris

Summary

Recently new research has appeared on using a co-movement measure to
construct the covariance matrix as part of the Modern Portfolio Theory (MPT)
style portfolio construction. Below is the abstract of the
Gerber, Markowith and Pujara (2015) paper whose methodology is also currently patent pending:

Markowitz’s mean-variance MPT has remained the cornerstone of portfolio selection methods after decades of research and debate. There is an extensive literature on MPT implementation, especially on estimation errors and expected return assumptions. However, covariance matrix estimation, an essential input, continues to be frequently based on historical correlations. There has been a recent new study that proposes replacing historical correlations with a robust co-movement measure called the Gerber Statistic.

In the research paper, it is stated that MPT using the Gerber Statistic outperformed portfolios using historical correlation as measured by ex-post returns under realistic investment constraints, including transaction costs and a broad range of investor types, for an investment universe of global stock indices, bonds and commodities for the period January 1994 to December 2013.

This post is to illustrate an implementation of the Gerber statistic. The focus is to compare the speed of computation for three different implementations with increasing performance

  • R
  • Rcpp
  • Rcpp with OpenMP for parallization

Implementation in R

<span class="n">gerber.correlation</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="k">function</span><span class="p">(</span><span class="n">hist.returns</span><span class="p">,</span><span class="w"> </span><span class="n">lookback</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="n">nrow</span><span class="p">(</span><span class="n">hist.returns</span><span class="p">),</span><span class="w"> </span><span class="n">threshold</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="m">0.5</span><span class="p">)</span><span class="w"> </span><span class="p">{</span><span class="w">
    </span><span class="n">n</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">ncol</span><span class="p">(</span><span class="n">hist.returns</span><span class="p">)</span><span class="w">
    </span><span class="n">nperiods</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">nrow</span><span class="p">(</span><span class="n">hist.returns</span><span class="p">)</span><span class="w">
  
    </span><span class="k">if</span><span class="w"> </span><span class="p">(</span><span class="n">lookback</span><span class="w"> </span><span class="o">></span><span class="w"> </span><span class="n">nperiods</span><span class="p">)</span><span class="w"> </span><span class="n">lookback</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">nperiods</span><span class="w">
    </span><span class="n">index</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="p">(</span><span class="n">nperiods</span><span class="w"> </span><span class="o">-</span><span class="w"> </span><span class="n">lookback</span><span class="w"> </span><span class="o">+</span><span class="w"> </span><span class="m">1</span><span class="p">)</span><span class="w"> </span><span class="o">:</span><span class="w"> </span><span class="n">nperiods</span><span class="w">
  
    </span><span class="n">standard.deviation</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">apply</span><span class="p">(</span><span class="n">hist.returns</span><span class="p">[</span><span class="n">index</span><span class="p">,,</span><span class="n">drop</span><span class="o">=</span><span class="nb">F</span><span class="p">],</span><span class="w"> </span><span class="m">2</span><span class="p">,</span><span class="w"> </span><span class="n">sd</span><span class="p">,</span><span class="w"> </span><span class="n">na.rm</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="nb">T</span><span class="p">)</span><span class="w">
    </span><span class="n">threshold</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">threshold</span><span class="w"> </span><span class="o">*</span><span class="w"> </span><span class="n">standard.deviation</span><span class="w">
  
    </span><span class="n">correlation</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">matrix</span><span class="p">(</span><span class="m">1</span><span class="p">,</span><span class="w"> </span><span class="n">n</span><span class="p">,</span><span class="w"> </span><span class="n">n</span><span class="p">)</span><span class="w">
    </span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="n">i</span><span class="w"> </span><span class="k">in</span><span class="w"> </span><span class="m">1</span><span class="o">:</span><span class="p">(</span><span class="n">n</span><span class="m">-1</span><span class="p">))</span><span class="w">
        </span><span class="k">for</span><span class="w"> </span><span class="p">(</span><span class="n">j</span><span class="w"> </span><span class="k">in</span><span class="w"> </span><span class="m">2</span><span class="o">:</span><span class="n">n</span><span class="p">)</span><span class="w"> </span><span class="p">{</span><span class="w">
            </span><span class="n">pos</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="nf">sum</span><span class="p">((</span><span class="n">hist.returns</span><span class="p">[,</span><span class="n">i</span><span class="p">]</span><span class="w"> </span><span class="o">>=</span><span class="w"> </span><span class="n">threshold</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="w"> </span><span class="o">&</span><span class="w"> </span><span class="n">hist.returns</span><span class="p">[,</span><span class="n">j</span><span class="p">]</span><span class="w"> </span><span class="o">>=</span><span class="w"> </span><span class="n">threshold</span><span class="p">[</span><span class="n">j</span><span class="p">])</span><span class="w"> </span><span class="o">|</span><span class="w">
                       </span><span class="p">(</span><span class="n">hist.returns</span><span class="p">[,</span><span class="n">i</span><span class="p">]</span><span class="w"> </span><span class="o"><=</span><span class="w"> </span><span class="o">-</span><span class="n">threshold</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="w"> </span><span class="o">&</span><span class="w"> </span><span class="n">hist.returns</span><span class="p">[,</span><span class="n">j</span><span class="p">]</span><span class="w"> </span><span class="o"><=</span><span class="w"> </span><span class="o">-</span><span class="n">threshold</span><span class="p">[</span><span class="n">j</span><span class="p">]),</span><span class="w"> </span><span class="n">na.rm</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="nb">T</span><span class="p">)</span><span class="w">
      
            </span><span class="n">neg</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="nf">sum</span><span class="p">((</span><span class="n">hist.returns</span><span class="p">[,</span><span class="n">i</span><span class="p">]</span><span class="w"> </span><span class="o">>=</span><span class="w"> </span><span class="n">threshold</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="w"> </span><span class="o">&</span><span class="w"> </span><span class="n">hist.returns</span><span class="p">[,</span><span class="n">j</span><span class="p">]</span><span class="w"> </span><span class="o"><=</span><span class="w"> </span><span class="o">-</span><span class="n">threshold</span><span class="p">[</span><span class="n">j</span><span class="p">])</span><span class="w"> </span><span class="o">|</span><span class="w">
                       </span><span class="p">(</span><span class="n">hist.returns</span><span class="p">[,</span><span class="n">i</span><span class="p">]</span><span class="w"> </span><span class="o"><=</span><span class="w"> </span><span class="o">-</span><span class="n">threshold</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="w"> </span><span class="o">&</span><span class="w"> </span><span class="n">hist.returns</span><span class="p">[,</span><span class="n">j</span><span class="p">]</span><span class="w"> </span><span class="o">>=</span><span class="w"> </span><span class="n">threshold</span><span class="p">[</span><span class="n">j</span><span class="p">]),</span><span class="w"> </span><span class="n">na.rm</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="nb">T</span><span class="p">)</span><span class="w">
      
            </span><span class="n">correlation</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="n">correlation</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">i</span><span class="p">]</span><span class="w"> </span><span class="o"><-</span><span class="w"> </span><span class="p">(</span><span class="n">pos</span><span class="w"> </span><span class="o">-</span><span class="w"> </span><span class="n">neg</span><span class="p">)</span><span class="w"> </span><span class="o">/</span><span class="w"> </span><span class="p">(</span><span class="n">pos</span><span class="w"> </span><span class="o">+</span><span class="w"> </span><span class="n">neg</span><span class="p">)</span><span class="w">
        </span><span class="p">}</span><span class="w">
    </span><span class="n">correlation</span><span class="w">
</span><span class="p">}</span><span class="w">
</span>

Implementation in Rcpp

<span class="c1">// [[Rcpp::depends(RcppArmadillo)]]
</span><span class="cp">#include <RcppArmadillo.h>
</span>
<span class="k">using</span> <span class="k">namespace</span> <span class="n">arma</span><span class="p">;</span> 
<span class="k">using</span> <span class="k">namespace</span> <span class="n">Rcpp</span><span class="p">;</span>

<span class="c1">// [[Rcpp::export]]
</span><span class="kt">double</span> <span class="nf">SD_RCPP</span><span class="p">(</span><span class="n">arma</span><span class="o">::</span><span class="n">vec</span> <span class="n">DATA_VEC</span><span class="p">)</span> <span class="p">{</span>
  
    <span class="kt">double</span> <span class="n">SD</span><span class="p">;</span>
    <span class="kt">double</span> <span class="n">MEAN</span><span class="p">;</span>
  
    <span class="n">MEAN</span> <span class="o">=</span> <span class="n">std</span><span class="o">::</span><span class="n">accumulate</span><span class="p">(</span><span class="n">DATA_VEC</span><span class="p">.</span><span class="n">begin</span><span class="p">(),</span> <span class="n">DATA_VEC</span><span class="p">.</span><span class="n">end</span><span class="p">(),</span><span class="mf">0.0</span><span class="p">)</span><span class="o">/</span><span class="n">DATA_VEC</span><span class="p">.</span><span class="n">size</span><span class="p">();</span>
    <span class="n">DATA_VEC</span> <span class="o">=</span> <span class="n">pow</span><span class="p">(</span><span class="n">DATA_VEC</span> <span class="o">-</span> <span class="n">MEAN</span><span class="p">,</span><span class="mi">2</span><span class="p">);</span>
    <span class="n">SD</span> <span class="o">=</span> <span class="n">pow</span><span class="p">(</span><span class="n">std</span><span class="o">::</span><span class="n">accumulate</span><span class="p">(</span><span class="n">DATA_VEC</span><span class="p">.</span><span class="n">begin</span><span class="p">(),</span> <span class="n">DATA_VEC</span><span class="p">.</span><span class="n">end</span><span class="p">(),</span><span class="mf">0.0</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="n">DATA_VEC</span><span class="p">.</span><span class="n">size</span><span class="p">()</span><span class="o">-</span><span class="mi">1</span><span class="p">),</span> <span class="mf">0.5</span><span class="p">);</span>
  
    <span class="k">return</span> <span class="n">SD</span><span class="p">;</span>
<span class="p">}</span>

<span class="c1">// [[Rcpp::export]]
</span><span class="n">mat</span> <span class="nf">GERBER_CORRELATION</span><span class="p">(</span><span class="n">arma</span><span class="o">::</span><span class="n">mat</span> <span class="n">HIST_RETURN</span><span class="p">,</span>
                       <span class="kt">int</span> <span class="n">LOOKBACK</span> <span class="o">=</span> <span class="mi">0</span><span class="p">,</span>
                       <span class="kt">double</span> <span class="n">THRESHOLD</span> <span class="o">=</span> <span class="mf">0.5</span><span class="p">,</span>
                       <span class="kt">bool</span> <span class="n">LOOKBACK_ORDER</span><span class="o">=</span> <span class="nb">false</span><span class="p">)</span> <span class="p">{</span>
  
    <span class="n">arma</span><span class="o">::</span><span class="n">mat</span> <span class="n">CORRELATION_MAT</span><span class="p">;</span>
    <span class="n">arma</span><span class="o">::</span><span class="n">mat</span> <span class="n">HIST_RETURN_SD</span><span class="p">;</span>
    <span class="n">arma</span><span class="o">::</span><span class="n">vec</span> <span class="n">SD_VEC</span><span class="p">;</span>
    <span class="n">arma</span><span class="o">::</span><span class="n">vec</span> <span class="n">THRESHOLD_VEC</span><span class="p">;</span>
  
    <span class="n">arma</span><span class="o">::</span><span class="n">uvec</span> <span class="n">TEMP_VEC_POS</span><span class="p">;</span>
    <span class="n">arma</span><span class="o">::</span><span class="n">uvec</span> <span class="n">TEMP_VEC_NEG</span><span class="p">;</span>
  
    <span class="n">arma</span><span class="o">::</span><span class="n">uvec</span> <span class="n">TEMP_UVEC_1</span><span class="p">;</span>
    <span class="n">arma</span><span class="o">::</span><span class="n">uvec</span> <span class="n">TEMP_UVEC_2</span><span class="p">;</span>
    <span class="n">arma</span><span class="o">::</span><span class="n">uvec</span> <span class="n">TEMP_UVEC_3</span><span class="p">;</span>
    <span class="n">arma</span><span class="o">::</span><span class="n">uvec</span> <span class="n">TEMP_UVEC_4</span><span class="p">;</span>
  
    <span class="kt">double</span> <span class="n">NCOL</span><span class="p">;</span>
    <span class="kt">double</span> <span class="n">NPERIODS</span><span class="p">;</span>
  
    <span class="kt">int</span> <span class="n">i</span><span class="p">;</span>
    <span class="kt">int</span> <span class="n">j</span><span class="p">;</span>
  
    <span class="kt">double</span> <span class="n">POS</span><span class="p">;</span>
    <span class="kt">double</span> <span class="n">NEG</span><span class="p">;</span>
  
    <span class="n">NCOL</span><span class="o">=</span> <span class="n">HIST_RETURN</span><span class="p">.</span><span class="n">n_cols</span><span class="p">;</span>
    <span class="n">NPERIODS</span> <span class="o">=</span> <span class="n">HIST_RETURN</span><span class="p">.</span><span class="n">n_rows</span><span class="p">;</span>
  
    <span class="k">if</span> <span class="p">(</span><span class="n">LOOKBACK</span> <span class="o">></span> <span class="mi">0</span><span class="p">)</span> <span class="p">{</span>
        <span class="k">if</span> <span class="p">(</span><span class="n">LOOKBACK_ORDER</span><span class="p">)</span> <span class="p">{</span>
            <span class="k">if</span> <span class="p">(</span><span class="n">LOOKBACK</span> <span class="o"><</span> <span class="n">NPERIODS</span> <span class="o">-</span> <span class="mi">1</span><span class="p">)</span> <span class="p">{</span> <span class="n">NPERIODS</span> <span class="o">=</span> <span class="n">LOOKBACK</span><span class="p">;</span> <span class="p">}</span>
            <span class="n">HIST_RETURN_SD</span> <span class="o">=</span> <span class="n">HIST_RETURN</span><span class="p">.</span><span class="n">rows</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="n">NPERIODS</span><span class="o">-</span><span class="mi">1</span><span class="p">);</span>
        <span class="p">}</span> <span class="k">else</span> <span class="p">{</span>
            <span class="k">if</span><span class="p">(</span><span class="n">LOOKBACK</span> <span class="o">></span> <span class="n">NPERIODS</span> <span class="o">-</span> <span class="mi">1</span><span class="p">){</span><span class="n">LOOKBACK</span> <span class="o">=</span> <span class="n">NPERIODS</span><span class="p">;}</span>
            <span class="n">HIST_RETURN_SD</span> <span class="o">=</span> <span class="n">HIST_RETURN</span><span class="p">.</span><span class="n">rows</span><span class="p">(</span><span class="n">NPERIODS</span> <span class="o">-</span> <span class="n">LOOKBACK</span><span class="p">,</span> <span class="n">NPERIODS</span><span class="o">-</span><span class="mi">1</span><span class="p">);</span>
        <span class="p">}</span>
    <span class="p">}</span>

    <span class="n">SD_VEC</span><span class="p">.</span><span class="n">set_size</span><span class="p">(</span><span class="n">NCOL</span><span class="p">);</span>
    <span class="k">for</span> <span class="p">(</span><span class="n">i</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span> <span class="n">i</span> <span class="o"><</span> <span class="n">NCOL</span><span class="p">;</span> <span class="n">i</span><span class="o">++</span><span class="p">)</span> <span class="p">{</span>
        <span class="n">SD_VEC</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">SD_RCPP</span><span class="p">(</span><span class="n">HIST_RETURN_SD</span><span class="p">.</span><span class="n">col</span><span class="p">(</span><span class="n">i</span><span class="p">));</span>
    <span class="p">}</span>

    <span class="n">THRESHOLD_VEC</span><span class="p">.</span><span class="n">set_size</span><span class="p">(</span><span class="n">NCOL</span><span class="p">);</span>
    <span class="n">THRESHOLD_VEC</span> <span class="o">=</span> <span class="n">THRESHOLD</span> <span class="o">*</span> <span class="n">SD_VEC</span><span class="p">;</span>

    <span class="n">CORRELATION_MAT</span><span class="p">.</span><span class="n">set_size</span><span class="p">(</span><span class="n">NCOL</span><span class="p">,</span> <span class="n">NCOL</span><span class="p">);</span>
    <span class="n">CORRELATION_MAT</span><span class="p">.</span><span class="n">fill</span><span class="p">(</span><span class="mi">1</span><span class="p">);</span>
  
    <span class="k">for</span> <span class="p">(</span><span class="n">i</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span> <span class="n">i</span> <span class="o"><</span> <span class="n">NCOL</span><span class="o">-</span><span class="mi">1</span><span class="p">;</span> <span class="n">i</span><span class="o">++</span><span class="p">)</span> <span class="p">{</span>
        <span class="k">for</span> <span class="p">(</span><span class="n">j</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span> <span class="n">j</span> <span class="o"><</span> <span class="n">i</span><span class="p">;</span> <span class="n">j</span><span class="o">++</span><span class="p">)</span> <span class="p">{</span>
            <span class="k">if</span> <span class="p">(</span><span class="n">i</span> <span class="o">==</span> <span class="n">j</span><span class="p">)</span> <span class="p">{</span>
                <span class="n">CORRELATION_MAT</span><span class="p">.</span><span class="n">at</span><span class="p">(</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">)</span> <span class="o">=</span> <span class="mi">1</span><span class="p">;</span>
            <span class="p">}</span> <span class="k">else</span> <span class="p">{</span>
                <span class="n">TEMP_UVEC_1</span> <span class="o">=</span> <span class="n">find</span><span class="p">(</span><span class="n">HIST_RETURN</span><span class="p">.</span><span class="n">col</span><span class="p">(</span><span class="n">i</span><span class="p">)</span><span class="o">></span><span class="n">THRESHOLD_VEC</span><span class="p">[</span><span class="n">i</span><span class="p">]);</span>
                <span class="n">TEMP_UVEC_2</span> <span class="o">=</span> <span class="n">find</span><span class="p">(</span><span class="n">HIST_RETURN</span><span class="p">.</span><span class="n">col</span><span class="p">(</span><span class="n">j</span><span class="p">)</span><span class="o">></span><span class="n">THRESHOLD_VEC</span><span class="p">[</span><span class="n">j</span><span class="p">]);</span>
                <span class="n">TEMP_UVEC_3</span> <span class="o">=</span> <span class="n">find</span><span class="p">(</span><span class="n">HIST_RETURN</span><span class="p">.</span><span class="n">col</span><span class="p">(</span><span class="n">i</span><span class="p">)</span><span class="o"><-</span><span class="n">THRESHOLD_VEC</span><span class="p">[</span><span class="n">i</span><span class="p">]);</span>
                <span class="n">TEMP_UVEC_4</span> <span class="o">=</span> <span class="n">find</span><span class="p">(</span><span class="n">HIST_RETURN</span><span class="p">.</span><span class="n">col</span><span class="p">(</span><span class="n">j</span><span class="p">)</span><span class="o"><-</span><span class="n">THRESHOLD_VEC</span><span class="p">[</span><span class="n">j</span><span class="p">]);</span>
        
                <span class="n">TEMP_VEC_POS</span> <span class="o">=</span> <span class="n">find_unique</span><span class="p">(</span><span class="n">join_cols</span><span class="p">(</span><span class="n">TEMP_UVEC_1</span><span class="p">,</span> <span class="n">TEMP_UVEC_2</span><span class="p">));</span>
                <span class="n">TEMP_VEC_NEG</span> <span class="o">=</span> <span class="n">find_unique</span><span class="p">(</span><span class="n">join_cols</span><span class="p">(</span><span class="n">TEMP_UVEC_3</span><span class="p">,</span> <span class="n">TEMP_UVEC_4</span><span class="p">));</span>
        
                <span class="n">POS</span> <span class="o">=</span> <span class="p">(</span><span class="n">TEMP_UVEC_1</span><span class="p">.</span><span class="n">size</span><span class="p">()</span> <span class="o">+</span> <span class="n">TEMP_UVEC_2</span><span class="p">.</span><span class="n">size</span><span class="p">()</span> <span class="o">-</span> <span class="n">TEMP_VEC_POS</span><span class="p">.</span><span class="n">size</span><span class="p">())</span> <span class="o">+</span> <span class="p">(</span><span class="n">TEMP_UVEC_3</span><span class="p">.</span><span class="n">size</span><span class="p">()</span> <span class="o">+</span> <span class="n">TEMP_UVEC_4</span><span class="p">.</span><span class="n">size</span><span class="p">()</span> <span class="o">-</span> <span class="n">TEMP_VEC_NEG</span><span class="p">.</span><span class="n">size</span><span class="p">());</span>
        
                <span class="n">TEMP_UVEC_1</span> <span class="o">=</span> <span class="n">find</span><span class="p">(</span><span class="n">HIST_RETURN</span><span class="p">.</span><span class="n">col</span><span class="p">(</span><span class="n">i</span><span class="p">)</span><span class="o">></span><span class="n">THRESHOLD_VEC</span><span class="p">[</span><span class="n">i</span><span class="p">]);</span>
                <span class="n">TEMP_UVEC_2</span> <span class="o">=</span> <span class="n">find</span><span class="p">(</span><span class="n">HIST_RETURN</span><span class="p">.</span><span class="n">col</span><span class="p">(</span><span class="n">j</span><span class="p">)</span><span class="o"><-</span><span class="n">THRESHOLD_VEC</span><span class="p">[</span><span class="n">j</span><span class="p">]);</span>
                <span class="n">TEMP_UVEC_3</span> <span class="o">=</span> <span class="n">find</span><span class="p">(</span><span class="n">HIST_RETURN</span><span class="p">.</span><span class="n">col</span><span class="p">(</span><span class="n">i</span><span class="p">)</span><span class="o"><-</span><span class="n">THRESHOLD_VEC</span><span class="p">[</span><span class="n">i</span><span class="p">]);</span>
                <span class="n">TEMP_UVEC_4</span> <span class="o">=</span> <span class="n">find</span><span class="p">(</span><span class="n">HIST_RETURN</span><span class="p">.</span><span class="n">col</span><span class="p">(</span><span class="n">j</span><span class="p">)</span><span class="o">></span><span class="n">THRESHOLD_VEC</span><span class="p">[</span><span class="n">j</span><span class="p">]);</span>
        
                <span class="n">TEMP_VEC_POS</span> <span class="o">=</span> <span class="n">find_unique</span><span class="p">(</span><span class="n">join_cols</span><span class="p">(</span><span class="n">TEMP_UVEC_1</span><span class="p">,</span> <span class="n">TEMP_UVEC_2</span><span class="p">));</span>
                <span class="n">TEMP_VEC_NEG</span> <span class="o">=</span> <span class="n">find_unique</span><span class="p">(</span><span class="n">join_cols</span><span class="p">(</span><span class="n">TEMP_UVEC_3</span><span class="p">,</span> <span class="n">TEMP_UVEC_4</span><span class="p">));</span>
        
                <span class="n">NEG</span> <span class="o">=</span> <span class="p">(</span><span class="n">TEMP_UVEC_1</span><span class="p">.</span><span class="n">size</span><span class="p">()</span> <span class="o">+</span> <span class="n">TEMP_UVEC_2</span><span class="p">.</span><span class="n">size</span><span class="p">()</span> <span class="o">-</span> <span class="n">TEMP_VEC_POS</span><span class="p">.</span><span class="n">size</span><span class="p">())</span> <span class="o">+</span> <span class="p">(</span><span class="n">TEMP_UVEC_3</span><span class="p">.</span><span class="n">size</span><span class="p">()</span> <span class="o">+</span> <span class="n">TEMP_UVEC_4</span><span class="p">.</span><span class="n">size</span><span class="p">()</span> <span class="o">-</span> <span class="n">TEMP_VEC_NEG</span><span class="p">.</span><span class="n">size</span><span class="p">());</span>
        
                <span class="n">CORRELATION_MAT</span><span class="p">.</span><span class="n">at</span><span class="p">(</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">)</span> <span class="o">=</span> <span class="p">(</span><span class="n">POS</span> <span class="o">-</span> <span class="n">NEG</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="n">POS</span> <span class="o">+</span> <span class="n">NEG</span><span class="p">);</span>
            <span class="p">}</span>
            <span class="n">CORRELATION_MAT</span><span class="p">.</span><span class="n">at</span><span class="p">(</span><span class="n">j</span><span class="p">,</span><span class="n">i</span><span class="p">)</span> <span class="o">=</span> <span class="n">CORRELATION_MAT</span><span class="p">.</span><span class="n">at</span><span class="p">(</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">);</span>
        <span class="p">}</span>
    <span class="p">}</span>
  
    <span class="k">return</span> <span class="n">CORRELATION_MAT</span><span class="p">;</span>
<span class="p">}</span>

Implementation in Rcpp with OpenMP

<span class="cp">#include <Rcpp.h>
#include <omp.h>
</span>
<span class="c1">// [[Rcpp::plugins(openmp)]]
</span>
<span class="k">using</span> <span class="k">namespace</span> <span class="n">Rcpp</span><span class="p">;</span>

<span class="c1">// [[Rcpp::export]]
</span><span class="n">NumericMatrix</span> <span class="nf">GERBER_CORRELATION_PARALLEL_OMP</span><span class="p">(</span><span class="n">NumericMatrix</span> <span class="n">HIST_RETURN_RAW</span><span class="p">,</span>
                                              <span class="kt">int</span> <span class="n">LOOKBACK</span> <span class="o">=</span> <span class="mi">0</span><span class="p">,</span>
                                              <span class="kt">double</span> <span class="n">THRESHOLD</span> <span class="o">=</span> <span class="mf">0.5</span><span class="p">,</span>
                                              <span class="kt">bool</span> <span class="n">LOOKBACK_ORDER</span> <span class="o">=</span> <span class="nb">false</span><span class="p">)</span> <span class="p">{</span>
    <span class="kt">double</span> <span class="n">NPERIODS</span><span class="p">;</span>
  
    <span class="kt">int</span> <span class="n">NROW_BEGIN</span><span class="p">;</span>
    <span class="kt">int</span> <span class="n">NROW_END</span><span class="p">;</span>
  
    <span class="kt">int</span> <span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="n">k</span><span class="p">;</span>
    <span class="kt">double</span> <span class="n">POS</span><span class="p">,</span> <span class="n">NEG</span><span class="p">;</span>
  
    <span class="n">NPERIODS</span> <span class="o">=</span> <span class="n">HIST_RETURN_RAW</span><span class="p">.</span><span class="n">nrow</span><span class="p">();</span>
    <span class="n">NROW_BEGIN</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span>
    <span class="n">NROW_END</span> <span class="o">=</span> <span class="n">NPERIODS</span><span class="p">;</span>
  
    <span class="n">NumericMatrix</span> <span class="n">HIST_RETURN</span><span class="p">;</span>
    <span class="k">if</span> <span class="p">(</span><span class="n">LOOKBACK</span> <span class="o">></span> <span class="mi">0</span><span class="p">)</span> <span class="p">{</span>
        <span class="k">if</span><span class="p">(</span><span class="n">LOOKBACK_ORDER</span><span class="p">)</span> <span class="p">{</span>
            <span class="k">if</span><span class="p">(</span><span class="n">LOOKBACK</span> <span class="o"><</span> <span class="n">NPERIODS</span><span class="p">)</span> <span class="p">{</span> <span class="n">NPERIODS</span> <span class="o">=</span> <span class="n">LOOKBACK</span><span class="p">;</span> <span class="p">}</span>
            <span class="n">NROW_BEGIN</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span>
            <span class="n">NROW_END</span> <span class="o">=</span> <span class="n">NPERIODS</span> <span class="o">-</span> <span class="mi">1</span><span class="p">;</span>
        <span class="p">}</span> <span class="k">else</span> <span class="p">{</span>
            <span class="k">if</span> <span class="p">(</span><span class="n">LOOKBACK</span> <span class="o">></span> <span class="n">NPERIODS</span><span class="p">)</span> <span class="p">{</span> <span class="n">LOOKBACK</span> <span class="o">=</span> <span class="n">NPERIODS</span><span class="p">;</span> <span class="p">}</span>
            <span class="n">NROW_BEGIN</span> <span class="o">=</span> <span class="n">NPERIODS</span> <span class="o">-</span> <span class="n">LOOKBACK</span><span class="p">;</span>
            <span class="n">NROW_END</span> <span class="o">=</span> <span class="n">NPERIODS</span> <span class="o">-</span> <span class="mi">1</span><span class="p">;</span>
        <span class="p">}</span>
        <span class="n">HIST_RETURN</span> <span class="o">=</span> <span class="n">HIST_RETURN_RAW</span><span class="p">(</span><span class="n">Range</span><span class="p">(</span><span class="n">NROW_BEGIN</span><span class="p">,</span><span class="n">NROW_END</span><span class="p">),</span><span class="n">Range</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="n">HIST_RETURN_RAW</span><span class="p">.</span><span class="n">ncol</span><span class="p">()</span><span class="o">-</span><span class="mi">1</span><span class="p">));</span>
    <span class="p">}</span> <span class="k">else</span> <span class="p">{</span>
        <span class="n">HIST_RETURN</span> <span class="o">=</span> <span class="n">clone</span><span class="p">(</span><span class="n">HIST_RETURN_RAW</span><span class="p">);</span>
    <span class="p">}</span>
  
    <span class="c1">//calculate standard deviation matrix
</span>    <span class="n">NumericMatrix</span> <span class="n">SD_MAT</span><span class="p">(</span><span class="n">HIST_RETURN</span><span class="p">.</span><span class="n">ncol</span><span class="p">(),</span><span class="mi">1</span><span class="p">);</span>
    <span class="n">std</span><span class="o">::</span><span class="n">fill</span><span class="p">(</span><span class="n">SD_MAT</span><span class="p">.</span><span class="n">begin</span><span class="p">(),</span> <span class="n">SD_MAT</span><span class="p">.</span><span class="n">end</span><span class="p">(),</span> <span class="mf">0.0</span><span class="p">);</span>
    <span class="k">for</span><span class="p">(</span><span class="n">j</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span> <span class="n">j</span> <span class="o"><</span> <span class="n">HIST_RETURN</span><span class="p">.</span><span class="n">ncol</span><span class="p">();</span> <span class="n">j</span><span class="o">++</span><span class="p">)</span> <span class="p">{</span>
        <span class="k">for</span><span class="p">(</span><span class="n">i</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span> <span class="n">i</span> <span class="o"><</span> <span class="n">HIST_RETURN</span><span class="p">.</span><span class="n">nrow</span><span class="p">();</span> <span class="n">i</span><span class="o">++</span><span class="p">)</span> <span class="p">{</span>
            <span class="n">SD_MAT</span><span class="p">(</span><span class="n">j</span><span class="p">,</span> <span class="mi">0</span><span class="p">)</span> <span class="o">=</span> <span class="n">SD_MAT</span><span class="p">(</span><span class="n">j</span><span class="p">,</span><span class="mi">0</span><span class="p">)</span> <span class="o">+</span> <span class="n">HIST_RETURN</span><span class="p">(</span><span class="n">i</span><span class="p">,</span> <span class="n">j</span><span class="p">);</span>
        <span class="p">}</span>
        <span class="n">SD_MAT</span><span class="p">(</span><span class="n">j</span><span class="p">,</span> <span class="mi">0</span><span class="p">)</span> <span class="o">=</span> <span class="n">SD_MAT</span><span class="p">(</span><span class="n">j</span><span class="p">,</span> <span class="mi">0</span><span class="p">)</span><span class="o">/</span><span class="n">HIST_RETURN</span><span class="p">.</span><span class="n">nrow</span><span class="p">();</span>
        <span class="k">for</span><span class="p">(</span><span class="n">i</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span> <span class="n">i</span> <span class="o"><</span> <span class="n">HIST_RETURN</span><span class="p">.</span><span class="n">nrow</span><span class="p">();</span> <span class="n">i</span><span class="o">++</span><span class="p">)</span> <span class="p">{</span>
            <span class="n">HIST_RETURN</span><span class="p">(</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">)</span> <span class="o">=</span> <span class="n">std</span><span class="o">::</span><span class="n">pow</span><span class="p">((</span><span class="n">HIST_RETURN</span><span class="p">(</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">)</span> <span class="o">-</span> <span class="n">SD_MAT</span><span class="p">(</span><span class="n">j</span><span class="p">,</span> <span class="mi">0</span><span class="p">)),</span> <span class="mi">2</span><span class="p">);</span>
        <span class="p">}</span>
        <span class="n">SD_MAT</span><span class="p">(</span><span class="n">j</span><span class="p">,</span> <span class="mi">0</span><span class="p">)</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span>
        <span class="k">for</span><span class="p">(</span><span class="n">i</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span> <span class="n">i</span> <span class="o"><</span> <span class="n">HIST_RETURN</span><span class="p">.</span><span class="n">nrow</span><span class="p">();</span> <span class="n">i</span><span class="o">++</span><span class="p">)</span> <span class="p">{</span>
            <span class="n">SD_MAT</span><span class="p">(</span><span class="n">j</span><span class="p">,</span> <span class="mi">0</span><span class="p">)</span> <span class="o">=</span> <span class="n">SD_MAT</span><span class="p">(</span><span class="n">j</span><span class="p">,</span><span class="mi">0</span><span class="p">)</span> <span class="o">+</span> <span class="n">HIST_RETURN</span><span class="p">(</span><span class="n">i</span><span class="p">,</span> <span class="n">j</span><span class="p">);</span>
        <span class="p">}</span>
        <span class="n">SD_MAT</span><span class="p">(</span><span class="n">j</span><span class="p">,</span> <span class="mi">0</span><span class="p">)</span> <span class="o">=</span> <span class="n">std</span><span class="o">::</span><span class="n">pow</span><span class="p">(</span><span class="n">SD_MAT</span><span class="p">(</span><span class="n">j</span><span class="p">,</span> <span class="mi">0</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="n">HIST_RETURN</span><span class="p">.</span><span class="n">nrow</span><span class="p">()</span><span class="o">-</span><span class="mi">1</span><span class="p">),</span> <span class="mf">0.5</span><span class="p">)</span> <span class="o">*</span> <span class="n">THRESHOLD</span><span class="p">;</span>
    <span class="p">}</span>
  
    <span class="c1">//calculate correlation matrix
</span>    <span class="n">NumericMatrix</span> <span class="n">CORRELATION_MAT</span><span class="p">(</span><span class="n">HIST_RETURN_RAW</span><span class="p">.</span><span class="n">ncol</span><span class="p">(),</span> <span class="n">HIST_RETURN_RAW</span><span class="p">.</span><span class="n">ncol</span><span class="p">());</span>
    <span class="n">std</span><span class="o">::</span><span class="n">fill</span><span class="p">(</span><span class="n">CORRELATION_MAT</span><span class="p">.</span><span class="n">begin</span><span class="p">(),</span> <span class="n">CORRELATION_MAT</span><span class="p">.</span><span class="n">end</span><span class="p">(),</span> <span class="mf">1.0</span><span class="p">);</span>
  
    <span class="cp">#pragma omp parallel for private(i, j, POS, NEG)
</span>    <span class="k">for</span> <span class="p">(</span><span class="n">i</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span> <span class="n">i</span> <span class="o"><</span> <span class="n">CORRELATION_MAT</span><span class="p">.</span><span class="n">nrow</span><span class="p">();</span> <span class="n">i</span><span class="o">++</span><span class="p">)</span> <span class="p">{</span>
        <span class="k">for</span><span class="p">(</span><span class="n">j</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span> <span class="n">j</span> <span class="o"><</span> <span class="n">i</span><span class="p">;</span> <span class="n">j</span><span class="o">++</span><span class="p">)</span> <span class="p">{</span>
            <span class="k">if</span><span class="p">(</span><span class="n">i</span> <span class="o">==</span> <span class="n">j</span><span class="p">)</span> <span class="p">{</span>
                <span class="n">CORRELATION_MAT</span><span class="p">(</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">)</span> <span class="o">=</span> <span class="mi">1</span><span class="p">;</span>
            <span class="p">}</span> <span class="k">else</span> <span class="p">{</span>
                <span class="n">POS</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span>
                <span class="n">NEG</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span>
        
                <span class="k">for</span> <span class="p">(</span><span class="n">k</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span> <span class="n">k</span> <span class="o"><</span> <span class="n">HIST_RETURN_RAW</span><span class="p">.</span><span class="n">nrow</span><span class="p">();</span> <span class="n">k</span><span class="o">++</span><span class="p">)</span> <span class="p">{</span>
                    <span class="k">if</span><span class="p">(((</span><span class="n">HIST_RETURN_RAW</span><span class="p">(</span><span class="n">k</span><span class="p">,</span><span class="n">i</span><span class="p">)</span> <span class="o">></span> <span class="n">SD_MAT</span><span class="p">(</span><span class="n">i</span><span class="p">,</span><span class="mi">0</span><span class="p">))</span> <span class="o">&</span> <span class="p">(</span><span class="n">HIST_RETURN_RAW</span><span class="p">(</span><span class="n">k</span><span class="p">,</span><span class="n">j</span><span class="p">)</span> <span class="o">></span> <span class="n">SD_MAT</span><span class="p">(</span><span class="n">j</span><span class="p">,</span><span class="mi">0</span><span class="p">)))</span> <span class="o">|</span>
                       <span class="p">((</span><span class="n">HIST_RETURN_RAW</span><span class="p">(</span><span class="n">k</span><span class="p">,</span><span class="n">i</span><span class="p">)</span> <span class="o"><</span> <span class="o">-</span><span class="mf">1.0</span><span class="o">*</span><span class="n">SD_MAT</span><span class="p">(</span><span class="n">i</span><span class="p">,</span><span class="mi">0</span><span class="p">))</span> <span class="o">&</span> <span class="p">(</span><span class="n">HIST_RETURN_RAW</span><span class="p">(</span><span class="n">k</span><span class="p">,</span><span class="n">j</span><span class="p">)</span> <span class="o"><</span> <span class="o">-</span><span class="mf">1.0</span><span class="o">*</span><span class="n">SD_MAT</span><span class="p">(</span><span class="n">j</span><span class="p">,</span><span class="mi">0</span><span class="p">))))</span> <span class="p">{</span>
                        <span class="n">POS</span><span class="o">++</span><span class="p">;</span>
                    <span class="p">}</span>
                    <span class="k">if</span><span class="p">(((</span><span class="n">HIST_RETURN_RAW</span><span class="p">(</span><span class="n">k</span><span class="p">,</span><span class="n">i</span><span class="p">)</span> <span class="o">></span> <span class="n">SD_MAT</span><span class="p">(</span><span class="n">i</span><span class="p">,</span><span class="mi">0</span><span class="p">))</span> <span class="o">&</span> <span class="p">(</span><span class="n">HIST_RETURN_RAW</span><span class="p">(</span><span class="n">k</span><span class="p">,</span><span class="n">j</span><span class="p">)</span> <span class="o"><</span> <span class="o">-</span><span class="mf">1.0</span><span class="o">*</span><span class="n">SD_MAT</span><span class="p">(</span><span class="n">j</span><span class="p">,</span><span class="mi">0</span><span class="p">)))</span> <span class="o">|</span>
                       <span class="p">((</span><span class="n">HIST_RETURN_RAW</span><span class="p">(</span><span class="n">k</span><span class="p">,</span><span class="n">i</span><span class="p">)</span> <span class="o"><</span> <span class="o">-</span><span class="mf">1.0</span><span class="o">*</span><span class="n">SD_MAT</span><span class="p">(</span><span class="n">i</span><span class="p">,</span><span class="mi">0</span><span class="p">))</span> <span class="o">&</span> <span class="p">(</span><span class="n">HIST_RETURN_RAW</span><span class="p">(</span><span class="n">k</span><span class="p">,</span><span class="n">j</span><span class="p">)</span> <span class="o">></span> <span class="n">SD_MAT</span><span class="p">(</span><span class="n">j</span><span class="p">,</span><span class="mi">0</span><span class="p">))))</span> <span class="p">{</span>
                        <span class="n">NEG</span><span class="o">++</span><span class="p">;</span>
                    <span class="p">}</span>
                <span class="p">}</span>
                <span class="n">CORRELATION_MAT</span><span class="p">(</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">)</span> <span class="o">=</span> <span class="p">(</span><span class="n">POS</span> <span class="o">-</span> <span class="n">NEG</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="n">POS</span> <span class="o">+</span> <span class="n">NEG</span><span class="p">);</span>
            <span class="p">}</span>
            <span class="n">CORRELATION_MAT</span><span class="p">(</span><span class="n">j</span><span class="p">,</span><span class="n">i</span><span class="p">)</span> <span class="o">=</span> <span class="n">CORRELATION_MAT</span><span class="p">(</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">);</span>
        <span class="p">}</span>
    <span class="p">}</span>
    
    <span class="k">return</span> <span class="n">CORRELATION_MAT</span><span class="p">;</span>
<span class="p">}</span>

Speed Comparison

Finally, let’s compare the speed gain result. The test data is based on a
return matrix of 30 securities with 2500 data points. It can be seen that the
OpenMP version of the calculation is clearly faster than the serial version
which itself is much faster than the R version.

Implementation replications elapsed relative
R Version 500 0.707 1.000
Rcpp Version 500 78.807 111.467
Rcpp + OpenMP Version 500 227.384 321.618

To leave a comment for the author, please follow the link and comment on their blog: Next Level Analytics' Blog.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Copyright © 2020 | MH Corporate basic by MH Themes

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)