How to generate correlated samples?

Question

The variation of the average price of a product over a longer period is generally lower than a shorter period. I am interested to capture both uncertainties as to the input of the stochastic programming problem. Let's say the average price of oil in a week has a mean of $\overline{\lambda}_{w}$ and std of $\sigma_{w}$. Thus, the N samples representing the weakly average price of oil can be shown by $\lambda_{wi}, \ wi=1...N$.

On the other hand, let's represent the daily variations of oil price regarding each above mentioned generated sample with M other samples, e.g. $\lambda_{wi}^{dj},\ j=1:M$. Specifically, for each $wi$, the uncertainty around the daily price, despite having the same mean as the weekly sample, i.e. $\overline{\lambda_{wi}}^{dj=1:M} = \lambda_{wi}$, has a larger standard deviation, $\text{i.e } \sigma_{d}$.

To better illustrate, lets say we have N=2 samples regarding the average oil price during one week with the mean of $\overline{\lambda}_{w} = 100 \\\$$ and $\sigma_{w} = 20 \\\$$. Also, for each weakly sample, we have M=3 samples discirbing the daily variation of the price with $\overline{\lambda_{wi}}^{dj=1:M} = \lambda_{wi}$ and $ \sigma_{d} = 30$.

$\text{for } \lambda_{w1} = 80, \ \lambda_{d1}^{w1} = 50,\lambda_{d1}^{w1} =80,\lambda_{d1}^{w1}=110 \\ \text{for } \lambda_{w2} = 120, \ \lambda_{d1}^{w2} = 90,\lambda_{d1}^{w2} = 120,\lambda_{d1}^{w2}=150$

Assume that the uncertainty in both time periods can be modeled as a Gaussian distribution function.

Q1) How can I generate N samples with a mean of $\overline{\lambda}_{w}$ and std of $\sigma_{w}$ as well as M other samples with mean of $\overline{\lambda_{wi}}^{dj=1:M} = \lambda_{wi}$ and $\sigma_{d} $

Q2) How can I generate the above samples while considering the std of the daily price variations as a function of weekly average sample: e.g. $\sigma_{di} = \lambda_{wi}/10$

If possible please give some hints for Matlab or Python implementation of such a sampling method.

Answer 1

The issue you are describing has to do with the necessity of accounting for both short- and long-term dynamics in a decision problem under uncertainty, or in general uncertainty at different levels of resolution. There are two issues here.

1. The practical implementation of a stochastic program lives on a scenario tree. So the first issue is how to arrange random data in a scenario tree accounting for different time resolutions. The authors of this work suggest a so called "multi-horizon" scenario tree, which allows you to model exactly the issue you are describing, that is, uncertainty at different levels of resolution. In practice, when you populate your scenario tree, you would sample daily prices conditioning on the realization of the weekly prices.

2. How to concretely sample from a multivariate Gaussian in order to populate you scenario tree (and your stochastic program) with data. In Python you can use numpy to sample from a multivariate Gaussian. It is explained here. For example:

Import the necessary stuff

from scipy.stats import multivariate_normal

Example data which you would have to populate based on your case

n_random_vars = 5
means = [10 for i in range(n_random_vars)] 
covs = [[0 for x in range(n_random_vars))] for y in range(n_random_vars)] 

Create a multivariate_normal object

mn = multivariate_normal(mean = means, cov = covs)

Draw random samples

mn.rvs()

Answer 2

How to generate a multivariate Gaussian? It must be answered somewhere on Cross Validated, but I cannot find it now, some comments at https://stats.stackexchange.com/questions/341805/are-mvrnorm-in-mass-r-package-and-rmvn-in-mgcv-r-package-equivalent/341808#341808.

Let $X \sim \mathcal{N}(\mu, \Sigma)$ and $\epsilon \sim \mathcal{N}(0,I)$. Then we can decompose the covariance matrix in different ways, say the Cholesky decomposition $\Sigma= C C^T$ or by the spectral theorem $\Sigma=U \Lambda U^T = (U\Lambda^{1/2}) (U\Lambda^{1/2})^T$ and then either $$ \mu + C\epsilon$$ or $$ \mu + U \Lambda^{1/2} \epsilon$$ does the trick.