10
$\begingroup$

As I mentioned in the question title, I want to generate specific uniformly distributed samples.

I need them to model a real world scenario. For my real data, I estimated a function, which approximates the auto-correlation function by a e-function. I also modeled a distribution model for my real data.

I can generate samples, which follow my distribution model with the inverse transform sampling.

If i think correct, instead of "standard" uniform distributed samples, which i feed into the inverse transform sampling, i need to adapt those uniform distributed samples, so that their auto-correlation function follows my estimated e-function.

Now the question is, if i am correct, how to adapt my uniform distributed samples, that they correspond to my requirements.

It would be great if, somebody could help.

EDIT 1:

Here is a plot of the auto-correlation function of my real data (different scenarios) and the modeled function (green). The ACF of my uniform distributed samples should follow the green (respectively blue) line.

Image 1

EDIT 2:

Thank you so much so far for the answers. To make my problem a bit clearer, i added a illustration of the problem. For a simulation, i need to generate samples, which follow my distribution model (with inverse transform sampling a easy task) AND have the same auto-correlation behavior like the real data.

My idea was, that i need to adapt the uniformly distributed samples, which i feed into the inverse transform sampling, so that their ACF follows my model.

The problem is, that i don't know if this idea is correct and that i don't know how to adapt the uniformly distributed samples.

Image 2

EDIT 3: Following the approach from jblood94, the solution is based on a reordering of the output samples after the inverse transform sampling. So, the underlying distribution stays the same AND the variables follow the desired acf.

I tried to translate the R code from jblood94 to matlab:

clear;
close all;
clc;


x = gamrnd(0.9, 1, 1e6, 1);
alpha = [1, sort(betarnd(1, 3, 1, 50), 'descend')];

tic
xx = acf_reorder(x, alpha);
toc

figure;
autocorr(xx, 'NumLags', length(alpha)-1);
hold on;
plot((1:length(alpha))-1, alpha);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function xx = acf_reorder(x, alpha)

tol = 1e-5;
maxIter = 10;
n = length(x);
xx = sort(x);
y = randn(n,1);
[w0, w, alpha1] = deal(alpha);
m = length(alpha);
m1 = 1:(m-1);
tol10 = tol/10;
iter = 0;
x = xx(tiedrank(filter(w, 1, y)));
SSE0 = inf;

acfInit = autocorr(x, 'NumLags', m-1)';
SSE = sum((acfInit(1:m)-alpha).^2);

while SSE > tol
    if(SSE < SSE0)
        SSE0 = SSE;
        w = w0;
        iter = iter + 1;
        if(iter == maxIter)
            break;
        end
    else
        break;
    end

    w1 = w0;
    a = 0;
    sse0 = inf;

    while (max(abs(alpha1 - a)) > tol10)
        a = [1, arrayfun(@(x) sum(head(w1, -x).*tail(w1, -x)), m1)/sum(w1.^2)];

        sse = sum((a - alpha1).^2);
        if (sse < sse0)
            sse0 = sse;
            w0 = w1;
        else
            % w0 failed to converge; try optim
            disp('w0 failed to converge; try optim');
            w1 = exp(fminsearch(@fun, log(w0), [], m1));

            a = [1, arrayfun(@(x) sum(head(w1, -x).*tail(w1, -x)), m1)/sum(w1.^2)];
            if (sum((a - alpha1).^2) < sse0)
                w0 = w1;
                break;
            end
        end

        w1 = (w1.*alpha1./a + w1)/2;
    end

    x = xx(tiedrank(filter(w0, 1, y)));
    acf = autocorr(x, 'NumLags', m-1);
    alpha1 = (alpha1.*alpha./acf(1:m) + alpha1)/2;

    SSE = sum((acf(1:m)-alpha).^2);
end

xx = xx(tiedrank(filter(w, 1, y)));

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function output = fun(ww, m1)

ww = exp(ww);
output = sum([1, arrayfun(@(x) sum(head(ww, -x).*tail(ww, -x)), m1)/sum(ww.^2)]);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function dataOut = head(dataIn, num)

dataIn = dataIn(:);

if num>0
    dataOut = dataIn(1:num, :);
elseif num<0
    dataOut = dataIn(1:(length(dataIn)+num));
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function dataOut = tail(dataIn, num)

dataIn = dataIn(:);

if num>0
    dataOut = dataIn((length(dataIn) - num + 1):length(dataIn),:);
elseif num<0
    dataOut = dataIn((-num+1):length(dataIn),:);
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

This yields SOMETIMES plots like this, and sometimes it's trapped in the inner while loop (I need to fix this). As you can see, the plot shows, that the result is not ideal.

enter image description here

EDIT 4:

Maybe I didn't formulate my problem correctly at first: I want to generate a time series that follows a defined distribution (in my case approximately Half-Gaussian) and whose autocorrelation function follows a given autocorrelation function (e.g. up to lag 30). So I have two properties that my time series must fulfil.

I tried the following approaches (without success):

1)

  • generation of uniformly distributed samples
  • invese transform sampling with desired distribution
  • Reordering the samples to match the ACF
  • Failed: it takes far too long for my time series (1e6 samples) + acf could not be reached.
  • generation of normally distributed samples with desired autocorrelation (with the approach from jblood94)
  • invese transform sampling with uniform distribution (based on this solution)
  • invese transform sampling with desired distribution
  • Failed: the autocorrelation of the transformed samples does not correspond to the desired Acf, which was correct before the transformation.

I think approach 2 might work. However, it must be prevented that the autocorrelation is changed by the transformation. Or the initial autocorrelation must be adjusted in such a way that it matches the desired Acf again, at least after the transformation.

Now the question is how the initial autocorrelation needs to be adjusted. Or is the approach correct? First generate normally distributed samples with autocorrelation, transform them into normally distributed samples, and then transform them into my desired distribution?

In addition, so far I have only managed to change individual lags of the normal distribution, but not that lag 1-30, for example, follow my function.

EDIT 5:

I tested approach 2 again. It seems, that the change in the autocorrelation is dependent from the acf itself. I generated samples, which follow now my distribution and my acf. I will post plots later.

This approach seems to work for me, but I can't guarantee, that it works in general.

$\endgroup$
8
  • 3
    $\begingroup$ Greetings! It may be helpful to show specifically what you have tried so far. This question is fairly vague and doesn't give much for people to work with I'm afraid. $\endgroup$ Feb 22, 2023 at 14:46
  • 1
    $\begingroup$ Hi! So far i do not have a entry point for my problem with the correlated uniform distributed samples. All references, which i found so far only explain how to generate random samples with a single correlation-coefficent, but not with a auto-correlation function. For my auto-correlation i have 30 lags with the corresponding value of my e-function which approximates the auto-correlation-function. $\endgroup$
    – Christian
    Feb 22, 2023 at 14:55
  • 3
    $\begingroup$ I vote to reopen. The answer linked to in the close vote does not answer this question. $\endgroup$ Feb 22, 2023 at 20:50
  • 1
    $\begingroup$ @JarleTufto the nature of the question was a bit changed due to the answers giving more straightforward approaches to generating correlated variables. But your answer explains nicely why this question is not a duplicate and how finding the Gaussian coppula with a matching correlation structure for the uniform variables is still a difficult task and seperate from the linked duplicate. $\endgroup$ Feb 22, 2023 at 23:50
  • 4
    $\begingroup$ Related literature Luc Devroye, Gérard Letac (2015) Copulas with prescribed correlation matrix $\endgroup$ Feb 23, 2023 at 0:04

8 Answers 8

3
$\begingroup$

This problem can be approached by first generating samples from the desired distribution and then reordering them to match the desired autocorrelation function.

The R code below demonstrates an approach based on this answer that can be modified for any desired ACF and distribution. The example generates $n=10^6$ samples from $\text{Gamma}(0.9,1)$ whose ACF follows 50 random samples from $\text{Beta}(1,3)$, sorted descending.

The process is as follows.

  1. Generate the desired number of samples, $X=\{x_1,...,x_n\}$, from the target distribution. Set $\alpha_0$ equal to the desired ACF. Initialize the target ACF, $\alpha$, as the desired ACF.
  2. Find a set of weights that, when passed to filter along with $n$ random normal variates, results in a series, $Y$, with $ACF=\alpha$ (see the answer linked above).
  3. Reorder $X$ so that its rank ordering matches that of $Y$. If $X$ is normally distributed, the resulting series should have the desired ACF; however, the more $X$ deviates from normality, the more the ACF will deviate from $\alpha_0$ (the example below has a target distribution of $\text{Gamma}(0.9,1)$, which is very "non-normal"). Update the target ACF, $\alpha$, according to $\alpha'=\frac{\alpha}{2}\Big(\frac{\alpha_0}{ACF}+1\Big)$ and repeat steps 1-3 until the ACF of the reordered $X$ converges.

The function that performs the reordering (it works only for positive values for alpha):

acf.reorder <- function(x, alpha) {
  tol <- 1e-5
  maxIter <- 10L
  n <- length(x)
  xx <- sort(x)
  y <- rnorm(n)
  w0 <- w <- alpha1 <- alpha
  m <- length(alpha)
  i1 <- sequence((m - 1):1)
  i2 <- sequence((m - 1):1, 2:m)
  i3 <- cumsum((m - 1):1)
  tol10 <- tol/10
  iter <- 0L
  x <- xx[rank(filter(y, w, circular = TRUE))]
  SSE0 <- Inf
  f <- function(ww) {
    sum((c(1, diff(c(0, cumsum(ww[i1]*(ww[i2]))[i3]))/sum(ww^2)) - alpha1)^2)
  }
  ACF <- function(x) acf(x, lag.max = m - 1, plot = FALSE)$acf[1:m]
  
  while ((SSE <- sum((ACF(x) - alpha)^2)) > tol) {
    if (SSE < SSE0) {
      SSE0 <- SSE
      w <- w0
    }
    if ((iter <- iter + 1L) == maxIter) break
    w1 <- w0
    a <- 0
    sse0 <- Inf
    
    while (max(abs(alpha1 - a)) > tol10) {
      a <- c(1, diff(c(0, cumsum(w1[i1]*(w1[i2]))[i3]))/sum(w1^2))
      
      if ((sse <- sum((a - alpha1)^2)) < sse0) {
        sse0 <- sse
        w0 <- w1
      } else {
        # w0 failed to converge; try optim
        w1 <- optim(w0, f, method = "L-BFGS-B")$par
        a <- c(1, diff(c(0, cumsum(w1[i1]*(w1[i2]))[i3]))/sum(w1^2))
        if (sum((a - alpha1)^2) < sse0) w0 <- w1
        break
      }
      
      w1 <- (w1*alpha1/a + w1)/2
    }
    
    x <- xx[rank(filter(y, w0, circular = TRUE))]
    alpha1 <- (alpha1*alpha/ACF(x) + alpha1)/2
  }
  
  xx[rank(filter(y, w, circular = TRUE))]
}

Generate samples from the target distribution and specify the desired ACF:

set.seed(1960841256)
x <- rgamma(1e6, 0.9, 1)
alpha <- c(1, sort(rbeta(50, 1, 3), TRUE))

Reorder x and plot its ACF against alpha:

system.time(x <- acf.reorder(x, alpha))
#>    user  system elapsed 
#>    7.13    0.41    7.53
acf(x, lag.max = length(alpha) - 1)
lines(seq_along(alpha) - 1, alpha, col = "green")

enter image description here

The resulting ACF is a good match with the target, and since x has simply been reordered, it is known to have the desired distribution.


Update: Worst of 100 iterations

The original code was breaking too early for some cases as noted by the OP in the updated question. The above code has been modified to correct the behavior. I also removed the restriction that the ACF be strictly positive.

I ran the same procedure 100 times, each with a new random vector for x and alpha. The longest-running iteration took 7.36 seconds, and the worst performing iteration (by maximum absolute difference of the achieved vs. desired ACF) had a maximum absolute error of 0.056. It is plotted below.

enter image description here

$\endgroup$
16
  • 1
    $\begingroup$ Hi jblodd94. Thanks for your answer. this seems to be the solution of my problem. After the sorting, the samples still follow the initial distribution AND also follow the desired acf. Do you have literature for your approach? I translated your R code to matlab, but till now i did not manage to get your algorithm working. $\endgroup$
    – Christian
    Mar 8, 2023 at 16:18
  • 1
    $\begingroup$ Good to hear you got it working in Matlab. Sorry, no literature to cite here. The approach was just off the top of my head without any literature search, so it is possible someone else has published something similar. $\endgroup$
    – jblood94
    Mar 8, 2023 at 16:50
  • $\begingroup$ As an aside, I'm fairly confident this approach can be extended to work with mixed-sign autocorrelation values. Your "real data" showed strictly positive values, so I didn't bother. $\endgroup$
    – jblood94
    Mar 8, 2023 at 16:55
  • 1
    $\begingroup$ +1 clever solution. But gadzooks! Using optim in a while statement might be a CPU smasher if $n$ is big and $\alpha$ is high. I wonder if you can't just calculate the fraction $p$ of the sample $U = U_1, \ldots U_n$ that should be replaced with the order statistics... $\endgroup$
    – AdamO
    Mar 8, 2023 at 17:25
  • $\begingroup$ @AdamO, it's actually optim in a nested while :). The break conditions seem to keep it manageable. I tried to demonstrate with a hard example (1M samples from a very non-normal distribution and a not-very-smooth ACF) and it took ~1 min. I just added the timing. $\endgroup$
    – jblood94
    Mar 8, 2023 at 17:41
9
$\begingroup$

Letting $U_t=\Phi(X_t)$ where $X_t$ is a zero-mean and unit-variance stationary Gaussian process with autocorrelation function $\rho_h$ and $\Phi$ is the standard normal cdf, it follows that each $U_t$ is marginally uniformly distibuted. It follows that the relation between $\rho_h$ and the autocovariance function of $U_t$ is \begin{align} \gamma_h&=\operatorname{Cov}(U_t,U_{t+h}) \\&=E(U_tU_{t+h})-E(U_t)E(U_{t+h}) \\&=E(\Phi(X_t)\Phi(X_{t+h}))-1/4 \\&=P(Z_1\le X_t\cap Z_2\le X_{t+h})-1/4 \\&=P\left(\frac{Z_1-X_1}{\sqrt{2}}\le0\cap \frac{Z_2-X_{t+h}}{\sqrt{2}}\le 0\right)-1/4 \\&=\Phi_2\left(\begin{bmatrix}0\\0\end{bmatrix};\frac{\rho_h}2\right)-1/4. \end{align} Here $Z_1$ and $Z_2$ are independent standard normal random variables, and $\Phi_2$ is the cdf of the standard bivariate normal distribution with correlation $\rho_h/2$.

Using this result, the above relation simplifies to $$ \gamma_h = \frac1{2\pi}\operatorname{arcsin}\frac{\rho_h}2. $$ Solving this equation, we find that the autocorrelation of $X_t$ needed to achieve a target autocovariance $\gamma_h$ of $U_t$ is $$ \rho_h = 2\sin(2\pi \gamma_h). $$

Even if the target autocovariance function $\gamma_h$ is positive semi-definite, the above construction may not be feasible. For example, if the target autocovariance function is that of a MA(1)-process with a unit root, $$ \gamma_h=\begin{cases} 1/12 &\text{for }h=0 \\ 1/24 &\text{for }h=1 \\ 0 &\text{for }h>1 \end{cases}, $$ this would imply that $$ \rho_h=\begin{cases} 1 &\text{for }h=0 \\ 0.518 &\text{for }h=1 \\ 0 &\text{for }h>1 \end{cases} $$ which is not a positive semi-definite autocorrelation function.

$\endgroup$
4
  • $\begingroup$ Hi Jarle Tufto, i edited my post, so that my problem should be clearer. As i understand you correct, its not feasible to generate uniformly distributed samples with a desired ACF or at least a approximation to the ACF of my real data? $\endgroup$
    – Christian
    Feb 23, 2023 at 11:12
  • $\begingroup$ @Christian My answer shows that how you can generate uniformly distributed $U_t$ with a given autocovariance function $\gamma_h$. The construction may or may not exists since $\rho_h$ needs to positive semi-definite. But in your EDIT2 you seem to want to further thansform $U_t$ into something that matches your empirical distribution function (via the inverse cdf method). This additional transformation will also further change the autocovariance function so you need to further tweak the autocorrelation function of $X_t$ to achieve your goal (maybe via numerical integration and root finding). $\endgroup$ Feb 23, 2023 at 11:24
  • $\begingroup$ Could you provide an example code for your approach? $\endgroup$
    – Christian
    Mar 14, 2023 at 16:24
  • $\begingroup$ @Christian In general, you would need to use mvtnorm::rmvnorm to simulate $X_t$ with a given autocovariance function $\rho_h$. The toeplitz function can be used to construct the required covariance matrix from a vector of autocorrelations. If $X_t$ is a moving average process ($\rho_h=0$ for lags $h>q$), then arima.sim would be much more efficient. $\endgroup$ Mar 15, 2023 at 8:57
7
$\begingroup$

Here's a pragmatic and easy approach with room to expand and establish proofs, with the focus on the main problem: how do you generate a correlated uniform sample?

Let $U_1, U_2$ be uniformly distributed and independent on the unit interval. Let $V_1 = U_1$. For a desired correlation $d$, let $G_{1,2}$ be yet another uniformly distributed random variable.

Let $$V_2 = \left\{ \begin{array}{ccc} U_1 & \text{if} & G_{1,2} < d \\ U_2 & \text{if} &G_{1,2} \ge d\end{array} \right.$$

I claim that:

  • One can show that $V_2$ is indeed uniformly distributed, with a correlation of $d$ with $V_1$
  • By way of induction, one can set an arbitrary sequence of uniform random variables and create a new sample having a desired covariance structure.

As usual, this site is always most compelled by code, so I can show a couple cases of the spherical AR-1 auto correlation, where the correlation between observations with a lag of 1 is set to $d$, but otherwise it's relatively straightforward to use any structure you want.

do.one <- function(n,N,d) {
  u <- matrix(runif(n*N), N, n)
  g <- matrix(runif((n-1)*N/2), N, n-1)

  v <- u
  for( i in 2:n) {
    v[g[, i-1] < d,i] <- v[g[, i-1] < d, i-1]
  }  
    
  acf <- sapply(2:n, function(i) cor(v[,i], v[, 1]))
  acf
}

set.seed(123)
ds <- c(0,0.25, 0.5, 0.95)
acfs <- sapply(ds, do.one, n=10, N=1000)
matplot(acfs,type='l', xlab='Lags')
legend('topright', title = 'AR-1 spherical correlation', legend = ds, lty=1:4, col=1:4)

enter image description here

$\endgroup$
5
  • 1
    $\begingroup$ Interesting. Basically a sort of random walk where you stay in the same place with probability $d$ and transition to a random position in the unit interval with probability $1-d$. $\endgroup$
    – Luca Citi
    Feb 23, 2023 at 0:09
  • $\begingroup$ What's also interesting is that the Gaussian version of your technique is indistinguishable from a standard AR(1) based on mean, variance and ACF alone, while the actual time series is very different. $\endgroup$
    – Luca Citi
    Feb 23, 2023 at 7:37
  • 1
    $\begingroup$ You might extend your answer to explain that your method gives ACF $\rho_h=d^{|h|}$. Should a more general ACF be required, one can increase the lag, for example with lag 2: $U_k=\begin{cases}U_{k-1}&\text{with probability $d_1$}\\ U_{k-2}&\text{with probability $d_2$}\\ \mathrm{Unif}(0,1)& \text{otherwise}\end{cases}$ $\endgroup$
    – Luca Citi
    Feb 23, 2023 at 7:47
  • $\begingroup$ It is also worth noting that some ACF profiles may not be achievable with the above approach. For example one cannot have $\rho_1>0$ and $\rho_2<\rho_1^2$. This can be easily overcome by allowing some of the lags to have terms $1-U_{k-\ell}$ rather than $U_{k-\ell}$. For example, the case I just mentioned can be obtained with: $U_k = \begin{cases} U_{k-1} & \text{with probability $d_1$}\\ 1 - U_{k-2} & \text{with probability $d_2$}\\ \mathrm{Unif}(0,1) & \mathrm{otherwise} \end{cases} $. $\endgroup$
    – Luca Citi
    Feb 24, 2023 at 22:04
  • $\begingroup$ @LucaCiti and AdamO: Out of interest I posted another answer spelling out some of details of what you propose. $\endgroup$ Feb 26, 2023 at 22:12
6
$\begingroup$

If I understand correctly your problem, you need to simulate random variables which marginally follow an uniform distribution while the joint distribution is a multivariate uniform distribution with some correlation between your marginals.

Formally:

$X_i\sim U(0,1), \forall_i, i = 1,2,...,K$ , where $X\sim U(a_i, b_i, Cor)$. In the method I am firstly proposing, your uniform random variables will be uniformly distributed between 0 and 1. However, I also provide a "Normal to Anything" kind of approach.

One trick is to use generate a Multivariate normal distribution, specifying a correlation matrix. Let's say $\Sigma$ is your auto-correlation structure for your $K$ random variables.

$$X \sim MVN(\mu, \Sigma)$$

$\mu$ is a vector of means of size $K$, while $\Sigma$ is a $K*K$ matrix

Then you have to transform your quantiles to probabilities from any normal distribution you want, for one $X$:

$$CDF_{Normal}(X_{i}, \mu_i, \Sigma_{ii}) \sim U(0,1)$$

Now if you assume that your variables are not uniformly distributed between a common support $[a,b], \forall K$ but a variable-dependent support instead. You just need to convert your probabilities obtained previously with an inverse uniform distribution specifying a and b. For a random variable,

$$\theta^{-1}(CDF_{Normal}(X_{i}, \mu_i, \Sigma_{ii}), a_i, b_i)$$

It is worth noting that you will respect globally your correlation structure, while under-estimate it a little bit. Some force-brute methods should be used to correct your final correlation structure.

I can provide script to illustrate the method.

Hope this helps.


set.seed(1234)

#Example function to generate auto-correlation matrix
autocorr.mat <- function(p = 100, rho = 0.9) {
  mat <- diag(p)
  autocorr <- sapply(seq_len(p), function(i) rho^i)
  mat[lower.tri(mat)] <- autocorr[unlist(sapply(seq.int(p-1, 1), function(i) seq_len(i)))]
  mat[upper.tri(mat)] <- autocorr[unlist(sapply(seq_len(p - 1), function(i) seq.int(i, 1)))]
  return(mat)
}


#Simulate normal data with some autocorrelation structure
autocorrelated.normal = MASS::mvrnorm(1000, mu= rep(0,100), autocorr.mat())

#Let's see the heatmap
heatmap(cor(autocorrelated.normal))

#If you are interested in simulating random uniform between 0 and 1 with your auto-correlation structure
pnorm(autocorrelated.normal)

#If you are interested in simulating random uniform between q and p, here [-10; 10]

qunif(pnorm(autocorrelated.normal), -10,10)

A brief comparison of the respective heatmaps give global concordance of initial correlation structure

heatmap(cor(autocorrelated.normal))
heatmap(cor(pnorm(autocorrelated.normal)))
heatmap(cor(qunif(pnorm(autocorrelated.normal), -10,10)))
$\endgroup$
13
  • 1
    $\begingroup$ +1 this is essentially a copula approach en.wikipedia.org/wiki/Copula_(probability_theory) $\endgroup$
    – AdamO
    Feb 22, 2023 at 18:49
  • $\begingroup$ Yes indeed, intuitive and easy to implement for this kind of problem $\endgroup$ Feb 22, 2023 at 19:04
  • 2
    $\begingroup$ More specifically, it is a Gaussian coppula. The sidenote is very important "It is worth noting that you will respect globally your correlation structure, while under-estimate it a little bit.". The correlation after transformation will not be the same. $\endgroup$ Feb 22, 2023 at 23:42
  • $\begingroup$ Hi Mangnier Loïc, i need to generate random samples for a simulation. The samples should have the probability distribution of my model AND should have the same ACF like my real data (the ACF is modeled by an e-function, which approximates the ACF of my real data quite well). Currently i generate the samples by feeding in uniformly distributed samples in a inverse transfom sampling. my generated samples follow now my desired distribution, BUT do not have the desired ACF. $\endgroup$
    – Christian
    Feb 23, 2023 at 11:22
  • 1
    $\begingroup$ @Christian If you are interested only in generating uniform RVs so my first solution is the way to go. However, if you are interested to simulate others RVs from others distributions, you will need to correct your correlation structure. I can provide additional details if you are interested. $\endgroup$ Feb 23, 2023 at 13:49
3
$\begingroup$

You can try implying AR(p) process coefficients $\phi_i$ from the given ACF $r(p)$. You could apply Yule Walker equations:

  • form a vector $r$ of ACF for lags $p$: $1, r_1, r_2,\dots, r_p$
  • construct a correlation matrix $R$ as described in the link above from $r_p$, e.g. the third row would be $(r_2,r_1,1,p_1,\dots,r_{p-2})$
  • calculate $\phi=R^{-1}r$

Use these coefficients to produce autocorrelated samples

$\endgroup$
2
  • 1
    $\begingroup$ This wouldn't generate samples that are uniform. $\endgroup$
    – AdamO
    Feb 22, 2023 at 18:32
  • $\begingroup$ @AdamO easy to transform them to uniform using CDF $\endgroup$
    – Aksakal
    Feb 22, 2023 at 20:05
2
$\begingroup$

The following spells out the details of the approach proposed in the other answer by @AdamO and in its comments by @LucaCiti.

For $i=1,2,\dots,\infty$, let $|\phi_i|$ denote the probability that $U_t$ takes a value identical to either $U_{t-i}$ or $1 - U_{t-i}$ and let these two possibilities be determined by the sign of $\phi_i$. Let the remaining fraction $$ \phi_0=1-\sum_{i=1}^\infty |\phi_i| $$ denote the probability that $U_t$ takes a uniformly distributed value independent of the history of the process. Clearly, we must have $$ 0\le \phi_0\le 1. \tag{1} $$ and $$ -1\le \phi_i\le 1 \tag{2} $$ for $i=1,2,\dots,\infty$.

Letting $V_t=U_t-\frac12$ denote the mean-centered process, and using the law of total expectation, we have \begin{align} E(V_t|V_{t-1},V_{t-2},\dots) &=|\phi_1|\operatorname{sgn}(\phi_1) V_{t-1} + |\phi_2|\operatorname{sgn}(\phi_2) V_{t-2} + \dots \\&=\phi_1 V_{t-1} + \phi_2 V_{t-2} + \dots. \end{align} Thus it is immedeately clear that $\phi_1,\phi_2,\dots$ are the coefficients in the $\operatorname{AR}(\infty)$ representation of the model.

Unlike ordinary Gaussian ARMA models, constraint (1) and (2) implies that not all positive semi-definite autocovariance functions are possible via this construction, however. For example, if the target autocovariance function is that of a MA(1) model with MA polynomial $1-\theta B$, the infinite AR polynomial would equal $$ \frac1{1-\theta B}=1+\theta B+\theta^2 B^2+\dots, $$ and we would have $$ \sum_{i=1}^\infty|\phi_i|=\sum_{i=1}^\infty |\theta^i|=\sum_{i=1}^\infty |\theta|^i = \frac{|\theta|}{1-|\theta|}. $$ Combined with (1) this limits possible values of $\theta$ to $$ |\theta|\le \frac12 $$ and the correlation at lag 1 to $$ -\frac25\le \rho_1=\frac{\theta}{1+\theta^2}\le\frac25 $$ In contrast, via the copula described in my other answer the correlation at lag 1 is limited to $|\rho_1|<0.4825837$ only . Semipositive definiteness in itself limits the same correlation to $|\rho_1|\le 1/2$.

$\endgroup$
1
$\begingroup$

Uniform distributed samples are the set of samples in which every element is distributed uniformly, if we place further constraints on this, it ceases to be uniform distributed samples. But anyway, if you meanе something else, we can easily upgrade my answer.

Let us just to generate samples with given autocorrelation function. Our idea is put all needed constraints on samples and let scipy.optimize to do everything for us.

While the key concept is simple, optimization and numerical problems can be arbitrary hard in real application and it can be necessary to adjust scipy solvers parameters here and/or optimize some computations.

The code implements this idea:

import numpy as np
import scipy.optimize as sco
import matplotlib.pyplot as plt

def cov(x,y):
    return np.sum(x*y)/len(x)

def cor(x,y):
    return cov(x,y)/(np.std(x)*np.std(y))

#define our autocorrelation
def autoMoment(x, moment, n):
    out = []
    for i in range(n):
        out = [moment(x,x)]
        for i in range(1,n):
            out.append(moment(x[i:],x[:-i]))
    return np.array(out)

#define generation of function for scipy.optimize
def generateWithSpecificAutomoment(moment, values): 
    def f(x):
        out = autoMoment(x, moment, len(values))
        return np.sum((np.array(out) - np.array(values))**2)
    return f

desiredAutocorr = [1,0.3,0.1,-0.23,0.03,0.07,-0.03]

initial = np.random.randn(50) #starting point for the solver
solution = sco.minimize(generateWithSpecificAutomoment(cor, desiredAutocorr), initial)

print(solution.success)
out = solution.x #samples obtained after optimization

#check autocorrelation on generated samples
realAutocorr = autoMoment(out, cor, len(desiredAutocorr))

#compare desired and real result
plt.title("Autocorrelations")
plt.plot(desiredAutocorr, linewidth=7, label="desired")
plt.plot(realAutocorr, label="real")
plt.legend()
plt.show()

Generated picture: Comparison desired and real results

$\endgroup$
4
  • 3
    $\begingroup$ "... if we place further constraints on this, it ceases to be uniform distributed samples." I think this statement needs to be qualified in some ways. For instance, you may set U1, U2 uniformly distributed. You can set U1=U2 or U1=-U2 to create a perfect counter example. $\endgroup$
    – AdamO
    Feb 22, 2023 at 18:07
  • $\begingroup$ Hi Kirill Erofeev, i need to generate random samples for a simulation. The samples should have the probability distribution of my model AND should have the same ACF like my real data (the ACF is modeled by an e-function, which approximates the ACF of my real data quite well). Currently i generate the samples by feeding in uniformly distributed samples in a inverse transfom sampling. my generated samples follow now my desired distribution, BUT do not have the desired ACF. $\endgroup$
    – Christian
    Feb 23, 2023 at 11:24
  • $\begingroup$ i am not sure, if its correct, that i simply adapt the uniformly distributed input-samples. But if so, then i am looking for an solution on how to modify the uniform distributed samples as mentioned above. $\endgroup$
    – Christian
    Feb 23, 2023 at 11:24
  • $\begingroup$ But with your approach i can solve at least the problem with the ACF of the input samples to my inverse transform sampling. Thank you!! $\endgroup$
    – Christian
    Feb 23, 2023 at 11:25
1
$\begingroup$

You can try 1d-OU Process:
$dX_t=-\alpha X_tdt+\sqrt{2\alpha}dW_t$
The autocorrelation function of $X_t$ is $\text{exp}(-\alpha t)$, where $dW_t\sim \mathcal{N}(0, \sqrt{dt})$
Assume that your exponential ACF fitting parameter is $0.08$ and $dt=0.1$, then the SDE is $dX_t = -0.08dX_t+0.4dW_t$.
Use memoryless transformations on simulation results $Y$ to get the results you need:
$Y = U^{-1}[\varPhi(X)] = a+\varPhi(X)\cdot (b-a)$
where $\varPhi(\cdot)$ is the CDF of the standard normal distribution with zero mean and unit variance; $U(a,b)$ is the CDF of the uniform distribution with lower bound $a$ and upper bound $b$.

ACF verification:

enter image description here

Histogram of $Y$:
$a=1, b=3$ in this case

enter image description here

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import ccf, acf
import scipy.stats as stats
import seaborn as sns

def oup(alpha, dt, num_steps, x0):
    sigma = np.sqrt(2*alpha)
    dW = np.sqrt(dt) * np.random.normal(size=num_steps)
    X = np.zeros(num_steps + 1)
    X[0] = x0
    for i in range(1, num_steps + 1):
        dX = - X[i-1] * alpha * dt + sigma * dW[i-1]
        X[i] = X[i-1] + dX
    return X

a = 0.08
dt = 0.1
T = 1e6
num_steps = int(T / dt)
lower = 1
upper = 3

X = oup(a, dt, num_steps, 0)
Y = lower + stats.norm.cdf(X, 0, 1) * (upper - lower)

num = int(30 / dt)
t = np.linspace(0, 30, num + 1)
acf_original = np.exp(-a*t)
X_acf = acf(X, nlags = num, fft = True)
Y_acf = acf(Y, nlags = num, fft = True)

plt.figure(figsize = (4, 3), dpi = 120)
plt.plot(t, acf_original, label = 'Original ACF', color = 'black')
plt.plot(t, X_acf, label = r'$X$', color = 'r', linestyle = '--')
plt.plot(t, Y_acf, label = r'$Y$', color = 'b', linestyle = '-')
plt.ylim(-0.1, 1.1)
plt.ylabel('ACF')
plt.xlabel('Lag [s]')
plt.legend()
plt.show()

sns.histplot(samples, stat = 'density', bins = 100)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.