Nothing
func.input <- function(x, formula, p, q, f=NULL, eps=1e-5, off.diag=TRUE){
n <- dim(x)[1]; T <- dim(x)[2]; r <- 2*max(p, q)
denom <- res <- tx <- matrix(0, nrow=n, ncol=T-r); h <- x*0
c.mat <- c()
for(i in 1:n){
gf <- garchFit(formula = formula, data = x[i, ], trace=FALSE, include.mean=FALSE)
c <- coef(gf)
#Sanity check to examine if a+b<1. If not then run a simple GARCH(0,1) model which is a EWMA equivalent.
if (sum(c[-1]) >= 1) {
print(paste("Component series",i,"returns a+b > 1: Re-estimating assuming GARCH(0,1)"))
c <- c*0
gf <- garchFit(formula = ~garch(1, 0), data = x[i, ], trace=FALSE, include.mean=FALSE)
c[1:2] <- coef(gf)
}
c.mat <- rbind(c.mat, c)
h[i, ] <- gf@h.t
res[i, ] <- x[i, -(1:r)]/sqrt(h[i, -(1:r)])
}
if(is.null(f)) f <- apply(c.mat, 1, function(z){max(1, min(.99, sum(z[-1]))/max(.01, 1-sum(z[-1])))})
for(t in (r+1):T){
denom[, t-r] <- c.mat[, 1] # apply(c.mat, 1, max)
if(p > 0) denom[, t-r] <- denom[, t-r] + apply(c.mat[, 1+1:p, drop=FALSE]*x[, (t-1):(t-p), drop=FALSE]^2, 1, sum)/f
if(q > 0) denom[, t-r] <- denom[, t-r] + apply(c.mat[, 1+p+1:q, drop=FALSE]*h[, (t-1):(t-q), drop=FALSE], 1, sum)/f
}
denom <- denom + eps*x[, -(1:r)]^2
tx <- x[, -(1:r)]/sqrt(denom)
ttx <- t(apply(tx^2, 1, function(z){z/mean(z)}))
if(off.diag){ sgn <- sign(tx%*%t(tx)); ttx <- rbind(ttx, t(func_input_off(tx, sgn, 0, 0)$input)) } else sgn <- c()
list(denom=denom, tx=tx, ttx=ttx, h=h, res=res, c.mat=c.mat, sgn=sgn, f=f)
}
make.tree <- function(y, op, gam, dw, rule=NULL){
len <- ncol(y)
if(is.null(rule)) rule <- round(log(len, 2)/2)
tree <- list(matrix(0, 6, 1))
mat <- c()
fd <- func_density(y, gam)
stat <- fd$res[, op]
test.stat <- max(stat[-c((1:dw), (len-dw):len)])
hat.chp <- min(which(stat==test.stat))
tree[[1]][1, 1] <- 1
tree[[1]][2, 1] <- 1
tree[[1]][3, 1] <- hat.chp
tree[[1]][4, 1] <- len
tree[[1]][5, 1] <- 0
tree[[1]][6, 1] <- test.stat
mat <- cbind(mat, c(tree[[1]][-5, ], 1, 1))
j <- 1
while(length(tree)==j & j < rule){
npc <- dim(tree[[j]])[2]
if(sum(tree[[j]][4,]-tree[[j]][2,]-rep(4*dw, npc)>0)){
ncc <- 0; i <- 1
while(i <= npc){
if(tree[[j]][3, i]-tree[[j]][2, i]+1>4*dw){
s <- tree[[j]][2, i]; e <- tree[[j]][3, i]
fd <- func_density(y[, s:e], gam)
stat <- fd$res[, op]
test.stat <- max(stat[-c((1:dw), (e-s+1-dw):(e-s+1))])
hat.chp <- s+min(which(stat==test.stat))-1
if(length(tree)==j) tree <- c(tree, list(matrix(0, 6, 0)))
ncc <- ncc+1
tree[[j+1]] <- matrix(c(tree[[j+1]], matrix(0, 6, 1)), 6, ncc)
tree[[j+1]][1, ncc] <- 2*tree[[j]][1, i]-1
tree[[j+1]][2, ncc] <- s
tree[[j+1]][3, ncc] <- hat.chp
tree[[j+1]][4, ncc] <- e
tree[[j+1]][5, ncc] <- 0
tree[[j+1]][6, ncc] <- test.stat
mat <- cbind(mat, c(tree[[j+1]][-5, ncc], j+1, ncc))
}
if(tree[[j]][4, i]-tree[[j]][3, i]>4*dw){
s <- tree[[j]][3, i]+1; e <- tree[[j]][4, i]
fd <- func_density(y[, s:e], gam)
stat <- fd$res[, op]
test.stat <- max(stat[-c((1:dw), (e-s+1-dw):(e-s+1))])
hat.chp <- s+min(which(stat==test.stat))-1
if(length(tree)==j) tree <- c(tree, list(matrix(0, 6, 0)))
ncc <- ncc+1
tree[[j+1]] <- matrix(c(tree[[j+1]], matrix(0, 6, 1)), 6, ncc)
tree[[j+1]][1, ncc] <- 2*tree[[j]][1, i]
tree[[j+1]][2, ncc] <- s
tree[[j+1]][3, ncc] <- hat.chp
tree[[j+1]][4, ncc] <- e
tree[[j+1]][5, ncc] <- 0
tree[[j+1]][6, ncc] <- test.stat
mat <- cbind(mat, c(tree[[j+1]][-5, ncc], j+1, ncc))
}
i <- i+1
}
j <- j+1
} else{
break
}
}
list(tree=tree, mat=mat)
}
tri.kern <- function(len, h){
filter <- rep(0, h+1)
i <- 0
while (i <= h) {
u <- i/h
if (u < 1/2)
filter[i+1] <- 1
if (u >= 1/2 & u < 1)
filter[i+1] <- 2 * (1 - u)
if (u > 1)
break
i <- i + 1
}
filter
}
get.gg <- function(z, M=NULL, C=2, max.K=5){
len <- length(z)
max.K <- max(max.K, sqrt(log(len)))
acv <- acf(z, type="covariance", lag.max=len-1, plot=FALSE)$acf[,,1]
if(is.null(M)){
l <- 1; ind <- 0
while(l < sqrt(len)){
if(abs(acv[l+1])/acv[1] < C*sqrt(log(len)/len)){
ind <- ind+1
} else{
if(ind>0) ind <- 0
}
if(ind==max.K) break
l <- l+1
}
lam <- max(1, l-max.K); M <- 2*lam
}
k <- tri.kern(len, M)
c(acv[1]+2*sum(k[-1]*acv[2:(M+1)]), 2*sum(k[-1]*(1:M)*acv[2:(M+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.