Nothing
#############################################################
# #
# wle.ar.ao function #
# Author: Claudio Agostinelli #
# E-mail: claudio@unive.it #
# Date: December, 20, 2003 #
# Version: 0.1-3 #
# #
# Copyright (C) 2003 Claudio Agostinelli #
# #
#############################################################
wle.ar.ao <- function(x, x.init, x.seasonal.init, coef, ncoef, ncoef.seasonal, period, sigma2, xreg=NULL, raf, smooth, w.level, verbose=FALSE, ao.list, population.size=20, population.choose=5, elements.random=4, num.max, approx.w=TRUE) {
nused <- length(x)
xx <- wle.ar.matrix(x=x, x.init=x.init, x.seasonal.init=x.seasonal.init, ncoef=ncoef, ncoef.seasonal=ncoef.seasonal, period=period, xreg=xreg)
resid <- x - c(xx%*%coef)
ww <- weights <- .Fortran("wlew",
as.double(resid),
as.integer(nused),
as.double(resid),
as.integer(nused),
as.integer(raf),
as.double(smooth),
as.double(sigma2),
totweight=double(1),
weights=double(nused),
PACKAGE="wle")$weights
ao.position <- 0
pos.temp <- 1:nused
pos.temp <- pos.temp[order(weights)]
weights.sort <- sort(weights)
ao.temp <- weights.sort <= w.level
pos.temp <- pos.temp[ao.temp]
ao <- rep(FALSE, nused)
if (length(pos.temp)) {
pos.temp <- pos.temp[1:min(length(pos.temp),num.max)]
ao[pos.temp] <- TRUE
}
pos <- (1:nused)[ao]
if (verbose) {
cat("We have the following observations under the w.level=",w.level,":\n",pos,"\n")
}
if (any(ao)) {
model.in <- vector(length=0)
for (i in 1:length(ao.list)) {
if (all(is.element(ao.list[[i]], pos))) {
temp <- vector(length=0)
for (j in 1:length(ao.list[[i]])) {
temp <- c(temp,(1:length(pos))[pos==ao.list[[i]][j]])
}
model.in <- c(model.in, sum(2^(temp-1)))
}
}
num.model <- max(length(model.in),population.size)
num.pos <- (2^sum(ao))-1
dim.dim <- floor(log(num.pos,2))+1
w.tilde <- rep(0,num.model)
model.in <- c(model.in, sample(x=(1:num.pos), size=(num.model-length(model.in)), replace=TRUE))
for (isearch in 1:num.model) {
pos.ao <- sort(pos[binary(model.in[isearch],dim.dim)$dicotomy])
num.ao <- length(pos.ao)
x.ao <- x
xx.ao <- xx
for (t in pos.ao) {
if (ncoef) {
x.ao[t] <- xx.ao[t,]%*%coef
for (tt in 1:ncoef) {
if ((t+tt)<=nused) {
xx.ao[t+tt,tt] <- x.ao[t]
}
}
}
if (ncoef.seasonal) {
for (tt in 1:ncoef.seasonal) {
if ((t+tt*period)<=nused) {
xx.ao[t+tt*period,tt+ncoef] <- x.ao[t]
}
}
}
}
resid.ao <- x - c(xx.ao%*%coef)
resid.ao <- resid.ao[-pos.ao]
if (approx.w) {
weights.temp <- approx(x=resid, y=ww, xout=resid.ao)$y
weights.temp[is.na(weights.temp)] <- 0
weights.temp[weights.temp > 1] <- 1
weights.temp[weights.temp < 0] <- 0
} else {
weights.temp <- wle.weights(x=resid.ao, smooth=smooth, sigma2=sigma2, raf=raf, location=TRUE)$weights
}
w.tilde[isearch] <- sum(weights.temp)/nused
}
model.in <- model.in[order(w.tilde)]
w.tilde <- sort(w.tilde)
while ((model.in[1]-model.in[num.model])!=0) {
num.model.sel <- population.choose
cum.wtilde <- cumsum(w.tilde)[num.model.sel:num.model]
pos.child <- vector(length=0)
while (length(pos.child)==0) {
temp <- runif(2,0,cum.wtilde[length(cum.wtilde)])
pos.aaa <- min((num.model.sel:num.model)[cum.wtilde > temp[1]])
pos.bbb <- min((num.model.sel:num.model)[cum.wtilde > temp[2]])
pos.aa <- pos[binary(model.in[pos.aaa],dim.dim)$dicotomy]
pos.bb <- pos[binary(model.in[pos.bbb],dim.dim)$dicotomy]
pos.child <- c(pos.aa,pos.bb,pos[sample(x=(1:length(pos)), size=elements.random, replace=TRUE)])
pos.child <- pos.child[as.logical(sample(x=c(0,1), size=length(pos.child), replace=TRUE))]
pos.child <- sort(unique(pos.child))
}
temp.child <- vector(length=0)
for (i in 1:length(pos.child)) {
temp.child <- c(temp.child,(1:length(pos))[pos==pos.child[i]])
}
model.child <- sum(2^(temp.child-1))
num.child <- length(pos.child)
x.ao <- x
xx.ao <- xx
for (t in pos.child) {
x.ao[t] <- xx.ao[t,]%*%coef
if (ncoef) {
for (tt in 1:ncoef) {
if ((t+tt)<=nused) {
xx.ao[t+tt,tt] <- x.ao[t]
}
}
}
if (ncoef.seasonal) {
for (tt in 1:ncoef.seasonal) {
if ((t+tt*period)<=nused) {
xx.ao[t+tt*period,tt+ncoef] <- x.ao[t]
}
}
}
}
resid.ao <- x - c(xx.ao%*%coef)
resid.ao <- resid.ao[-pos.child]
if (approx.w) {
weights.temp <- approx(x=resid, y=ww, xout=resid.ao)$y
weights.temp[is.na(weights.temp)] <- 0
weights.temp[weights.temp > 1] <- 1
weights.temp[weights.temp < 0] <- 0
} else {
weights.temp <- wle.weights(x=resid.ao, smooth=smooth, sigma2=sigma2, raf=raf, location=TRUE)$weights
}
w.tilde.child <- sum(weights.temp)/nused
w.tilde <- c(w.tilde,w.tilde.child)
model.in <- c(model.in,model.child)
model.in <- model.in[order(w.tilde)][-1]
w.tilde <- sort(w.tilde)[-1]
}
if (max(w.tilde)<(sum(weights)/nused)) {
ao.position <- NULL
} else {
ao.position <- sort(pos[binary(model.in[1],dim.dim)$dicotomy])
}
} else {
ao.position <- NULL
}
x.ao <- x
xx.ao <- xx
for (t in ao.position) {
x.ao[t] <- xx.ao[t,]%*%coef
if (ncoef) {
for (tt in 1:ncoef) {
if ((t+tt)<=nused) {
xx.ao[t+tt,tt] <- x.ao[t]
}
}
}
if (ncoef.seasonal) {
for (tt in 1:ncoef.seasonal) {
if ((t+tt*period)<=nused) {
xx.ao[t+tt*period,tt+ncoef] <- x.ao[t]
}
}
}
}
resid.ao <- x - c(xx.ao%*%coef)
w.temp <- wle.weights(x=resid.ao, smooth=smooth, sigma2=sigma2, raf=raf, location=TRUE)
resid.ao <- resid.ao - w.temp$location
if (verbose) {
cat("Additive outliers: \n", ao.position, "\n")
}
return(list(x.ao=x.ao, resid.ao=resid.ao, ao.position=ao.position))
}
#############################################################
# #
# wle.ar function #
# Author: Claudio Agostinelli #
# E-mail: claudio@unive.it #
# Date: January, 5, 2004 #
# Version: 0.1-5 #
# #
# Copyright (C) 2004 Claudio Agostinelli #
# #
#############################################################
wle.ar <- function(x, order=c(1,0), seasonal=list(order=c(0,0), period=NA), group, group.start, group.step=group.start, xreg=NULL, include.mean=TRUE, na.action=na.fail, tol=10^(-6), tol.step=tol, equal=10^(-3), equal.step=equal, raf="HD", smooth=0.0031, smooth.ao=smooth, boot=10, boot.start=10, boot.step=boot.start, num.sol=1, x.init=0, x.seasonal.init=0, max.iter.out=20, max.iter.in=50, max.iter.start=200, max.iter.step=500, verbose=FALSE, w.level=0.4, min.weights=0.5, population.size=10, population.choose=5, elements.random=2, wle.start=FALSE, init.values=NULL, num.max=NULL, num.sol.step=2, approx.w=TRUE) {
if (length(order)!=2) stop("order must have two components")
if (length(seasonal$order)!=2) stop("seasonal$order must have two components")
ncoef <- order[1]
ncoef.seasonal <- seasonal$order[1]
raf <- switch(raf,
HD = 1,
NED = 2,
SCHI2 = 3,
-1)
if (raf==-1) stop("Please, choose the RAF: HD=Hellinger Disparity, NED=Negative Exponential Disparity, SCHI2=Symmetric Chi-squares Disparity")
result <- list()
series <- deparse(substitute(x))
if(NCOL(x) > 1)
stop("only implemented for univariate time series")
x <- na.action(as.ts(x))
n <- length(x)
if(is.null(seasonal$period) || is.na(seasonal$period)
|| seasonal$period == 0) seasonal$period <- frequency(x)
period <- seasonal$period
arma <- c(ncoef,0,ncoef.seasonal,0,period,order[2],seasonal$order[2])
if (!(length(x.init)==1 | length(x.init)==ncoef)) stop("x.init must have one or order[1] elements\n")
if (!(length(x.seasonal.init)==1 | length(x.seasonal.init)==ncoef.seasonal*period)) stop("x.seasonal.init must have one or seasonal$order[1]*period elements\n")
d <- order[2]
d.s <- seasonal$order[2]
#if (d | d.s) {
# if (any(.packages(all.available=TRUE)=="ts")) {
# library(ts)
# } else {
# stop("For Integrated model you need function diff in package ts")
# }
#}
if(d) x <- diff(x, 1, d)
if(d.s) x <- diff(x, seasonal$period, d.s)
xtsp <- tsp(x)
tsp(x) <- NULL
nd <- d+d.s
nused <- length(x)
if(is.null(xreg)) {
ncxreg <- 0
} else {
if(NROW(xreg) != nused) stop("lengths of x and xreg do not match")
ncxreg <- NCOL(xreg)
}
class(xreg) <- NULL
if (include.mean && (nd==0)) {
if (is.matrix(xreg) && is.null(colnames(xreg)))
colnames(xreg) <- paste("xreg", 1:ncxreg, sep="")
xreg <- cbind(intercept=rep(1, nused), xreg=xreg)
ncxreg <- ncxreg + 1
}
if (ncxreg) {
xreg <- as.matrix(xreg)
if (qr(xreg)$rank < ncol(xreg)) stop("xreg is collinear")
}
if (ncxreg) {
if (d) xreg <- diff(xreg, 1, d)
if (d.s) xreg <- diff(xreg, seasonal$period, d.s)
xreg <- as.matrix(xreg)
if (qr(xreg)$rank < ncol(xreg)) stop("xreg is collinear")
}
if (is.null(num.max)) num.max <- nused
# start bootstrap iteration
first.time <- TRUE
iboot <- 1
tot.sol <- 0
not.conv <- 0
while(iboot<=boot & tot.sol<num.sol) {
pos.iboot <- round(runif(1,(group+1),nused))
x.boot <- x[(pos.iboot-group+1):pos.iboot]
if (is.matrix(xreg))
xreg.boot <- xreg[(pos.iboot-group+1):pos.iboot,]
else
xreg.boot <- xreg[(pos.iboot-group+1):pos.iboot]
if (!is.null(init.values)) {
temp <- list(conv=TRUE)
temp.n <- length(init.values)
temp$coef <- init.values[1:(temp.n-1)]
temp$sigma2 <- init.values[temp.n]
} else {
temp <- wle.ar.start(x=x.boot, x.init=x.init, x.seasonal.init=x.seasonal.init, ncoef=ncoef, ncoef.seasonal=ncoef.seasonal, period=period, xreg=xreg.boot, raf=raf, smooth=smooth, group=group.start, boot=boot.start, max.iter=max.iter.start, wle.start=wle.start, verbose=verbose)
}
if (temp$conv) {
coef.init <- temp$coef
sigma2.init <- temp$sigma2
xx <- wle.ar.matrix(x=x, x.init=x.init, x.seasonal.init=x.seasonal.init, ncoef=ncoef, ncoef.seasonal=ncoef.seasonal, period=period, xreg=xreg)
resid.init <- x - c(xx%*%coef.init)
if (verbose) {
cat("Initial values from the subsample ",iboot,": \n parameters: ", coef.init,"\n sigma2: ",sigma2.init," \n")
}
weights <- .Fortran("wlew",
as.double(resid.init),
as.integer(nused),
as.double(resid.init),
as.integer(nused),
as.integer(raf),
as.double(smooth.ao),
as.double(sigma2.init),
totweights=double(1),
weights=double(nused),
PACKAGE="wle")$weights
if (sum(weights)/nused >= min.weights) {
ao.list <- list(0)
wres <- wle.ar.ao(x=x, x.init=x.init, x.seasonal.init=x.seasonal.init, coef=coef.init, ncoef=ncoef, ncoef.seasonal=ncoef.seasonal, period=period, sigma2=sigma2.init, xreg=xreg, raf=raf, smooth=smooth.ao, w.level=w.level, verbose=verbose, ao.list=ao.list, population.size=population.size, population.choose=population.choose, elements.random=elements.random, num.max=num.max, approx.w=approx.w)
resid.ao <- wres$resid.ao
x.ao <- wres$x.ao
ao.position <- wres$ao.position
if (!is.null(ao.position)) {
ao.list <- c(ao.list,list(ao.position))
}
ao.position.old <- c(ao.position,0)
conv <- TRUE
iter <- 0
while (!setequal(ao.position,ao.position.old) & conv) {
iter <- iter + 1
max.tol <- tol + 1
ao.position.old <- ao.position
iter.int <- 0
while(max.tol>tol & conv) {
iter.int <- iter.int + 1
coef.old <- coef.init
res <- wle.ar.step(coef=coef.init, ncoef=ncoef, ncoef.seasonal=ncoef.seasonal, period=period, x=x, xreg=xreg, x.init=x.init, x.seasonal.init=x.seasonal.init, raf=raf, smooth=smooth, sigma2=sigma2.init, num.sol=num.sol.step, ao.position=ao.position, group=group.step, boot=boot.step, max.iter=max.iter.step, verbose=verbose, tol=tol.step, equal=equal.step)
coef.init <- res$coef
sigma2.init <- res$sigma2
max.tol <- max(abs(coef.old-coef.init))
if (iter.int > max.iter.in) {
if (verbose) cat("Convergence problem: maximum iteration number reached in the inner loop\n")
conv <- FALSE
}
if(any(!is.finite(coef.init)) | any(!is.finite(sigma2.init)) | (sum(res$weights)/nused < min.weights)) {
if (verbose) {
cat("Convergence problem: some values are not finite, bad starting point or sum of the weights less than min.weights\n")
cat("Parameters: ",coef.init,"\n")
cat("Sigma2: ",sigma2.init,"\n")
cat("Sum of weights/size: ",sum(res$weights)/nused,"\n")
}
conv <- FALSE
}
}
# end while(max.tol>tol & conv)
if (conv) {
if(verbose) {
cat("iteration: ",iter," \n parameters: ",coef.init," \n sigma2: ",sigma2.init," \n")
}
wres <- wle.ar.ao(x=x, x.init=x.init, x.seasonal.init=x.seasonal.init, coef=coef.init, ncoef=ncoef, ncoef.seasonal=ncoef.seasonal, period=period, sigma2=sigma2.init, xreg=xreg, raf=raf, smooth=smooth.ao, w.level=w.level, verbose=verbose, ao.list=ao.list, population.size=population.size, population.choose=population.choose, elements.random=elements.random, num.max=num.max, approx.w=approx.w)
resid.ao <- wres$resid.ao
x.ao <- wres$x.ao
ao.position <- wres$ao.position
if (!is.null(ao.position)) {
ao.list <- c(ao.list,list(ao.position))
}
if ((lao <- length(ao.list))>2) {
for (i.ao in 1:(lao-2)) {
if (setequal(ao.list[[lao]],ao.list[[i.ao]])) {
conv <- FALSE
}
}
}
if(iter > max.iter.out) {
if (verbose) cat("Convergence problem: maximum iteration number reached in the outer loop\n")
conv <- FALSE
}
}
}
# end while (!setequal(ao.position,ao.position.old) & conv)
if (conv) {
xx <- wle.ar.matrix(x=x, x.init=x.init, x.seasonal.init=x.seasonal.init, ncoef=ncoef, ncoef.seasonal=ncoef.seasonal, period=period, xreg=xreg)
resid.init <- c(x - c(xx%*%coef.init))
resid.init <- ts(resid.init, start=start(x), end=end(x), frequency=frequency(x))
class(resid.init) <- "ts"
weights.with.ao <- wle.weights(x=resid.init, smooth=smooth, sigma2=res$sigma2, raf=raf, tol=tol, location=TRUE)$weights
resid <- ts(res$resid, start=start(x), end=end(x), frequency=frequency(x))
class(resid) <- "ts"
xx.ao <- wle.ar.matrix(x=x.ao, x.init=x.init, x.seasonal.init=x.seasonal.init, ncoef=ncoef, ncoef.seasonal=ncoef.seasonal, period=period, xreg=xreg)
resid.ao <- c(x - c(xx.ao%*%coef.init))
resid.ao <- ts(resid.ao, start=start(x), end=end(x), frequency=frequency(x))
weights.without.ao <- wle.weights(x=resid.ao, smooth=smooth, sigma2=res$sigma2, raf=raf, tol=tol, location=TRUE)$weights
class(resid.ao) <- "ts"
x.ao <- ts(res$x.ao, start=start(x), end=end(x), frequency=frequency(x))
class(x.ao) <- "ts"
if(first.time) {
coef.final <- res$coef
sigma2.coef.final <- res$sigma2.coef
weights.final <- res$weights
weights.final.with.ao <- weights.with.ao
weights.final.without.ao <- weights.without.ao
sigma2.final <- res$sigma2
resid.final <- resid
resid.ao.final <- resid.ao
x.ao.final <- x.ao
resid.init.final <- resid.init
ao.position.final <- list(wres$ao.position)
first.time <- FALSE
tot.sol <- 1
} else {
if(min(abs(coef.final-res$coef))>equal) {
tot.sol <- tot.sol+1
coef.final <- rbind(coef.final,res$coef)
sigma2.coef.final <- rbind(sigma2.coef.final,res$sigma2.coef)
weights.final <- rbind(weights.final,res$weights)
weights.final.with.ao <- rbind(weights.final.with.ao, weights.with.ao)
weights.final.without.ao <- rbind(weights.final.without.ao, weights.without.ao)
sigma2.final <- c(sigma2.final,res$sigma2)
resid.final <- rbind(resid.final,resid)
ao.position.final <- c(ao.position.final,list(wres$ao.position))
resid.ao.final <- rbind(resid.ao.final,resid.ao)
x.ao.final <- rbind(x.ao.final,x.ao)
resid.init.final <- rbind(resid.init.final,resid.init)
}
}
} else {
not.conv <- not.conv+1
}
} else {
not.conv <- not.conv+1
}
} else {
not.conv <- not.conv+1
}
iboot <- iboot+1
}
# end bootstrap iteration
if(tot.sol==0) {
result$coef <- NULL
result$sigma2.coef <- NULL
result$sigma2 <- NULL
result$arma <- arma
result$resid <- NULL
result$resid.with.ao <- NULL
result$resid.without.ao <- NULL
result$x.ao <- NULL
result$call <- match.call()
result$series <- series
result$weights <- NULL
result$weights.with.ao <- NULL
result$weights.without.ao <- NULL
result$tot.sol <- tot.sol
result$not.conv <- not.conv
result$ao.position <- NULL
} else {
nm <- NULL
if(arma[1] > 0) nm <- c(nm, paste("ar", 1:arma[1], sep=""))
if(arma[2] > 0) nm <- c(nm, paste("ma", 1:arma[2], sep=""))
if(arma[3] > 0) nm <- c(nm, paste("sar", 1:arma[3], sep=""))
if(arma[4] > 0) nm <- c(nm, paste("sma", 1:arma[4], sep=""))
if(ncxreg > 0)
if(!is.null(cn <- colnames(xreg))) nm <- c(nm, cn)
else nm <- c(nm, paste("xreg", 1:ncxreg, sep=""))
if(tot.sol==1) {
names(coef.final) <- nm
names(sigma2.coef.final) <- nm
} else {
colnames(coef.final) <- nm
rownames(coef.final) <- paste("root", 1:tot.sol, sep=" ")
colnames(sigma2.coef.final) <- nm
rownames(sigma2.coef.final) <- paste("root", 1:tot.sol, sep=" ")
rownames(weights.final) <- paste("root", 1:tot.sol, sep=" ")
rownames(resid.init.final) <- paste("root", 1:tot.sol, sep=" ")
rownames(resid.final) <- paste("root", 1:tot.sol, sep=" ")
rownames(resid.ao.final) <- paste("root", 1:tot.sol, sep=" ")
rownames(x.ao.final) <- paste("root", 1:tot.sol, sep=" ")
}
names(arma) <- c("ar", "ma", "sar", "sma", "period", "diff", "sdiff")
result$coef <- coef.final
result$sigma2.coef <- sigma2.coef.final
result$sigma2 <- sigma2.final
result$arma <- arma
result$resid <- resid.final
result$resid.without.ao <- resid.ao.final
result$resid.with.ao <- resid.init.final
result$x <- x
result$x.ao <- x.ao.final
result$call <- match.call()
result$series <- series
result$weights <- weights.final
result$weights.with.ao <- weights.final.with.ao
result$weights.without.ao <- weights.final.without.ao
result$tot.sol <- tot.sol
result$not.conv <- not.conv
result$ao.position <- ao.position.final
result$w.level <- w.level
}
class(result) <- "wle.arima"
return(result)
}
#############################################################
# #
# wle.ar.step function #
# Author: Claudio Agostinelli #
# E-mail: claudio@unive.it #
# Date: December, 30, 2003 #
# Version: 0.1-3 #
# #
# Copyright (C) 2003 Claudio Agostinelli #
# #
#############################################################
wle.ar.step <- function(coef, ncoef, ncoef.seasonal, period, x, xreg, x.init, x.seasonal.init, raf, smooth, sigma2, ao.position, group, boot, max.iter, verbose, num.sol=2, tol, equal) {
result <- list()
nused <- length(x)
if (is.null(xreg)) {
ncxreg <- 0
} else {
ncxreg <- NCOL(xreg)
}
x.ao <- x
xx.ao <- wle.ar.matrix(x=x.ao, x.init=x.init, x.seasonal.init=x.seasonal.init, ncoef=ncoef, ncoef.seasonal=ncoef.seasonal, period=period, xreg=xreg)
for (t in ao.position) {
x.ao[t] <- xx.ao[t,]%*%coef
if (ncoef) {
for (tt in 1:ncoef) {
if ((t+tt)<=nused) {
xx.ao[t+tt,tt] <- x.ao[t]
}
}
}
if (ncoef.seasonal) {
for (tt in 1:ncoef.seasonal) {
if ((t+tt*period)<=nused) {
xx.ao[t+tt*period,tt+ncoef] <- x.ao[t]
}
}
}
}
# print(coef)
# print(cbind(ao.position,x.ao[ao.position]))
if (qr(xx.ao)$rank==NCOL(xx.ao)) {
temp.wle <- wle.lm(x.ao~xx.ao -1, boot=boot, smooth=smooth, num.sol=num.sol, group=group, max.iter=max.iter, tol=tol, equal=equal)
} else {
if (verbose) cat("wle.ar.step: the matrix is not full rank\n")
temp.wle <- list()
temp.wle$tot.sol==0
}
if (temp.wle$tot.sol!=0) {
if (verbose) {
cat("Number of solutions: ",temp.wle$tot.sol," found on ",num.sol,"\n")
cat("Parameters: ",temp.wle$coefficients,"\n")
cat("Sigma2: ",temp.wle$scale^2,"\n")
}
if (temp.wle$tot.sol>1) {
ccc <- c(coef,sigma2)
dist <- rep(0,temp.wle$tot.sol)
for (k in 1:temp.wle$tot.sol) {
dist[k] <- sum((ccc-c(temp.wle$coefficients[k,],temp.wle$scale[k]^2))^2)
}
root <- (1:temp.wle$tot.sol)[dist==min(dist)]
if (verbose) cat("We use root: ",root,"\n")
ttt <- list()
ttt$coefficients <- temp.wle$coefficients[root,]
ttt$residuals <- temp.wle$residuals[root,]
ttt$scale <- temp.wle$scale[root]
ttt$weights <- temp.wle$weights[root,]
temp.wle <- ttt
}
result$coef <- temp.wle$coefficients
result$resid <- temp.wle$residuals
result$sigma2 <- temp.wle$scale^2
result$weights <- temp.wle$weights
result$weights[ao.position] <- 0
result$sigma2.coef <- diag(result$sigma2*solve(t(xx.ao)%*%diag(result$weights)%*%xx.ao, tol=1e-10))
result$resid.ao <- temp.wle$residuals
result$x.ao <- x.ao
result$sol <- 1
} else {
result$sol <- 0
result$coef <- rep(NA,ncol(xx.ao))
result$sigma2 <- NA
}
return(result)
}
#############################################################
# #
# wle.ar.matrix function #
# Author: Claudio Agostinelli #
# E-mail: claudio@unive.it #
# Date: September, 26, 2001 #
# Version: 0.1 #
# #
# Copyright (C) 2001 Claudio Agostinelli #
# #
#############################################################
wle.ar.matrix <- function(x, x.init=0, x.seasonal.init=0, ncoef, ncoef.seasonal, period, xreg=NULL) {
nused <- length(x)
if(is.null(xreg)) {
ncxreg <- 0
xreg <- NULL
} else {
ncxreg <- NCOL(xreg)
}
xx <- vector(length=0)
if (length(x.init)==ncoef) {
x.temp <- c(x.init,x)
} else {
x.temp <- c(rep(x.init,ncoef),x)
}
for (i in 1:ncoef) {
xx <- cbind(xx,x.temp[(ncoef-i+1):(nused+ncoef-i)])
}
if (ncoef.seasonal) {
if (length(x.seasonal.init)==ncoef.seasonal*period) {
x.temp <- c(x.seasonal.init,x)
} else {
x.temp <- c(rep(x.seasonal.init,ncoef.seasonal*period),x)
}
for (i in 1:ncoef.seasonal) {
xx <- cbind(xx,x.temp[((ncoef.seasonal-i)*period+1):(nused+(ncoef.seasonal-i)*period)])
}
}
if (ncxreg) xx <- cbind(xx,xreg)
return(xx)
}
#############################################################
# #
# wle.ar.start function #
# Author: Claudio Agostinelli #
# E-mail: claudio@unive.it #
# Date: September, 26, 2001 #
# Version: 0.1 #
# #
# Copyright (C) 2001 Claudio Agostinelli #
# #
#############################################################
wle.ar.start <- function(x, x.init, x.seasonal.init, ncoef, ncoef.seasonal, period, xreg, raf, smooth, group, boot, max.iter, wle.start, verbose=FALSE) {
result <- list()
nused <- length(x)
xx <- wle.ar.matrix(x=x, x.init=x.init, x.seasonal.init=x.seasonal.init, ncoef=ncoef, ncoef.seasonal=ncoef.seasonal, period=period, xreg=xreg)
if (qr(xx)$rank==NCOL(xx)) {
if (wle.start) {
temp <- wle.lm(x~xx -1, raf=raf, smooth=smooth, group=group, boot=boot, max.iter=max.iter, verbose=verbose, num.sol=1)
tot.sol <- temp$tot.sol
} else {
temp <- lm(x~xx -1)
tot.sol <- 1
}
if (tot.sol) {
resid <- temp$resid
result$coef <- temp$coef
result$resid <- temp$resid
if (wle.start) {
result$sigma2 <- temp$scale^2
} else {
result$sigma2 <- (summary(temp)$sigma)^2
}
result$conv <- TRUE
} else {
result$coef <- rep(NA,ncoef)
result$resid <- rep(NA,nused)
result$sigma2 <- NA
result$conv <- FALSE
}
} else {
result$coef <- rep(NA,ncoef)
result$resid <- rep(NA,nused)
result$sigma2 <- NA
result$conv <- FALSE
}
return(result)
}
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.