# Standard deviation in generated data with specified autocorrelation - correction factor?

I want to generate time-series data with specified AR(1) and standard deviation. I am using arima.sim to generate univariate time series. (Ultimately, I want to use VAR.sim in tsDyn package for multivariate data, but I will keep it simple here with arima.sim for demonstration in a univariate situation).

I noticed that when AR is higher, the standard error will also be higher. A correction factor is needed to estimate the standard deviation.

This was discussed in both Section 1.1 of Beran (1994) and Bence (1995).

In Beran (1994)'s discussion on correction factor c(a) with a denoting AR(1) rho and n denoting number of observations,

$$\operatorname{var}(\bar{X}) =\sigma^{2}\left[1+\delta_{n}(a)\right] n^{-1}=\sigma^{2} c_{n}(a) n^{-1},$$ with $$\delta_{n}=\frac{2 a}{1-a}\left[1-n^{-1} \frac{1}{1-a}+n^{-1} \frac{a^{n}}{1-a}\right] .$$

From these equations, I would expect I can get the corrected variance by dividing the calculated variance with c(a). Or, dividing calculated standard error by square-root of c(a) to get the corrected estimation of standard deviation.

Putting these into my data-generating R codes, I noticed 2 things and have 2 questions:

1. The square-root of the correction factor has to be applied. Why is that so? (See the code - I needed to apply sqrt(ca) to get the corrected variance. I couldn't figure out why from the equations stated in the references)
2. The performance of the correction factor seems to be the worst when rho approaches 0.5 (or -0.5). Why? And is there some further adjustment I can do?

Sample R codes below. I applied the correction factor according to estimated AR(1) from acf instead of the specified rho because it is always a bit smaller than what I specified.

SD <- 1
n <- 100 # number of observations in one time series
nsim <- 500 # number of repetitions in simulation

colSdApply <- function(x, ...)apply(X=x, MARGIN=2, FUN=sd, ...)
decimal <- function(x, k) trimws(format(round(x, k), nsmall=k))

ca <- function(x){
# see section 1.1 of Beran (1994)
x<- abs(x)
cf <- 1+2*x/(1-x) # eq. 1.14
cf.smallN <- cf*(1-(1/(1-x)/n+ x^n/(1-x)/n)) # eq 1.12
# see correction factor k in https://ljmartin.github.io/technical-notes/stats/estimators-autocorrelated/
delta <- ((n-1)*x - n*x^2 + x^(n+1))/(1-x)^2
corrk <- sqrt((1+2*delta/n)/(1-2*delta/(n*(n-1))))
cf.smallN # or return corrk
}

test_arimasim = function(AR){
ts <- arima.sim(list(order=c(1,0,0), ar=AR), n=n, sd=SD)
resacf <- acf(ts,plot = FALSE)
return(list(ts,resacf\$acf[2]))
}

df.VAR <-  data.frame(AR = unlist(t(x)[,2]))
df.VAR$$SD <- colSdApply(do.call(cbind, x[1,])) df.VAR$$ca <- ca(df.VAR$$AR) df.VAR$$adjvar <- df.VAR$$SD^2 / sqrt(df.VAR$$ca)
df.VAR$$adjSD <- sqrt(df.VAR$$adjvar)
paste0("AR(empirical)=",decimal(mean(df.VAR$$AR),2), " | Naive SD=",decimal(mean(df.VAR$$SD),2),
" | Corrected SD=",decimal(mean(df.VAR$$adjSD),2), " | Correction Factor=",decimal(mean(df.VAR$$ca),2))
}

sdresult <- function(AR){
sim.arima <- (replicate(nsim,test_arimasim(AR),simplify = TRUE))