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 )
}
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)
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)
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)
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 )
}
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.