Nothing
#data(dfs)
#data(const.Q)
#data(var.n)
mscarm.filter <- function(time.series,
# input arguments of the univariate SCARM
right.width=30,
min.left.width=right.width,
min.width=floor(right.width/3),
max.width=200,
sign.level=0.001,
bound.noise.sd=0.01,
rtr=TRUE,
autocorrelations="automatic",
# input arguments of the SSM
c.bound=3,
r.bound=0
){
#######################################
# Stopping and rules and preparations #
#######################################
if(missing(time.series))
stop("The input data is missing with no default.\n")
if(!is.numeric(time.series))
stop("Data must be numeric.\n")
if(!is.matrix(time.series))
stop("'time.series' must be a matrix.\n")
length.series <- dim(time.series)[1]
dim.series <- dim(time.series)[2]
if(dim.series<2)
stop("'time.series' must be a multivariate time series of dimension >=2.\n")
for(i in 1:dim.series){
if(all(is.na(time.series[,i])))
stop("At least one variable has only missing values.\n")
}
if(length.series < dim.series){
time.series <- t(time.series)
length.series <- dim(time.series)[1]
dim.series <- dim(time.series)[2]
}
if(is.matrix(time.series)==FALSE){
time.series <- as.matrix(time.series)
}
r <- round(right.width)
ell.min <- round(min.left.width)
n.min <- round(min.width)
n.max <- round(max.width)
if(length.series < n.min)
stop("Length of time series must be at least 'n.min' =", n.min, ".\n")
if(n.min < 5){
n.min <- 5
warning("'n.min' must be an integer >=5; 'n.min' is set to 5.\n")
}
if(r < 5){
r <- 5
warning("'r' must be an integer >=5; 'r' is set to 5.\n")
}
if(ell.min < 5){
ell.min <- 5
warning("'ell.min' must be an integer >=5; 'ell.min' is set to 5.\n")
}
if(r > ell.min){
ell.min <- r
warning("'ell.min' must not be smaller than 'r'; 'ell.min' is set to 'r' = ", r,".\n")
}
if(sign.level <= 0 | sign.level > 0.5)
stop("'sign.level' must be a value in (0,0.5)")
if(n.max < n.min)
stop("'n.max' must not be smaller than 'n.min'.\n")
if(n.max < r + ell.min | n.max%%1!=0)
stop("'n.max' must not be smaller than 'r'+'ell.min' = ", r+ell.min,".\n")
if(bound.noise.sd <= 0){
bound.noise.sd <- 0.01
warning("'bound.noise.sd' must be a value >0; 'bound.noise.sd' is set to 0.01.\n")
}
if(c.bound <= 0){
c.bound <- 3
warning("'c.bound' must be a value >0; 'c.bound' is set to 3.\n")
}
if(r.bound < 0 | r.bound > 1){
r.bound <- 0
warning("'r.bound' must be a value in [0,1]; 'r.bound' is set to 0.\n")
}
if(!is.logical(rtr))
stop("'rtr' must be either TRUE or FALSE.\n")
if(all(autocorrelations!=c("high.positive", "moderate.positive", "small.positive", "no", "small.negative", "moderate.negative", "high.negative", "automatic")))
stop("'autocorrelations' must be either 'high.positive', 'moderate.positive', 'small.positive', 'no', 'small.negative', 'moderate.negative', 'high.negative' or 'automatic'.\n")
######################
# required functions #
######################
# Q scale estimator with bias-correction c.q depending on chosen value for 'autocorrelations'
if(autocorrelations=="high.negative"){
c.q.sim <- const.Q[,1]
c.q.model <- function(n){
c.q.sim[300]
}
emp.var <- var.n[,1]
v.model <- function(n){
out <- 0.0003514981 + (-0.0629949278/n) + (1.6446684583/(n^2)) + (-3.1859306183/(n^3))
return(out)
}
}
if(autocorrelations=="moderate.negative"){
c.q.sim <- const.Q[,2]
c.q.model <- function(n){
c.q.sim[300]
}
emp.var <- var.n[,2]
v.model <- function(n){
out <- 0.0002560433 + (-0.0434806450/n) + (0.9224130380/(n^2)) + (5.4903584072/(n^3))
return(out)
}
}
if(autocorrelations=="small.negative"){
c.q.sim <- const.Q[,3]
c.q.model <- function(n){
c.q.sim[300]
}
emp.var <- var.n[,3]
v.model <- function(n){
out <- 0.0002385982 + (-0.0395730416/n) + (0.7798496595/(n^2)) + (9.7595002645/(n^3))
return(out)
}
}
if(autocorrelations=="no"){
c.q.sim <- const.Q[,4]
c.q.model <- function(n){
1.21 * (n/(n+0.44))
}
emp.var <- var.n[,4]
v.model <- function(n){
out <- 4.77e-07 + (17.71*(1/n^3))
return(out)
}
}
if(autocorrelations=="small.positive"){
c.q.sim <- const.Q[,5]
c.q.model <- function(n){
c.q.sim[300]
}
emp.var <- var.n[,5]
v.model <- function(n){
out <- 0.000473727 + (-0.086723762/n) + (2.442400804/(n^2)) + (10.578915504/(n^3))
return(out)
}
}
if(autocorrelations=="moderate.positive"){
c.q.sim <- const.Q[,6]
c.q.model <- function(n){
c.q.sim[300]
}
emp.var <- var.n[,6]
v.model <- function(n){
out <- 0.0007355669 + (-0.1485876684/n) + (5.5086877566/(n^2)) + (-5.3249305456/(n^3))
return(out)
}
}
if(autocorrelations=="high.positive"){
c.q.sim <- const.Q[,7]
c.q.model <- function(n){
c.q.sim[300]
}
emp.var <- var.n[,7]
v.model <- function(n){
out <- -0.0003318919 + (0.0476670483/n) + (2.2422001348/(n^2)) + (-5.8731008014/(n^3))
return(out)
}
}
get.Q.adj <- function(y, phi.est){
n <- length(y)
m <- length(na.omit(y))
if(m < 3)
stop("data vector must contain at least 3 observations")
x <- which(!is.na(y))
y <- y[x]
if(autocorrelations=="automatic"){
autocor.type <- which.min(abs(phi.est-c(-0.9,-0.6,-0.3,0,0.3,0.6,0.9)))
c.q.sim <- const.Q[,autocor.type]
c.q.model <- function(n){
c.q.sim[300]
}
}
if(m > 4){
if(m <= 150){
c.q <- c.q.sim[m]
} else {
c.q <- c.q.model(m)
}
} else {
c.q <- 1
}
h <- numeric(m-2)
for(j in 1:(m-2)){
h[j] <- abs(y[j+1] - y[j] - (x[j+1] - x[j])*((y[j+2]-y[j])/(x[j+2]-x[j])))
}
h <- sort(h)
out <- c.q * sort(h)[max(1,floor(0.5*(m-2)))]
return(max(out,bound.noise.sd))
}
# covariance estimator OGK based on scale estimator Q.adj
get.Ogk <- function(x, noise.sd.est.s){
D.Ogk <- diag(noise.sd.est.s)
Y.Ogk <- x%*%solve(D.Ogk)
R.Ogk <- diag(1, nrow=ncol(x))
Rij.Ogk <- function(x){
((get.Q.adj(apply(x,1,sum), phi.est=0))^2 - (get.Q.adj(apply(x,1,diff), phi.est=0))^2)/4
}
comparisons <- t(combn(1:ncol(x),2))
for(i in 1:nrow(comparisons)){
column1 <- comparisons[i,1]
column2 <- comparisons[i,2]
R.Ogk[column1,column2] <- R.Ogk[column2,column1] <- Rij.Ogk(cbind(Y.Ogk[,column1], Y.Ogk[,column2]))
}
Eigen.R.Ogk <- eigen(R.Ogk)
E.Ogk <- Eigen.R.Ogk$vectors
Lambda.Ogk <- Eigen.R.Ogk$values
A.Ogk <- D.Ogk%*% E.Ogk
pre.Gamma.Ogk <- x%*%solve(t(A.Ogk))
Gamma.Ogk <- diag(apply(pre.Gamma.Ogk, 2, get.Q.adj, phi.est=0)^2)
Ogk.est <- A.Ogk %*% Gamma.Ogk %*% t(A.Ogk)
return(Ogk.est)
}
# function which gives the matrix of the slopes between each pair of observations in a window sample,
# required for fast update of RM slope estimation
get.B.t <- function(y){
n <- length(y)
j <- 1:n
B.t <- matrix(NA,n,n)
for (i in j) {
B.t[j[j<i],i] <- B.t[i,j[j<i]] <- (y[i] - y[j[j<i]]) / (i - j[j<i])
}
return(B.t)
}
# function which gives the RM slope estimation based on the output of function "get.B.t"
get.beta.RM <- function(B.t){
beta.i <- apply(B.t,1,median, na.rm=T)
beta.RM <- median(beta.i, na.rm=T)
return(beta.RM)
}
# function which gives the RM level estimation based on the RM slope estimation given by the function "get.beta.RM"
get.mu.RM <- function(y, beta.RM){
n <- length(y)
j <- 1:n
mu.RM <- median(y - (beta.RM*(j-n)), na.rm=T)
return(mu.RM)
}
# degrees of freedom for t-distribution to obtain critical values for the SCARM test statistic
get.critval <- function(ell, r, sign.level){
if(ell <= 100){
ell <- ell - ell%%5
r <- r - r%%5
degree.of.freedom <- dfs[ell/5,r/5]
critval <- qt(p=1-(sign.level/2), df=degree.of.freedom)
} else {
critval <- qnorm(p=1-(sign.level/2))
}
return(critval)
}
# obtain the scarm test statistic based on the outputs "B" and "noise.sd" of the functions "get.B.t" and "get.Q.adj"
scarm.test.statistic <- function(y, B, ell, r, noise.sd, phi.est){
n <- ell+r
left.sample <- y[1:ell]
right.sample <- y[(ell+1):n]
B.t.left <- B[1:ell,1:ell]
B.t.right <- B[(ell+1):n,(ell+1):n]
beta.left <- get.beta.RM(B.t.left)
beta.right <- get.beta.RM(B.t.right)
d.t <- beta.left - beta.right
if(autocorrelations=="automatic"){
autocor.type <- which.min(abs(phi.est-c(-0.9,-0.6,-0.3,0,0.3,0.6,0.9)))
emp.var <- var.n[,autocor.type]
v.model <- function(n){
emp.var[300]
}
}
if(ell<=300){
v.left <- emp.var[ell]
} else {
v.left <- v.model(ell)
}
if(r<=300){
v.right <- emp.var[r]
} else {
v.right <- v.model(r)
}
var.t <- noise.sd^2 * (v.left+v.right)
test.statistic <- d.t/(sqrt(var.t))
return(list(test.statistic=test.statistic, beta.left=beta.left, beta.right=beta.right))
}
# model for the correlation of two RM-slopes; depends on cross-correlation of the errors and the window widths
corr.RM.slope.model <- function(ratio,rho){
out <- (rho*0.08986669) + ((rho^3)*0.01312361) - (rho*ratio*1.17617121) + (rho*(ratio^2)*3.57618663) - (rho*(ratio^3)*1.72785517) - ((rho^3)*ratio*0.18067041)+ ((rho^3)*(ratio^2)*0.58289707) - ((rho^3)*(ratio^3)*0.19580252)
return(out)
}
# model for the variance of the RM-slope for a linear signal with Gaussian error in a window sample of length "width"
var.RM.slope.model <- function(width){
out <- 4.77e-07 + (17.71*(1/width^3))
return(out)
}
# function for block building based on the comparison statistics in "S.est"
build.blocks <- function(S){
S <- abs(S)
S[is.na(S)] <- c.bound+1
K <- 0.5 + sqrt(2*length(S)+0.25)
j <- 1
R <- 1:K
B <- list()
combs <- combn(R,2)
repeat{
if(all(S >= c.bound)){
if(length(R)>0){
for(i in 1:length(R)){
B[[j+i-1]] <- R[i]
}
}
K.t <- length(B)
break
}
if(all(S < c.bound)){
B[[j]] <- R
K.t <- length(B)
break
}
if(any(S < c.bound)){
pos <- which.min(S)
B[[j]] <- combs[,pos]
R <- R[is.na(match(R,B[[j]]))]
S[pos] <- Inf
combs[,pos] <- 0
}
repeat{
if(all(S >= c.bound)){
K.t <- length(B)
break
}
aligned <- integer()
aligned.pos <- integer()
for(i in 1:length(R)){
k <- R[i]
check1 <- which(apply(combs,2,function(x){any(!is.na(match(x,k))) & any(!is.na(match(x,B[[j]])))}))
check2 <- which(apply(combs,2,function(x){any(!is.na(match(x,k))) & any(!is.na(match(x,B[[j]])))}) & (S < c.bound))
if(length(check1)==length(check2) & length(check1)>0 & length(check2)>0 ){
aligned <- c(aligned,k)
aligned.pos <- c(aligned.pos,check1)
}
}
if(length(aligned)==0){
for(i in 1:dim(combs)[2]){
if(any(B[[j]]==combs[1,i]) | any(B[[j]]==combs[2,i])){
combs[,i] <- 0
S[i] <- Inf
}
}
break
}
best.aligned.pos <- min(which(S==min(S[aligned.pos])))
if(any(R==combs[1,best.aligned.pos])){
k <- combs[1,best.aligned.pos]
} else {
k <- combs[2,best.aligned.pos]
}
B[[j]] <- sort(c(B[[j]],k))
R <- R[-which(R==k)]
for(i in 1:length(B[[j]])){
k <- B[[j]][i]
remove.pos <- apply(combs,2,function(x){any(x==k)})
combs[,remove.pos] <- 0
S[remove.pos] <- Inf
}
}
if(length(R)==0){
K.t <- length(B)
break
}
j <- j+1
}
blocks <- rep(0,K)
for(i in 1:K.t){
blocks[B[[i]]] <- i
}
return(blocks)
}
# function to allocate a color to each block
block.col <- function(x){
cols <- x
color <- 1
for(i in 1:length(unique(x))){
if(length(which(x==i))==1){
cols[x==i]<- 0
} else {
cols[x==i] <- color
color <- color + 1
}
}
return(cols)
}
# function to sort the colors which are allocated to the blocks
sort.colors <- function(x){
if(all(is.na(x)) | all(x==0)){
return(x)
} else {
n <- length(x)
not.zero <- which(x!=0)
y <- x[not.zero]
K.t <- max(x)
if(K.t>1){
y <- -y
m <- length(y)
z <- c()
i <- 1
repeat{
block.i <- which(y==y[1])
z <- c(z,rep(i,length(block.i)))
y <- y[-block.i]
if(length(y)==0) break
i <- i+1
}
x[not.zero] <- z
}
return(x)
}
}
#################
# Main function #
#################
calc.steps <- dim.series*(dim.series-1)/2
combs <- combn(1:dim.series,2)
# empty vectors, matrices and lists for results
scarm.width.mat <- scarm.signal.mat <- slope.mat <- noise.sd.mat <- scarm.statistic <- scarm.critval <- acf.lag.one.mat <- matrix(NA, ncol=dim.series, nrow=length.series)
blocks <- block.colors <- matrix(NA, ncol=dim.series, nrow=length.series)
signal.mat <- width.mat <- matrix(NA, ncol=dim.series, nrow=length.series)
width.mat[n.min-1,] <- n.min-1
S.est <- S.est2 <- matrix(NA, ncol=calc.steps, nrow=length.series)
Cov.est <- Corr.est <- NULL
K.t <- rep(NA, length.series)
width.mat[n.min-1,] <- n.min-1
# matrix of the pairwise slopes, obtained by "get.B.t"; required for the very first signal estimation
B.matrices <- NULL
for(k in 1:dim.series){
B.matrices[[k]] <- get.B.t(time.series[1:n.min,k])
}
# begin loop over time points i
for(i in n.min:length.series){
# 1.) apply SCARM for each univariate component of the multivariate time series
for(k in 1:dim.series){
# compute autocorrelation at lag 1 for each variable at time point i
phi.est.sample.data <- time.series[(max(i-n.max+1,1)):i,k]
phi.est.sample.signal <- signal.mat[(max(i-n.max+1,1)):i,k]
phi.est.sample <- na.omit(phi.est.sample.data-phi.est.sample.signal)
if(length(phi.est.sample)>=n.min){
phi.est <- acf(phi.est.sample, plot=F)$acf[2]
if(is.na(phi.est)){
phi.est <- 0
}
} else {
phi.est <- 0
}
if(is.na(slope.mat[i-1,k])){
n <- min(width.mat[i-1,k]+1, ell.min+r)
} else {
n <- min(width.mat[i-1,k]+1,n.max)
}
y <- time.series[(i-n+1):i,k]
B <- B.matrices[[k]]
if(all(is.na(y[(n-n.min+1):n])) || length(which(!is.na(y[(n-(min(r,n))+1):n]))) < n.min){
# there are not enough recent observations!
mu.est <- NA
beta.est <- NA
noise.sd.est <- NA
# Update step
m <- min(n,ell.min+r)
#adapted.width[i,k] <- m
if(i != length.series){
b.new <- (time.series[i+1,k] - time.series[i+c((-m + 1):0),k]) / (i+1 - (i+c((-m + 1):0)))
B <- B[(ncol(B)-m+1):ncol(B),(ncol(B)-m+1):ncol(B)]
B <- cbind(rbind(B, b.new), c(b.new, NA))
#y <- time.series[(i - m + 1):(i + 1)]
#n <- m + 1
}
} else {
if(n < ell.min+r){
# n is not large enough for test application!
noise.sd.est <- get.Q.adj(y, phi.est=phi.est)
beta.est <- get.beta.RM(B)
mu.est <- get.mu.RM(y=y, beta.RM=beta.est)
if(rtr==TRUE){
rtr.sample <- time.series[(i-n.min+1):i,k]
min.rtr.sample <- min(rtr.sample, na.rm=T)
max.rtr.sample <- max(rtr.sample, na.rm=T)
mu.est <- max(min.rtr.sample, min(mu.est, max.rtr.sample))
}
# Update step
if(i != length.series){
if(n == n.max){
# is window width maximal?
b.new <- (time.series[i+1,k] - time.series[i+c((-n + 2):0),k]) / (i+1 - (i+c((-n + 2):0)))
B <- cbind(rbind(B[-1, -1], b.new), c(b.new,NA))
} else {
b.new <- (time.series[i+1,k] - time.series[i+c((-n + 1):0),k]) / (i+1 - (i+c((-n + 1):0)))
B <- cbind(rbind(B, b.new), c(b.new, NA))
}
}
} else {
# is n large enough for test application?
ell <- n - r
if(length(which(!is.na(y[1:ell])))<ell.min/2 || length(which(!is.na(y[(ell+1):n])))<r/2){
# there are not enough observations for testing!
B <- B[(ncol(B)-ell.min-r+1):ncol(B),(ncol(B)-ell.min-r+1):ncol(B)]
noise.sd.est <- get.Q.adj(y, phi.est=phi.est)
beta.est <- get.beta.RM(B)
mu.est <- get.mu.RM(beta.RM=beta.est, y=y[(n-ell.min-r):n])
n <- ell.min+r
if(rtr==TRUE){
rtr.sample <- time.series[(i-n.min+1):i,k]
min.rtr.sample <- min(rtr.sample, na.rm=T)
max.rtr.sample <- max(rtr.sample, na.rm=T)
mu.est <- max(min.rtr.sample, min(mu.est, max.rtr.sample))
}
if(i != length.series){
b.new <- (time.series[i+1,k] - time.series[i+c((-n + 1):0),k]) / (i+1 - (i+c((-n + 1):0)))
B <- cbind(rbind(B, b.new), c(b.new,NA))
}
} else {
# are there enough observations for testing?
noise.sd.est <- get.Q.adj(y, phi.est=phi.est)
scarm.critval[i,k] <- get.critval(ell=ell, r=r, sign.level=sign.level)
scarm.test <- scarm.test.statistic(y=y, B=B, ell=ell, r=r, noise.sd=noise.sd.est, phi.est=phi.est)
scarm.statistic[i,k] <- scarm.test$test.statistic
slope.diff <- scarm.test$beta.left - scarm.test$beta.right
if(abs(scarm.statistic[i,k]) > scarm.critval[i,k]){
# does test reject the null hypothesis?
B <- B[(ncol(B)-n.min+1):ncol(B),(ncol(B)-n.min+1):ncol(B)]
beta.est <- get.beta.RM(B)
mu.est <- get.mu.RM(beta.RM=beta.est, y=y[(n-n.min+1):n])
n <- n.min
if(rtr==TRUE){
rtr.sample <- time.series[(i-n.min+1):i,k]
min.rtr.sample <- min(rtr.sample, na.rm=T)
max.rtr.sample <- max(rtr.sample, na.rm=T)
mu.est <- max(min.rtr.sample, min(mu.est, max.rtr.sample))
}
if(i != length.series) {
b.new <- (time.series[i+1,k] - time.series[i+c((-n + 1):0),k]) / (i+1 - (i+c((-n + 1):0)))
B <- cbind(rbind(B, b.new), c(b.new,NA))
}
} else {
# difference of RM slopes is NOT too large => do not decrease time window
beta.est <- get.beta.RM(B)
mu.est <- get.mu.RM(beta.RM=beta.est, y=y)
if(i != length.series){
if(n == n.max){
# is window width maximal?
b.new <- (time.series[i+1,k] - time.series[i+c((-n + 2):0),k]) / (i+1 - (i+c((-n + 2):0)))
B <- cbind(rbind(B[-1, -1], b.new), c(b.new,NA))
} else {
b.new <- (time.series[i+1,k] - time.series[i+c((-n + 1):0),k]) / (i+1 - (i+c((-n + 1):0)))
B <- cbind(rbind(B, b.new), c(b.new, NA))
}
}
if(rtr==TRUE){
rtr.sample <- time.series[(i-n.min+1):i,k]
min.rtr.sample <- min(rtr.sample, na.rm=T)
max.rtr.sample <- max(rtr.sample, na.rm=T)
mu.est <- max(min.rtr.sample, min(mu.est, max.rtr.sample))
}
}
}
}
}
# outputs of the SCARM
B.matrices[[k]] <- B
scarm.width.mat[i,k] <- n
scarm.signal.mat[i,k] <- mu.est
slope.mat[i,k] <- beta.est
noise.sd.mat[i,k] <- noise.sd.est
acf.lag.one.mat[i,k] <- phi.est
width.mat[i,] <- scarm.width.mat[i,]
}
# treatment of missing values
non.NA.pos <- which(!is.na(slope.mat[i,]))
dim.i <- length(non.NA.pos)
if(dim.i==1){
signal.mat[i,non.NA.pos] <- scarm.signal.mat[i,non.NA.pos]
width.mat[i,non.NA.pos] <- scarm.width.mat[i,non.NA.pos]
}
if(dim.i > 1){
replace.pos <- apply(combs, 2 , function(x){any(x[1]==non.NA.pos) & any(x[2]==non.NA.pos)})
scarm.width.i <- scarm.width.mat[i,non.NA.pos]
slope.i <- slope.mat[i,non.NA.pos]
noise.sd.i <- noise.sd.mat[i,non.NA.pos]
combs.i <- combn(1:dim.i,2)
calc.steps.i <- dim.i*(dim.i-1)/2
korr.beta <- corr.ogk.est <- rep(NA,calc.steps.i)
# 2.) obtain the comparison statistics in "S.est" for each time point i;
# use SCARM-outputs "slope.mat", "noise.sd.mat" and "scarm.width.mat" from step 1.)
# compute the ration "width.ratio" for each pair of components (smaller width is divided by larger width)
max.w <- apply(combs.i, 2, function(x){apply(cbind(scarm.width.i[x[1]],scarm.width.i[x[2]]),1,max)})
min.w <- apply(combs.i, 2, function(x){apply(cbind(scarm.width.i[x[1]],scarm.width.i[x[2]]),1,min)})
width.ratio <- min.w/max.w
# estimate the matrix "Cov.est" and "Corr.est" of the error covariances and error correlations
Cov.est[[i]] <- Corr.est[[i]] <- matrix(NA,dim.series,dim.series)
n <- max(scarm.width.i)
sample.matrix <- time.series[(i-n+1):i,non.NA.pos]
if(any(is.na(sample.matrix))){
noise.sd.est <- as.vector(noise.sd.mat[i,non.NA.pos])
RM.slopes <- as.vector(slope.mat[i,non.NA.pos])
RM.levels <- as.vector(scarm.signal.mat[i,non.NA.pos])
RM.lines <- matrix(NA, nrow=n, ncol=length(non.NA.pos))
for(k in 1:length(non.NA.pos)){
RM.lines[,k] <- (RM.levels[k] + (RM.slopes[k]*((1-n):0)))
na.pos <- which(is.na(sample.matrix[,k]))
#set.seed(i*k)
sample.matrix[na.pos,k] <- RM.lines[na.pos,k]# + rnorm(n=length(na.pos), sd=noise.sd.est[k])
}
}
Ogk.est <- get.Ogk(sample.matrix, noise.sd.est.s=noise.sd.mat[i,non.NA.pos])
Cov.est[[i]][non.NA.pos,non.NA.pos] <- Ogk.est
Corr.est[[i]][non.NA.pos,non.NA.pos] <- cov2cor(Ogk.est)
# estimate the correlation "korr.beta" between each pair of RM-slopes
# by using the estimated error correlations and the model "corr.RM.slope.model"
for(m in 1:calc.steps.i){
korr.beta[m] <- corr.RM.slope.model(ratio=width.ratio[m], rho=cov2cor(Ogk.est)[combs.i[1,m],combs.i[2,m]])
}
# estimate the variances of the RM-slopes ("var.slopes")
# by using the estimated error variance and the model "var.RM.slope.model"
var.slopes <- apply(matrix(1:dim.i), 1, function(x){noise.sd.i[x]^2 * var.RM.slope.model(scarm.width.i[x])})
# estimate the variance of the difference between each pair of RM-slopes
var.slope.diffs <- apply(combs.i, 2, function(x){var.slopes[x[1]] + var.slopes[x[2]] - 2*(korr.beta[which(combs.i[1,]==x[1]&combs.i[2,]==x[2])]*sqrt(var.slopes[x[1]])*sqrt(var.slopes[x[2]]))})
# estimate the comparison statistics "S.est"
slope.diffs <- apply(combs.i, 2, function(x){slope.i[x[1]] - slope.i[x[2]]})
S.est.i <- slope.diffs/sqrt(var.slope.diffs)
S.est[i,replace.pos] <- S.est.i
if(r.bound > 0){
r.bound.Index.no <- width.ratio < r.bound
r.bound.Index.yes <- width.ratio >= r.bound
S.est.i2 <- (c.bound+1)*r.bound.Index.no + r.bound.Index.yes*S.est.i
} else {
S.est.i2 <- S.est.i
}
S.est2[i,replace.pos] <- S.est.i2
# 3.) build blocks based on "S.est"
blocks[i,] <- build.blocks(S.est2[i,])
block.colors[i,] <- block.col(blocks[i,])
# 4.) if a block contains two or more components, the multivariate TRM-LS is applied for signal estimation
# using the estimated error covariance matrix "Cov.est" from step 2.)
# blocks and number of blocks at time point i
B <- blocks[i,]
K.t[i] <- length(unique(B))
# loop over all blocks
for(j in 1:K.t[i]){
block.comps <- which(B==j)
block.size <- length(block.comps)
if(block.size==1){
signal.mat[i,block.comps] <- scarm.signal.mat[i,block.comps]
width.mat[i,block.comps] <- scarm.width.mat[i,block.comps]
}
if(block.size>1){
# obtain the RM-residuals within the window of width "block.width"
indiv.widths <- as.vector(scarm.width.mat[i,block.comps])
block.width <- min(indiv.widths)
RM.slopes <- as.vector(slope.mat[i,block.comps])
RM.levels <- as.vector(scarm.signal.mat[i,block.comps])
RM.lines <- matrix(NA, nrow=block.width, ncol=block.size)
for(k in 1:block.size){
n <- indiv.widths[k]
RM.lines[,k] <- (RM.levels[k] + (RM.slopes[k]*((1-n):0)))[(n-block.width+1):n]
}
block.sample <- time.series[(i-block.width+1):i,block.comps]
block.sample[is.na(block.sample)] <- RM.lines[is.na(block.sample)]
res <- block.sample - RM.lines
# detect outlying observation vectors
Cov.est.j <- Cov.est[[i]][block.comps,block.comps]
mahala.dist <- mahalanobis(res, center=rep(0,block.size), cov=Cov.est.j, tol=2.220446e-16)
cutoff <- qchisq(p=0.99, df=block.size)*median(mahala.dist)/qchisq(p=0.5, df=block.size)
I.j <- which(mahala.dist <= cutoff)
if(length(I.j) < block.width/2){
I.j <- 1:block.width
}
# LS-Regression based on the trimmed window sample
design.points <- 1:block.width
LS.cov <- cov(cbind(design.points,block.sample)[I.j,])
LS.slopes <- LS.cov[1,2:(block.size+1)] / LS.cov[1,1]
LS.levels <- apply(block.sample[I.j,], 2, mean) - (LS.slopes * mean(design.points))
LS.signals <- LS.levels + (LS.slopes*block.width)
# apply restrict-to-range-rule
if(rtr==TRUE){
for(k in 1:block.size){
LS.signal <- LS.signals[k]
rtr.sample <- time.series[(i-n.min+1):i,block.comps[k]]
min.rtr <- min(rtr.sample, na.rm=T)
max.rtr <- max(rtr.sample, na.rm=T)
LS.signals[k] <- min(max(LS.signal, min.rtr),max.rtr)
}
}
# store TRM-LS signal estimations and window widths
signal.mat[i,block.comps] <- LS.signals
width.mat[i,block.comps] <- block.width
}
}
# the block building may induce that some individual window widths are decreased;
# in this case the matrices of the pairwise slopes have to be adjusted
if(any(width.mat[i,]!=scarm.width.mat[i,])){
new.widths <- which(width.mat[i,]!=scarm.width.mat[i,])
for(k in 1:length(new.widths)){
B <- B.matrices[[new.widths[k]]]
old.dim.B <- dim(B)[1]
new.dim.B <- width.mat[i,new.widths[k]]+1
B.matrices[[new.widths[k]]] <- B[(old.dim.B-new.dim.B+1):old.dim.B,(old.dim.B-new.dim.B+1):old.dim.B]
}
}
}
} # end of loop over time points i
width.mat[n.min-1,] <- NA
block.colors <- t(apply(block.colors, 1, sort.colors))
###########
# Outputs #
###########
result <- list(signal.est=signal.mat,
slope.est=slope.mat,
adapted.width=width.mat,
noise.sd.est=noise.sd.mat,
#
scarm.signal.est=scarm.signal.mat,
scarm.width=scarm.width.mat,
scarm.statistic=scarm.statistic,
scarm.critval=scarm.critval,
#
ssm.statistic=S.est,
blocks=blocks,
acf.lag.one=acf.lag.one.mat,
#
dim.series=dim.series,
length.series=length.series,
S.est2=S.est2,
block.colors=block.colors,
# all (maybe adjusted) inputs are given as outputs
time.series=time.series,
right.width=r,
min.left.width=ell.min,
min.width=n.min,
max.width=n.max,
sign.level=sign.level,
bound.noise.sd=bound.noise.sd,
rtr=rtr,
autocorrelations=autocorrelations,
c.bound=c.bound,
r.bound=r.bound
)
return(structure(result, class = "mscarm.filter"))
}
##################
# Default output #
##################
print.mscarm.filter <- function(x, ...){
N <- dim(x$signal.est)[1]
cat('$signal.est \n')
if(N <= 100){
print(x$signal.est, ...)
} else {
print(x$signal.est[1:100,])
cat('Only the first 100 signal estimations are printed.\n')
}
}
########
# Plot #
########
plot.mscarm.filter <- function(x, info=FALSE,...){
if(info==TRUE){
par(mfrow=c(2,1))
} else {
par(mfrow=c(1,1))
}
y.lim.min <- min(x$time.series,na.rm=T)
y.lim.max <- max(x$time.series,na.rm=T)
y.lim.diff <- y.lim.max-y.lim.min
plot(x$time.series[,1], ylim=c(y.lim.min-(y.lim.diff/5),y.lim.max+(y.lim.diff/5)), ylab="", main="Multivariate time series and mSCARM signal extraction", col="grey", xlab="t", type="l")
points(x$time.series[,1], pch=20, cex=0.5, col="grey")
lines(x$signal.est[,1])
legend("topleft", c("Data", "Signal extraction"), col=c("grey","black"), lty=c(1,1), lwd=c(1,1), bty="n")
#if(info==TRUE){
# legend("topright","Detected signal changes", col=1, pch=4, lty=NA, lwd=2, bty="n")
#}
for(i in 1:x$dim.series){
text(paste("V",i,sep=""), x=1, y=na.omit(x$time.series[,i])[1], pos=2)
}
for(j in 2:x$dim.series){
lines(x$time.series[,j], col="grey")
points(x$time.series[,j], pch=20, cex=0.5, col="grey")
lines(x$signal.est[,j])
}
if(info==TRUE){
#for(i in 1:x$dim.series){
# breaks.i <- which(abs(x$scarm.statistic)[,i] > x$scarm.critval[,i]) - floor(x$min.width/3)
# points(x=breaks.i, y=x$signal.est[breaks.i,i], pch=4, cex=1, lwd=2)
#}
comparisons <- (x$dim.series)*(x$dim.series-1)/2
col.names <- c()
j <- 1
for(i in seq(1,comparisons*2,2)){
col.names[j] <- paste("V",combn(x$dim.series,2)[i],"-V",combn(x$dim.series,2)[i+1], sep="")
j <- j+1
}
plot(rep(NA,x$length.series), type="l", ylim=c(1,comparisons+1), xlab="t", ylab="", yaxt="n",
main=paste("Pairwise slope comparison between the variables V1,...,V",x$dim.series, sep=""))
axis(2, at=1:comparisons, labels=col.names, las=1)
for(i in 1:comparisons){
lines(x=1:(x$length.series),y=rep(i,x$length.series), col="lightgrey")
crosses <- rep(NA,x$length.series)
crosses[which(abs(x$S.est2[,i]) <= x$c.bound)] <- i
points(crosses, pch=3)
}
legend("topleft", "Similar slope detected", lty=NA, lwd=NA, pch=3, bty="n", col=1)
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.