ARMAacf <- function(ar = numeric(0), ma = numeric(0), lag.max = r,
                    pacf = FALSE)
{
    p <- length(ar)
    q <- length(ma)
    if(!p && !q) stop("empty model supplied")
    r <- max(p, q + 1)
    if(p > 0) {
        if(p > 1) {
            x <- c(1, -ar)
            A <- matrix(0, p + 1, 2 * p + 1)
            ind <- as.matrix(expand.grid(1:(p + 1), 1:(p+1)))[, 2:1]
            ind[, 2] <- ind[, 1] + ind[, 2] - 1
            A[ind] <- x
            A[,  1:p] <- A[, 1:p] + A[, (2 * p + 1):(p + 2)]
            rhs <- c(1, rep(0,p))
            if(q > 0) {
                psi <- c(1, ARMAtoMA(ar, ma, q))
                theta <- c(1, ma, rep(0, q+1))
                for(k in 1 + 0:q) rhs[k] <- sum(psi * theta[k + 0:q])
            }
            ind <- (p+1):1
            Acf <- solve(A[ind, ind], rhs)
	    Acf <- Acf[-1]/Acf[1]
        } else Acf <- ar
        if(lag.max > p) {
            xx <- rep(0, lag.max - p)
            Acf <- c(Acf, filter(xx, ar, "recursive", init = Acf))
        }
        Acf <- c(1, Acf)
    } else if(q > 0) {
        x <- c(1, ma)
        Acf <- filter(c(x, rep(0, q)), rev(x), sides=1)[-(1:q)]
        if(lag.max > q) Acf <- c(Acf, rep(0, lag.max - q))
        Acf <- Acf/Acf[1]
    }
    names(Acf) <- 0:max(p, lag.max)
    if(pacf) .C("uni_pacf", as.double(Acf), pacf = double(lag.max),
                as.integer(lag.max), PACKAGE = "ts")$pacf
    else Acf
}

acf2AR <- function(acf)
{
    r <- as.double(drop(acf))
    order.max <- length(r) - 1
    if(order.max <= 0) stop("`acf' must be of length two or more")
    z <- .Fortran("eureka", as.integer(order.max), r, r,
                  coefs = double(order.max^2), vars = double(order.max),
                  double(order.max), PACKAGE = "ts")
    nm <- paste("ar(",1:order.max, ")", sep="")
    matrix(z$coefs, order.max, order.max, dimnames=list(nm, 1:order.max))
}

ARMAtoMA <- function(ar = numeric(0), ma = numeric(0), lag.max)
    .Call("ARMAtoMA", as.double(ar), as.double(ma), as.integer(lag.max),
          PACKAGE ="ts")
HoltWinters <-
function (x,

          # smoothing parameters
          alpha    = NULL, # level
          beta     = NULL, # trend
          gamma    = NULL, # seasonal component
          seasonal = c("additive", "multiplicative"),
          start.periods = 3,

          # starting values
          l.start  = NULL, # level
          b.start  = NULL, # trend
          s.start  = NULL  # seasonal components vector of length `period'
          )
{
    x <- as.ts(x)
    seasonal <- match.arg(seasonal)
    f <- frequency(x)

    if (!is.null(alpha) && alpha==0)
        stop ("cannot fit models without level (alpha must not be 0).")
    if (!all(is.null(c(alpha, beta, gamma))) &&
        any(c(alpha, beta, gamma) < 0 || c(alpha, beta, gamma) > 1))
        stop ("alpha, beta and gamma must be within the unit interval.")
    if ((is.null(gamma) || gamma > 0)) {
        if (seasonal == "multiplicative" && any(x <= 0))
            stop ("data must be strictly non-negative for multiplicative Holt-Winters")
        if (start.periods < 3)
            stop ("Need at least 3 periods to compute seasonal start values")
    }

    ## initialization
    if (!is.null(gamma) && gamma == 0) {
        ## non-seasonal Holt-Winters
        if (is.null(l.start)) l.start <- x[2] #FIXME: use x[1] in case of beta=0?
        if (is.null(b.start))
            if (is.null(beta) || beta > 0) b.start <- x[2] - x[1]
        start.time <- 3
        s.start    <- 0
    } else {
        ## seasonal Holt-Winters
        start.time <- f + 1
        wind       <- start.periods * f

        ## decompose series
        st <- decompose(ts(x[1:wind], start = start(x), frequency = f),
                        seasonal)

        ## level & intercept
        m  <- lm (na.omit(st$trend) ~ c(1:(wind - f + 1)))

        if (is.null(l.start)) l.start <- as.vector(coef(m)[1])
        if (is.null(b.start)) b.start <- as.vector(coef(m)[2])
        if (is.null(s.start)) s.start <- st$figure
    }

    ## Call to filtering loop
    hw <- function(alpha, beta, gamma)
        .C ("HoltWinters",
            as.double(x),
            as.integer(length(x)),
            as.double(alpha),
            as.double(beta),
            as.double(gamma),
            as.integer(start.time),
            as.integer(! + (seasonal == "multiplicative")),
            as.integer(f),

            a = as.double(l.start),
            b = as.double(b.start),
            s = as.double(s.start),
            SSE = as.double(0),
            xhat = double(length(x) - start.time + 1),
            PACKAGE = "ts"
            )

    ## if alpha and/or beta and/or gamma are omitted, use optim to find the
    ## values minimizing the squared prediction error
    if (is.null(gamma)) {
        ## Optimize seasonal parameter
        if(is.null(alpha) & is.null(beta))
        {
            error <- function (p) hw(p[1], p[2], p[3])$SSE
            sol   <- optim(c(0.3, 0.1, 0.1), error, method = "L-BFGS-B",
                           lower = c(0, 0, 0), upper = c(1, 1, 1))
            alpha <- sol$par[1]
            beta  <- sol$par[2]
            gamma <- sol$par[3]
        }
        else if(is.null(alpha) & !is.null(beta))
        {
            error <- function (p) hw(p[1], beta, p[2])$SSE
            sol   <- optim(c(0.3, 0.1), error, method = "L-BFGS-B",
                           lower = c(0, 0), upper = c(1, 1))
            alpha <- sol$par[1]
            gamma <- sol$par[2]
        }
        else if(is.null(beta))
        {
            error <- function (p) hw(alpha, p[1], p[2])$SSE
            sol   <- optim(c(0.1, 0.1), error, method = "L-BFGS-B",
                           lower = c(0, 0), upper = c(1, 1))
            beta  <- sol$par[1]
            gamma <- sol$par[2]
        }
    } else {
        ## Seasonal parameter fixed or 0
        if(is.null(alpha) & is.null(beta))
        {
            error <- function (p) hw(p[1], p[2], gamma)$SSE
            sol   <- optim(c(0.3, 0.1), error, method = "L-BFGS-B",
                           lower = c(0, 0), upper = c(1, 1))
            alpha <- sol$par[1]
            beta  <- sol$par[2]
        }
        else if(is.null(alpha) & !is.null(beta))
        {
            error <- function (p) hw(p, beta, gamma)$SSE
            sol   <- optim(0.3, error, method = "L-BFGS-B",
                           lower = 0, upper = 1)
            alpha <- sol$par
        }
        else if(is.null(beta))
        {
            error <- function (p) hw(alpha, p, gamma)$SSE
            sol   <- optim(0.1, error, method = "L-BFGS-B",
                           lower = 0, upper = 1)
            beta  <- sol$par
        }
    }

    # get (final) results
    final.fit <- hw(alpha, beta, gamma)

    ## return fitted values and estimated coefficients along with parameters used
    structure(list(
                   fitted    = ts(final.fit$xhat,
                                  start = start(lag(x, k = 1 - start.time)),
                                  freq  = frequency(x)),
                   x         = x,
                   alpha     = alpha,
                   beta      = beta,
                   gamma     = gamma,
                   coefficients = c(a = final.fit$a,
                                    b = if (beta > 0) final.fit$b,
                                    s = if (gamma > 0) final.fit$s),
                   seasonal  = seasonal,
                   SSE       = final.fit$SSE,
                   call      = match.call()
                   ),
              class = "HoltWinters"
              )
}

## Predictions, optionally with prediction intervals
predict.HoltWinters <-
    function (object, n.ahead = 1, prediction.interval = FALSE,
              quantile = qnorm(0.975), ...)
{
    f <- frequency(object$x)

    vars <- function(h) {
        psi <- function(j)
            object$alpha * (1 + j * object$beta) +
                (j %% f == 0) * object$gamma * (1 - object$alpha)
        var(residuals(object)) * if (object$seasonal == "additive")
            sum(1, (h > 1) * sapply(1:(h-1), function(j) crossprod(psi(j))))
        else {
            rel <- 1 + (h - 1) %% f
            sum(sapply(0:(h-1), function(j) crossprod (psi(j) * object$coefficients[2 + rel] / object$coefficients[2 + (rel - j) %% f])))
        }
    }

    ## compute predictions
    # level
    fit <- rep(as.vector(object$coefficients[1]),n.ahead)
    # trend
    if (object$beta > 0)
        fit <- fit + as.vector((1:n.ahead)*object$coefficients[2])
        # seasonal component
    if (object$gamma > 0)
        if (object$seasonal == "additive")
            fit <- fit + rep(object$coefficients[-(1:(1+(object$beta>0)))],
                             length.out=length(fit))
        else
            fit <- fit * rep(object$coefficients[-(1:(1+(object$beta>0)))],
                             length.out=length(fit))

    ## compute prediction intervals
    if (prediction.interval) int <- quantile*sqrt(sapply(1:n.ahead,vars))
    ts(
       cbind(fit = fit,
             upr = if(prediction.interval) fit + int,
             lwr = if(prediction.interval) fit - int
             ),
       start = end(lag(object$fitted, k = -1)),
       freq  = frequency(object$fitted)
       )
}

residuals.HoltWinters <- function (object, ...) object$x - object$fitted


plot.HoltWinters <-
    function (x, predicted.values = NA, intervals = TRUE,
              separator = TRUE, col = 1, col.predicted = 2, col.intervals = 4,
              lty.separator = 3, ylab = "Observed / Fitted",
              main = "Holt-Winters filtering", ...)
{
    ## plot fitted/predicted values
    plot(ts(c(x$fitted, if(!is.na(predicted.values)) predicted.values[,1]),
            start = start(fitted(x)),
            frequency = frequency(fitted(x))),
         col = col.predicted,
         ylim = range(na.omit(c(x$fitted,x$x,predicted.values))),
         ylab = ylab, main = main,
         ...
         )

    ## plot prediction interval
    if(!is.na(predicted.values) && intervals && ncol(predicted.values) > 1) {
        lines(predicted.values[,2], col = col.intervals)
        lines(predicted.values[,3], col = col.intervals)
    }

    ## plot observed values
    lines(x$x, col = col)

    ## plot separator
    if (separator && !is.na(predicted.values))
        abline (v = time(x$x)[length(x$x)], lty = lty.separator)
}

## print function
print.HoltWinters <- function (x, ...)
{
    cat ("Holt-Winters exponential smoothing",
         if (x$beta == 0) "without" else "with", "trend and",
         if (x$gamma == 0) "without" else
         paste(if (x$beta==0) "with ", x$seasonal, sep=""),
         "seasonal componenent.\n")
    cat ("\nCall:\n", deparse (x$call), "\n\n")
    cat ("Smoothing parameters:\n")
    cat (" alpha: ", x$alpha, "\n")
    cat (" beta : ", x$beta, "\n")
    cat (" gamma: ", x$gamma, "\n\n")

    cat ("Coefficients:\n")
    print(t(t(x$coefficients)))
}

## decompose additive/multiplicative series into trend/seasonal figures/noise
decompose <- function (x, type = c("additive", "multiplicative"))
{
    type <- match.arg(type)
    l <- length(x)
    f <- frequency(x)
    if (f == 1) stop ("Time series has no period")
    if (l < 3*f) stop ("Need at least 3 periods")

    ## filter out seasonal components
    trend <- filter (x, rep(1, f)/f)

    ## compute seasonal components
    season <- if (type == "additive") x - trend else x / trend

    ## remove incomplete seasons at beginning/end
    season <- window(season, start(x) + c(1, 0), c(end(x)[1] - 1, f))

    ## average seasonal figures
    periods <- l %/% f
    index <- c(0, cumsum(rep (f, periods - 3)))
    figure <- numeric(f)
    for (i in 1:f) figure[i] <- mean(season[index + i])

    ## normalize figure
    figure <- if (type == "additive") figure - mean(figure)
    else figure / mean(figure)

    ## return values
    seasonal <- ts(rep(figure, periods), start = start(x), frequency = f)

    structure(
              list(seasonal = seasonal,
                   trend    = trend,
                   random   = x -
                   if (type == "additive") seasonal + trend
                   else seasonal * trend,
                   figure   = figure,
                   type     = type
                   ),
              class = "decomposed.ts"
              )
}

plot.decomposed.ts <- function(x, ...)
{
    plot(cbind(
               observed = x$random +
               if (x$type == "additive") x$trend + x$seasonal
               else x$trend * x$seasonal,
               trend    = x$trend,
               seasonal = x$seasonal,
               random   = x$random
               ),
         main = paste("Decomposition of", x$type, "time series"),
         ...)
}


KalmanLike <- function(y, mod, nit = 0)
{
    ## next call changes objects a, P, Pn: beware!
    x <- .Call("KalmanLike", y, mod$Z, mod$a, mod$P, mod$T, mod$V, mod$h,
               mod$Pn, as.integer(nit), FALSE, PACKAGE = "ts")
    names(x) <- c("ssq", "sumlog")
    s2 <- x[1]/length(y)
    list(Lik = 0.5*(log(x[1]/length(y)) + x[2]/length(y)), s2 = s2)
}

KalmanRun <- function(y, mod, nit = 0)
{
    ## next call changes objects a, P, Pn: beware!
    z <- .Call("KalmanLike", y, mod$Z, mod$a, mod$P, mod$T, mod$V, mod$h,
               mod$Pn, as.integer(nit), TRUE, PACKAGE = "ts")
    names(z) <- c("values", "resid", "states")
    x <- z$values
    s2 <- x[1]/length(y)
    z[[1]] <- c(Lik = 0.5*(log(x[1]/length(y)) + x[2]/length(y)), s2 = s2)
    z
}

KalmanForecast <- function(n.ahead = 10, mod)
{
    a <- an <- numeric(p <- length(mod$a))
    P <- Pn <- matrix(0, p, p)
    a[] <- mod$a
    P[] <- mod$P
    ## next call changes objects a, P
    x <- .Call("KalmanFore", as.integer(n.ahead), mod$Z, a, P,
               mod$T, mod$V, mod$h, PACKAGE = "ts")
    names(x) <- c("pred", "var")
    x
}

KalmanSmooth <- function(y, mod, nit = 0)
{
    z <- .Call("KalmanSmooth", y, mod$Z, mod$a, mod$P, mod$T, mod$V, mod$h,
               mod$Pn, as.integer(nit), PACKAGE = "ts")
    names(z) <- c("smooth", "var")
    dn <- dim(z$smooth)
    dim(z$var) <- dn[c(1, 2, 2)]
    z
}
StructTS <- function(x, type = c("level", "trend", "BSM"),
                     init = NULL, fixed = NULL, optim.control = NULL)
{
    KalmanLike2 <- function (y, mod, nit = 0)
    {
        x <- .Call("KalmanLike", y, mod$Z, mod$a, mod$P, mod$T, mod$V,
                   mod$h, mod$Pn, as.integer(nit), FALSE, PACKAGE = "ts")
        0.5 * sum(x)/length(y)
    }
    makeLevel <- function(x)
    {
        T <- matrix(1, 1, 1)
        Z <- 1
        a <- x[1]
        P <- Pn <- matrix(0, 1, 1)
        h <- 1
        V <- diag(1)
        return(Z, a, P, T, V, h, Pn)
    }
    makeTrend <- function(x)
    {
        T <- matrix(c(1,0,1,1), 2, 2)
        Z <- c(1, 0)
        a <- c(x[1], 0)
        P <- Pn <- matrix(0, 2, 2)
        h <- 1
        V <- diag(2)
        return(Z, a, P, T, V, h, Pn)
    }
    makeBSM <- function(x, nf)
    {
        if(nf <= 1) stop("frequency must be a positive integer for BSM")
        T <- matrix(0, nf + 1, nf + 1)
        T[1:2, 1:2] <- c(1, 0, 1, 1)
        T[3, ] <- c(0, 0, rep(-1, nf - 1))
        ind <- 3:nf
        T[cbind(ind+1, ind)] <- 1
        Z <- c(1, 0, 1, rep(0, nf - 2))
        a <- c(x[1], rep(0, nf))
        P <- Pn <- matrix(0, nf+1, nf+1)
        h <- 1
        V <- diag(c(1, 1, 1, rep(0, nf-2)))
        return(Z, a, P, T, V, h, Pn)
    }
    getLike <- function(par)
    {
        p <- cf
        p[mask] <- par
        if(all(p == 0)) return(1000)
        Z$V[cbind(1:np, 1:np)] <- p[-(np+1)]*vx
        Z$h <- p[np+1]*vx
        Z$P[] <- 1e6*vx
        Z$a <- a0
        res <- KalmanLike2(y, Z, -1)
#        print(c(res, p))
        res
    }

    series <- deparse(substitute(x))
    if(NCOL(x) > 1)
        stop("only implemented for univariate time series")
    x <- as.ts(x)
    type <- if(missing(type)) if(frequency(x) > 1) "BSM" else "trend"
    else match.arg(type)
    dim(x) <- NULL
    n <- length(x)
    xtsp <- tsp(x)
    nf <- frequency(x)
    Z <- switch(type,
                "level" = makeLevel(x),
                "trend" = makeTrend(x),
                "BSM" = makeBSM(x, nf)
                )
    a0 <- Z$a
    vx <- var(x, na.rm=TRUE)/100
    np <- switch(type, "level" = 1, "trend" = 2, "BSM" = 3)
    if (is.null(fixed)) fixed <- rep(NA, np+1)
    mask <- is.na(fixed)
    if(!any(mask)) stop("all parameters were fixed")
    cf <- fixed/vx
    if(is.null(init)) init <- rep(1, np+1) else init <- init/vx

    y <- x
    res <- optim(init[mask], getLike, method = "L-BFGS-B",
                 lower = rep(0, np+1), upper = rep(Inf, np+1),
                 control = optim.control)
        if(res$convergence > 0)
            warning(paste("possible convergence problem: optim gave code=",
                          res$convergence, res$message))
    coef <- cf
    coef[mask] <- res$par
    Z$V[cbind(1:np, 1:np)] <- coef[1:np]*vx
    Z$h <- coef[np+1]*vx
    Z$P[] <- 1e6*vx
    Z$a <- a0
    z <- KalmanRun(y, Z, -1)
    resid <- ts(z$resid)
    tsp(resid) <- xtsp
    Z0 <- Z; Z0$P[] <- 1e6*vx; Z0$a <- a0

    cn <- switch(type,
                 "level" = c("level"),
                 "trend" = c("level", "slope"),
                 "BSM" = c("level", "slope", "sea")
                 )
    states <- z$states
    if(type == "BSM") states <- states[, 1:3]
    dimnames(states) <- list(time(x), cn)
    states <- ts(states, start = xtsp[1], frequency = nf)

    coef <- coef*vx
    names(coef) <- switch(type,
                          "level" = c("level", "epsilon"),
                          "trend" = c("level", "slope", "epsilon"),
                          "BSM" = c("level", "slope", "seas", "epsilon")
                          )
    loglik <- -length(y) * res$value + length(y) * log(2 * pi)
    res <- list(coef = coef, loglik = loglik, data = y,
                residuals = resid, fitted = states,
                call = match.call(), series = series,
                code = res$convergence, model = Z, model0 = Z0, xtsp = xtsp)
    class(res) <- "StructTS"
    res
}

print.StructTS <- function(x, digits = max(3, getOption("digits") - 3), ...)
{
    cat("\nCall:", deparse(x$call, width = 75), "", sep = "\n")
    cat("Variances:\n")
    print.default(x$coef, print.gap = 2, digits=digits)
    invisible(x)
}

predict.StructTS <- function(object, n.ahead = 1, se.fit = TRUE, ...)
{
    xtsp <- object$xtsp
    z <- KalmanForecast(n.ahead, object$model)
    pred <- ts(z[[1]], start = xtsp[2] + 1/xtsp[3], frequency = xtsp[3])
    if (se.fit) {
        se <- ts(sqrt(z[[2]]), start = xtsp[2] + 1/xtsp[3],
                 frequency = xtsp[3])
        return(pred, se)
    }
    else return(pred)
}

tsdiag.StructTS <- function(object, gof.lag = 10, ...)
{
    ## plot standardized residuals, acf of residuals, Ljung-Box p-values
    oldpar<- par(mfrow = c(3, 1))
    on.exit(par(oldpar))
    rs <- object$resid
    stdres <- rs
    plot(stdres, type = "h", main = "Standardized Residuals", ylab = "")
    abline(h = 0)
    acf(object$resid, plot = TRUE, main = "ACF of Residuals",
        na.action = na.pass)
    nlag <- gof.lag
    pval <- numeric(nlag)
    for(i in 1:nlag) pval[i] <- Box.test(rs, i, type = "Ljung-Box")$p.value
    plot(1:nlag, pval, xlab = "lag", ylab = "p value", ylim = c(0,1),
         main = "p values for Ljung-Box statistic")
    abline(h = 0.05, lty = 2, col = "blue")
}


tsSmooth <- function(object, ...) UseMethod("tsSmooth")

tsSmooth.StructTS <- function(object, ...)
{
    res <- KalmanSmooth(object$data, object$model0, -1)$smooth
    dn <- dim(fitted(object))
    res <- res[, 1:dn[2], drop = FALSE]
    dimnames(res) <- dimnames(fitted(object))
    ts(res, start = object$xtsp[1], frequency = object$xtsp[3])
}
acf <-
    function (x, lag.max = NULL,
              type = c("correlation", "covariance", "partial"),
              plot = TRUE, na.action = na.fail, demean= TRUE, ...)
{
    type <- match.arg(type)
    if(type == "partial") {
        m <- match.call()
        m[[1]] <- as.name("pacf")
        m$type <- NULL
        return(eval(m, parent.frame()))
    }
    series <- deparse(substitute(x))
    x <- na.action(as.ts(x))
    x.freq <- frequency(x)
    x <- as.matrix(x)
    sampleT <- nrow(x)
    nser <- ncol(x)
    if (is.null(lag.max))
        lag.max <- floor(10 * (log10(sampleT) - log10(nser)))
    lag.max <- min(lag.max, sampleT - 1)
    if (lag.max < 1) stop("lag.max must be at least 1")
    if(demean) x <- sweep(x, 2, colMeans(x, na.rm = TRUE))
    lag <- matrix(1, nser, nser)
    lag[lower.tri(lag)] <- -1
    acf <- array(.C("acf",
                    as.double(x), as.integer(sampleT), as.integer(nser),
                    as.integer(lag.max), as.integer(type=="correlation"),
                    acf=double((lag.max+1) * nser * nser), NAOK = TRUE,
                    PACKAGE = "ts")$acf, c(lag.max + 1, nser, nser))
    lag <- outer(0:lag.max, lag/x.freq)
    acf.out <- structure(.Data = list(acf = acf, type = type,
        n.used = sampleT, lag = lag, series = series, snames = colnames(x)),
        class = "acf")
    if (plot) {
        plot.acf(acf.out, ...)
        return(invisible(acf.out))
    } else return(acf.out)
}

pacf <- function(x, lag.max, plot, na.action, ...) UseMethod("pacf")

pacf.default <- function(x, lag.max = NULL, plot = TRUE,
                         na.action = na.fail, ...)
{
    series <- deparse(substitute(x))
    if(is.matrix(x)) {
        m <- match.call()
        m[[1]] <- as.name("pacf.mts")
        return(eval(m, parent.frame()))
    }
    x <- na.action(as.ts(x))
    x.freq <- frequency(x)
    if(is.matrix(x))
        if(ncol(x) > 1) stop("univariate ts method")
        else x <- drop(x)
    sampleT <- length(x)
    if (is.null(lag.max))
        lag.max <- floor(10 * (log10(sampleT)))
    lag.max <- min(lag.max, sampleT - 1)
    if (lag.max < 1) stop("lag.max must be at least 1")
    x <- scale(x, TRUE, FALSE)
    acf <- drop(acf(x, lag.max = lag.max, plot = FALSE,
                    na.action = na.action)$acf)
    pacf <- array(.C("uni_pacf",
               as.double(acf),
               pacf = double(lag.max),
               as.integer(lag.max), PACKAGE="ts")$pacf, dim=c(lag.max,1,1))
    acf.out <- structure(.Data = list(acf = pacf, type = "partial",
                         n.used = sampleT,
                         lag = array((1:lag.max)/x.freq, dim=c(lag.max,1,1)),
                         series = series, snames = NULL),
                         class = "acf")
    if (plot) {
        plot.acf(acf.out, ...)
        return(invisible(acf.out))
    } else return(acf.out)
}

pacf.mts <- function(x, lag.max = NULL, plot = TRUE, na.action = na.fail, ...)
{
    series <- deparse(substitute(x))
    x <- na.action(as.ts(x))
    x.freq <- frequency(x)
    x <- as.matrix(x)
    if(any(is.na(x))) stop("NAs in x")
    sampleT <- nrow(x)
    nser <- ncol(x)
    if (is.null(lag.max))
        lag.max <- floor(10 * (log10(sampleT) - log10(nser)))
    lag.max <- min(lag.max, sampleT - 1)
    if (lag.max < 1) stop("lag.max must be at least 1")
    x <- sweep(x, 2, colMeans(x))
    lag <- matrix(1, nser, nser)
    lag[lower.tri(lag)] <- -1
    acf <- ar.yw(x, order.max = lag.max)$partialacf
    lag <- outer(1:lag.max, lag/x.freq)
    acf.out <- structure(.Data = list(acf = acf, type = "partial",
                         n.used = sampleT, lag = lag, series = series,
                         snames = colnames(x)),
                         class = "acf")
    if (plot) {
        plot.acf(acf.out, ...)
        return(invisible(acf.out))
    } else return(acf.out)
}

plot.acf <-
    function (x, ci = 0.95, type = "h", xlab = "Lag", ylab = NULL,
              ylim = NULL, main = NULL, ci.col="blue",
              ci.type = c("white", "ma"),
              max.mfrow = 6,
              ask = Npgs > 1 && dev.interactive(),
              mar = if(nser > 2) c(3,2,2,0.8) else par("mar"),
              oma = if(nser > 2) c(1,1.2,1,1) else par("oma"),
              mgp = if(nser > 2) c(1.5,0.6,0) else par("mgp"),
              xpd = par("xpd"),
              cex.main = if(nser > 2) 1 else par("cex.main"),
              verbose = getOption("verbose"),
              ...)
{
    ci.type <- match.arg(ci.type)
    if((nser <- ncol(x$lag)) < 1) stop("x$lag must have at least 1 column")
    if (is.null(ylab))
        ylab <- switch(x$type,
                       correlation = "ACF",
                       covariance = "ACF (cov)",
                       partial = "Partial ACF")
    if (is.null(snames <- x$snames))
        snames <- paste("Series ", if (nser == 1) x$series else 1:nser)

    with.ci <- ci > 0 && x$type != "covariance"
    with.ci.ma <- with.ci && ci.type == "ma" && x$type == "correlation"
    clim0 <- if (with.ci) qnorm((1 + ci)/2)/sqrt(x$n.used) else c(0, 0)

    Npgs <- 1 ## we will do [ Npgs x Npgs ] pages !
    nr <- nser
    if(nser > 1) { ## at most m x m (m := max.mfrow)  panels per page
        sn.abbr <- if(nser > 2) abbreviate(snames) else snames

        if(nser > max.mfrow) {
            ##  We need more than one page: The plots are laid out
            ##  such that we can manually paste the paper pages and get a
            ##  nice square layout with diagonal !
            ## NB: The same applies to pairs() where we'd want several pages
            Npgs <- ceiling(nser / max.mfrow)
            nr <- ceiling(nser / Npgs)  # <= max.mfrow
        }
        opar <- par(mfrow = rep(nr, 2), mar = mar, oma = oma, mgp = mgp,
                    ask = ask, xpd = xpd, cex.main = cex.main)
        on.exit(par(opar))
        if(verbose) {
            cat("par(*) : ")
            str(par("mfrow","cex", "cex.main","cex.axis","cex.lab","cex.sub"))
        }
    }

    for (I in 1:Npgs) for (J in 1:Npgs) {
        ## Page [ I , J ] : Now do   nr x nr  `panels' on this page
        iind <- (I-1)*nr + 1:nr
        jind <- (J-1)*nr + 1:nr
        if(verbose)
            cat("Page [",I,",",J,"]: i =",
                paste(iind,collapse=","),"; j =",
                paste(jind,collapse=","),"\n")
        for (i in iind) for (j in jind)
            if(max(i,j) > nser) {
                frame(); box(col = "light gray")
                ## the above is EXTREMELY UGLY; should have a version
                ## of frame() that really does advance a frame !!
            }
            else {
                clim <- if (with.ci.ma && i == j)
                    clim0 * sqrt(cumsum(c(1, 2*x$acf[-1, i, j]^2))) else clim0
                if (is.null(ylim)) {
                    ymin <- min(c(x$acf[, i, j], -clim))
                    ymax <- max(c(x$acf[, i, j], clim))
                    ylim <- c(ymin, ymax)
                }
                plot(x$lag[, i, j], x$acf[, i, j], type = type, xlab = xlab,
                     ylab = if(j==1) ylab else "", ylim = ylim, ...)
                abline(h = 0)
                if (with.ci && ci.type == "white")
                    abline(h = c(clim, -clim), col = ci.col, lty = 2)
                else if (with.ci.ma && i == j) {
                    lines(x$lag[, i, j], clim, col = ci.col, lty = 2)
                    lines(x$lag[, i, j], -clim, col = ci.col, lty = 2)
                }
                title(if (!is.null(main)) main else
                      if (i == j) snames[i]
                      else paste(sn.abbr[i], "&", sn.abbr[j]),
                      line = if(nser > 2) 1 else 2)
            }
        if(Npgs > 1) {                  # label the page
            mtext(paste("[",I,",",J,"]"), side=1, line = -0.2, adj=1,
                  col = "dark gray", cex = 1, outer = TRUE)
        }
    }
}

ccf <- function(x, y, lag.max = NULL,
                type = c("correlation", "covariance"),
                plot = TRUE, na.action = na.fail, ...)
{
    type <- match.arg(type)
    if(is.matrix(x) || is.matrix(y))
        stop("univariate time series only")
    X <- na.action(ts.union(x, y))
    colnames(X) <- c(deparse(substitute(x)), deparse(substitute(y)))
    acf.out <- acf(X, lag.max = lag.max, plot = FALSE, type = type)
    lag <- c(rev(acf.out$lag[-1,2,1]), 0, acf.out$lag[,1,2])
    y   <- c(rev(acf.out$acf[-1,2,1]), 0, acf.out$acf[,1,2])
    acf.out$acf <- array(y, dim=c(length(y),1,1))
    acf.out$lag <- array(lag, dim=c(length(y),1,1))
    acf.out$snames <- paste(acf.out$snames, collapse = " & ")
    if (plot) {
        plot.acf(acf.out, ...)
        return(invisible(acf.out))
    } else return(acf.out)
}
## based on, especially multivariate case, code by Martyn Plummer
ar <-
    function (x, aic = TRUE, order.max = NULL,
              method=c("yule-walker","burg", "ols", "mle", "yw",),
              na.action = na.fail, series = deparse(substitute(x)), ...)
{
    res <- switch(match.arg(method),
        "yule-walker" = ar.yw(x, aic=aic, order.max=order.max,
                  na.action = na.action, series=series, ...),
	"burg" = ar.burg(x, aic=aic, order.max=order.max,
                              na.action = na.action, series=series, ...),
	"ols" = ar.ols(x, aic=aic, order.max=order.max,
                              na.action = na.action, series=series, ...),
 	"mle" = ar.mle(x, aic=aic, order.max=order.max,
                              na.action = na.action, series=series, ...),
        "yw" = ar.yw(x, aic=aic, order.max=order.max,
                  na.action = na.action, series=series, ...)
   )
    res$call <- match.call()
    res
}

ar.yw <- function(x, ...) UseMethod("ar.yw")

ar.yw.default <-
    function (x, aic = TRUE, order.max = NULL, na.action = na.fail,
              demean = TRUE, series = NULL, ...)
{
    if(is.null(series)) series <- deparse(substitute(x))
    ists <- is.ts(x)
    x <- na.action(as.ts(x))
    if(ists)  xtsp <- tsp(x)
    xfreq <- frequency(x)
    x <- as.matrix(x)
    if(any(is.na(x))) stop("NAs in x")
    nser <- ncol(x)
    if (demean) {
        xm <- colMeans(x)
        x <- sweep(x, 2, xm)
    } else xm <- rep(0, nser)
    n.used <- nrow(x)
    order.max <- if (is.null(order.max)) floor(10 * log10(n.used))
                 else round(order.max)
    if (order.max < 1) stop("order.max must be >= 1")
    xacf <- acf(x, type = "covariance", lag.max = order.max, plot = FALSE,
                demean = demean)$acf
    if(nser > 1) {
        ## multivariate case
        snames <- colnames(x)
        A <- B <- array(0, dim = c(order.max + 1, nser, nser))
        A[1, , ] <- B[1, , ] <- diag(nser)
        EA <- EB <- xacf[1, , , drop = TRUE]
        partialacf <- array(dim = c(order.max, nser, nser))
        xaic <- numeric(order.max + 1)
        solve.yw <- function(m) {
            # Solve Yule-Walker equations with Whittle's
            # generalization of the Levinson(-Durbin) algorithm
            betaA <- betaB <- 0
            for (i in 0:m) {
                betaA <- betaA + A[i + 1, , ] %*% xacf[m + 2 - i, , ]
                betaB <- betaB + B[i + 1, , ] %*% t(xacf[m + 2 - i, , ])
            }
            KA <- -t(qr.solve(t(EB), t(betaA)))
            KB <- -t(qr.solve(t(EA), t(betaB)))
            EB <<- (diag(nser) - KB %*% KA) %*% EB
            EA <<- (diag(nser) - KA %*% KB) %*% EA
            Aold <- A
            Bold <- B
            for (i in 1:(m + 1)) {
                A[i + 1, , ] <<- Aold[i + 1, , ] + KA %*% Bold[m + 2 - i, , ]
                B[i + 1, , ] <<- Bold[i + 1, , ] + KB %*% Aold[m + 2 - i, , ]
            }
        }
        cal.aic <- function() { # omits mean params, that is constant adj
            det <- abs(prod(diag(qr(EA)$qr)))
            return(n.used * log(det) + 2 * m * nser * nser)
        }
        cal.resid <- function() {
            resid <- array(0, dim = c(n.used - order, nser))
            for (i in 0:order) {
                resid <- resid + x[(order - i + 1):(n.used - i),
                                   , drop = FALSE] %*% t(ar[i + 1, , ])
            }
            return(rbind(matrix(NA, order, nser), resid))
        }
        order <- 0
        for (m in 0:order.max) {
            xaic[m + 1] <- cal.aic()
            if (!aic || xaic[m + 1] == min(xaic[1:(m + 1)])) {
                ar <- A
                order <- m
                var.pred <- EA * n.used/(n.used - nser * (m + 1))
            }
            if (m < order.max) {
                solve.yw(m)
                partialacf[m + 1, , ] <- -A[m + 2, , ]
            }
        }
        xaic <- xaic - min(xaic)
        names(xaic) <- 0:order.max
        resid <- cal.resid()
        if(order > 0 ) {
            ar <- -ar[2:(order + 1), , , drop = FALSE]
            dimnames(ar) <- list(1:order, snames, snames)
        } else ar <- array(0, dim=c(0, nser, nser),
                           dimnames=list(NULL, snames, snames))
        dimnames(var.pred) <- list(snames, snames)
        dimnames(partialacf) <- list(1:order.max, snames, snames)
        colnames(resid) <- colnames(x)
    } else {
        ## univariate case
        r <- as.double(drop(xacf))
        z <- .Fortran("eureka",
                      as.integer(order.max),
                      r, r,
                      coefs=double(order.max^2),
                      vars=double(order.max),
                      double(order.max), PACKAGE="ts")
        coefs <- matrix(z$coefs, order.max, order.max)
        partialacf <- array(diag(coefs), dim=c(order.max, 1, 1))
        var.pred <- c(r[1], z$vars)
        xaic <- n.used * log(var.pred) + 2 * (0:order.max) + 2 * demean
        xaic <- xaic - min(xaic)
        names(xaic) <- 0:order.max
        order <- if (aic) (0:order.max)[xaic == 0] else order.max
        ar <- if (order > 0) coefs[order, 1:order] else numeric(0)
        var.pred <- var.pred[order+1]
        ## Splus compatibility fix
        var.pred <- var.pred * n.used/(n.used - (order + 1))
        if(order > 0)
            resid <- c(rep(NA, order), embed(x, order+1) %*% c(1, -ar))
        else resid <- as.vector(x)
        if(ists) {
            attr(resid, "tsp") <- xtsp
            attr(resid, "class") <- "ts"
        }
    }
    res <- list(order=order, ar=ar, var.pred=var.pred, x.mean = drop(xm),
                aic = xaic, n.used=n.used, order.max=order.max,
                partialacf=partialacf, resid=resid, method = "Yule-Walker",
                series=series, frequency=xfreq, call=match.call())
    if(nser == 1 && order > 0)
        res$asy.var.coef <-
            solve(toeplitz(drop(xacf)[seq(length=order)]))*var.pred/n.used
    class(res) <- "ar"
    res
}

print.ar <- function(x, digits = max(3, getOption("digits") - 3), ...)
{
    cat("\nCall:\n", deparse(x$call), "\n\n", sep = "")
    nser <- NCOL(x$var.pred)
    if(nser > 1) {
        if(!is.null(x$x.intercept))
            res <- x[c("ar", "x.intercept", "var.pred")]
        else res <- x[c("ar", "var.pred")]
        res$ar <- aperm(res$ar, c(2,3,1))
        print(res, digits=digits)
    } else {
        if(x$order > 0) {
            cat("Coefficients:\n")
            coef <- drop(round(x$ar, digits = digits))
            names(coef) <- seq(length=x$order)
            print.default(coef, print.gap = 2)
        }
        if(!is.null(xint <- x$x.intercept) && !is.na(xint))
            cat("\nIntercept: ", format(xint, digits = digits),
                " (", format(x$asy.se.coef$x.mean, digits = digits),
                ") ", "\n", sep="")
        cat("\nOrder selected", x$order, " sigma^2 estimated as ",
            format(x$var.pred, digits = digits),"\n")

    }
    invisible(x)
}

predict.ar <- function(object, newdata, n.ahead = 1, se.fit=TRUE, ...)
{
    if(missing(newdata)) {
        newdata <- eval(parse(text=object$series))
        if (!is.null(nas <- object$call$na.action))
            newdata <- eval(call(nas, newdata))
    }
    nser <- NCOL(newdata)
    ar <- object$ar
    p <- object$order
    st <- tsp(as.ts(newdata))[2]
    dt <- deltat(newdata)
    xfreq <- frequency(newdata)
    tsp(newdata) <- NULL
    class(newdata) <- NULL
    if(NCOL(ar) != nser)
        stop("number of series in fit and newdata do not match")
    n <- NROW(newdata)
    if(nser > 1) {
        if(is.null(object$x.intercept)) xint <- rep(0, nser)
        else xint <- object$x.intercept
        x <- rbind(sweep(newdata, 2, object$x.mean),
                   matrix(rep(0, nser), n.ahead, nser, byrow=TRUE))
        if(p > 0) {
            for(i in 1:n.ahead) {
                x[n+i,] <- ar[1,,] %*% x[n+i-1,] + xint
                if(p > 1) for(j in 2:p)
                    x[n+i,] <- x[n+i,] + ar[j,,] %*% x[n+i-j,]
            }
            pred <- x[n+(1:n.ahead), ]
        } else {
            pred <- matrix(xint, n.ahead, nser, byrow=TRUE)
        }
        pred <- pred + matrix(object$x.mean, n.ahead, nser, byrow=TRUE)
        colnames(pred) <- colnames(object$var.pred)
        if(se.fit) {
            warning("se.fit not yet implemented for multivariate models")
            se <- matrix(NA, n.ahead, nser)
        }
    } else {
        if(is.null(object$x.intercept)) xint <- 0
        else xint <- object$x.intercept
        x <- c(newdata - object$x.mean, rep(0, n.ahead))
        if(p > 0) {
            for(i in 1:n.ahead) {
                x[n+i] <- sum(ar * x[n+i - (1:p)]) + xint
            }
            pred <- x[n+(1:n.ahead)]
            if(se.fit) {
                npsi <- n.ahead - 1
                psi <- .C("artoma",
                        as.integer(object$order), as.double(ar),
                        psi = double(npsi+object$order+1),
                        as.integer(npsi), PACKAGE="ts")$psi[1:npsi]
                vars <- cumsum(c(1, psi^2))
                se <- sqrt(object$var.pred*vars)[1:n.ahead]
            }
        } else {
            pred <- rep(xint, n.ahead)
            if (se.fit) se <- rep(sqrt(object$var.pred), n.ahead)
        }
        pred <- pred + rep(object$x.mean, n.ahead)
    }
    pred <- ts(pred, start = st + dt, frequency=xfreq)
    if(se.fit) se <- ts(se, start = st + dt, frequency=xfreq)
    if(se.fit) return(pred, se) else return(pred)
}

## ar.burg by B.D. Ripley based on R version by Martyn Plummer
ar.burg <- function(x, ...) UseMethod("ar.burg")
ar.burg.default <-
    function (x, aic = TRUE, order.max = NULL, na.action = na.fail,
                   demean = TRUE, series = NULL, var.method = 1, ...)
{
    if(is.null(series)) series <- deparse(substitute(x))
    if (!is.null(dim(x)))
        stop("Burg's algorithm only implemented for univariate series")
    if (ists <- is.ts(x)) xtsp <- tsp(x)
    x <- na.action(as.ts(x))
    if(any(is.na(x))) stop("NAs in x")
    if (ists)  xtsp <- tsp(x)
    xfreq <- frequency(x)
    x <- as.vector(x)
    if (demean) {
        x.mean <- mean(x)
        x <- x - x.mean
    } else x.mean <- 0
    n.used <- length(x)
    order.max <- if (is.null(order.max)) floor(10 * log10(n.used))
    else floor(order.max)
    if (order.max < 1) stop("order.max must be >= 1")
    xaic <- numeric(order.max + 1)
    z <- .C("burg",
            as.integer(n.used),
            as.double(x),
            as.integer(order.max),
            coefs=double(order.max^2),
            var1=double(1+order.max),
            var2=double(1+order.max), PACKAGE="ts"
            )
    coefs <- matrix(z$coefs, order.max, order.max)
    partialacf <- array(diag(coefs), dim = c(order.max, 1, 1))
    var.pred <- if(var.method == 1) z$var1 else z$var2
    xaic <- n.used * log(var.pred) + 2 * (0:order.max) + 2 * demean
    xaic <- xaic - min(xaic)
    names(xaic) <- 0:order.max
    order <- if (aic) (0:order.max)[xaic == 0] else order.max
    ar <- if (order > 0) coefs[order, 1:order] else numeric(0)
    var.pred <- var.pred[order+1]
    if(order > 0)
        resid <- c(rep(NA, order), embed(x, order+1) %*% c(1, -ar))
    else resid <- as.vector(x)
    if(ists) {
        attr(resid, "tsp") <- xtsp
        attr(resid, "class") <- "ts"
    }
    res <- list(order = order, ar = ar, var.pred = var.pred, x.mean = x.mean,
                aic = xaic, n.used = n.used, order.max = order.max,
                partialacf = partialacf, resid = resid,
                method = ifelse(var.method==1,"Burg","Burg2"),
                series = series, frequency = xfreq, call = match.call())
    xacf <- acf(x, type = "covariance", lag.max = order, plot=FALSE)$acf
    if(order > 0) res$asy.var.coef <- solve(toeplitz(drop(xacf)[seq(length=order)]))*var.pred/n.used
    class(res) <- "ar"
    return(res)
}

"ar.burg.mts" <-
function (x, aic = TRUE, order.max = NULL, na.action = na.fail,
    demean = TRUE, series = NULL, var.method = 1, ...)
{
    if (is.null(series))
        series <- deparse(substitute(x))
    if (ists <- is.ts(x))
        xtsp <- tsp(x)
    x <- na.action(as.ts(x))
    if (any(is.na(x)))
        stop("NAs in x")
    if (ists)
        xtsp <- tsp(x)
    xfreq <- frequency(x)
    x <- as.matrix(x)
    nser <- ncol(x)
    n.used <- nrow(x)
    if (demean) {
        x.mean <- colMeans(x)
        x <- sweep(x, 2, x.mean)
    }
    else x.mean <- rep(0, nser)
    order.max <- if (is.null(order.max))
        floor(10 * log10(n.used))
    else floor(order.max)
    xaic <- numeric(order.max + 1)
    z <- .C("multi_burg",
            as.integer(n.used),
            resid = as.double(x),
            as.integer(order.max),
            as.integer(nser),
            coefs = double((1 + order.max) * nser * nser),
            pacf = double((1 + order.max) * nser * nser),
            var = double((1 + order.max) * nser * nser),
            aic = double(1 + order.max),
            order = integer(1),
            as.integer(aic),
            as.integer(var.method),
            PACKAGE = "ts")
    partialacf <- aperm(array(z$pacf, dim = c(nser, nser, order.max +
        1)), c(3, 2, 1))[-1, , , drop = FALSE]
    var.pred <- aperm(array(z$var, dim = c(nser, nser, order.max +
        1)), c(3, 2, 1))
    xaic <- z$aic - min(z$aic)
    names(xaic) <- 0:order.max
    order <- z$order
    ar <- if (order > 0)
        coefs <- -aperm(array(z$coefs, dim = c(nser, nser, order.max +
            1)), c(3, 2, 1))[2:(order + 1), , , drop = FALSE]
    else array(dim = c(0, nser, nser))
    var.pred <- var.pred[order + 1, , , drop = TRUE]
    resid <- matrix(z$resid, nrow = n.used, ncol = nser)
    if (order > 0)
        resid[1:order, ] <- NA
    if (ists) {
        attr(resid, "tsp") <- xtsp
        attr(resid, "class") <- "mts"
    }
    snames <- colnames(x)
    colnames(resid) <- snames
    dimnames(ar) <- list(seq(length=order), snames, snames)
    dimnames(var.pred) <- list(snames, snames)
    dimnames(partialacf) <- list(1:order.max, snames, snames)
    res <- list(order = order, ar = ar, var.pred = var.pred,
        x.mean = x.mean, aic = xaic, n.used = n.used, order.max = order.max,
        partialacf = partialacf, resid = resid, method = ifelse(var.method ==
            1, "Burg", "Burg2"), series = series, frequency = xfreq,
        call = match.call())
    class(res) <- "ar"
    return(res)
}
ar.mle <- function (x, aic = TRUE, order.max = NULL, na.action = na.fail,
                    demean = TRUE, series = NULL, ...)
{
    if(is.null(series)) series <- deparse(substitute(x))
    ists <- is.ts(x)
    if (!is.null(dim(x)))
        stop("MLE only implemented for univariate series")
    x <- na.action(as.ts(x))
    if(any(is.na(x))) stop("NAs in x")
    if(ists)  xtsp <- tsp(x)
    xfreq <- frequency(x)
    x <- as.vector(x)
    n.used <- length(x)
    order.max <- if (is.null(order.max)) min(12, floor(10 * log10(n.used)))
    else round(order.max)

    if (order.max < 0) stop ("order.max must be >= 0")
    if (aic) {
        coefs <- matrix(NA, order.max+1, order.max+1)
        var.pred <- numeric(order.max+1)
        xaic <- numeric(order.max+1)
        xm <- if(demean) mean(x) else 0
        coefs[1, 1] <- xm
        var0 <- sum((x-xm)^2)/n.used
        var.pred[1] <- var0
        xaic[1] <- n.used * log(var0) + 2 * demean + 2 + n.used + n.used * log(2 * pi)
        for(i in 1:order.max) {
            fit <- arima0(x, order=c(i, 0, 0), include.mean=demean)
            coefs[i+1, 1:(i+demean)] <- fit$coef[1:(i+demean)]
            xaic[i+1] <- fit$aic
            var.pred[i+1] <- fit$sigma2
        }
        xaic <- xaic - min(xaic)
        names(xaic) <- 0:order.max
        order <- (0:order.max)[xaic == 0]
        ar <- coefs[order+1, 1:order]
        x.mean <- coefs[order+1, order+1]
        var.pred <- var.pred[order+1]
    } else {
        order <- order.max
        fit <- arima0(x, order=c(order, 0, 0), include.mean=demean)
        coefs <- fit$coef
        if(demean) {
            ar <- coefs[-length(coefs)]
            x.mean <- coefs[length(coefs)]
        } else {
            ar <- coefs
            x.mean <- 0
        }
        var.pred <- fit$sigma2
        xaic <- structure(0, names=order)
    }
    if(order > 0)
        resid <- c(rep(NA, order), embed(x - x.mean, order+1) %*% c(1, -ar))
    else resid <- as.vector(x) - x.mean
    if(ists) {
        attr(resid, "tsp") <- xtsp
        attr(resid, "class") <- "ts"
    }
    res <- list(order = order, ar = ar, var.pred = var.pred,
                x.mean = x.mean, aic = xaic,
                n.used = n.used, order.max = order.max,
                partialacf=NULL, resid=resid, method = "MLE",
                series = series, frequency = xfreq, call = match.call())
    if(order > 0) {
        xacf <- acf(x, type = "covariance", lag.max = order, plot=FALSE)$acf
        res$asy.var.coef <- solve(toeplitz(drop(xacf)[seq(length=order)])) *
            var.pred/n.used
    }
    class(res) <- "ar"
    res
}
## code by Adrian Trapletti
ar.ols <- function (x, aic = TRUE, order.max = NULL, na.action = na.fail,
                    demean = TRUE, intercept = demean, series = NULL, ...)
{
    if(is.null(series)) series <- deparse(substitute(x))
    rescale <- TRUE
    ists <- is.ts(x)
    x <- na.action(as.ts(x))
    xfreq <- frequency(x)
    if(any(is.na(x))) stop("NAs in x")
    if(ists)  xtsp <- tsp(x)
    x <- as.matrix(x)
    n.used <- nrow(x)
    nser <- ncol(x)
    if(rescale) {
        sc <- sqrt(drop(apply(x, 2, var)))
        x <- x/rep(sc, rep(n.used, nser))
    } else sc <- rep(1, nser)
    order.max <- if (is.null(order.max)) floor(10 * log10(n.used))
    else round(order.max)

    if (order.max < 0) stop ("order.max must be >= 0")
    if (aic) order.min <- 0
    else order.min <- order.max
    A <- vector("list", order.max - order.min + 1)
    varE <- vector("list", order.max - order.min + 1)
    seA <- vector("list", order.max - order.min + 1)
    aic <- rep(Inf, order.max - order.min + 1)

    det <- function(x) { prod(diag(qr(x)$qr))*(-1)^(ncol(x)-1) }

    ## remove means for conditioning
    if(demean) {
        xm <- colMeans(x)
        x <- sweep(x, 2, xm)
    } else xm <- rep(0, nser)
    ## Fit models of increasing order

    for (m in order.min:order.max)
    {
        y <- embed(x, m+1)
        if(intercept) {
            if (m > 0) X <- cbind(rep(1,nrow(y)), y[, (nser+1):ncol(y)])
            else X <- as.matrix(rep(1, nrow(y)))
        } else {
            if (m > 0) X <- y[, (nser+1):ncol(y)]
            else X <- matrix(0, nrow(y), 0)
        }
        Y <- t(y[, 1:nser])
        N <- ncol(Y)
        XX <- t(X)%*%X
        rank <- qr(XX)$rank
        if (rank != nrow(XX))
        {
            warning (paste("Model order", m))
            warning ("Singularities in the computation of the projection matrix")
            warning (paste("Results are only valid up to model order", m - 1))
            break
        }
        P <- if(ncol(XX) > 0) solve(XX) else XX
        A[[m - order.min + 1]] <- Y %*% X %*% P
        YH <- A[[m - order.min + 1]] %*% t(X)
        E <- (Y - YH)
        varE[[m - order.min + 1]] <- E %*% t(E)/N
        varA <- P %x% (varE[[m - order.min + 1]])
        seA[[m - order.min+1]] <- if(ncol(varA) > 0) sqrt(diag(varA))
        else numeric(0)
        aic[m - order.min+1] <-
            n.used*log(det(varE[[m-order.min+1]]))+2*nser*(nser*m+intercept)
    }

    m <- which(aic==min(aic)) + order.min - 1 # Determine best model

    ## Recalculate residuals of best model

    y <- embed(x, m+1)
    AA <- A[[m - order.min + 1]]
    if(intercept) {
        xint <- AA[, 1]
        ar <- AA[, -1]
        if (m > 0) X <- cbind(rep(1,nrow(y)), y[, (nser+1):ncol(y)])
        else X <- as.matrix(rep(1, nrow(y)))
    } else {
        if (m > 0) X <- y[, (nser+1):ncol(y)]
        else X <- matrix(0, nrow(y), 0)
        xint <- NULL
        ar <- AA
    }
    Y <- t(y[, 1:nser, drop=FALSE])
    YH <- AA %*% t(X)
    E <- drop(rbind(matrix(NA, m, nser), t(Y - YH)))

    aic <- aic - min(aic)
    names(aic) <- order.min:order.max
    dim(ar) <- c(nser, nser, m)
    ar <- aperm(ar, c(3,1,2))
    ses <- seA[[m - order.min + 1]]
    if(intercept) {
        sem <- ses[1:nser]
        ses <- ses[-(1:nser)]
    } else sem <- rep(0, nser)
    dim(ses) <- c(nser, nser, m)
    ses <- aperm(ses, c(3,1,2))
    var.pred <- varE[[m - order.min + 1]]
    if(nser > 1) {
        snames <- colnames(x)
        dimnames(ses) <- dimnames(ar) <- list(seq(length=m), snames, snames)
        dimnames(var.pred) <- list(snames, snames)
        names(sem) <- colnames(E) <- snames
    }
    if(ists) {
        attr(E, "tsp") <- xtsp
        attr(E, "class") <- "ts"
    }
    if(rescale) {
        xm <- xm * sc
        if(!is.null(xint)) xint <- xint * sc
        aa <- outer(sc, 1/sc)
        if(nser > 1 && m > 0)
            for(i in 1:m) ar[i,,] <- ar[i,,]*aa
        var.pred <- var.pred * outer(sc, sc)
        E <- E * rep(sc, rep(NROW(E), nser))
        sem <- sem*sc
        if(m > 0)
            for(i in 1:m) ses[i,,] <- ses[i,,]*aa
    }
    res <- list(order = m, ar = ar, var.pred = var.pred,
                x.mean = xm, x.intercept = xint, aic = aic,
                n.used = n.used, order.max = order.max,
                partialacf = NULL, resid = E, method = "Unconstrained LS",
                series = series, frequency = xfreq, call = match.call(),
                asy.se.coef = list(x.mean = sem, ar=drop(ses)))
    class(res) <- "ar"
    res
}
"ar.yw.mts" <-
function (x, aic = TRUE, order.max = NULL, na.action = na.fail,
    demean = TRUE, series = NULL, var.method = 1, ...)
{
    if (is.null(series))
        series <- deparse(substitute(x))
    if (ists <- is.ts(x))
        xtsp <- tsp(x)
    x <- na.action(as.ts(x))
    if (any(is.na(x)))
        stop("NAs in x")
    if (ists)
        xtsp <- tsp(x)
    xfreq <- frequency(x)
    x <- as.matrix(x)
    nser <- ncol(x)
    n.used <- nrow(x)
    if (demean) {
        x.mean <- colMeans(x)
        x <- sweep(x, 2, x.mean)
    }
    else x.mean <- rep(0, nser)
    order.max <- if (is.null(order.max))
        floor(10 * log10(n.used))
    else floor(order.max)
    if (order.max < 1)
        stop("order.max must be >= 1")
    xacf <- acf(x, type = "cov", plot = FALSE, lag.max = order.max)$acf
    z <- .C("multi_yw",
            aperm(xacf, c(3, 2, 1)),
            as.integer(n.used),
            as.integer(order.max),
            as.integer(nser),
            coefs = double((1 + order.max) * nser * nser),
            pacf = double((1 + order.max) * nser * nser),
            var = double((1 + order.max) * nser * nser),
            aic = double(1 + order.max),
            order = integer(1),
            as.integer(aic),
####        as.integer(var.method),
            PACKAGE = "ts")
    partialacf <- aperm(array(z$pacf, dim = c(nser, nser, order.max +
        1)), c(3, 2, 1))[-1, , , drop = FALSE]
    var.pred <- aperm(array(z$var, dim = c(nser, nser, order.max +
        1)), c(3, 2, 1))
    xaic <- z$aic - min(z$aic)
    names(xaic) <- 0:order.max
    order <- z$order
    resid <- x
    if (order > 0) {
        ar <- -aperm(array(z$coefs, dim = c(nser, nser, order.max +
            1)), c(3, 2, 1))[2:(order + 1), , , drop = FALSE]
        for (i in 1:order) resid[-(1:order), ] <- resid[-(1:order),
            ] - x[(order - i + 1):(n.used - i), ] %*% t(ar[i,
            , ])
        resid[1:order, ] <- NA
    }
    else ar <- array(dim = c(0, nser, nser))
    var.pred <- var.pred[order + 1, , , drop = TRUE] * n.used/(n.used -
        nser * (demean + order))
    if (ists) {
        attr(resid, "tsp") <- xtsp
        attr(resid, "class") <- c("mts", "ts")
    }
    snames <- colnames(x)
    colnames(resid) <- snames
    dimnames(ar) <- list(seq(length=order), snames, snames)
    dimnames(var.pred) <- list(snames, snames)
    dimnames(partialacf) <- list(1:order.max, snames, snames)
    res <- list(order = order, ar = ar, var.pred = var.pred,
        x.mean = x.mean, aic = xaic, n.used = n.used, order.max = order.max,
        partialacf = partialacf, resid = resid, method = "Yule-Walker",
        series = series, frequency = xfreq, call = match.call())
    class(res) <- "ar"
    return(res)
}
arima <- function(x, order = c(0, 0, 0),
                  seasonal = list(order = c(0, 0, 0), period = NA),
                  xreg = NULL, include.mean = TRUE,
                  transform.pars = TRUE, fixed = NULL, init = NULL,
                  method = c("CSS-ML", "ML", "CSS"), n.cond,
                  optim.control = list(), kappa = 1e6)
{
    "%+%" <- function(a, b) .Call("convolve", a, b, PACKAGE = "ts")

    upARIMA <- function(mod, phi, theta)
    {
        p <- length(phi); q <- length(theta)
        mod$phi <- phi; mod$theta <- theta
        r <- max(p, q + 1)
        if(p > 0) mod$T[1:p, 1] <- phi
        if(r > 1)
            mod$Pn[1:r, 1:r] <- .Call("getQ0", phi, theta, PACKAGE = "ts")
        else mod$Pn[1, 1] <- 1/(1 - phi^2)
        mod$a[] <- 0
        mod
    }

    arimaSS <- function(y, mod)
    {
        ## next call changes objects a, P, Pn so beware!
        .Call("ARIMA_Like", y, mod$phi, mod$theta, mod$Delta,
              mod$a, mod$P, mod$Pn, as.integer(0), TRUE, PACKAGE = "ts")
    }

    armafn <- function(p, trans)
    {
        par <- coef
        par[mask] <- p
        trarma <- .Call("ARIMA_transPars", par, arma, trans, PACKAGE = "ts")
        Z <- upARIMA(mod, trarma[[1]], trarma[[2]])
        if(ncxreg > 0) x <- x - xreg %*% par[narma + (1:ncxreg)]
        ## next call changes objects a, P, Pn so beware!
        res <- .Call("ARIMA_Like", x, Z$phi, Z$theta, Z$Delta,
                     Z$a, Z$P, Z$Pn, as.integer(0), FALSE, PACKAGE = "ts")
        s2 <- res[1]/res[3]
        0.5*(log(s2) + res[2]/res[3])
    }

    armaCSS <- function(p)
    {
        par <- as.double(fixed)
        par[mask] <- p
        trarma <- .Call("ARIMA_transPars", par, arma, FALSE, PACKAGE = "ts")
        if(ncxreg > 0) x <- x - xreg %*% par[narma + (1:ncxreg)]
        res <- .Call("ARIMA_CSS", x, arma, trarma[[1]], trarma[[2]],
                     as.integer(ncond), FALSE, PACKAGE = "ts")
        0.5 * log(res)
    }

    arCheck <- function(ar)
    {
        p <- max(which(c(1, -ar) != 0)) - 1
        if(!p) return(TRUE)
        all(Mod(polyroot(c(1, -ar[1:p]))) > 1)
    }

    maInvert <- function(ma)
    {
        ## polyroot can't cope with leading zero.
        q <- length(ma)
        q0 <- max(which(c(1,ma) != 0)) - 1
        if(!q0) return(ma)
        roots <- polyroot(c(1, ma[1:q0]))
        ind <- Mod(roots) < 1
        if(all(!ind)) return(ma)
        if(q0 == 1) return(c(1/ma[1], rep(0, q-q0)))
        roots[ind] <- 1/roots[ind]
        x <- 1
        for(r in roots) x <- c(x, 0) - c(0, x)/r
        c(Re(x[-1]), rep(0, q-q0))
    }

    series <- deparse(substitute(x))
    if(NCOL(x) > 1)
        stop("only implemented for univariate time series")
    method <- match.arg(method)

    x <- as.ts(x)
    dim(x) <- NULL
    n <- length(x)

    if(!missing(order))
        if(!is.numeric(order) || length(order) != 3 || any(order < 0))
            stop("`order' must be a non-negative numeric vector of length 3")
    if(!missing(seasonal))
        if(is.list(seasonal)) {
            if(is.null(seasonal$order))
                stop("`seasonal' must be a list with component `order'")
            if(!is.numeric(seasonal$order) || length(seasonal$order) != 3
               || any(seasonal$order < 0))
                stop("`seasonal$order' must be a non-negative numeric vector of length 3")
        } else if(is.numeric(order)) {
            if(length(order) == 3) seasonal <- list(order=seasonal)
            else ("`seasonal' is of the wrong length")
        } else stop("`seasonal' must be a list with component `order'")

    if (is.null(seasonal$period) || is.na(seasonal$period)
        ||seasonal$period == 0) seasonal$period <- frequency(x)
    arma <- as.integer(c(order[-2], seasonal$order[-2], seasonal$period,
                         order[2], seasonal$order[2]))
    narma <- sum(arma[1:4])
    if(narma == 0) stop("no AR nor MA terms to be fitted")

    xtsp <- tsp(x)
    tsp(x) <- NULL
    Delta <- 1
    for(i in seq(len = order[2])) Delta <- Delta %+% c(1, -1)
    for(i in seq(len = seasonal$order[2]))
        Delta <- Delta %+% c(1, rep(0, seasonal$period-1), -1)
    Delta <- - Delta[-1]
    nd <- order[2] + seasonal$order[2]
    n.used <- sum(!is.na(x)) - length(Delta)
    if (is.null(xreg)) {
        ncxreg <- 0
    } else {
        nmxreg <- deparse(substitute(xreg))
        if (NROW(xreg) != n) stop("lengths of x and xreg do not match")
        ncxreg <- NCOL(xreg)
        xreg <- as.matrix(xreg)
    }
    class(xreg) <- NULL
    if (ncxreg > 0 && is.null(colnames(xreg)))
        colnames(xreg) <-
            if(ncxreg == 1) nmxreg else paste(nmxreg, 1:ncxreg, sep = "")
    if (include.mean && (nd == 0)) {
        xreg <- cbind(intercept = rep(1, n), xreg = xreg)
        ncxreg <- ncxreg + 1
    }
    if(method == "CSS-ML") {
        anyna <- any(is.na(x))
        if(ncxreg) anyna <- anyna || any(is.na(xreg))
        if(anyna) method <- "ML"
    }

    if (method == "CSS" || method == "CSS-ML") {
        ncond <- order[2] + seasonal$order[2] * seasonal$period
        ncond1 <- order[1] + seasonal$period * seasonal$order[1]
        ncond <- if (!missing(n.cond)) ncond + max(n.cond, ncond1)
        else ncond + ncond1
    } else ncond <- 0

    if (is.null(fixed)) fixed <- rep(NA, narma + ncxreg)
    else if(length(fixed) != narma + ncxreg) stop("wrong length for fixed")
    mask <- is.na(fixed)
    if(!any(mask)) stop("all parameters were fixed")
    init0 <- rep(0, narma)
    parscale <- rep(1, narma)
    if (ncxreg) {
        cn <- colnames(xreg)
        orig.xreg <- (ncxreg == 1) || any(!mask[narma + 1:ncxreg] )
        if(!orig.xreg) {
            S <- svd(na.omit(xreg))
            xreg <- xreg %*% S$v
        }
        fit <- lm(x ~ xreg - 1, na.action = na.omit)
        n.used <- sum(!is.na(resid(fit))) - length(Delta)
        init0 <- c(init0, coef(fit))
        ses <- summary(fit)$coef[,2]
        parscale <- c(parscale, 10*ses)
    }
    if (n.used <= 0) stop("too few non-missing observations")

    if(!is.null(init)) {
        if(length(init) != length(init0))
            stop("`init' is of the wrong length")
        if(any(ind <- is.na(init))) init[ind] <- init0[ind]
        if(method == "ML") {
            ## check stationarity
            if(arma[1] > 0)
                if(!arCheck(init[1:arma[1]]))
                    stop("non-stationary AR part")
            if(arma[3] > 0)
                if(!arCheck(init[sum(arma[1:2]) + 1:arma[3]]))
                    stop("non-stationary seasonal AR part")
            if(transform.pars)
                init <- .Call("ARIMA_Invtrans", as.double(init), arma,
                              PACKAGE = "ts")
        }
    } else init <- init0

    coef <- as.double(fixed)
    if(!("parscale" %in% names(optim.control)))
       optim.control$parscale <- parscale[mask]

    if(method == "CSS") {
        res <- optim(init[mask], armaCSS,  method = "BFGS", hessian = TRUE,
                     control = optim.control)
        if(res$convergence > 0)
            warning(paste("possible convergence problem: optim gave code=",
                          res$convergence))
        coef[mask] <- res$par
        ## set model for predictions
        trarma <- .Call("ARIMA_transPars", coef, arma, FALSE,
                        PACKAGE = "ts")
        mod <- makeARIMA(trarma[[1]], trarma[[2]], Delta, kappa)
        if(ncxreg > 0) x <- x - xreg %*% coef[narma + (1:ncxreg)]
        arimaSS(x, mod)
        val <- .Call("ARIMA_CSS", x, arma, trarma[[1]], trarma[[2]],
                     as.integer(ncond), TRUE, PACKAGE = "ts")
        sigma2 <- val[[1]]
        var <- solve(res$hessian * n.used)
    } else {
        if(method == "CSS-ML") {
            res <- optim(init[mask], armaCSS,  method = "BFGS",
                         hessian = FALSE, control = optim.control)
            if(res$convergence == 0) init[mask] <- res$par
            ## check stationarity
            if(arma[1] > 0)
                if(!arCheck(init[1:arma[1]]))
                    stop("non-stationary AR part from CSS")
            if(arma[3] > 0)
                if(!arCheck(init[sum(arma[1:2]) + 1:arma[3]]))
                    stop("non-stationary seasonal AR part from CSS")
            ncond <- 0
        }
        if(transform.pars) {
            init <- .Call("ARIMA_Invtrans", init, arma, PACKAGE = "ts")
            ## enforce invertibility
            if(arma[2] > 0) {
                ind <- arma[1] + 1:arma[2]
                init[ind] <- maInvert(init[ind])
            }
            if(arma[4] > 0) {
                ind <- sum(arma[1:3]) + 1:arma[4]
                init[ind] <- maInvert(init[ind])
            }
        }
        trarma <- .Call("ARIMA_transPars", init, arma, transform.pars,
                        PACKAGE = "ts")
        mod <- makeARIMA(trarma[[1]], trarma[[2]], Delta, kappa)
        res <- optim(init[mask], armafn,  method = "BFGS",
                     hessian = TRUE, control = optim.control,
                     trans = as.logical(transform.pars))
        if(res$convergence > 0)
            warning(paste("possible convergence problem: optim gave code=",
                          res$convergence))
        coef[mask] <- res$par
        if(transform.pars) {
            ## enforce invertibility
            if(arma[2] > 0) {
                ind <- arma[1] + 1:arma[2]
                if(all(mask[ind]))
                    coef[ind] <- maInvert(coef[ind])
            }
            if(arma[4] > 0) {
                ind <- sum(arma[1:3]) + 1:arma[4]
                if(all(mask[ind]))
                    coef[ind] <- maInvert(coef[ind])
            }
            if(any(coef[mask] != res$par))  {  # need to re-fit
                res <- optim(coef[mask], armafn, method = "BFGS",
                             hessian = TRUE,
                             control = list(maxit = 0,
                             parscale = optim.control$parscale),
                             trans = TRUE)
                coef[mask] <- res$par
            }
            ## do it this way to ensure hessian was computed inside
            ## stationarity region
            A <- .Call("ARIMA_Gradtrans", as.double(coef), arma,
                       PACKAGE = "ts")
            A <- A[mask, mask]
            var <- t(A) %*% solve(res$hessian * n.used) %*% A
            coef <- .Call("ARIMA_undoPars", coef, arma, PACKAGE = "ts")
        } else var <- solve(res$hessian * n.used)
        trarma <- .Call("ARIMA_transPars", coef, arma, FALSE,
                        PACKAGE = "ts")
        mod <- makeARIMA(trarma[[1]], trarma[[2]], Delta, kappa)
        val <- if(ncxreg > 0)
            arimaSS(x - xreg %*% coef[narma + (1:ncxreg)], mod)
        else arimaSS(x, mod)
        sigma2 <- val[[1]][1]/n.used
    }
    value <- 2 * n.used * res$value + n.used + n.used * log(2 * pi)
    aic <- if(method != "CSS") value + 2*sum(mask) + 2 else NA
    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) {
        nm <- c(nm, cn)
        if(!orig.xreg) {
            ind <- narma + 1:ncxreg
            coef[ind] <- S$v %*% coef[ind]
            A <- diag(narma + ncxreg)
            A[ind, ind] <- S$v
            A <- A[mask, mask]
            var <- A %*% var %*% t(A)
        }
    }
    names(coef) <- nm
    dimnames(var) <- list(nm[mask], nm[mask])
    resid <- val[[2]]
    tsp(resid) <- xtsp
    class(resid) <- "ts"
    res <- list(coef = coef, sigma2 = sigma2, var.coef = var, mask = mask,
                loglik = -0.5 * value, aic = aic, arma = arma,
                residuals = resid, call = match.call(), series = series,
                code = res$convergence, n.cond = ncond, model = mod)
    class(res) <- "Arima"
    res
}


print.Arima <-
    function (x, digits = max(3, getOption("digits") - 3), se = TRUE, ...)
{
    cat("\nCall:", deparse(x$call, width = 75), "", sep = "\n")
    cat("Coefficients:\n")
    coef <- round(x$coef, digits = digits)
    if (se && nrow(x$var.coef)) {
        ses <- rep(0, length(coef))
        ses[x$mask] <- round(sqrt(diag(x$var.coef)), digits = digits)
        coef <- matrix(coef, 1, dimnames = list(NULL, names(coef)))
        coef <- rbind(coef, s.e. = ses)
    }
    print.default(coef, print.gap = 2)
    cm <- x$call$method
    if(is.null(cm) || cm != "CSS")
        cat("\nsigma^2 estimated as ", format(x$sigma2, digits = digits),
            ":  log likelihood = ", format(round(x$loglik, 2)),
            ",  aic = ", format(round(x$aic, 2)), "\n", sep = "")
    else
        cat("\nsigma^2 estimated as ",
            format(x$sigma2, digits = digits),
            ":  part log likelihood = ", format(round(x$loglik,2)),
            "\n", sep="")
    invisible(x)
}


predict.Arima <-
    function (object, n.ahead = 1, newxreg = NULL, se.fit = TRUE, ...)
{
    myNCOL <- function(x) if (is.null(x)) 0 else NCOL(x)
    rsd <- object$residuals
    xr <- object$call$xreg
    xreg <- if (!is.null(xr)) eval.parent(xr) else NULL
    ncxreg <- myNCOL(xreg)
    if (myNCOL(newxreg) != ncxreg)
        stop("xreg and newxreg have different numbers of columns")
    class(xreg) <- NULL
    xtsp <- tsp(rsd)
    n <- length(rsd)
    arma <- object$arma
    coefs <- object$coef
    narma <- sum(arma[1:4])
    if (length(coefs) > narma) {
        if (names(coefs)[narma + 1] == "intercept") {
            xreg <- cbind(intercept = rep(1, n), xreg)
            newxreg <- cbind(intercept = rep(1, n.ahead), newxreg)
            ncxreg <- ncxreg + 1
        }
        xm <- drop(as.matrix(newxreg) %*% coefs[-(1:narma)])
    }
    else xm <- 0
    if (arma[2] > 0) {
        ma <- coefs[arma[1] + 1:arma[2]]
        if (any(Mod(polyroot(c(1, ma))) < 1))
            warning("ma part of model is not invertible")
    }
    if (arma[4] > 0) {
        ma <- coefs[sum(arma[1:3]) + 1:arma[4]]
        if (any(Mod(polyroot(c(1, ma))) < 1))
            warning("seasonal ma part of model is not invertible")
    }
    z <- KalmanForecast(n.ahead, object$mod)
    pred <- ts(z[[1]] + xm, start = xtsp[2] + deltat(rsd),
               frequency = xtsp[3])
    if (se.fit) {
        se <- ts(sqrt(z[[2]] * object$sigma2),
                 start = xtsp[2] + deltat(rsd),
                 frequency = xtsp[3])
        return(pred, se)
    }
    else return(pred)
}

makeARIMA <- function(phi, theta, Delta, kappa = 1e6)
{
    p <- length(phi); q <- length(theta)
    r <- max(p, q + 1); d <- length(Delta)
    rd <- r + d
    Z <- c(1, rep(0, r-1), Delta)
    T <- matrix(0, rd, rd)
    if(p > 0) T[1:p, 1] <- phi
    if(r > 1) {
        ind <- 2:r
        T[cbind(ind-1, ind)]<- 1
    }
    if(d > 0) {
        T[r+1, ] <- Z
        if(d > 1) {
            ind <- r + 2:d
            T[cbind(ind, ind-1)]<- 1
        }
    }
    if(q < r - 1) theta <- c(theta, rep(0, r-1-q))
    R <- c(1, theta, rep(0, d))
    V <- R %o% R
    h <- 0
    a <- rep(0, rd)
    Pn <- P <- matrix(0, rd, rd)
    if(r > 1) Pn[1:r, 1:r] <- .Call("getQ0", phi, theta, PACKAGE = "ts")
    else Pn[1, 1] <- 1/(1 - phi^2)
    if(d > 0) Pn[cbind(r+1:d, r+1:d)] <- kappa
    return(phi, theta, Delta, Z, a, P, T, V, h, Pn)
}

arima0 <- function(x, order = c(0, 0, 0),
                   seasonal = list(order = c(0, 0, 0), period = NA),
                   xreg = NULL, include.mean = TRUE, delta = 0.01,
                   transform.pars = TRUE, fixed = NULL, init = NULL,
                   method = c("ML", "CSS"), n.cond,
                   optim.control = list())
{
    arma0f <- function(p)
    {
        par <- as.double(fixed)
        par[mask] <- p
        .Call("arma0fa", G, par, PACKAGE = "ts")
    }

    arCheck <- function(ar)
    {
        p <- max(which(c(1, -ar) != 0)) - 1
        if(!p) return(TRUE)
        all(Mod(polyroot(c(1, -ar[1:p]))) > 1)
    }

    maInvert <- function(ma)
    {
        ## polyroot can't cope with leading zero.
        q <- length(ma)
        q0 <- max(which(c(1,ma) != 0)) - 1
        if(!q0) return(ma)
        roots <- polyroot(c(1, ma[1:q0]))
        ind <- Mod(roots) < 1
        if(all(!ind)) return(ma)
        warning("converting non-invertible initial MA values")
        if(q0 == 1) return(c(1/ma[1], rep(0, q-q0)))
        roots[ind] <- 1/roots[ind]
        x <- 1
        for(r in roots) x <- c(x, 0) - c(0, x)/r
        c(Re(x[-1]), rep(0, q-q0))
    }

    series <- deparse(substitute(x))
    if(NCOL(x) > 1)
        stop("only implemented for univariate time series")
    method <- match.arg(method)
    x <- as.ts(x)
    dim(x) <- NULL
    n <- length(x)

    if(!missing(order))
        if(!is.numeric(order) || length(order) != 3 || any(order < 0))
            stop("`order' must be a non-negative numeric vector of length 3")
    if(!missing(seasonal))
        if(is.list(seasonal)) {
            if(is.null(seasonal$order))
                stop("`seasonal' must be a list with component `order'")
            if(!is.numeric(seasonal$order) || length(seasonal$order) != 3
               || any(seasonal$order < 0))
                stop("`seasonal$order' must be a non-negative numeric vector of length 3")
        } else if(is.numeric(order)) {
            if(length(order) == 3) seasonal <- list(order=seasonal)
            else ("`seasonal' is of the wrong length")
        } else stop("`seasonal' must be a list with component `order'")

    if(is.null(seasonal$period) || is.na(seasonal$period)
       || seasonal$period == 0) seasonal$period <- frequency(x)
    arma <- c(order[-2], seasonal$order[-2], seasonal$period,
              order[2], seasonal$order[2])
    narma <- sum(arma[1:4])
    if(d <- order[2]) x <- diff(x, 1, d)
    if(d <- seasonal$order[2]) x <- diff(x, seasonal$period, d)
    xtsp <- tsp(x)
    tsp(x) <- NULL
    nd <- order[2] + seasonal$order[2]
    n.used <- length(x)
    ncond <- n - n.used
    if(method == "CSS") {
        ncond1 <- order[1] + seasonal$period * seasonal$order[1]
        ncond <- if(!missing(n.cond)) ncond + max(n.cond, ncond1)
        else ncond + ncond1
    }
    if(is.null(xreg)) {
        ncxreg <- 0
    } else {
        if(NROW(xreg) != n) 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, n), xreg = xreg)
        ncxreg <- ncxreg + 1
    }

    if (is.null(fixed)) fixed <- rep(NA, narma + ncxreg)
    else if(length(fixed) != narma + ncxreg) stop("wrong length for fixed")
    mask <- is.na(fixed)
    if(!any(mask)) stop("all parameters were fixed")
    if(transform.pars && any(!mask[1:narma])) {
        warning("some ARMA parameters were fixed: setting transform.pars = FALSE")
        transform.pars <- FALSE
    }

    if(ncxreg) {
        if(d <- order[2]) xreg <- diff(xreg, 1, d)
        if(d <- seasonal$order[2]) xreg <- diff(xreg, seasonal$period, d)
        xreg <- as.matrix(xreg)
        if(qr(na.omit(xreg))$rank < ncol(xreg)) stop("xreg is collinear")
        if(is.null(cn <- colnames(xreg)))
            cn <- paste("xreg", 1:ncxreg, sep = "")
    }
    if(any(is.na(x)) || (ncxreg && any(is.na(xreg))))
        ## only exact recursions handle NAs
        if(method == "ML" && delta >= 0) {
            warning("NAs present: setting delta to -1")
            delta <- -1
        }

    init0 <- rep(0, narma)
    parscale <- rep(1, narma)
    if (ncxreg) {
        orig.xreg <- (ncxreg == 1) || any(!mask[narma + 1:ncxreg])
        if(!orig.xreg) {
            S <- svd(na.omit(xreg))
            xreg <- xreg %*% S$v
        }
        fit <- lm(x ~ xreg - 1, na.action = na.omit)
        init0 <- c(init0, coef(fit))
        ses <- summary(fit)$coef[,2]
        parscale <- c(parscale, ses)
    }

    storage.mode(x) <- storage.mode(xreg) <- "double"
    if(method == "CSS") transform.pars <- 0
    G <- .Call("setup_starma", as.integer(arma), x, n.used, xreg,
               ncxreg, delta, transform.pars > 0,
               ncond - (n - n.used), PACKAGE = "ts")
    on.exit(.Call("free_starma", G, PACKAGE = "ts"))

    if(!is.null(init)) {
        if(length(init) != length(init0))
            stop("`init' is of the wrong length")
        if(any(ind <- is.na(init))) init[ind] <- init0[ind]
        if(transform.pars) {
            if(any(!mask[1:narma]))
                warning("transformed ARMA parameters were fixed")
            ## check stationarity
            if(arma[1] > 0)
                if(!arCheck(init[1:arma[1]]))
                    stop("non-stationary AR part")
            if(arma[3] > 0)
                if(!arCheck(init[sum(arma[1:2]) + 1:arma[3]]))
                    stop("non-stationary seasonal AR part")
            ## enforce invertibility
            if(arma[2] > 0) {
                ind <- arma[1] + 1:arma[2]
                init[ind] <- maInvert(init[ind])
            }
            if(arma[4] > 0) {
                ind <- sum(arma[1:3]) + 1:arma[4]
                init[ind] <- maInvert(init[ind])
            }
            init <- .Call("Invtrans", G, as.double(init), PACKAGE = "ts")
        }
    } else init <- init0


    .Call("Starma_method", G, method == "CSS", PACKAGE = "ts")
    if(!("parscale" %in% names(optim.control)))
       optim.control$parscale <- parscale[mask]
    res <- optim(init[mask], arma0f, method = "BFGS",
                 hessian = TRUE, control = optim.control)
    if((code <- res$convergence) > 0)
        warning(paste("possible convergence problem: optim gave code=",
                      code))
    coef <- res$par

    if(transform.pars) {
        cf <- fixed
        cf[mask] <- coef
        ## do it this way to ensure hessian was computed inside
        ## stationarity region
        A <- .Call("Gradtrans", G, as.double(cf), PACKAGE = "ts")[mask, mask]
        var <- t(A) %*% solve(res$hessian*length(x)) %*% A
        coef <- .Call("Dotrans", G, as.double(cf), PACKAGE = "ts")[mask]
        .Call("set_trans", G, 0, PACKAGE = "ts")
    } else var <- solve(res$hessian*length(x))
    arma0f(coef)  # reset pars
    sigma2 <- .Call("get_s2", G, PACKAGE = "ts")
    resid <- .Call("get_resid", G, PACKAGE = "ts")
    tsp(resid) <- xtsp
    class(resid) <- "ts"
    n.used <- sum(!is.na(resid))
    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 = ""))
    fixed[mask] <- coef
    if(ncxreg > 0) {
        nm <- c(nm, cn)
        if(!orig.xreg) {
            ind <- narma + 1:ncxreg
            fixed[ind] <- S$v %*% fixed[ind]
            A <- diag(narma + ncxreg)
            A[ind, ind] <- S$v
            A <- A[mask, mask]
            var <- A %*% var %*% t(A)
        }
    }
    names(fixed) <- nm
    names(arma) <- c("ar", "ma", "sar", "sma", "period", "diff", "sdiff")
    dimnames(var) <- list(nm[mask], nm[mask])
    value <- 2 * n.used * res$value + n.used + n.used*log(2*pi)
    aic <- if(method != "CSS") value + 2*length(coef) + 2 else NA
    res <- list(coef = fixed, sigma2 = sigma2, var.coef = var, mask = mask,
                loglik = -0.5*value, aic = aic, arma = arma,
                residuals = resid,
                call = match.call(), series = series,
                code = code, n.cond = ncond)
    class(res) <- "arima0"
    res
}

print.arima0 <- function(x, digits = max(3, getOption("digits") - 3),
                         se = TRUE, ...)
{
    cat("\nCall:", deparse(x$call, width = 75), "", sep = "\n")
    cat("Coefficients:\n")
    coef <- round(x$coef, digits = digits)
    if (se && nrow(x$var.coef)) {
        ses <- rep(0, length(coef))
        ses[x$mask] <- round(sqrt(diag(x$var.coef)), digits = digits)
        coef <- matrix(coef, 1, dimnames = list(NULL, names(coef)))
        coef <- rbind(coef, s.e. = ses)
    }
    print.default(coef, print.gap = 2)
    cm <- x$call$method
    if(is.null(cm) || cm != "CSS")
        cat("\nsigma^2 estimated as ",
            format(x$sigma2, digits = digits),
            ":  log likelihood = ", format(round(x$loglik,2)),
            ",  aic = ", format(round(x$aic,2)),
            "\n", sep="")
    else
        cat("\nsigma^2 estimated as ",
            format(x$sigma2, digits = digits),
            ":  part log likelihood = ", format(round(x$loglik,2)),
            "\n", sep="")
    invisible(x)
}

predict.arima0 <-
    function(object, n.ahead = 1, newxreg = NULL, se.fit=TRUE, ...)
{
    myNCOL <- function(x) if(is.null(x)) 0 else NCOL(x)
    data <- eval.parent(parse(text = object$series))
    xr <- object$call$xreg
    xreg <- if(!is.null(xr)) eval.parent(xr) else NULL
    ncxreg <- myNCOL(xreg)
    if(myNCOL(newxreg) != ncxreg)
        stop("xreg and newxreg have different numbers of columns")
    class(xreg) <- NULL
    xtsp <- tsp(object$residuals)
    n <- length(data)
    arma <- object$arma
    coefs <- object$coef
    narma <- sum(arma[1:4])
    if(length(coefs) > narma) {
        if(names(coefs)[narma+1] == "intercept") {
            xreg <- cbind(intercept = rep(1, n), xreg)
            newxreg <- cbind(intercept = rep(1, n.ahead), newxreg)
            ncxreg <- ncxreg+1
        }
        data <- data - as.matrix(xreg) %*% coefs[-(1:narma)]
        xm <- drop(as.matrix(newxreg) %*% coefs[-(1:narma)])
    } else xm <- 0
    ## check invertibility of MA part(s)
    if(arma[2] > 0) {
        ma <- coefs[arma[1] + 1:arma[2]]
        if(any(Mod(polyroot(c(1, ma))) < 1))
            warning("ma part of model is not invertible")
    }
    if(arma[4] > 0) {
        ma <- coefs[sum(arma[1:3]) + 1:arma[4]]
        if(any(Mod(polyroot(c(1, ma))) < 1))
            warning("seasonal ma part of model is not invertible")
    }
    storage.mode(data) <- "double"
    G <- .Call("setup_starma", as.integer(arma), data, n, rep(0, n),
               0, -1, 0, 0, PACKAGE = "ts")
    on.exit(.Call("free_starma", G, PACKAGE = "ts"))
    .Call("Starma_method", G, TRUE, PACKAGE = "ts")
    .Call("arma0fa", G, as.double(coefs), PACKAGE = "ts")
    z <- .Call("arma0_kfore", G, arma[6], arma[7], n.ahead, PACKAGE = "ts")
    pred <- ts(z[[1]] + xm, start = xtsp[2] + deltat(data),
               frequency = xtsp[3])
    if(se.fit) {
        se <- ts(sqrt(z[[2]]),
                 start = xtsp[2] + deltat(data), frequency = xtsp[3])
        return(pred, se)
    } else return(pred)
}

arima0.diag <- function(object, gof.lag = 10)
{
    .Deprecated("tsdiag")
    ## plot standardized residuals, acf of residuals, Box-Pierce p-values
    oldpar<- par(mfrow = c(3, 1))
    on.exit(par(oldpar))
    rs <- object$resid
    stdres <- rs/sqrt(object$sigma2)
    plot(stdres, type = "h", main = "Standardized Residuals", ylab = "")
    abline(h = 0)
    acf(object$resid, plot = TRUE, main = "ACF of Residuals")
    nlag <- gof.lag
    pval <- numeric(nlag)
    for(i in 1:nlag) pval[i] <- Box.test(rs, i)$p.value
    plot(1:nlag, pval, xlab = "lag", ylab = "p value", ylim = c(0,1),
         main = "p values for Box-Pierce statistic")
    abline(h = 0.05, lty = 2, col = "blue")
}

tsdiag.Arima <- tsdiag.arima0 <- function(object, gof.lag = 10, ...)
{
    ## plot standardized residuals, acf of residuals, Ljung-Box p-values
    oldpar<- par(mfrow = c(3, 1))
    on.exit(par(oldpar))
    rs <- object$resid
    stdres <- rs/sqrt(object$sigma2)
    plot(stdres, type = "h", main = "Standardized Residuals", ylab = "")
    abline(h = 0)
    acf(object$resid, plot = TRUE, main = "ACF of Residuals",
        na.action = na.pass)
    nlag <- gof.lag
    pval <- numeric(nlag)
    for(i in 1:nlag) pval[i] <- Box.test(rs, i, type="Ljung-Box")$p.value
    plot(1:nlag, pval, xlab = "lag", ylab = "p value", ylim = c(0,1),
         main = "p values for Ljung-Box statistic")
    abline(h = 0.05, lty = 2, col = "blue")
}

tsdiag <- function(object, gof.lag, ...) UseMethod("tsdiag")
## from MASS library: (C) 1994-9 W. N. Venables and B. D. Ripley

cpgram <- function(ts, taper=0.1,
		   main=paste("Series: ", deparse(substitute(ts))), ci.col="blue")
{
    main
    if(NCOL(ts) > 1)
	stop("only implemented for univariate time series")
    x <- as.vector(ts)
    x <- x[!is.na(x)]
    x <- spec.taper(scale(x, TRUE, FALSE), p=taper)
    y <- Mod(fft(x))^2/length(x)
    y[1] <- 0
    n <- length(x)
    x <- (0:(n/2))*frequency(ts)/n
    if(length(x)%%2==0) {
	n <- length(x)-1
	y <- y[1:n]
	x <- x[1:n]
    } else y <- y[seq(along=x)]
    xm <- frequency(ts)/2
    mp <- length(x)-1
    crit <- 1.358/(sqrt(mp)+0.12+0.11/sqrt(mp))
    oldpty <- par(pty ="s")
    on.exit(par(oldpty))
    plot(x, cumsum(y)/sum(y), type="s", xlim=c(0, xm),
	 ylim=c(0, 1), xaxs="i", yaxs="i", xlab="frequency",
	 ylab="")
    lines(c(0, xm*(1-crit)), c(crit, 1), col = ci.col, lty = 2)
    lines(c(xm*crit, xm), c(0, 1-crit), col = ci.col, lty = 2)
    title(main = main)
    invisible()
}
# Copyright (C) 1997-1999  Adrian Trapletti
#

diffinv <- function (x, ...) { UseMethod("diffinv") }

diffinv.vector <- function (x, lag = 1, differences = 1,
                            xi = rep(0.0,lag*differences), ...)
{
    if (!is.vector(x)) stop ("x is not a vector")
    if (lag < 1 | differences < 1) stop ("Bad value for lag or differences")
    if (length(xi) != lag*differences) stop ("xi has not the right length")
    if (differences == 1) {
        n <- length(x)
        x <- as.vector(x,mode="double")
        y <- as.vector(numeric(n+lag))
        xi <- as.vector(xi,mode="double")
        for (i in 1:lag) y[i] <- xi[i]
        res <- .C("R_intgrt_vec", x, y=y, as.integer(lag), as.integer(n),
                  PACKAGE="ts")
        res$y
    }
    else
        diffinv.vector(diffinv.vector(x, lag, differences-1,
                                      diff(xi, lag=lag, differences=1)),
                       lag, 1, xi[1:lag])
}

diffinv.default <-
    function (x, lag = 1, differences = 1,
              xi = rep(0.0,lag*differences*dim(as.matrix(x))[2]), ...)
{
    if (is.matrix(x)) {
        n <- nrow(x)
        m <- ncol(x)
        y <- matrix(0, nr=n+lag*differences, nc=m)
        dim(xi) <- c(lag*differences, m)
        for (i in 1:m)
            y[,i] <- diffinv.vector(as.vector(x[,i]), lag, differences,
                                    as.vector(xi[,i]))
    }
    else if (is.vector(x))
        y <- diffinv.vector(x, lag, differences, xi)
    else
        stop ("x is not a vector or matrix")
    y
}

diffinv.ts <- function (x, lag = 1, differences = 1,
                       xi = rep(0.0, lag*differences*NCOL(x)), ...)
{
    if (is.ts(x) & is.null(dim(x)))
        y <- diffinv.default(as.vector(x), lag, differences, xi)
    else
        y <- diffinv.default(as.matrix(x), lag, differences, xi)
    ts(y, frequency = frequency(x), end = end(x))
}

toeplitz <- function (x)
{
    if (!is.vector(x)) stop ("x is not a vector")
    n <- length (x)
    A <- matrix (0, n, n)
    matrix (x[abs(col(A) - row(A)) + 1], n, n)
}





# Copyright (C) 1997-1999  Adrian Trapletti
#
# Rewritten to use R indexing (C) 1999 R Core Development Team
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# Library General Public License for more details.
#
# You should have received a copy of the GNU Library General Public
# License along with this library; if not, write to the Free
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

embed <- function (x, dimension = 1)
{
    if (is.matrix(x)) {
        n <- nrow(x)
        m <- ncol(x)
        if ((dimension < 1) | (dimension > n))
            stop ("wrong embedding dimension")
        y <- matrix(0.0, n - dimension + 1, dimension * m)
        for (i in (1:m))
            y[, seq(i, by = m,length = dimension)] <-
                Recall (as.vector(x[,i]), dimension)
        return (y)
    } else if (is.vector(x) || is.ts(x)) {
        n <- length (x)
        if ((dimension < 1) | (dimension > n))
            stop ("wrong embedding dimension")
        m <- n - dimension + 1
        return(matrix(x[1:m + rep(dimension:1, rep(m, dimension)) - 1], m))
    } else
        stop ("x is not a vector or matrix")
}
filter <- function(x, filter, method = c("convolution", "recursive"),
                   sides = 2, circular = FALSE, init=NULL)
{
    method <- match.arg(method)
    x <- as.ts(x)
    xtsp <- tsp(x)
    x <- as.matrix(x)
    n <- nrow(x)
    nser <- ncol(x)
    series <- colnames(x)
    nfilt <- length(filter)
    if(any(is.na(filter))) stop("missing values in filter")
    if(nfilt > n) stop("filter is longer than time series")
    y <- matrix(NA, n, nser)
    if(method == "convolution") {
        if(sides != 1 && sides != 2)
            stop("argument sides must be 1 or 2")
        for (i in 1:nser)
            y[, i] <- .C("filter1",
                         as.double(x[,i]),
                         as.integer(n),
                         as.double(filter),
                         as.integer(nfilt),
                         as.integer(sides),
                         as.integer(circular),
                         out=double(n), NAOK=TRUE, PACKAGE="ts")$out
    } else {
        if(missing(init)) {
            init <- matrix(0, nfilt, nser)
        } else {
            ni <- NROW(init)
            if(ni != nfilt)
                stop("length of init must equal length of filter")
            if(NCOL(init) != 1 && NCOL(init) != nser)
                stop(paste("init must have 1 or", nser, "cols"))
            if(!is.matrix(init)) init <- matrix(init, nfilt, nser)
        }
        for (i in 1:nser)
            y[, i] <- .C("filter2",
                         as.double(x[,i]),
                         as.integer(n),
                         as.double(filter),
                         as.integer(nfilt),
                         out=as.double(c(init[, i], double(n))),
                         NAOK=TRUE, PACKAGE="ts")$out[-(1:nfilt)]
    }
    y <- drop(y)
    tsp(y) <- xtsp
    class(y) <- if(nser > 1) c("mts", "ts") else "ts"
    y
}

# Copyright (C) 1997-1999  Adrian Trapletti
#

kernel <- function (coef, m = length(coef)+1, r, name="unknown")
{
    modified.daniell.kernel <- function (m)
    {
        if (any(m) < 0) stop ("m is negative")
        if(length(m) == 1)
            return (kernel(c(rep(1, m), 0.5)/(2*m), m,
                           name=paste("mDaniell(",m,")",sep="")))
        else {
            k <- Recall(m[1])
            for(i in 2:length(m)) k <- kernapply(k,  Recall(m[i]))
        }
        k
    }

    daniell.kernel <- function (m)
    {
        if (any(m) < 0) stop ("m is negative")
        if(length(m) == 1)
            return (kernel(rep(1/(2*m+1),m+1),m,
                           name=paste("Daniell(",m,")",sep="")))
        else {
            k <- Recall(m[1])
            for(i in 2:length(m)) k <- kernapply(k,  Recall(m[i]))
        }
        k
    }

    fejer.kernel <- function (m, r)
    {
        if (r < 1) stop ("r is less than 1")
        if (m < 1) stop ("m is less than 1")
        n <- 2*m+1
        wn <- double(m+1)
        for (j in (1:m))
        {
            wj <- 2*pi*j/n
            wn[j+1] <- sin(r*wj/2)^2/sin(wj/2)^2
        }
        wn <- wn/(n*r)
        wn[1] <- r/n
        wn <- wn/sum(c(rev(wn[2:(m+1)]),wn))
        kernel(wn, m, name=paste("Fejer(",m,",",r,")", sep=""))
    }

    dirichlet.kernel <- function (m, r)
    {
        if (r < 0) stop ("r is less than 0")
        if (m < 1) stop ("m is less than 1")
        n <- 2*m+1
        wn <- double(m+1)
        for (j in (1:m))
        {
            wj <- 2*pi*j/n
            wn[j+1] <- sin((r+0.5)*wj)/sin(wj/2)
        }
        wn <- wn/n
        wn[1] <- (2*r+1)/n
        wn <- wn/sum(c(rev(wn[2:(m+1)]),wn))
        return (kernel(wn, m, name=paste("Dirichlet(",m,",",r,")",sep="")))
    }

    if(is.character(coef)) {
        switch(coef,
               daniell = daniell.kernel(m),
               dirichlet = dirichlet.kernel(m, r),
               fejer = fejer.kernel(m, r),
               modified.daniell = modified.daniell.kernel(m),
               stop("unknown named kernel"))
    } else {
        if (!is.vector(coef))
            stop ("coef must be a vector")
        if ((length(coef) != m+1) | (length(coef) <= 0))
            stop ("coef has not the correct length")
        kernel <- list (coef=coef, m=m)
        attr(kernel, "name") <- name
        class(kernel) <- "tskernel"
        eps <- getOption("ts.eps")
        if ((sum(kernel[-m:m]) > 1.0+eps) || (sum(kernel[-m:m]) < 1.0-eps))
            stop ("coefficients do not add to 1")
        kernel
    }
}

print.tskernel <- function (x, digits = max(3,getOption("digits")-3), ...)
{
    y <- c(rev(x$coef[2:(x$m + 1)]), x$coef)
    i <- -x$m:x$m
    cat(attr(x, "name"), "\n")
    cat(paste("coef[", format(i), "] = ", format(y, digits = digits),
              sep = ""), sep = "\n")
}

plot.tskernel <- function (x, ...)
{
    y <- c(rev(x$coef[2:(x$m+1)]),x$coef)
    plot ((-x$m:x$m), y, xlab="k", ylab="W[k]", type="h", main=attr(x,"name"))
}

df.kernel <- function (k)
{
    2/sum(k[-k$m:k$m]^2)
}

bandwidth.kernel <- function (k)
{
    sqrt(sum((1/12 + (-k$m:k$m)^2) * k[-k$m:k$m]))
}


"[.tskernel" <- function (k, i)
{
    y <- c(rev(k$coef[2:(k$m+1)]), k$coef)
    y[i+k$m+1]
}

is.tskernel <- function (k)
{
    inherits(k, "tskernel")
}

kernapply <- function (x, ...)
{
    UseMethod("kernapply")
}

kernapply.vector <- function (x, k, circular = FALSE, ...)
{
    if (!is.vector(x)) stop ("x is not a vector")
    if (!is.tskernel(k)) stop ("k is not a kernel")
    if (length(x) <= 2*k$m)
        stop ("x is shorter than kernel k")
    if (k$m == 0)
        return (x)
    else
    {
        n <- length(x)
        w <- c(k[0:k$m], rep(0,n-2*k$m-1), k[-k$m:-1])
        y <- fft(fft(x)*fft(w), inverse = TRUE)/n
        if (is.numeric(x)) y <- Re(y)
        if (circular)
            return (y)
        else
            return (y[(1+k$m):(n-k$m)])
    }
}

kernapply.default <- function (x, k, circular = FALSE, ...)
{
    if (is.vector(x))
        return (kernapply.vector(x, k, circular=circular))
    else if (is.matrix(x))
        return (apply(x, MARGIN=2, FUN=kernapply, k, circular=circular))
    else
        stop ("kernapply is not available for object x")
}

kernapply.ts <- function (x, k, circular = FALSE, ...)
{
    if (!is.matrix(x))
        y <- kernapply.vector(as.vector(x), k, circular=circular)
    else
        y <- apply(x, MARGIN=2, FUN=kernapply, k, circular=circular)
    ts (y, end=end(x), frequency=frequency(x))
}

kernapply.tskernel <- function (k1, k2)
{
    if (!is.tskernel(k1))
        stop ("k1 is not a kernel")
    if (!is.tskernel(k2))
        stop ("k2 is not a kernel")
    n <- k2$m
    x <- c(rep(0,n), k1[-k1$m:k1$m], rep(0,n))
    coef <- kernapply(x, k2, circular = TRUE)
    m <- length(coef)%/%2
    kernel(coef[(m+1):length(coef)],m,
           paste("Composite(", attr(k1, "name"),",",attr(k2, "name"),")",sep=""))
}







lag <- function(x, ...) UseMethod("lag")

lag.default <- function(x, k = 1, ...)
{
    if(k != round(k)) {
        k <- round(k)
        warning("k is not an integer")
    }
    x <- hasTsp(x)
    p <- tsp(x)
    tsp(x) <- p - (k/p[3]) * c(1, 1, 0)
    x
}
## Function exists in S-plus

## Differences:
## 1) R has `type = "p"' argument
##    Idea: use "b" for n <= 10, else "p" as default, allow "text" / "labels" !
## 2) R uses `main', not `head' {consistency!}
## 3) R has `oma' and `...' args
## 4) R has  ask = par("ask") where S-plus has  ask = FALSE,
## ....

lag.plot <- function(x, lags = 1, layout = NULL, set.lags = 1:lags,
                     main = NULL, asp = 1,
                     font.main = par("font.main"), cex.main = par("cex.main"),  
                     diag = TRUE, diag.col = "gray", type = "p", oma = NULL, 
                     ask = NULL, do.lines = n <= 150, labels = do.lines, ...)
{
    xnam <- deparse(substitute(x))
    is.mat <- !is.null(ncol(x))
    nser <- ncol(x <- as.ts(as.matrix(x)))
    n <- nrow(x)
    
    if(is.null(oma)) {
        oma <- rep(2, 4)
        if (!is.null(main)) oma[3] <- oma[3] + 3*cex.main
    }
    if(missing(lags) && !missing(set.lags))
        lags <- length(set.lags <- as.integer(set.lags))
    tot.lags <- nser * lags

    if(is.null(ask))
        ask <-
            if(is.null(layout)) par("ask")## FALSE, since will have big layout
            else (interactive() && .Device != "postscript" &&
                  prod(layout) < tot.lags)
    if(is.null(layout))
        layout <-
            if(prod(pmf <- par("mfrow")) >= tot.lags) pmf
            else n2mfrow(tot.lags)

    ## Plotting
    opar <- par(mfrow = layout, mar = c(1.1,1.1, .5,.5) + is.mat*c(0,.5,0,.5),
                oma = oma, ask = ask)
    on.exit(par(opar))
    nR <- layout[1]
    nC <- layout[2]

    ii <- jj <- 0 ## current row and column in the layout
    for(i in 1:nser) {
        X <-  x[,i]
        xl <- range(X)
        nam <- if(is.mat) dimnames(x)[[2]][i] else xnam
        newX <- is.mat

        for (ll in set.lags) {
            jj <- 1 + jj %% nC
            if(jj == 1) #  new row
                ii <- 1 + ii %% nR
            ##  plot.ts(x,y) *does* a lag plot -> text, ...
            plot(lag(X,ll), X, xlim = xl, ylim = xl, asp = asp,
                 xlab = paste("lag",ll), ylab = nam, mgp = c(0,0,0),
                 axes = FALSE, type = type,
                 xy.lines = do.lines, xy.labels = labels,
                 col.lab = if(newX) "red",
                 font.lab = if(newX) 2, ...)
            box()
            if(diag) abline(c(0,1), lty = 2, col = diag.col)

            if (jj ==  1 && ii %% 2 == 1 && !newX)
                axis(2, xpd=NA)
            if (ii ==  1 && jj %% 2 == 1)
                axis(3, xpd=NA)

            do.4 <- (ii %% 2 == 0 && (jj == nC ||
                           ## very last one:
                           (i == nser && ll == set.lags[lags])))
            if (do.4) axis(4, xpd=NA)
            if (jj %% 2 == 0 && ii == nR)
                axis(1, xpd=NA)
            
            if(newX) {
                newX <- FALSE
                if(!do.4) axis(4, xpd = NA, mgp = c(0,.6,0))
            }

            if (!is.null(main) &&
                (jj == nC && ii == nR)  || ll == set.lags[lags])
                mtext(main, 3, 3, outer = TRUE, at = 0.5,
                      cex = cex.main, font = font.main)
        }
    }
    invisible(NULL)
}
monthplot <- function(x, ...) UseMethod("monthplot", x)

monthplot.StructTS <-
    function (x, labels = NULL, ylab = choice, choice = "sea", ...)
    monthplot(fitted(x)[, choice], labels = labels, ylab = ylab, ...)

monthplot.stl <-
    function (x, labels = NULL, ylab = choice, choice = "seasonal", ...)
    monthplot(x$time.series[, choice], labels = labels, ylab = ylab, ...)

monthplot.ts <-
    function (x, labels = NULL, times = time(x), phase = cycle(x),
              ylab = deparse(substitute(x)), ...)
{
    if (is.null(labels) & !missing(phase))
        return(monthplot.default(x, times = times, phase = phase,
                                 ylab = ylab, ...))
    if (is.null(labels)) {
        if (missing(phase)) {
            f <- frequency(x)
            if (f == 4) labels <- paste("Q", 1:4, sep = "")
            else if (f == 12)
                labels <- c("J", "F", "M", "A", "M", "J", "J",
                  "A", "S", "O", "N", "D")
            else labels <- 1:f
        }
    }
    monthplot.default(x, labels = labels, times = times, phase = phase,
                      ylab = ylab, ...)
}

monthplot.default <-
    function (x, labels = 1:12,
              times = 1:length(x),
              phase = (times - 1)%%length(labels) + 1,
              base = mean,
              xlim = c(0.55, f + 0.45), ylim = range(x, na.rm = TRUE),
              axes = TRUE, xlab = "", ylab = deparse(substitute(x)),
              type = "l", box = TRUE, add = FALSE, ...)
{
    if (is.null(labels) || (missing(labels) && !missing(phase))) {
        labels <- unique(phase)
        phase <- match(phase, labels)
    }
    f <- length(labels)
    if (!is.null(base))
        means <- tapply(x, phase, base)
    if (!add) {
        plot(NA, NA, axes = FALSE, xlim = xlim, ylim = ylim,
             xlab = xlab, ylab = ylab, ...)
        if (box)
            box()
        if (axes) {
            axis(1, at = 1:f, labels = labels, ...)
            axis(2, ...)
        }
        if (!is.null(base))
            segments(1:f - 0.45, means, 1:f + 0.45, means)
    }
    y <- as.numeric(times)
    scale <- 1 / diff(range(y, na.rm = TRUE)) * 0.9
    for (i in 1:f) {
        sub <- phase == i
        if (type != "h")
            lines((y[sub] - min(y)) * scale - 0.45 + i, x[sub],
                  type = type, ...)
        else segments((y[sub] - min(y)) * scale - 0.45 + i, means[i],
                      (y[sub] - min(y)) * scale - 0.45 + i, x[sub], ...)
    }
}
na.contiguous <- function(frame)
{
    tm <- time(frame)
    xfreq <- frequency(frame)
    ## use (first) maximal contiguous length of non-NAs
    if(is.matrix(frame))
        good <- apply(!is.na(frame), 1, all)
    else  good <- !is.na(frame)
    if(!sum(good)) stop("all times contain an NA")
    tt <- cumsum(!good)
    ln <- sapply(0:max(tt), function(i) sum(tt==i))
    seg <- (seq(along=ln)[ln==max(ln)])[1] - 1
    keep <- (tt == seg)
    st <- min(which(keep))
    if(!good[st]) st <- st + 1
    en <- max(which(keep))
    omit <- integer(0)
    n <- NROW(frame)
    if(st > 1) omit <- c(omit, 1:(st-1))
    if(en < n) omit <- c(omit, (en+1):n)
    cl <- class(frame)
    if(length(omit)) {
        frame <- if(is.matrix(frame)) frame[st:en,] else frame[st:en]
        attr(omit, "class") <- "omit"
        attr(frame, "na.action") <- omit
        tsp(frame) <- c(tm[st], tm[en], xfreq)
        if(!is.null(cl)) class(frame) <- cl
    }
    frame
}
## based on code by Martyn Plummer, plus kernel code by Adrian Trapletti
spectrum <- function (..., method = c("pgram", "ar"))
{
    switch(match.arg(method),
	   pgram = spec.pgram(...),
	   ar	 = spec.ar(...)
	   )
}

## spec.taper based on code by Kurt Hornik
spec.taper <- function (x, p = 0.1)
{
    if (any(p < 0) || any(p > 0.5))
        stop("p must be between 0 and 0.5")
    a <- attributes(x)
    x <- as.matrix(x)
    nc <- ncol(x)
    if (length(p) == 1)
        p <- rep(p, nc)
    else if (length(p) != nc)
        stop("length of p must be 1 or equal the number of columns of x")
    nr <- nrow(x)
    for (i in 1:nc) {
        m <- floor(nr * p[i])
        if(m == 0) next
        w <- 0.5 * (1 - cos(pi * seq(1, 2 * m - 1, by = 2)/(2 * m)))
        x[, i] <- c(w, rep(1, nr - 2 * m), rev(w)) * x[, i]
    }
    attributes(x) <- a
    x
}

spec.ar <- function(x, n.freq, order = NULL, plot = TRUE,
                    na.action = na.fail, method = "yule-walker", ...)
{
    ## can be called with a ts or a result of an AR fit.
    if(!is.list(x)) {
        series <- deparse(substitute(x))
        x <- na.action(as.ts(x))
        xfreq <- frequency(x)
        n <- NROW(x)
        nser <- NCOL(x)
        x <- ar(x, is.null(order), order, na.action=na.action, method=method)
    } else {
        cn <- match(c("ar", "var.pred", "order"), names(x))
        if(any(is.na(cn)))
            stop("x must be a time series or an ar() fit")
        series <- x$series
        xfreq <- x$frequency
        if(is.array(x$ar)) nser <- dim(x$ar)[2] else nser <- 1
    }
    n <- x$n.used
    order <- x$order
    if(missing(n.freq)) n.freq <- 500
    freq <- seq(0, 0.5, length = n.freq)
    if (nser == 1) {
        coh <- phase <- NULL
        if(order >= 1) {
            cs <- outer(freq, 1:order, function(x, y) cos(2*pi*x*y)) %*% x$ar
            sn <- outer(freq, 1:order, function(x, y) sin(2*pi*x*y)) %*% x$ar
            spec <- x$var.pred/(xfreq*((1 - cs)^2 + sn^2))
        } else
            spec <- rep(x$var.pred/(xfreq), length(freq))
    } else .NotYetImplemented()
    spg.out <- list(freq = freq, spec = spec, coh = coh, phase = phase,
                    n.used = nrow(x), series = series,
                    method = paste("AR (", order, ") spectrum ", sep="")
                    )
    class(spg.out) <- "spec"
    if(plot) {
	plot(spg.out, ci = 0, ...)
        return(invisible(spg.out))
    } else return(spg.out)
}

spec.pgram <-
    function (x, spans = NULL, kernel = NULL, taper = 0.1,
              pad = 0, fast = TRUE,
              demean = FALSE, detrend = TRUE,
              plot = TRUE, na.action = na.fail, ...)
{
    ## Estimate spectral density from (smoothed) periodogram.
    series <- deparse(substitute(x))
    x <- na.action(as.ts(x))
    xfreq <- frequency(x)
    x <- as.matrix(x)
    N <- nrow(x)
    nser <- ncol(x)
    if(!is.null(spans)) # allow user to mistake order of args
        if(is.tskernel(spans)) kernel <- spans
        else kernel <- kernel("modified.daniell", spans %/% 2)
    if(!is.null(kernel) && !is.tskernel(kernel))
        stop("must specify spans or a valid kernel")
    if (detrend) {
        t <- 1:N - (N + 1)/2
        sumt2 <- N * (N^2 - 1)/12
        for (i in 1:ncol(x))
            x[, i] <- x[, i] - mean(x[, i]) - sum(x[, i] * t) * t/sumt2
    }
    else if (demean) {
        x <- sweep(x, 2, colMeans(x))
    }
    x <- spec.taper(x, taper)
    ## to correct for tapering: Bloomfield (1976, p. 194)
    ## Total taper is taper*2
    u2 <- (1 - (5/8)*taper*2)
    u4 <- (1 - (93/128)*taper*2)
    if (pad > 0) {
        x <- rbind(x, matrix(0, nrow = N * pad, ncol = ncol(x)))
        N <- nrow(x)
    }
    NewN <- if(fast) nextn(N) else N
    x <- rbind(x, matrix(0, nrow = (NewN - N), ncol = ncol(x)))
    N <- nrow(x)
    Nspec <- floor(N/2)
    freq <- seq(from = xfreq/N, by = xfreq/N, length = Nspec)
    xfft <- mvfft(x)
    pgram <- array(NA, dim = c(N, ncol(x), ncol(x)))
    for (i in 1:ncol(x)) {
        for (j in 1:ncol(x)) {
            pgram[, i, j] <- xfft[, i] * Conj(xfft[, j])/(N*xfreq)
        }
    }
    ## value at zero is invalid as mean has been removed, so interpolate
    pgram[1, i, j] <- 0.5*(pgram[2, i, j] + pgram[N, i, j])
    if(!is.null(kernel)) {
        for (i in 1:ncol(x)) for (j in 1:ncol(x))
                pgram[, i, j] <- kernapply(pgram[, i, j], kernel,
                                           circular = TRUE)
        df <- df.kernel(kernel)/(u4/u2^2)
        bandwidth <- bandwidth.kernel(kernel) * xfreq/N
    } else {
        df <- 2/(u4/u2^2)
        bandwidth <- sqrt(1/12) * xfreq/N
    }
    pgram <- pgram[1+(1:Nspec),,, drop=FALSE]
    spec <- matrix(NA, nrow = Nspec, ncol = nser)
    for (i in 1:nser) spec[, i] <- Re(pgram[1:Nspec, i, i])
    if (nser == 1) {
        coh <- phase <- NULL
    } else {
        coh <- phase <- matrix(NA, nrow = Nspec, ncol = nser * (nser - 1)/2)
        for (i in 1:(nser - 1)) {
            for (j in (i + 1):nser) {
                coh[, i + (j - 1) * (j - 2)/2] <-
                    Mod(pgram[, i, j])^2/(spec[, i] * spec[, j])
                phase[, i + (j - 1) * (j - 2)/2] <- Arg(pgram[, i, j])
            }
        }
    }
    ## correct for tapering
    for (i in 1:nser) spec[, i] <- spec[, i]/u2
    spec <- drop(spec)
    spg.out <-
        list(freq = freq, spec = spec, coh = coh, phase = phase,
             kernel = kernel, df = df,
             bandwidth = bandwidth, n.used = nrow(x),
             series = series, snames = colnames(x),
             method = ifelse(!is.null(kernel), "Smoothed Periodogram",
                             "Raw Periodogram"),
             taper = taper, pad = pad, detrend = detrend, demean = demean)
    class(spg.out) <- "spec"
    if(plot) {
	plot(spg.out, ...)
        return(invisible(spg.out))
    } else return(spg.out)
}

plot.spec <-
    function (x, add = FALSE, ci = 0.95, log = c("yes", "dB", "no"),
              xlab = "frequency", ylab = NULL,
              type = "l", ci.col="blue", main = NULL, sub = NULL,
              plot.type = c("marginal", "coherency", "phase"), ...)
{
    spec.ci <- function (spec.obj, coverage = 0.95)
    {
        ## A utility function for plot.spec which calculates the confidence
        ## interval (centred around zero). We use a conditional argument to
        ## ensure that the ci always contains zero.

        if (coverage < 0 || coverage >= 1)
            stop("coverage probability out of range [0,1)")
        tail <- (1 - coverage)
        df <- spec.obj$df
        upper.quantile <- 1 - tail * pchisq(df, df, lower.tail = FALSE)
        lower.quantile <- tail * pchisq(df, df)
        1/(qchisq(c(upper.quantile, lower.quantile), df)/df)
    }

    plot.type <- match.arg(plot.type)
    m <- match.call()
    if(plot.type == "coherency") {
        m[[1]] <- as.name("plot.spec.coherency")
        m$plot.type <- m$log <- m$add <- NULL
        return(eval(m, parent.frame()))
    }
    if(plot.type == "phase") {
        m[[1]] <- as.name("plot.spec.phase")
        m$plot.type <- m$log <- m$add <- NULL
        return(eval(m, parent.frame()))
    }
    if(is.null(ylab))
        ylab <- if(log == "dB") "spectrum (dB)" else "spectrum"
    if(is.logical(log))
        log <- if(log) "yes" else "no"
    if(missing(log) && getOption("ts.S.compat")) log <- "dB"
    log <- match.arg(log)
    ylog <- ""
    if(log=="dB") x$spec <- 10 * log10(x$spec)
    if(log=="yes") ylog <- "y"
    if(add) {
        matplot(x$freq, x$spec, type = type, add=TRUE, ...)
    } else {
        matplot(x$freq, x$spec, xlab = xlab, ylab = ylab, type = type,
                log = ylog, ...)
        is.ar <- !is.na(pmatch("AR", x$method))
        if (ci <= 0 || log == "no" || is.ar) {
            ## No confidence limits
            ci.text <- ""
        } else {
            ## The position of the error bar has no meaning: only the width
            ## and height. It is positioned in the top right hand corner.
            ##
            conf.lim <- spec.ci(x, coverage = ci)
            if(log=="dB") {
                conf.lim <- 10*log10(conf.lim)
                conf.y <- max(x$spec) - conf.lim[2]
                conf.x <- max(x$freq) - x$bandwidth
                lines(rep(conf.x, 2), conf.y + conf.lim, col=ci.col)
                lines(conf.x + c(-0.5, 0.5) * x$bandwidth, rep(conf.y, 2),
                      col=ci.col)
                ci.text <- paste(", ", round(100*ci, 2),  "% C.I. is (",
                                 paste(format(conf.lim, digits = 3),
                                       collapse = ","), ")dB", sep="")
            } else {
                ci.text <- ""
                conf.y <- max(x$spec) / conf.lim[2]
                conf.x <- max(x$freq) - x$bandwidth
                lines(rep(conf.x, 2), conf.y * conf.lim, col=ci.col)
                lines(conf.x + c(-0.5, 0.5) * x$bandwidth, rep(conf.y, 2),
                      col=ci.col)
            }
        }
        if (is.null(main))
            main <- paste(paste("Series:", x$series), x$method, sep = "\n")
        if (is.null(sub) && !is.ar)
             sub <- paste("bandwidth = ", format(x$bandwidth, digits = 3),
                         ci.text, sep="")
        title(main = main, sub = sub)
    }
    invisible(x)
}

## based on code in Venables & Ripley
plot.spec.coherency <-
    function(x, ci = 0.95,
             xlab = "frequency", ylab = "squared coherency", ylim=c(0,1),
             type = "l", main = NULL, ci.lty = 3, ci.col="blue", ...)
{
    nser <- NCOL(x$spec)
    ## Formulae from Bloomfield (1976, p.225)
    gg <- 2/x$df
    se <- sqrt(gg/2)
    z <- -qnorm((1-ci)/2)
    if (is.null(main))
        main <- paste(paste("Series:", x$series),
                      "Squared Coherency", sep = " --  ")
    if(nser == 2) {
        plot(x$freq, x$coh, type=type, xlab=xlab, ylab=ylab, ylim=ylim, ...)
        coh <- pmin(0.99999, sqrt(x$coh))
        lines(x$freq, (tanh(atanh(coh) + z*se))^2, lty=ci.lty, col=ci.col)
        lines(x$freq, (pmax(0, tanh(atanh(coh) - z*se)))^2,
              lty=ci.lty, col=ci.col)
        title(main)
    } else {
        opar <- par(mfrow = c(nser-1, nser-1), mar = c(1.5, 1.5, 0.5, 0.5),
                    oma = c(4, 4, 6, 4))
        on.exit(par(opar))
        plot.new()
        for (j in 2:nser) for (i in 1:(j-1)) {
            par(mfg=c(j-1,i, nser-1, nser-1))
            ind <- i + (j - 1) * (j - 2)/2
            plot(x$freq, x$coh[, ind], type=type, ylim=ylim, axes=FALSE,
                 xlab="", ylab="", ...)
            coh <- pmin(0.99999, sqrt(x$coh[, ind]))
            lines(x$freq, (tanh(atanh(coh) + z*se))^2, lty=ci.lty, col=ci.col)
            lines(x$freq, (pmax(0, tanh(atanh(coh) - z*se)))^2,
                  lty=ci.lty, col=ci.col)
            box()
            if (i == 1) {
                axis(2, xpd = NA)
                title(ylab=x$snames[j], xpd = NA)
            }
            if (j == nser) {
                axis(1, xpd = NA)
                title(xlab=x$snames[i], xpd = NA)
            }
            mtext(main, 3, 3, TRUE, 0.5,
                  cex = par("cex.main"), font = par("font.main"))
        }
    }
}

plot.spec.phase <-
    function(x, ci = 0.95,
             xlab = "frequency", ylab = "phase", ylim=c(-pi, pi),
             type = "l", main = NULL, ci.lty = 3, ci.col="blue", ...)
{
    nser <- NCOL(x$spec)
    ## Formulae from Bloomfield (1976, p.225)
    gg <- 2/x$df
    if (is.null(main))
        main <- paste(paste("Series:", x$series),
                      "Phase spectrum", sep = "  -- ")
    if(nser == 2) {
        plot(x$freq, x$phase, type=type, xlab=xlab, ylab=ylab, ylim=ylim, ...)
        coh <- sqrt(x$coh)
        cl <- asin( pmin( 0.9999, qt(ci, 2/gg-2)*
                         sqrt(gg*(coh^{-2} - 1)/(2*(1-gg)) ) ) )
        lines(x$freq, x$phase + cl, lty=ci.lty, col=ci.col)
        lines(x$freq, x$phase - cl, lty=ci.lty, col=ci.col)
        title(main)
    } else {
        opar <- par(mfrow = c(nser-1, nser-1), mar = c(1.5, 1.5, 0.5, 0.5),
                    oma = c(4, 4, 6, 4))
        on.exit(par(opar))
        plot.new()
        for (j in 2:nser) for (i in 1:(j-1)) {
            par(mfg=c(j-1,i, nser-1, nser-1))
            ind <- i + (j - 1) * (j - 2)/2
            plot(x$freq, x$phase[, ind], type=type, ylim=ylim, axes=FALSE,
                 xlab="", ylab="", ...)
            coh <- sqrt(x$coh[, ind])
            cl <- asin( pmin( 0.9999, qt(ci, 2/gg-2)*
                             sqrt(gg*(coh^{-2} - 1)/(2*(1-gg)) ) ) )
            lines(x$freq, x$phase[, ind] + cl, lty=ci.lty, col=ci.col)
            lines(x$freq, x$phase[, ind] - cl, lty=ci.lty, col=ci.col)
            box()
            if (i == 1) {
                axis(2, xpd = NA)
                title(ylab=x$snames[j], xpd = NA)
            }
            if (j == nser) {
                axis(1, xpd = NA)
                title(xlab=x$snames[i], xpd = NA)
            }
            mtext(main, 3, 3, TRUE, 0.5,
                  cex = par("cex.main"), font = par("font.main"))
        }
    }
}
stl <- function(x, s.window,
		s.degree = 0,
		t.window = NULL, t.degree = 1,
		l.window = nextodd(period), l.degree = t.degree,
		s.jump = ceiling(s.window/10),
		t.jump = ceiling(t.window/10),
		l.jump = ceiling(l.window/10),
		robust = FALSE,
		inner = if(robust)  1 else 2,
		outer = if(robust) 15 else 0,
		na.action = na.fail)
{
    nextodd <- function(x){
	x <- round(x)
	if(x%%2==0) x <- x+1
	as.integer(x)
    }
    deg.check <- function(deg) {
	degname <- deparse(substitute(deg))
	deg <- as.integer(deg)
	if(deg < 0 || deg > 1) stop(paste(degname, "must be 0 or 1"))
	deg
    }
    x <- na.action(as.ts(x))
    if(is.matrix(x)) stop("only univariate series are allowed")
    n <- length(x)
    period <- frequency(x)
    if(period < 2 || n <= 2 * period)
	stop("series is not periodic or has less than two periods")
    periodic <- FALSE
    if(is.character(s.window)) {
	if(is.na(pmatch(s.window, "periodic")))
	    stop("unknown string value for s.window")
	else {
	    periodic <- TRUE
	    s.window <- 10 * n + 1
	    s.degree <- 0
	}
    }
    s.degree <- deg.check(s.degree)
    t.degree <- deg.check(t.degree)
    l.degree <- deg.check(l.degree)
    if(is.null(t.window))
	t.window <- nextodd(ceiling( 1.5 * period / (1- 1.5/s.window)))
    z <- .Fortran("stl",
		  as.double(x),
		  as.integer(n),
		  as.integer(period),
		  as.integer(s.window),
		  as.integer(t.window),
		  as.integer(l.window),
		  s.degree, t.degree, l.degree,
		  nsjump = as.integer(s.jump),
		  ntjump = as.integer(t.jump),
		  nljump = as.integer(l.jump),
		  ni = as.integer(inner),
		  no = as.integer(outer),
		  weights = double(n),
		  seasonal = double(n),
		  trend = double(n),
		  double((n+2*period)*5),
		  PACKAGE="ts")
    if(periodic) {
	## make seasonal part exactly periodic
	which.cycle <- cycle(x)
	z$seasonal <- tapply(z$seasonal, which.cycle, mean)[which.cycle]
    }
    remainder <- as.vector(x) - z$seasonal - z$trend
    y <- cbind(seasonal=z$seasonal, trend=z$trend, remainder=remainder)
    res <- list(time.series = ts(y, start=start(x), frequency = period),
		weights=z$weights, call=match.call(),
		win = c(s = s.window, t = t.window, l = l.window),
		deg = c(s = s.degree, t = t.degree, l = l.degree),
		jump= c(s = s.jump,   t = t.jump,   l = l.jump),
		inner = z$ni, outer = z$no)
    class(res) <- "stl"
    res
}

print.stl <- function(x, ...)
{
    cat(" Call:\n ")
    dput(x$call)
    cat("\nComponents\n")
    print(x$time.series, ...)
    invisible(x)
}

summary.stl <- function(object, digits = getOption("digits"), ...)
{
    cat(" Call:\n ")
    dput(object$call)
    cat("\n Time.series components:\n")
    print(summary(object$time.series, digits = digits, ...))
    cat(" IQR:\n")
    iqr <- apply(cbind(STL = object$time.series,
                       data = object$time.series %*% rep(1,3)),
		 2, IQR)
    print(rbind(format(iqr, digits = max(2, digits - 3)),
		"   %"= format(round(100 * iqr / iqr["data"], 1))),
	  quote = FALSE)
    cat("\n Weights:")
    if(all(object$weights == 1)) cat(" all == 1\n")
    else { cat("\n"); print(summary(object$weights, digits = digits, ...)) }
    cat("\n Other components: ")
    str(object[-(1:3)], give.attr = FALSE)
    invisible(object)
}

plot.stl <- function(x, labels = colnames(X),
		     set.pars = list(mar = c(0, 6, 0, 6), oma = c(6, 0, 4, 0),
				     tck = -0.01, mfrow = c(nplot, 1)),
		     main = NULL, range.bars = TRUE, ...)
{
    sers <- x$time.series
    ncomp <- ncol(sers)
    data <- drop(sers %*% rep(1, ncomp))
    X <- cbind(data, sers)
    colnames(X) <- c("data", colnames(sers))
    nplot <- ncomp + 1
    if(range.bars)
	mx <- min(apply(rx <- apply(X,2, range), 2, diff))
    if(length(set.pars)) {
	oldpar <- do.call("par", as.list(names(set.pars)))
	on.exit(par(oldpar))
	do.call("par", set.pars)
    }
    for(i in 1:nplot) {
	plot(X[, i], type = if(i < nplot) "l" else "h",
	     xlab = "", ylab = "", axes = FALSE, ...)
	if(range.bars) {
	    dx <- 1/64 * diff(ux <- par("usr")[1:2])
	    y <- mean(rx[,i])
	    rect(ux[2] - dx, y + mx/2, ux[2] - 0.4*dx, y - mx/2,
		 col = "light gray", xpd = TRUE)
	}
	if(i == 1 && !is.null(main))
	    title(main, line = 2, outer = par("oma")[3] > 0)
	if(i == nplot) abline(h=0)
	box()
	right <- i %% 2 == 0
	axis(2, labels = !right)
	axis(4, labels = right)
	axis(1, labels = i == nplot)
	mtext(labels[i], side = 2, 3)
    }
    mtext("time", side = 1, line = 3)
    invisible()
}
Box.test <- function (x, lag = 1, type=c("Box-Pierce", "Ljung-Box"))
{
    if (NCOL(x) > 1)
        stop ("x is not a vector or univariate time series")
    DNAME <- deparse(substitute(x))
    type <- match.arg(type)
    cor <- acf (x, lag.max = lag, plot = FALSE, na.action = na.pass)
    n <- sum(!is.na(x))
    PARAMETER <- lag
    obs <- cor$acf[2:(lag+1)]
    if (type=="Box-Pierce")
    {
        METHOD <- "Box-Pierce test"
        STATISTIC <- n*sum(obs^2)
        PVAL <- 1-pchisq(STATISTIC,lag)
    }
    else
    {
        METHOD <- "Box-Ljung test"
        STATISTIC <- n*(n+2)*sum(1/seq(n-1,n-lag)*obs^2)
        PVAL <- 1-pchisq(STATISTIC,lag)
    }
    names(STATISTIC) <- "X-squared"
    names(PARAMETER) <- "df"
    structure(list(statistic = STATISTIC,
                   parameter = PARAMETER,
                   p.value = PVAL,
                   method = METHOD,
                   data.name = DNAME),
              class = "htest")
}

PP.test <- function (x, lshort = TRUE)
{
    if (NCOL(x) > 1)
        stop ("x is not a vector or univariate time series")
    DNAME <- deparse(substitute(x))
    z <- embed (x, 2)
    yt <- z[,1]
    yt1 <- z[,2]
    n <- length (yt)
    tt <- (1:n)-n/2
    res <- lm (yt~1+tt+yt1)
    if (res$rank < 3)
        stop ("Singularities in regression")
    res.sum <- summary (res)
    tstat <- (res.sum$coefficients[3,1]-1)/res.sum$coefficients[3,2]
    u <- residuals (res)
    ssqru <- sum(u^2)/n
    if (lshort)
        l <- trunc(4*(n/100)^0.25)
    else
        l <- trunc(12*(n/100)^0.25)
    ssqrtl <- .C ("R_pp_sum", as.vector(u,mode="double"), as.integer(n),
                  as.integer(l), trm=as.double(ssqru), PACKAGE="ts")
    ssqrtl <- ssqrtl$trm
    n2 <- n^2
    trm1 <- n2*(n2-1)*sum(yt1^2)/12
    trm2 <- n*sum(yt1*(1:n))^2
    trm3 <- n*(n+1)*sum(yt1*(1:n))*sum(yt1)
    trm4 <- (n*(n+1)*(2*n+1)*sum(yt1)^2)/6
    Dx <- trm1-trm2+trm3-trm4
    STAT <- sqrt(ssqru)/sqrt(ssqrtl)*tstat-(n^3)/(4*sqrt(3)*sqrt(Dx)*sqrt(ssqrtl))*(ssqrtl-ssqru)
    table <- cbind(c(4.38,4.15,4.04,3.99,3.98,3.96),
                   c(3.95,3.80,3.73,3.69,3.68,3.66),
                   c(3.60,3.50,3.45,3.43,3.42,3.41),
                   c(3.24,3.18,3.15,3.13,3.13,3.12),
                   c(1.14,1.19,1.22,1.23,1.24,1.25),
                   c(0.80,0.87,0.90,0.92,0.93,0.94),
                   c(0.50,0.58,0.62,0.64,0.65,0.66),
                   c(0.15,0.24,0.28,0.31,0.32,0.33))
    table <- -table
    tablen <- dim(table)[2]
    tableT <- c(25,50,100,250,500,100000)
    tablep <- c(0.01,0.025,0.05,0.10,0.90,0.95,0.975,0.99)
    tableipl <- numeric(tablen)
    for (i in (1:tablen))
        tableipl[i] <- approx (tableT,table[,i],n,rule=2)$y
    PVAL <- approx (tableipl,tablep,STAT,rule=2)$y
    PARAMETER <- l
    METHOD <- "Phillips-Perron Unit Root Test"
    names(STAT) <- "Dickey-Fuller"
    names(PARAMETER) <- "Truncation lag parameter"
    structure(list(statistic = STAT, parameter = PARAMETER,
                   p.value = PVAL, method = METHOD, data.name = DNAME),
              class = "htest")
}
ts.plot <- function(..., gpars = list())
{
    dots <- list(...)
    pars <- c("xlab", "ylab", "xlim", "ylim", "col", "lty", "lwd",
              "type", "main", "sub", "log")
    m <- names(dots) %in% pars
    if(length(m)) {
        gpars <- c(gpars, dots[m])
        dots <- dots[!m]
    }
    sers <- do.call("ts.union", dots)
    if(is.null(gpars$ylab))
        gpars$ylab <- if(NCOL(sers) > 1) "" else deparse(substitute(...))
    do.call("plot.ts", c(list(sers, plot.type = "single"), gpars))
}

arima.sim <- function(model, n, rand.gen = rnorm,
                      innov = rand.gen(n, ...), n.start = NA, ...)
{
    if(!is.list(model)) stop("`model' must be list")
    nm <- names(model)
    p <- length(model$ar)
    if(p) {
        minroots <- min(Mod(polyroot(c(1, -model$ar))))
        if(minroots <= 1) stop("ar part of model is not stationary")
    }
    q <- length(model$ma)
    if(is.na(n.start)) n.start <- p + q +
        ifelse(p > 0, ceiling(6/log(minroots)), 0)
    if(n.start < p + q) stop("burn-in must be as long as ar + ma")
    d <- 0
    if(!is.null(ord <- model$order)) {
        if(length(ord) != 3) stop("`model$order' must be of length 3")
        d <- ord[2]
        if(d != round(d) || d < 0)
            stop("number of differences must be a positive integer")
    }
    x <- ts(c(rand.gen(n.start, ...), innov[1:n]), start = 1 - n.start)
    if(length(model$ma)) x <- filter(x, c(1, model$ma), sides = 1)
    if(length(model$ar)) x <- filter(x, model$ar, method = "recursive")
    x <- x[-(1:n.start)]
    if(d > 0) x <- diffinv(x, differences = d)
    else as.ts(x)
}
ts.union <- function(..., dframe = FALSE) {
    makeNames <- function(...) {
        l <- as.list(substitute(list(...)))[-1]
        nm <- names(l)
        fixup <- if(is.null(nm)) seq(along = l) else nm == ""
        dep <- sapply(l[fixup], function(x) deparse(x)[1])
        if(is.null(nm)) return(dep)
        if(any(fixup)) nm[fixup] <- dep
        nm
    }
    .cbind.ts(list(...), makeNames(...), dframe = dframe, union = TRUE)
}

ts.intersect <- function(..., dframe = FALSE) {
    makeNames <- function(...) {
        l <- as.list(substitute(list(...)))[-1]
        nm <- names(l)
        fixup <- if(is.null(nm)) seq(along = l) else nm == ""
        dep <- sapply(l[fixup], function(x) deparse(x)[1])
        if(is.null(nm)) return(dep)
        if(any(fixup)) nm[fixup] <- dep
        nm
    }
    .cbind.ts(list(...), makeNames(...), dframe = dframe, union = FALSE)
}
.First.lib <- function(lib, pkg)
{
    library.dynam("ts", pkg, lib)
#    if(interactive() || getOption("verbose"))
#	cat("\n	   This is a preliminary time series package for R\n",
#	    "	   See `library(help=ts)' for details.\n\n")
}
options(ts.S.compat = FALSE)
