September 19, 2017 / by / In /

Supplemental Material: A Study of the Allan Variance for Constant-Mean Non-Stationary Processes

In the paper A Study of the Allan Variance for Constant-Mean Non-Stationary Processes, we provide a new formula to compute the theoretical Allan Variance (AVar) for constant-mean non-stationary processes and discuss its practical implications. This formula is provided for both the Maximum Overlapping AVar (MOAV) and the Non Overlapping AVar (NOAV) and a few examples are given for the former (MOAV) using the following processes

  • Non-Stationary White Noise

  • Bias-Instability

  • AR1 Blocks

These processes are used to compare the empirical MOAV, the theoretical MOAV for stationary processes and the new form of the theoretical non-stationary MOAV. In the following paragraphs, we present R code snippets that allow to reproduce the examples presented in the paper.

Installing and Loading Required Packages

Several R packages are required to run the code presented in this document. These packages can be installed as follows:

install.packages(c("MASS","devtools","knitr","rmarkdown"))
devtools::install_github("smac-group/simts")
devtools::install_github("smac-group/wv")

Once the packages are installed, they can easily be “loaded” as follows:

library(simts)
library(wv)
library(MASS)

Examples

Non-Stationary White Noise

The first process we study in the paper is the non-stationary white noise by which we intend, without loss of generality, all those zero-mean processes whose variance changes with time. The goal of studying the AVar for these types of processes is to understand whether it is able to detect such a structure in a time series and if it can be distinguished from a stationary white noise process. To do so, we perform some simulation studies (see below) to compare the empirical MOAV and the two forms of theoretical MOAV (stationary and non-stationary). Since the true process is non-stationary (i.e. non-constant variance), the stationary theoretical form was computed using the average variance over all times considered.

# set up the process
n_total = 1000       # length of ts in total
n_block = 10         # length of ts in each block
n_sims = 50          # number of simulations for empirical AVar
nb_scales = floor(log2(n_total)) - 1

sigma2 = (1+n_total)/2  # theoretical sigma2


# plot empirical AVar
for (i in (1:n_sims)){

  # generate the process
  wn = gen_nswn(n_total = n_total, seed = i + 1)

  # calculate empirical allan variance
  prac_av_wn = avar(wn, type="mo")$allan

  # plot
  if (i == 1){

  # background plot
  name.x = parse(text = paste("2^", 1:nb_scales, sep = ""))
  name.y = c(0.1,   0.2,   0.4,   0.8,   1.6,   3.2,   
             6.4,  12.8,  25.6,  51.2, 102.4)
  plot(NA, ylim = c(0.1, 300), xlim = c(1, nb_scales), log = "y", 
       xaxt = "n",  yaxt = "n", xlab = bquote("Scale("*tau*")"), 
       ylab = "Allan Variance", cex.lab=1.2, cex.axis=1, cex.main=1.2, 
       cex.sub=1.2)
  grid()

  # draw an axis on the left
  axis(1, at=1:nb_scales,labels = name.x, cex.lab=1, cex.axis=1, cex.sub=1)
  axis(2, at=name.y,labels = name.y, cex.lab=1, cex.axis=1, cex.sub=1)

  # assign different colors
  coleur = hcl(h = seq(15, 375, length = 3), 
               l = 65, c = 100, alpha = 0.5)[1:2]
  coleur2 = hcl(h = seq(15, 375, length = 4), 
               l = 65, c = 100, alpha = 1)[1:3]
  lines(1:8, prac_av_wn, col = coleur[2])

  }else{
  lines(1:8, prac_av_wn, col = coleur[2])
  }

  # calculate theoretical stationary allan variance
  theo_av_wn = av_wn(sigma2 = sigma2 ,
                               n = avar(wn, type="mo")$clusters)

  # calculate theoretical non-stationary allan variance
  covmat = covmat_nswn(sigma2 = sigma2, n_total = n_total)
  theo_av_wn2 = sapply(1:(floor(log2(n_total))-1), 
                       function(y){MOAV(2^y, covmat)})

  # plot the theoretical ones

  lines(theo_av_wn, col = coleur2[1], lwd = 2)
  points(theo_av_wn, col = coleur2[1], pch = 17, cex = 2.5)

  lines(theo_av_wn2, col = "black", lwd = 3)
  points(theo_av_wn2, col = "black", pch = 16, cex = 1.2)

  legend("bottomleft", c("Empirical MOAV","Theoretical Non-Stationary MOAV",
                         "Theoretical MOAV for Stationary WN"), 
         col = c(coleur[2], "black", coleur[1]), lwd = c(1,2,2),
  pch = c(NA, 16, 17), bty = "n", cex = 1  )

}

Figure 1: Logarithm of the MOAV of the non-stationary white noise process for scales 2^n. Estimated MOAV (light-blue lines); theoretical non-stationary MOAV (black line with dots) and theoretical stationary AV based on the average variance (red line with triangles).

From Figure 1, we can see how both theoretical forms correspond and the estimated MOAV closely follows these quantities. This example confirms that the AVar is unable to distinguish between a stationary white noise process and a white noise process whose second-order behavior is non-stationary.

Bias-Instability

The second process we study in the paper is a bias-instability process which is a commonly known process in the engineering domain, specifically for inertial sensor calibration for navigation. The characteristic of this process is that it consists in different concatenated sequences (blocks) where, within each block, the realization of a random variable is repeated (i.e. constant). One realization of the bias-instability process is illustrated in Figure 2.

# set up the process
n_total = 2500       # length of ts in total
n_block = 10         # length of ts in each block
n_sims = 50          # number of simulations for empirical AVar
nb_scales = floor(log2(n_total)) - 1

sigma2 = 1

# realization
bi_real = gen_bi(sigma2 = sigma2, n_total = n_total, n_block = n_block)

plot(bi_real, type="l", xlab = "time", ylab = expression(X[t]), cex.lab=1.2, cex.axis=1, cex.main=1.2, cex.sub=1.2)

Figure 2: Realization of a bias-instability process with variance 1. We use a length of each block to be 10, and a total of 250 blocks.

As the theoretical form of the AVar for this process is not known exactly, it is often approximated by the AVar of a First-Order AutoRegressive (AR1) process. As before, some simulations are given below in which we study and compare the custom AR1 approximation (whose parameters are estimated using maximum likelihood) and the theoretical form of the MOAV given in the paper.

# set up the process
n_total = 2500       # length of ts in total
n_block = 10         # length of ts in each block
n_sims = 50          # number of simulations for empirical AVar
nb_scales = floor(log2(n_total)) - 1

sigma2 = 1

# plot empirical AVar

theta_matrix = matrix(0, nrow = n_sims, ncol = 2)
av_matrix = matrix(0, nrow = n_sims, ncol = (floor(log2(n_total))-1))

for (i in (1:n_sims)){
  
  # generate the process
  bi = gen_bi(sigma2 = sigma2, n_total = n_total, 
              n_block = n_block, seed = i + 1)

  # calculate empirical AVar
  prac_av_bi = avar(bi, type = "mo")$allan

  # plot empirical AVar
  if (i == 1){
  # background plot
  name.x = parse(text = paste("2^", 1:nb_scales, sep = ""))
  name.y = c(0.01,  0.02,  0.04,  0.08,  
             0.16,  0.32,  0.64,  1.28,  2.56,  5.12)
  plot(NA, ylim = c(0.01,1), xlim = c(1,nb_scales), 
       log = "y", xaxt = "n",  yaxt = "n", xlab = bquote("Scale("*tau*")"), 
       ylab = "Allan Variance", cex.lab=1.2, cex.axis=1, cex.main=1.2, cex.sub=1.2)
  grid()

  # draw an axis on the left
  axis(1, at=1:nb_scales,labels = name.x, cex.lab=1, cex.axis=1, cex.sub=1)
  axis(2, at=name.y,labels = name.y, cex.lab=1, cex.axis=1, cex.sub=1)

  # assign different colors
  coleur = hcl(h = seq(15, 375, length = 3), 
               l = 65, c = 100, alpha = 0.5)[1:2]
  coleur2 = hcl(h = seq(15, 375, length = 4), 
               l = 65, c = 100, alpha = 1)[1:3]
  lines(1:10, prac_av_bi, col = coleur[2])

  }else{
  lines(1:10, prac_av_bi, col = coleur[2])
  }


  # store values for approximate stationary
  fit = arima(bi, order = c(1,0,0))  # fit AR1 model to the bi data
  theta_matrix[i,] = c(as.numeric(fit$coef[1]), fit$sigma2)  # phi, sigma2
  av_matrix[i,] = av_ar1(n = avar(bi, type = "mo")$clusters,
                                 phi = theta_matrix[i,1],
                                 sigma2 = theta_matrix[i,2])

  legend("bottomleft", c("Empirical MOAV","Theoretical Non-Stationary MOAV",
                         "Theoretical MOAV for Closest Stationary AR1"), 
         col = c(coleur2[3], "black", coleur[1]), lwd = c(1,2,2),
  pch = c(NA, 16, 17), bty = "n", cex = 1  )

}

# calculate approximate stationary
theo_av_bi = colMeans(av_matrix, na.rm = FALSE, dims = 1)


# calculate theoretical non-stationary AVar
covmat = covmat_bi(sigma2 = sigma2, n_total = n_total, n_block = n_block)
theo_av_bi2 = sapply(1:(floor(log2(n_total))-1), 
                     function(y){MOAV(2^y, covmat)})

# plot the theoretical ones
lines(theo_av_bi, col = coleur2[1], lwd = 2)
points(theo_av_bi, col = coleur2[1], pch = 17, cex = 1.2)

lines(theo_av_bi2, col = "black", lwd = 3)
points(theo_av_bi2, col = "black", pch = 16, cex = 1.2)

Figure 3: Logarithm of the MOAV of the bias-instability process for scales 2^n. Estimated MOAV (light-blue lines); theoretical non-stationary MOAV (black line with dots) and theoretical stationary MOAV of an AR1 approximating the bias-instability MOAV (red line with triangles).

In Figure 3, the red line represents the AVar of a stationary AR1 process which is supposed to approximate the true AVar of a bias-instability process. Although close over some scales, it is clear that this approximation is not precise enough to be meaningful in practice. Therefore, knowing the exact form of the AVar for this process would allow to better interpret the signals characterized by bias-instability.

AR1 Blocks

The last process we study in the paper is a block-structure AR1 process. It is defined as a process whose parameters are fixed but is made by concatenated time periods (blocks) where observations within each block are generated independently from those in the other blocks. This setting is typical, for example, within longitudinal studies. One realization of the AR1 blocks process is shown in Figure 4.

# set up the process
n_total = 1000       # length of ts in total
n_block = 10         # length of ts in each block
n_sims = 50          # number of simulations for empirical AVar
nb_scales = floor(log2(n_total)) - 1

phi = 0.9
sigma2 = 1
scale = 100

# realization
ar1_real = gen_ar1blocks(phi = phi, sigma2 = sigma2, 
                         n_total = n_total, n_block = n_block, scale = 10)

plot(ar1_real, type="l", xlab = "time", ylab = expression(X[t]), cex.lab=1.2, cex.axis=1, cex.main=1.2, cex.sub=1.2)

Figure 4: Realization of a block-structure autoregressive process with phi = 0.9 and variance = 1. We use a length of each block to be 10, and a total of 100 blocks.

In order to study and compare the AVar when the block structure is either ignored or taken into consideration, the following simulations were carried out.

# set up the process
n_total = 1000       # length of ts in total
n_block = 10         # length of ts in each block
n_sims = 50          # number of simulations for empirical AVar
nb_scales = floor(log2(n_total)) - 1

phi = 0.9
sigma2 = 1
scale = 100

# plot empirical AVar
for (i in (1:n_sims)){

  # generate AR1 blocks process
  ar = gen_ar1blocks(phi=phi, sigma2=sigma2, n_total=n_total, 
                     n_block=n_block, scale=scale, seed = i + 1)

  # calculate empirical AVar
  prac_av_ar = avar(ar, type = "mo")$allan

  # plot
  if (i == 1){
  # Background plot
    name.x = parse(text = paste("2^", 1:nb_scales, sep = ""))
    name.y = c(0.01,  0.02,  0.04,  0.08,  0.16,  
             0.32,  0.64,  1.28,  2.56,  5.12)
    plot(NA, ylim = c(0.01,8), xlim = c(1,nb_scales), log = "y", 
         xaxt = "n",  yaxt = "n", xlab = bquote("Scale("*tau*")"), 
         ylab = "Allan Variance", cex.lab=1.2, cex.axis=1,
         cex.main=1.2, cex.sub=1.2)
    grid()

    # Draw an axis on the left
    axis(1, at=1:nb_scales,labels = name.x, cex.lab=1, 
         cex.axis=1, cex.sub=1)
    axis(2, at=name.y,labels = name.y, cex.lab=1, 
         cex.axis=1, cex.sub=1)

    # Assign different colors
    coleur = hcl(h = seq(15, 375, length = 3), 
               l = 65, c = 100, alpha = 0.5)[1:2]
    coleur2 = hcl(h = seq(15, 375, length = 4), 
                l = 65, c = 100, alpha = 1)[1:3]
    lines(1:8, prac_av_ar, col = coleur[2])

  }else{
    lines(1:8, prac_av_ar, col = coleur[2])
  }

  # theoretical covmat
  covmat = covmat_ar1blocks(n_total=n_total, n_block=n_block, 
                               phi=phi, sigma2=sigma2)

  # calculate theoretical stationary AVar
  theo_av_ar = av_ar1(n = avar(ar, type = "mo")$clusters,
                               phi = phi, sigma2 = sigma2)

  # calculate theoretical non-stationary AVar
  theo_av_ar2 = sapply(1:(floor(log2(n_total))-1), 
                       function(y){MOAV(2^y, covmat)})

  # plot the theoretical ones
  lines(theo_av_ar, col = coleur2[1], lwd = 2)
  points(theo_av_ar, col = coleur2[1], pch = 17, cex = 1.2)

  lines(theo_av_ar2, col = "black", lwd = 3)
  points(theo_av_ar2, col = "black", pch = 16, cex = 1.2)

  legend("bottomleft", c("Empirical MOAV","Theoretical Non-Stationary MOAV",
                         "Theoretical MOAV for Stationary AR1"), 
         col = c(coleur2[3], "black", coleur2[1]), lwd = c(1,2,2), 
         pch = c(NA, 16, 17), bty = "n", cex = 1 )
}

Figure 5: Logarithm of the MOAV of the block-structure firstorder autoregressive process for scales 2^n. Estimated MOAV (light-blue lines); theoretical non-stationary MOAV (black line with dots) and theoretical stationary MOAV assuming no block structure (red line with triangles).

In Figure 5, it can be observed how the stationary form of the AVar (that does not consider the block structure) is not close to the estimated AVar while the non-stationary form provided in the paper adequately represents this process and can therefore allow to distinguish between a stationary autoregressive process and a block-structure one.