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
 A: 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

A: Ok, I think have might have an "analytic" answer.
I'm a bit unsure about it, so if anyone fancies casting an eye over it that'd be much appreciated.

*

*Fix $q_x$ and $q_y$.


*Define indicator functions $I$ and $J$. $I$ is 1 when $x$ exceeds $q_x$ and zero otherwise. Similarly for $J$, but for $y$.


*Define $\hat{p}_x=\frac{\sum I}{N}$, $\hat{q}_y=\frac{\sum J}{N}$.
These are unbiased estimators for the probability of exceeding $q_x$ and $q_y$. Write those probabilities as $p_x$ and $p_y$.


*Show that $V(\hat{p}_x)=\frac{p_x(1-p_x)}{N}$, similarly for $V(\hat{p}_y)$


*Show that $cov(\hat{p}_x,\hat{p}_y)=\frac{cov(I,J)}{N}$


*Estimate $cov(I,J)$ using $\frac{\sum I J}{N}-\hat{p}_x\hat{p}_y$


*Assume that the CDF is locally linear for both $x$ and $y$. Write those linear models as $p_x=\alpha_x + \beta_x q_x$ and $p_y=\alpha_y + \beta_y q_y$


*Then $cov(\hat{p}_x,\hat{p}_y)=\beta_x \beta_y cov(\hat{q}_x,\hat{q}_y)$


*Rearrange and substitute in and that seems to give the answer


*I'm going to code this and Kozolovska's solution up and see if they give similar results.


*Not sure I'm really using the notation correctly.
