How to estimate correlations between the standard errors on empirical quantile estimates from two correlated series? I have two very long series of samples from a pair of correlated random variables, $x$ and $y$.
(they are fat tailed distributions, roughly Pareto, if that helps).
(if you're curious...$x$ and $y$ are annual simulated hurricane losses: $x$ is without climate change, and $y$ is with climate change. They are correlated because of how they are simulated).
I calculate the 90% quantiles from each series (often known as the 10 year return period losses).  Let's call them $q_x$ and $q_y$.
I also calculate the standard errors around these quantile estimates, using the binomial method. Let's call them $\sigma_{qx}$ and $\sigma_{qy}$.
Then I calculate the difference in the 90% quantiles, $d=q_y-q_x$, and the fractional change $f=(q_y-q_x)/q_x$. These are the things I'm really interested in (the change in the quantile due to climate change). So far so good.
But I'd also like to estimate the standard error around this difference $d$ and the fractional change $f$. Let's call them $\sigma_d$ and $\sigma_f$. That's where I get stuck.
If $x$ and $y$ were independent, I think I'd be ok, because I think $q_x$ and $q_y$ would be independent. I think I could say $\sigma_d^2=\sigma_{qx}^2+\sigma_{qy}^2$.
But $x$ and $y$ are not independent, and so the errors around the quantile estimates are not independent. I think I can say $\sigma_d^2=\sigma_{qx}^2+\sigma_{qy}^2-2cov(q_x,q_y)$.
My question is: is there any way I could calculate, or at least estimate, the covariance term?
My only idea so far is that I'm able to calculate whatever tail correlation metrics I need from the original series, which seems like it may be relevant to me.
thanks

You can try and use bootstrap. The idea to create many examples of the quantiles and then take the difference between them. If you sample them together you do not break the dependency between them.
Here is some code example:
# parameters 
n <- 10000 
b <- 1000
# generate data 
x <- rnorm(n, 0, 1)
y <- x * 0.5 + rnorm(n, 0, sqrt(0.75))
mat <- cbind(x, y)

# bootstrap 
all_quants <- NULL 

for (i in 1:b) {
  sample_mat <- mat[sample(n, n, replace=TRUE), ]
  iter_quant <- apply(sample_mat, 2, quantile, 0.9)
  all_quants <- rbind(all_quants, iter_quant)
}

# estimate covariance 
cov(all_quants)

# estimate the difference variance directly 
sd(all_quants[ ,1] - all_quants[ ,2])

Results
> # estimate covariance 
> cov(all_quants)
             x            y
x 3.756014e-04 9.541534e-05
y 9.541534e-05 2.621710e-04
> 
> # estimate the difference variance directly 
> sd(all_quants[ ,1] - all_quants[ ,2])
[1] 0.021141

