Nothing
##### TIMSAC78 #####
blocar <- function (y, max.order = NULL, span, plot = TRUE)
{
n <- length(y)
if (is.null(max.order))
max.order <- as.integer(2*sqrt(n))
morder <- max.order
if (span < 0)
span <- n
ns <- as.integer((n-morder+span-1)/span)
z <- .Fortran(C_blocarf,
as.double(y),
as.integer(n),
as.integer(morder),
as.integer(span),
as.integer(ns),
mean = double(1),
var = double(1),
aic = double(ns*ns),
bw = double(ns*ns),
b = double(ns*morder),
a = double(ns*morder),
v = double(ns),
ks = integer(ns),
ke = integer(ns),
pxx = double(ns*121))
aic <- list()
bweight <- list()
pacoef <- list()
arcoef <- list()
aic[[1]] <- NA
bweight[[1]] <- NA
for(i in 2:ns) {
j <- ((i-1)*ns+1):((i-1)*ns+i)
aic[[i]] <- z$aic[j]
bweight[[i]] <- z$bw[j]
}
for(i in 1:ns)
pacoef[[i]] <- z$b[((i-1)*morder+1):(i*morder)]
for(i in 1:ns)
arcoef[[i]] <- z$a[((i-1)*morder+1):(i*morder)]
pspec <- array(z$pxx, dim = c(121, ns))
if (plot == TRUE) {
x <- rep(0,121)
for(i in 1:121)
x[i] <- (i - 1) / 240
ymin <- min(pspec)
ymax <- max(pspec)
par(mfrow = c(ns, 1))
nk <- ns
if (ns == 4) {
par(mfrow = c(2, 2))
nk <- 4
}
if (ns > 4) {
par(mfrow = c(3, 2))
nk <- 6
}
if (ns > 6) {
par(mfrow = c(3, 3))
nk <- 9
}
nn <- 0
for(i in 1:ns) {
if (nn == nk) {
par(ask = TRUE)
nn <- 0
}
plot(x, pspec[,i], type = "l", ylim = c(ymin,ymax),
main = paste("y(", z$ks[i], "),...,y(", z$ke[i], ")"),
xlab = "Frequency", ylab = "Power Spectrum")
nn <- nn + 1
}
par(mfrow = c(1, 1))
}
blocar.out <- list(mean = z$mean, var = z$var, aic = aic,
bweight = bweight, pacoef = pacoef, arcoef = arcoef,
v = z$v, init = z$ks, end = z$ke, pspec = pspec)
class(blocar.out) <- "blocar"
return(blocar.out)
}
print.blocar <- function(x, ...)
{
cat(sprintf("\n Mean\t%f\n", x$mean))
cat(sprintf(" Variance\t%f\n", x$var))
ns <- length(x$bweight)
for(i in 1:ns) {
if (i != 1) {
cat("\n\nAR-Model fitted to\t! Bayesian weights\t")
cat("! AIC with respect to the present data\n")
cat("----------------------------------------")
cat("----------------------------------------\n")
cat(sprintf("Current block\t\t! %f\t\t! %f\n", x$bweight[[i]][1],
x$aic[[i]][1]))
for(k in 1:(i-1))
cat(sprintf("%i period former block\t! %f\t\t! %f\n", k,
x$bweight[[i]][k+1], x$aic[[i]][k+1]))
}
morder <- length(x$arcoef[[i]])
cat("\n\n.......... Current model (Average by the Bayesian weight)")
cat(" ..........\n\n")
cat(sprintf(" This model was fitted to the data y( %i ),...,y( %i )\n",
x$init[i], x$end[i]))
cat(sprintf(" Innovation variance = %e\n\n", x$v[i]))
cat(sprintf(" AR coefficients ( order %i ) \n", morder ))
print(x$arcoef[[i]])
}
cat("\n")
}
blomar <- function (y, max.order = NULL, span)
{
n <- nrow(y)
d <- ncol(y)
# upper limit of the order of AR-model
if (is.null(max.order))
max.order <- as.integer(2*sqrt(n))
icflag <- 0
if (max.order > n/(2*d))
icflag <- -1
max.order <- as.integer(min(max.order, n/(2*d)))
morder <- max.order
calb <- rep(1, d) # calibration constant for channel j (j=1,d)
if (span < 1)
span <- n
ns <- as.integer((n-morder+span-1)/span)
z <- .Fortran(C_blomarf,
as.double(y),
as.integer(n),
as.integer(d),
as.double(calb),
as.integer(morder),
as.integer(span),
as.integer(ns),
mean = double(d),
var = double(d),
bw = double(ns*ns),
raic = double(ns*ns),
a = double(d*d*morder*ns),
e = double(d*d*ns),
aic = double(ns),
ks = integer(ns),
ke = integer(ns),
nns = integer(1))
nns <- z$nns
bw <- array(z$bw, dim = c(ns, ns))
bweight <- list()
bweight[[1]] <- NA
if (nns > 1)
for(i in 2:nns)
bweight[[i]] <- bw[1:i, i]
raic <- array(z$raic, dim = c(ns, ns))
aic <- list()
aic[[1]] <- NA
if (nns > 1)
for(i in 2:nns)
aic[[i]] <- raic[1:i, i]
a <- array(z$a, dim = c(d, d, morder, ns))
arcoef <- list()
for(i in 1:nns)
arcoef[[i]] <- array(a[, , , i], dim = c(d, d, morder))
e <- array(z$e, dim = c(d, d, ns))
v <- list()
for(i in 1:nns)
v[[i]] <- array(e[, , i], dim = c(d, d))
aicbay <- z$aic[1:nns]
start <- z$ks[1:nns]
end <- z$ke[1:nns]
blomar.out <- list(mean = z$mean, var = z$var, bweight = bweight,
aic = aic, arcoef = arcoef, v = v, eaic = aicbay,
init = start, end = end)
class(blomar.out) <- "blomar"
if (icflag == -1)
warning(gettextf("max.order is corrected n/(2*d) = %d\n\n", max.order),
domain=NA)
return(blomar.out)
}
print.blomar <- function(x, ...)
{
id <- length(x$mean)
cat("\n\n Mean ")
for(i in 1:id)
cat(sprintf(" %f", x$mean[i]))
cat("\n Variance ")
for(i in 1:id)
cat(sprintf(" %f", x$var[i]))
cat("\n")
ns <- length(x$bweight)
for(i in 1:ns) {
if (i != 1) {
cat("\n\nAR-model fitted to\t! Bayesian weights\t")
cat("! AIC with respect to the present data\n")
cat("----------------------------------------")
cat("----------------------------------------\n")
cat(sprintf("Current block\t\t! %f\t\t! %f\n", x$bweight[[i]][1],
x$aic[[i]][1]))
for(k in 1:(i-1))
cat(sprintf("%i Period former block\t! %f\t\t! %f\n", k,
x$bweight[[i]][k+1], x$aic[[i]][k+1]))
}
mf <- dim(x$arcoef[[i]])[3]
cat("\n\n.......... Current model ..........\n\n")
cat(sprintf(" This model was fitted to the data y( %i ),...,y( %i )\n",
x$init[i], x$end[i]))
cat(sprintf(" aic = %f\n", x$eaic[i]))
cat("\n Innovation variance matrix\n")
print(x$v[[i]])
cat(sprintf("\n AR coefficient matrix ( order %i )\n", mf))
print(x$arcoef[[i]])
}
}
bsubst <- function (y, mtype, lag = NULL, nreg = NULL, reg = NULL,
term.lag = NULL, cstep = 5, plot = TRUE)
{
if (mtype<0 || mtype>4)
stop("allowed model type are
1 : autoregressive model,
2 : polynomial type non-linear model (lag's read in),
3 : polynomial type non-linear model (lag's automatically set),
4 : AR-model with polynomial mean value function\n")
# specification of number of regressors (mtype=1,2,4)
if (is.null(nreg))
if (mtype != 3)
stop("number of regressors is not specified")
# specification of maximum time lag (mtype=2)
if (is.null(reg))
if (mtype == 2)
stop("regressors are not specified")
# specification of regressor (mtype=3)
if (is.null(term.lag))
if (mtype == 3)
stop("maximum time lags are not specified")
n <- length(y)
if (is.null(lag))
lag <- as.integer(2*sqrt(n)) # maximum time lag use in the model
# number of regressors
if (mtype == 3) {
k <- 0
lag5 <- term.lag[5]
for(i in 1:lag5)
k <- k + i * (i + 1) / 2 - 1
nreg <- term.lag[1] + term.lag[2] + term.lag[3] + term.lag[4] + k
}
if (is.null(reg))
reg <- array(0, dim = c(3, nreg))
if (is.null(term.lag))
term.lag <- rep(0, 5)
# f <- "" # specification of regressor (mtype=5)
# cnst <- 0 # constant value (mtype=6)
z <- .Fortran(C_bsubstf,
as.double(y),
as.integer(n),
as.integer(mtype),
as.integer(lag),
as.integer(nreg),
as.integer(cstep),
as.integer(reg),
as.integer(term.lag),
ymean = double(1),
yvar = double(1),
m = integer(1),
aicm = double(1),
vm = double(1),
a1 = double(nreg),
v = double(nreg+1),
aic = double(nreg+1),
daic = double(nreg+1),
aicb = double(1),
vb = double(1),
ek = double(1),
a2 = double(nreg),
ind = integer(nreg),
c = double(nreg),
c1 = double(nreg+1),
c2 = double(nreg),
b = double(nreg),
eicmin = double(1),
esum = double(nreg+1),
npm = double(1),
npmnreg = double(1),
e = double(n*cstep),
mean = double(cstep),
var = double(cstep),
skew = double(cstep),
peak = double(cstep),
cov = double(101),
pxx = double(121))
mo <- z$m
arcoefm=z$a1[1:mo]
if (mtype != 4) {
nps <- lag + 1
perr <- array(z$e, dim = c(n, cstep))
perr <- perr[nps:n, ]
lagh <- 100
nn <- n - nps - 1
if ((lagh > nn) || (lagh == nn))
lagh <- nn - 1
lag1 <- lagh + 1
autcor1=z$cov[1:lag1]
}
if (plot == TRUE) {
if (mtype != 4) {
mm <- 2
if (mtype == 1)
mm <- mm + 1
nc <- (mm + cstep + 1) / 2
par(mfcol = c(nc, 2))
for(i in 1:cstep) {
sig <- sqrt(z$var[i]) * 0.5
hist(perr[,i], breaks = c(-sig*10:1, 0, sig*1:10),
main = "", xlab = paste(i, "-step ahead prediction error"))
}
}
if (mtype == 1) {
k <- lag
k1 <- k + 1
plot((0:k), z$daic[1:k1], ylim = c(0,40), type = "l", xlab = "Lag",
ylab = "AIC(m)-aicmin")
abline(h=0, lty=1)
} else {
plot((0:nreg), z$daic, ylim = c(0,40), type = "l", xlab = "Lag",
ylab = "AIC(m)-aicmin")
abline(h=0, lty=1)
}
if (mtype == 1) {
x <- rep(0, 121)
for(i in 1:121)
x[i] <- (i - 1) / 240
plot(x, z$pxx, type = "l", xlab = "Frequency",
ylab = paste("Power Spectrum"))
}
if (mtype != 4) {
plot(order <- c(0:lagh), autcor1, ylim = c(-1, 1), bty = "l",
type = "h",
main = "Autocorrelation of\n1-step ahead prediction error",
xlab = "Lines show +/-2sd\n ( sd = sqrt(1/n) )",
ylab = "peautcor")
abline(h = 0, lty = 1)
abline(h = 2*sqrt(1/(n-lag)), lty = 3)
abline(h = -2*sqrt(1/(n-lag)), lty = 3)
}
par(mfrow = c(1, 1))
}
if (mtype == 1) {
k <- lag
k1 <- k + 1
bsubst.out <- list(ymean = z$ymean, yvar = z$yvar, v = z$v[1:k1],
aic = z$aic[1:k1], aicmin = z$aicm, daic = z$daic[1:k1],
order.maice = mo, v.maice = z$vm, arcoef.maice = arcoefm,
v.bay = z$vb, aic.bay = z$aicb, np.bay = z$ek,
arcoef.bay = z$a2[1:k], ind.c = z$ind[1:k],
parcor2 = z$c[1:k], damp = z$c1[1:k1],
bweight = z$c2[1:k], parcor.bay = z$b[1:k],
eicmin = z$eicmin, esum = z$esum[1:k1], npmean = z$npm,
npmean.nreg = z$npmnreg, perr = perr, mean = z$mean,
var = z$var, skew = z$skew, peak = z$peak,
peautcor = autcor1, pspec = z$pxx)
}
if (mtype==2 || mtype==3)
bsubst.out <- list(ymean = z$ymean, yvar = z$yvar, v = z$v,
aic = z$aic, aicmin = z$aicm, daic = z$daic,
order.maice = mo, v.maice = z$vm, arcoef.maice = arcoefm,
v.bay = z$vb, aic.bay = z$aicb, np.bay = z$ek,
arcoef.bay = z$a2, ind.c = z$ind, parcor2 = z$c,
damp = z$c1, bweight = z$c2, parcor.bay = z$b,
eicmin = z$eicmin, esum = z$esum, npmean = z$npm,
npmean.nreg = z$npmnreg, perr = perr, mean = z$mean,
var = z$var, skew = z$skew, peak = z$peak,
peautcor = autcor1)
if (mtype == 4)
bsubst.out <- list(ymean = z$ymean, yvar = z$yvar, v = z$v,
aic = z$aic, aicmin = z$aicm, daic = z$daic,
order.maice = mo, v.maice = z$vm, arcoef.maice = arcoefm,
v.bay = z$vb, aic.bay = z$aicb, np.bay = z$ek,
arcoef.bay = z$a2, ind.c = z$ind, parcor2 = z$c,
damp = z$c1, bweight = z$c2, parcor.bay = z$b,
eicmin = z$eicmin, esum = z$esum, npmean = z$npm,
npmean.nreg = z$npmnreg)
return(bsubst.out)
}
exsar <- function (y, max.order = NULL, plot = FALSE)
{
n <- length(y)
if (is.null(max.order))
max.order <- as.integer(2*sqrt(n))
morder <- max.order
morder1 <- morder + 1
z <- .Fortran(C_exsarf,
as.double(y),
as.integer(n),
as.integer(morder),
mean = double(1),
var = double(1),
v = double(morder1),
aic = double(morder1),
daic = double(morder1),
m = integer(1),
aicm = double(1),
sdm1 = double(1),
a1 = double(morder),
sdm2 = double(1),
a2 = double(morder),
ier = integer(1))
if (z$ier != 0)
stop("in FUNCT : SD is less than or equal to 0")
if (plot == TRUE) {
plot((0:morder), z$daic, ylim = c(0,40), type = "l", xlab = "Lag",
ylab = "AIC(m)-aicmin (Truncated at 40.0)")
abline(h = 0, lty = 1)
}
m <- z$m
exsar.out <- list(mean = z$mean, var = z$var, v = z$v, aic = z$aic,
aicmin = z$aicm, daic = z$daic, order.maice = m,
v.maice = z$sdm1, arcoef.maice = z$a1[1:m],
v.mle = z$sdm2, arcoef.mle = z$a2[1:m])
return(exsar.out)
}
mlocar <- function (y, max.order = NULL, span, const = 0, plot = TRUE)
{
n <- length(y)
if (is.null(max.order))
max.order <- as.integer(2*sqrt(n))
morder <- max.order
if (span < 1)
span <- n
ns <- as.integer((n-morder+span-1)/span)
k <- morder + const
z <- .Fortran(C_mlocarf,
as.double(y),
as.integer(n),
as.integer(morder),
as.integer(span),
as.integer(const),
as.integer(ns),
mean = double(1),
var = double(1),
a = double(k*ns),
mf = integer(ns),
sdf = double(ns),
ks = integer(ns),
ke = integer(ns),
pxx = double(121*ns),
ld1 = integer(ns),
ld2 = integer(ns),
ms = integer(ns),
sdms = double(ns),
aics = double(ns),
mp = integer(ns),
sdmp = double(ns),
aicp = double(ns))
a <- array(z$a, dim = c(k, ns))
mf <- z$mf
mj.max <- 0
for(i in 1:ns)
if (mf[i] != 0) {
mj.max <- i
}
ns <- mj.max
arcoef <- list()
for(i in 1:ns) {
arcoef[[i]] <- a[1:mf[i], i]
}
pspec <- array(z$pxx, dim=c(121,ns))
pspec <- pspec[, 1:ns]
npre <- z$ld1
order.const <- z$mp
v.const <- z$sdmp
aic.const <- z$aicp
npre[1] <- NA
order.const[1] <- NA
v.const[1] <- NA
aic.const[1] <- NA
if (plot == TRUE) {
oldpar <- par(no.readonly = TRUE)
x <- rep(0, 121)
for(i in 1:121)
x[i] <- (i - 1) / 240
ymin <- min(pspec)
ymax <- max(pspec)
par(mfrow=c(ns, 1))
nk <- ns
if (ns == 4) {
par(mfrow = c(2, 2))
nk <- 4
}
if (ns > 4) {
par(mfrow = c(3, 2))
nk <- 6
}
if (ns > 6) {
par(mfrow = c(3, 3))
nk <- 9
}
nn <- 0
for(i in 1:ns) {
if (nn == nk) {
par(ask = TRUE)
nn <- 0
}
plot(x, pspec[, i], type = "l", ylim = c(ymin, ymax),
main = paste("y(", z$ks[i], "),...,y(", z$ke[i], ")"),
xlab = "Frequency", ylab = "Power Spectrum")
nn <- nn + 1
}
par(oldpar)
}
mlocar.out <- list(mean = z$mean, var = z$var, ns = ns, order = mf,
arcoef = arcoef, v = z$sdf, init = z$ks,
end = z$ke, pspec = pspec, npre = npre,
nnew = z$ld2, order.mov = z$ms, v.mov = z$sdms,
aic.mov = z$aics, order.const = order.const,
v.const = v.const, aic.const = aic.const)
class(mlocar.out) <- "mlocar"
return(mlocar.out)
}
print.mlocar <- function(x, ...)
{
cat(sprintf("\n\n Mean\t%f\n", x$mean))
cat(sprintf(" Variance\t%f\n", x$var))
ns <- x$ns
for(i in 1:ns) {
if (i == 1) {
cat(sprintf("\n Initial local model: (nnew = %i)", x$nnew[i]))
cat(sprintf("\tvariance = %e\taic = %f\n", x$v.mov[i],
x$aic.mov[i]))
}
if (i != 1) {
cat("\n\n >>> The following two models are compared <<<\n")
np <- x$npre[i] + x$nnew[i]
cat(sprintf(" Moving model: (npre = %i, nnew = %i)\t", x$npre[i],
x$nnew[i]))
cat(sprintf("variance = %e\taic = %f\n", x$v.mov[i], x$aic.mov[i]))
cat(sprintf(" Constant model: (npre+nnew = %i)\tvariance = %e",
np, x$v.const[i]))
cat(sprintf("\taic = %f\n", x$aic.const[i]))
if (x$aic.mov[i] < x$aic.const[i]) {
cat(" ---> New model adopted\n")
} else {
cat(" ---> Constant model adopted\n")
}
}
cat("\n\n.......... Current model ..........\n\n")
cat(sprintf(" This model was fitted to the data y( %i ),...,y( %i )\n",
x$init[i], x$end[i]))
cat(sprintf(" Innovation variance = %e\n\n", x$v[i]))
cat(sprintf(" AR coefficients ( order %i ) \n", x$order[i]))
print(x$arcoef[[i]])
}
}
mlomar <- function (y, max.order = NULL, span, const = 0)
{
n <- nrow(y)
d <- ncol(y)
if (is.null(max.order))
max.order <- as.integer(2*sqrt(n))
icflag <- 0
if (max.order > n/(2*d))
icflag <- -1
max.order <- as.integer(min(max.order, n/(2*d)))
morder <- max.order
calb <- rep(1, d) # calibration for channel j (j=1,d)
if (span < 1)
span <- n
ns <- as.integer((n-morder+span-1)/span)
z <- .Fortran(C_mlomarf,
as.double(y),
as.integer(n),
as.integer(d),
as.double(calb),
as.integer(morder),
as.integer(span),
as.integer(const),
as.integer(ns),
mean = double(d),
var = double(d),
ld1 = integer(ns),
ld2 = integer(ns),
ms = integer(ns),
aicm = double(ns),
mp = integer(ns),
aicc = double(ns),
mf = integer(ns),
aic = double(ns),
a = double(d*d*morder*ns),
e = double(d*d*ns),
ks = integer(ns),
ke = integer(ns),
nns = integer(1))
nns <- z$nns
npre <- z$ld1[1:nns]
nnew <- z$ld2[1:nns]
order.mov <- z$ms[1:nns]
aic.mov <- z$aicm[1:nns]
order.const <- z$mp[1:nns]
aic.const <- z$aicc[1:nns]
order <- z$mf[1:nns]
aic <- z$aic[1:nns]
start <- z$ks[1:nns]
end <- z$ke[1:nns]
npre[1] <- NA
order.mov[1] <- NA
aic.mov[1] <- NA
order.const[1] <- NA
aic.const[1] <- NA
a <- array(z$a, dim = c(d, d, morder, ns))
arcoef <- list()
for(i in 1:nns)
arcoef[[i]] <- array(a[, , , i], dim = c(d, d, z$mf[i]))
e <- array(z$e, dim = c(d, d, ns))
v <- list()
for(i in 1:nns)
v[[i]] <- array(e[, , i], dim = c(d, d))
mlomar.out <- list(mean = z$mean, var = z$var, ns = nns, order = order,
aic = aic, arcoef = arcoef, v = v, init = start,
end = end, npre = npre, nnew = nnew,
order.mov = order.mov, aic.mov = aic.mov,
order.const = order.const, aic.const = aic.const)
class(mlomar.out) <- "mlomar"
if (icflag == -1)
warning(gettextf("max.order is corrected n/(2*d) = %d\n\n", max.order),
domain=NA)
return(mlomar.out)
}
print.mlomar <- function(x, ...)
{
id <- length(x$mean)
cat("\n\n Mean ")
for(i in 1:id)
cat(sprintf(" %f", x$mean[i]))
cat("\n Variance ")
for(i in 1:id)
cat(sprintf(" %f", x$var[i]))
ns <- x$ns
for(i in 1:ns) {
if (i == 1)
cat(sprintf("\n\n Initial local model: (nnew = %i), aic = %f\n",
x$nnew[i], x$aic[i]))
if (i != 1) {
cat("\n\n >>> The following two models are compared <<<\n")
np <- x$npre[i] + x$nnew[i]
cat(sprintf(" Moving model: (npre = %i, nnew = %i)\taic = %f\n",
x$npre[i], x$nnew[i], x$aic.mov[i]))
cat(sprintf(" Constant model: (npre+nnew = %i)\taic = %f\n",
np,x$aic.const[i]))
if (x$aic.mov[i] < x$aic.const[i]) {
cat(" ---> New model adopted \n")
} else {
cat(" ---> Constant model adopted \n")
}
}
cat("\n\n.......... Current model ..........\n\n")
cat(sprintf(" This model was fitted to the data y( %i ),...,y( %i )\n",
x$init[i], x$end[i]))
cat(sprintf(" aic = %f\n", x$aic[i]))
cat("\n Innovation variance matrix\n")
print(x$v[[i]])
cat(sprintf("\n AR coefficient matrix ( order %i )\n", x$order[i]))
print(x$arcoef[[i]])
}
}
mulbar <- function (y, max.order = NULL, plot = FALSE)
{
n <- nrow(y)
d <- ncol(y)
calb <- rep(1, d) # calibration of channel i (i=1,d)
if (is.null(max.order)) # upper limit of the order of AR-model
max.order <- as.integer(2*sqrt(n))
icflag <- 0
if (max.order > n/(2*d))
icflag <- -1
max.order <- as.integer(min(max.order, n/(2*d)))
morder <- max.order
morder1 <- morder + 1
d2 <- d * d
z <- .Fortran(C_mulbarf,
as.double(y),
as.integer(n),
as.integer(d),
as.double(calb),
as.integer(morder),
mean = double(d),
var = double(d),
v = double(morder1),
aic = double(morder1),
daic = double(morder1),
m = integer(1),
aicm = double(1),
vm = double(1),
w1 = double(morder1),
w2 = double(morder),
a = double(d2*morder),
b = double(d2*morder),
g = double(d2*morder),
h = double(d2*morder),
e = double(d2),
aicb = double(1))
if (plot == TRUE) {
plot((0:morder), z$daic, ylim = c(0,40), type = "l", xlab = "Lag",
ylab = "AIC - aicmin (Truncated at 40.0)")
abline(h = 0, lty = 1)
}
mulbar.out <- list(mean = z$mean, var = z$var, v = z$v,
aic = z$aic, aicmin = z$aicm, daic = z$daic,
order.maice = z$m, v.maice = z$vm,
bweight = z$w1, integra.bweight = z$w2,
arcoef.for = array(z$a, dim = c(d, d, morder)),
arcoef.back = array(z$b, dim = c(d, d, morder)),
pacoef.for = array(z$g, dim = c(d, d, morder)),
pacoef.back = array(z$h, dim = c(d, d, morder)),
v.bay = array(z$e, dim = c(d, d)),
aic.bay = z$aicb)
if (icflag == -1)
warning(gettextf("max.order is corrected n/(2*d) = %d\n\n", max.order),
domain=NA)
return(mulbar.out)
}
mulmar <- function (y, max.order = NULL, plot = FALSE)
{
n <- nrow(y)
d <- ncol(y)
calb <- rep(1, d)
if (is.null(max.order))
max.order <- as.integer(2*sqrt(n))
icflag <- 0
if (max.order > n/(2*d))
icflag <- -1
max.order <- as.integer(min(max.order, n/(2*d)))
lag <- max.order
lag1 <- lag + 1
z <- .Fortran(C_mulmarf,
as.double(y),
as.integer(n),
as.integer(d),
as.double(calb),
as.integer(max.order),
mean = double(d),
var = double(d),
v = double(lag1*d),
aic = double(lag1*d),
daic = double(lag1*d),
m = integer(d),
aicm = double(d),
vm = double(d),
npr = integer(d),
jnd = integer(lag1*d*d),
a = double(lag1*d*d),
rv = double(d),
aicf = double(d),
ei = double(d*d),
bi = double(d*d*lag),
matv = double(d*d),
arcoef = double(d*d*lag),
morder = integer(1),
aics = double(1))
v <- list()
aic <- list()
daic <- list()
for(i in 1:d) {
j <- ((i - 1) * lag1 + 1):(i * lag1)
v[[i]] <- z$v[j]
aic[[i]] <- z$aic[j]
daic[[i]] <- z$daic[j]
}
jnd <- list()
subregcoef <- list()
ind <- array(z$jnd, dim = c(lag1*d, d))
a <- array(z$a, dim = c(lag1*d, d))
for(i in 1:d)
jnd[[i]] <- ind[(1:z$npr[[i]]), i]
for(i in 1:d)
subregcoef[[i]] <- a[(1:z$npr[[i]]), i]
if (plot == TRUE) {
par(mfrow = c(d, 1))
for(i in 1:d) {
plot((0:max.order), daic[[i]], ylim = c(0, 40), type = "l",
main = paste(" d=", i), xlab = "Lag",
ylab = "AIC(m)-aicmin (Truncated at 40.0)")
abline(h = 0, lty = 1)
}
par(mfrow = c(1, 1))
}
morder <- z$morder
mulmar.out <- list(mean = z$mean, var = z$var, v = v, aic = aic,
aicmin = z$aicm, daic = daic, order.maice = z$m,
v.maice = z$vm, np = z$npr, jnd = jnd,
subregcoef = subregcoef, rvar = z$rv, aicf = z$aicf,
respns = array(z$ei, dim = c(d, d)),
regcoef = array(z$bi, dim = c(d, d, morder)),
matv = array(z$matv, dim = c(d, d)), morder = morder,
arcoef = array(z$arcoef, dim = c(d, d, morder)),
aicsum = z$aics)
if (icflag == -1)
warning(gettextf("max.order is corrected n/(2*d) = %d\n\n", max.order),
domain=NA)
return(mulmar.out)
}
perars <- function (y, ni, lag = NULL, ksw = 0)
{
n <- length(y)
if (is.null(lag))
lag <- as.integer(2*sqrt(ni))
lag1 <- lag + 1
k <- (lag1*ni + ksw)*ni
z <- .Fortran(C_perarsf,
as.double(y),
as.integer(n),
as.integer(ni),
as.integer(lag),
as.integer(ksw),
mean = double(1),
var = double(1),
np = integer(ni),
jnd = integer(k),
a = double(k),
aic = double(ni),
b = double(ni*ni*lag),
v = double(ni*ni),
c = double(ni),
osd = double(ni),
mo = integer(1))
npara <- z$np
ind <- array(z$jnd, dim = c(lag1*ni+ksw, ni))
jnd <- list()
for(i in 1:ni)
jnd[[i]] <- ind[1:npara[i], i]
a <- array(z$a, dim = c(lag1*ni+ksw, ni))
regcoef <- list()
for(i in 1:ni)
regcoef[[i]] <- a[1:npara[i], i]
perars.out <- list(mean = z$mean, var = z$var, subset = jnd,
regcoef = regcoef, rvar = z$osd, np = npara,
aic = z$aic, v = array(z$v, dim = c(ni, ni)),
arcoef = array(z$b, dim = c(ni, ni, z$morder)),
const = z$c, morder = z$mo)
class(perars.out) <- "perars"
return(perars.out)
}
print.perars <- function(x, ...)
{
cat(sprintf("\n\n Mean = %f\n", x$mean))
cat(sprintf(" Variance = %f\n", x$var))
ni <- nrow(x$v)
for(i in 1:ni) {
cat("\n ........ Regression Model for the regressand")
cat(sprintf(" i = %i ........\n", i))
cat("\n subset")
sorder <- length(x$subset[[i]])
for(j in 1:sorder)
cat(sprintf("\t\t%i", x$subset[[i]][j]))
cat("\n regcoef")
for(j in 1:sorder)
cat(sprintf("\t%f", x$regcoef[[i]][j]))
cat(sprintf("\n\n residual variance (rvar) = %f\n", x$rvar[i]))
cat(sprintf(" number of parameter (np) = %i\n", x$np[i]))
cat(sprintf(" AIC = n*log(rvar) + 2*np = %f\n\n", x$aic[i]))
}
cat("\n\nmatrix of regression coefficients\n")
print(x$v)
cat("\n\nregression coefficients within the present period\n")
print(x$arcoef)
cat("\n\nconstants within the regression models\n")
print(x$const)
}
unibar <- function (y, ar.order = NULL, plot = TRUE)
{
n <- length(y)
if (is.null(ar.order))
ar.order <- as.integer(2*sqrt(n))
ar.order1 <- ar.order + 1
z <- .Fortran(C_unibarf,
as.double(y),
as.integer(n),
as.integer(ar.order),
mean = double(1),
var = double(1),
v = double(ar.order1),
aic = double(ar.order1),
daic = double(ar.order1),
m = integer(1),
aicm = double(1),
vm = double(1),
pa = double(ar.order),
bw = double(ar.order1),
sbw = double(ar.order),
pab = double(ar.order),
aicb = double(1),
vb = double(1),
pn = double(1),
a = double(ar.order),
pxx = double(121))
if (plot == TRUE) {
oldpar <- par(no.readonly = TRUE)
par(mfrow = c(3, 1))
plot((0:ar.order), z$daic, ylim = c(0,40), type = "l", xlab = "Lag",
ylab = "AIC(m)-aicmin (Truncated at 40.0)")
abline(h = 0, lty = 1)
plot(z$pa, type = "h", xlab = "Lag", ylab = "Partial autocorrelation")
abline(h = 0, lty = 1)
# abline(h = 1/sqrt(n), lty = 3)
# abline(h = -1/sqrt(n),lty = 3)
abline(h = 2/sqrt(n), lty = 3)
abline(h = -2/sqrt(n), lty = 3)
x <- rep(0, 121)
for(i in 1:121)
x[i] <- (i - 1) / 240
plot(x, z$pxx, type = "l", xlab = "Frequency",
ylab = "Power Spectral Density")
par(oldpar)
unibar.out <- list(mean = z$mean, var = z$var, v = z$v,
aic = z$aic, aicmin = z$aicm, order.maice = z$m,
v.maice = z$vm, bweight = z$bw[2:(ar.order1)],
integra.bweight = z$sbw, v.bay = z$vb,
aic.bay = z$aicb, np = z$pn,
pacoef.bay = z$pab, arcoef = z$a)
} else {
unibar.out <- list(mean = z$mean, var = z$var, v = z$v,
aic = z$aic, aicmin = z$aicm, daic = z$daic,
order.maice = z$m, v.maice = z$vm,
pacoef = z$pa, bweight = z$bw[2:(ar.order1)],
integra.bweight = z$sbw, v.bay = z$vb,
aic.bay = z$aicb, np = z$pn,
pacoef.bay = z$pab, arcoef = z$a,
pspec = z$pxx)
}
return(unibar.out)
}
unimar <- function (y, max.order = NULL, plot = FALSE)
{
n <- length(y)
if (is.null(max.order))
max.order <- as.integer(2*sqrt(n)) # upper limit of AR-order
morder <- max.order
morder1 <- morder + 1
z <- .Fortran(C_unimarf,
as.double(y),
as.integer(n),
as.integer(morder),
mean = double(1),
var = double(1),
v = double(morder1),
aic = double(morder1),
daic = double(morder1),
m = integer(1),
aicm = double(1),
vm = double(1),
a = double(morder))
m <- z$m
if (plot == TRUE) {
plot((0:morder), z$aicm, ylim = c(0,40), type = "l", xlab = "Lag",
ylab = "AIC(m)-aicmin (Truncated at 40.0)")
abline(h = 0, lty = 1)
unimar.out <- list(mean = z$mean, var = z$var, v = z$v,
aic = z$aic, aicmin = z$aicm, order.maice = m,
v.maice = z$vm, arcoef = z$a[1:m])
} else {
unimar.out <- list(mean = z$mean, var = z$var, v = z$v,
aic = z$aic, aicmin = z$aicm, daic = z$daic,
order.maice = m, v.maice = z$vm,
arcoef = z$a[1:m])
}
return(unimar.out)
}
xsarma <- function (y, arcoefi, macoefi)
{
n <- length(y)
arcoefi <- -arcoefi # Initial estimates of AR coefficients
macoefi <- -macoefi # Initial estimates of MA coefficients
p <- length(arcoefi) # AR-ORDER
q <- length(macoefi) # MA-ORDER
p01 <- c(arcoefi, macoefi) # INITIAL ESTIMATES OF AR- AND MA-COEFFICIENTS
z <- .Fortran(C_xsarmaf,
as.double(y),
as.integer(n),
as.integer(p),
as.integer(q),
as.double(p01),
g1 = double(p+q),
tl1 = double(1),
p02 = double(p+q),
g2 = double(p+q),
alpar = double(p),
alpma = double(q),
tl2 = double(1),
sigma2 = double(1))
coef <- -z$p02
xsarma.out <- list(gradi = z$g1, lkhoodi = z$tl1, arcoef = coef[1:p],
macoef = coef[(p+1):(p+q)], grad = z$g2,
alph.ar = z$alpar, alph.ma = z$alpma, lkhood = z$tl2,
wnoise.var = z$sigma2)
return(xsarma.out)
}
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.