Quick estimation for empirical permutation and correlation test in R
Using apply and replicate, you can get a slight speedup (below). Setting a seed, you can check that it gives the same result as your version. In my benchmark test, it's about 10% faster. For further speedup, you could consider using clusterApply and parLapply from the parallel package. 

# combining matrices into one to use apply
combinedMat &lt;- cbind(a[rep(1:nrow(a), each=3), ], b[rep(1:nrow(b), times=3), ])
combinedMat &lt;- cbind(combinedMat, obs = as.numeric(t(rcorr.observed)))
# using apply instead of loop
Pval_estim &lt;- t(matrix(apply(combinedMat, 1, function(combined) {
  obs &lt;- combined['obs']
  if (is.na(obs)) return(NA)
  combined &lt;- as.numeric(combined[names(combined) != 'obs'])
  # using replicate instead of loop
  diff.random &lt;- replicate(number_of_permutations,
  {
    shuffled = sample(combined)
    # you only use the correlation, so there is no need to use rcorr. 
    diff.random = cor(shuffled[1 : 6], 
                      shuffled[7 : length(combined)], 
                      method="spearman")
  })
  sum(abs(diff.random) &gt;= abs(obs)) / number_of_permutations
}), nrow = 3, ncol = 3))