6

I want to create a random variable with a given autocorrelation in R.

The target autocorrelation is defined by: $$acf_{target}=(lag+1)^{(-b)}$$ with $b=1.41519$ which I derived from a natural dataset. You can see the function in the next plot generated with this code:

lag <- seq(0,52)
b <- 1.41519
acf.target <- (lag+1)^(-b)
plot(lag, acf.target, col="green", type="l", xaxs="i", yaxs="i", las=1, lwd=2)

target autocorrelation function

Using the R approach shown in the answer to this question I created a random variable and then filtered the variable with my target autocorrelation function:

n <- 100000
var <- rnorm(n)
var_autocor <- filter(var, filter=acf.target, circular = T)

But if I calculate the autocorrelation of var_autocor I get a difference to the target autocorrelation function:

acf(var_autocor)
lines(lag, acf.target, col="green")

enter image description here

I calculate that the acf of the generated variable (var_autocor) has an $b_{result}=1.215306$ which slightly deviates from $b_{target}=1.41519$.

This deviation is more than I want and I don't understand why it appears. Do I use the filter function wrong? Do I have to give additional arguments? Or does someone have a totally different approach to create a random variable with the target autocorrelation?

Community
  • 1
nnn
  • 161

1 Answers1

2

In this case, each value returned by filter is a weighted sum of $m$ out of $n$ normal random variates $\{x_1,...,x_n\}$ with the $m$ weights, $w$, set equal to acf.target in the example provided by the OP. The $i^\text{th}$ term is:

$$y_i=\sum_{j=1}^m{w_j{}x_{\text{mod}(i-j+\lfloor(m-1)/2\rfloor,n)+1}}$$

The $i^{th}$ term returned by the acf function is:

$$a_i=\text{cor}(\{y_1,...,y_{n-i+1}\},\{y_i,...,y_{n}\})$$ $$=\frac{\text{cov}(\{y_1,...,y_{n-i+1}\},\{y_i,...,y_{n}\})}{\sum_{j=1}^mw_j^2}$$ $$\lim_{n\to\infty}a_i=\frac{\sum_{j=1}^{m-i+1}w_jw_{j+i-1}}{\sum_{j=1}^mw_j^2}\neq{w_i}$$

x <- rnorm(1e5)
b <- 1.41519
acf.target <- (1:53)^(-b)
w <- acf.target
n <- 1:(length(w) - 1)

y <- filter(x, w, circular = TRUE)
acf(y)
a <- c(1, sapply(n, function(i) sum(head(w, -i)*tail(w, -i)))/sum(w^2))
lines(seq_along(a) - 1, a, col = "green")

enter image description here

If $\{\alpha_i,...\alpha_m\}$ are the desired autocorrelation values, a naive approach would be to iteratively update $w_i$:

$$w_i'=w_i\frac{\alpha_i}{a_i}$$

This works for the OP's example:

tol <- sqrt(.Machine$double.eps)

while (max(abs(acf.target - a)) > tol) {
  a <- c(1, sapply(n, function(i) sum(head(w, -i)tail(w, -i)))/sum(w^2))
  w <- wacf.target/a
}


y <- filter(x, w, circular = TRUE)
acf(y)
lines(seq_along(acf.target) - 1, acf.target, col = "green")

enter image description here

jblood94
  • 1,459