Fréchet Hoeffding bounds for symmetric random variables (Edited to clarify the question).
The Hakan & Demirtas (2012 doi: 10.1198/tast.2011.10090) approach to approximating Pearson correlation bounds uses the concept of the Fréchet-Hoeffding bounds by simulating a (very large) sample of the desired marginal distributions, sorting them so that either both variables are ordered from smallest to largest (positive correlation) or only one of them from smallest to largest and the other from largest to smallest (negative correlation). This ends up revealing which pairs of marginal distributions do not fully span the [-1, +1] correlation range.
Classic example: If X and Y are independent, standard normally distributed, take $A=e^{X}$ and $B=e^{Y}$ so that both A and Y are log-normally distributed. It's not hard to show that the correlation range between these log-normal distributed variables is $[-1/e, +1]$. The Hakan & Demirtas method would work as follows for the case of the lower bound:
n <- 10000000

z1 <- rnorm(n, 0, 1)
z2 <- rnorm(n, 0, 1)

x1 <- exp(z1)
x2 <- exp(z2)

cor(sort(x1, decreasing=T), sort(x2, decreasing=F), method=pearson)
[1] -0.3678432

My own simulations have shown me that the range $[-1, +1]$ is only restricted in the negative side of the range if the marginal distributions are skewed in in the same direction OR the positive side of the range if the variables have opposite skews. I have yet to find a case where the range $[-1,+1]$ is restricted if the variables are symmetric. So my question is:
Would it be the general case that, with symmetric marginals, for continuous, non-degenerate distributions, the correlation range fully spans $[-1,+1]$ when sorted and correlated? Or is there a counter-example where the marginals are symmetric but sorting them and correlating them yields a reduction of range?
 A: I'll try to address the multiple questions that are packed into this.


*

*Strictly speaking, the Frechet-Hoeffding bounds do not place bounds on the correlation.  The upper bound $M(u,v)$ represents the dependence structure when two random variables are perfectly co-monotonic, while the lower bound $W(u,v)$ represents the dependence structure when two random variables are perfectly counter-monotonic.  The dependence structure says nothing about the marginal distributions, they can be whatever you want them to.

*Correlation (or more generally, measures of association) attempt to reduce the full dependence structure, which is represented by the copula function, into one number.  The correlation coefficient is a measure of linear association, while measures of concordance (such as Kendall's $\tau$ or Spearman's $\rho$) are measures of concordance (but limited to monotonic dependence structures).  Non-monotonic measures of association include distance correlation (dCor) and the Randomized Dependence Coefficient (RDC).  

*You need to think of these tools to measure correlation as estimators of the "true" underlying correlation.  So, you may estimate Kendall's $\tau$ or Spearman's $\rho$ using the available estimators in routines such as https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.kendalltau.html or https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.spearmanr.html#scipy.stats.spearmanr.  Depending on the skew of the data, your estimator may have some bias associated with its measurement, which is probably what you are experiencing.  Just because the estimate is less than one, it doesn't necessarily mean that your underlying dependence structure is not the Frechet-Hoeffding bounds.

*You can easily create a bivariate distribution which has a correlation less than 1, where the marginal distributions are symmetric.  What about the multivariate Gaussian, with a correlation coefficient of 0.6 for example?  In this case, your dependence structure is captured by the Gaussian copula, and your marginal distributions are Gaussian as well.  Or, you can go with any copula such as a Gaussian or Archimedean, where the dependence parameter is less than perfect dependence, and go with any marginal distribution you like.  A simple simulation of this could be done in Matlab as follows (here, we have a linear dependence structure captured by the Gaussian copula, with T distributions as the marginal, with 4 and 2 degrees of freedom respectively:

U = copularnd('Gaussian', 0.6, 1000); X = [tinv(U(:,1),4),
  tinv(U(:,2),2)];

