I want to sample from the empirical distribution of returns. To do so, I do not want to make the preliminary assumption of which distribution the returns follow, rather I would like to sample from the empirical unknown distribution of returns. The sampled value will help me in a Montecarlo simulation.
I thought of applying Metropolis-Hastings or Gibbs sampler having the empirical distribution as target distribution.
Do you know if any statistical technique already exists to accommodate this?
Financial returns exhibit well known time-dependancy in its higher conditional moments. For starters, just about no matter how you produce a time series of conditional volatility, it will be exhibit clustering patterns and almost always a high degree of persistence. So, regardless of what you want to do here, avoid sampling from the unconditional distribution of returns.
One option for doing what you want would be to take a clue from Rosenberg and Engel (2002). We're going to start by specifying a simple GARCH model: \begin{align} r_{t+1} - r_{ft} = \mu + \lambda h_{t+1} + \sqrt{h_{t+1}} \epsilon_t \\ h_{t+1} = w_h + \sum_{i=1}^2 \left( b_i h_{t+1-i} + a_i(\epsilon_t - \gamma_i \sqrt{h_{t+1-i}})^2 \right). \end{align} Assume $\epsilon_{t+1} \sim N(0,1)$ and estimate this by maximum likelihood. That's a GARCH(2,2) version of the Heston and Nandi (2000) option pricing model. I chose this one because it can be shown that it is equivalent to a component GARCH model à la Engel and Lee (1993) and consequently allows higher persistance than a GARCH(1,1) AND because it allows for asymmetry to be built into through parameters $(\gamma_i)_{i=1}^2$.
The idea here is that you just need to estimate the model because you need some kind of filter for estimates of conditional variance and shocks $(\hat{h}_t, \hat{\epsilon_t})_{t=1}^T$. You can see the above as a very convenient way to get a filter. Even if the shocks aren't normal, you can still have a good filter and your estimated shocks will likely have fat tails and skewness built into them. The only real downside here is that you then have the problem of working from a conditional density, meaning you need one such density per horizon -- and the simplest way to get it would be to use a parametric model for the conditional volatility such as the one you just used as a filter, assuming this proxies for all required changes over time.
Another course of action would involve using generative models from the machine learning literature. You could train a Generative Adverserial Network to spit out appropriate returns -- if you want sequences of length $\tau$, train it to mimmick those; if you just want the final returns at date $t+\tau$, train it to spit out this number.
