.packageName <- "tseries"
## Copyright (C) 1997-2000  Adrian Trapletti
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2, or (at your option)
## any later version.
##
## This program 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
## General Public License for more details.
##
## A copy of the GNU General Public License is available via WWW at
## http://www.gnu.org/copyleft/gpl.html.  You can also obtain it by
## writing to the Free Software Foundation, Inc., 59 Temple Place,
## Suite 330, Boston, MA  02111-1307  USA.

##
## ARMA class
##

arma <-
function(x, order = c(1, 1), lag = NULL, coef = NULL,
         include.intercept = TRUE, series = NULL, qr.tol = 1e-07, ...)
{
    seqN <- function(N) {
        if(0==length(N)) NULL else if(N<=0) NULL else seq(N)
    }
  
    err <- function(coef) {
        u <- double(n)
        u[seqN(max.order)] <- 0
        u <- .C("arma",
                as.vector(x, mode = "double"),
                u = as.vector(u),
                as.vector(coef, mode = "double"),
                as.integer(lag$ar),
                as.integer(lag$ma),
                as.integer(ar.l),
                as.integer(ma.l),
                as.integer(max.order),
                as.integer(n),
                as.integer(include.intercept),
                PACKAGE="tseries")$u
        return(sum(u^2))
    }
  
    resid <- function(coef) {
        u <- double(n)
        u[seqN(max.order)] <- 0
        u <- .C("arma",
                as.vector(x, mode = "double"),
                u = as.vector(u),
                as.vector(coef, mode = "double"),
                as.integer(lag$ar),
                as.integer(lag$ma),
                as.integer(ar.l),
                as.integer(ma.l),
                as.integer(max.order),
                as.integer(n),
                as.integer(include.intercept),
                PACKAGE="tseries")$u
        return(u)
    }
  
    arma.init <- function() {
        k <- round(1.1*log(n))
        e <- na.omit(drop(ar.ols(x, order.max = k, aic = FALSE,
                                 demean = FALSE,
                                 intercept = include.intercept)$resid))
        ee <- embed(e, max.order+1)
        xx <- embed(x[-(1:k)], max.order+1)
        if(include.intercept == TRUE) {
            if(is.null(lag$ar))
                coef <- lm(xx[,1]~ee[,lag$ma+1])$coef
            else if(is.null(lag$ma))
                coef <- lm(xx[,1]~xx[,lag$ar+1])$coef
            else
                coef <- lm(xx[,1]~xx[,lag$ar+1]+ee[,lag$ma+1])$coef
            coef <- c(coef[-1], coef[1])
        } 
        else {
            if(is.null(lag$ar))
                coef <- lm(xx[,1]~ee[,lag$ma+1]-1)$coef
            else if(is.null(lag$ma))
                coef <- lm(xx[,1]~xx[,lag$ar+1]-1)$coef
            else
                coef <- lm(xx[,1]~xx[,lag$ar+1]+ee[,lag$ma+1]-1)$coef
        }
        return(coef) 
    }
  
    if(!is.null(order) & !is.null(lag))
        warning("order is ignored")
    if(is.null(order) & is.null(lag))
        stop("order or lag must be given")
    if(is.null(lag) & !is.null(order))
        lag <- list(ar=seqN(order[1]), ma=seqN(order[2]))
    lag$ar <- unique(lag$ar)
    lag$ma <- unique(lag$ma)
    max.order <- max(unlist(lag),0)
    ar.l <- length(lag$ar)
    ma.l <- length(lag$ma)
    if(NCOL(x) > 1)
        stop("x is not a vector or univariate time series")
    if(is.null(series)) series <- deparse(substitute(x))
    ists <- is.ts(x)
    x <- as.ts(x)
    xfreq <- frequency(x)
    if(any(is.na(x))) stop("NAs in x")
    if(ists) xtsp <- tsp(x)
    n <- length(x)
    if(!is.null(unlist(lag)))
        if(min(unlist(lag)) < 1 | max(unlist(lag)) > (n-1))
            stop("invalid lag")
    ncoef <- length(unlist(lag))+as.numeric(include.intercept)
    if(is.null(coef)) {
        if(!is.null(unlist(lag)))
            coef <- arma.init()
        else
            coef <- 0
    }
    if(length(coef) != ncoef) stop("invalid coef")
    md <- optim(coef, err, gr=NULL, hessian=TRUE, ...)
    coef <- md$par
    rank <- qr(md$hessian, qr.tol)$rank
    if(rank != ncoef) {
        se <- rep(NA, ncoef)
        cat("Warning: singular Hessian\n")
    }
    else {
        di <- diag(2*md$value/n*solve(md$hessian))
        if(any(di < 0)) 
            cat("Warning: Hessian negative-semidefinite\n")
        se <- sqrt(di)
    }
    e <- resid(coef)
    e[seqN(max.order)] <- NA
    f <- x-e
    if(ists) {
        attr(e, "tsp") <- xtsp
        attr(e, "class") <- "ts"
        attr(f, "tsp") <- xtsp
        attr(f, "class") <- "ts"
    }
    nam.ar <- if(!is.null(lag$ar))
        paste("ar", lag$ar, sep = "")
    else
        NULL
    nam.ma <- if(!is.null(lag$ma))
        paste("ma", lag$ma, sep = "")
    else
        NULL
    nam.int <- if(include.intercept) "intercept" else NULL
    nam.coef <- c(nam.ar, nam.ma, nam.int)
    names(coef) <- nam.coef
    names(se) <- nam.coef
    arma <- list(coef = coef,
                 css = md$value,
                 n.used = n,
                 residuals = e,
                 fitted.values = f,
                 series = series,
                 frequency = xfreq,
                 call = match.call(),
                 asy.se.coef = se,
                 lag = lag,
                 convergence = md$convergence,
                 include.intercept = include.intercept)
    class(arma) <- "arma"
    return(arma)
}

coef.arma <-
function(object, ...)
{
    if(!inherits(object, "arma"))
        stop("method is only for arma objects")
    return(object$coef)
}

residuals.arma <-
function(object, ...)
{
    if(!inherits(object, "arma"))
        stop("method is only for arma objects")
    return(object$residuals)
}

fitted.arma <-
function(object, ...)
{
    if(!inherits(object, "arma"))
        stop("method is only for arma objects")
    return(object$fitted.values)
}

print.arma <-
function(x, digits = max(3, getOption("digits") - 3), ...)
{
    if(!inherits(x, "arma"))
        stop("method is only for arma objects")
    cat("\nCall:\n", deparse(x$call), "\n\n", sep = "")
    cat("Coefficient(s):\n")
    print.default(format(coef(x), digits = digits),
                  print.gap = 2, quote = FALSE)
    cat("\n")
    invisible(x)
}

summary.arma <-
function(object, ...)
{
    if(!inherits(object, "arma"))
        stop("method is only for arma objects")
    ans <- NULL
    ans$residuals <- na.remove(object$residuals)
    tval <- object$coef / object$asy.se.coef
    ans$coef <- cbind(object$coef, object$asy.se.coef, tval,
                      2 * (1-pnorm(abs(tval))))
    dimnames(ans$coef) <-
        list(names(object$coef),
             c(" Estimate"," Std. Error"," t value","Pr(>|t|)"))
    ans$call <- object$call
    ans$nn <- object$nn
    ans$css <- object$css
    ans$var <- var(ans$residuals)
    ans$aic <- (object$n.used * (1+log(2*pi)) + object$n.used *
                log(ans$var) + 2 * length(object$coef))
    ans$p <- max(object$lag$ar, 0)
    ans$q <- max(object$lag$ma, 0)
    class(ans) <- "summary.arma"
    return(ans)
}

print.summary.arma <-
function(x, digits = max(3, getOption("digits") - 3),
         signif.stars = getOption("show.signif.stars"), ...)
{
    if(!inherits(x, "summary.arma"))
        stop("method is only for summary.arma objects")
    cat("\nCall:\n", deparse(x$call), "\n", sep = "")
    cat("\nModel:\nARMA(",x$p,",",x$q,")\n", sep = "")
    cat("\nResiduals:\n")
    rq <- structure(quantile(x$residuals),
                    names = c("Min","1Q","Median","3Q","Max"))
    print(rq, digits = digits, ...)
    cat("\nCoefficient(s):\n")
    printCoefmat(x$coef, digits = digits,
                 signif.stars = signif.stars, ...)
    cat("\nFit:\n")
    cat("sigma^2 estimated as ", format(x$var, digits = digits), 
        ",  Conditional Sum-of-Squares = ", format(round(x$css, 2)), 
        ",  AIC = ", format(round(x$aic, 2)), "\n", sep = "")
    cat("\n")
    invisible(x)
}

plot.arma <-
function(x, ask = interactive(), ...)
{
    if(!inherits(x, "arma"))
        stop("method is only for arma objects")
    op <- par()
    par(ask = ask, mfrow = c(2, 1))
    data <- eval.parent(parse(text = x$series))
    if(any(is.na(data))) stop(paste("NAs in", x$series))
    plot(data, main = x$series, ylab = "Series")
    plot(x$residuals, main = "Residuals", ylab = "Series")
    acf(data, main = paste("ACF of", x$series))
    acf(x$residuals, main = "ACF of Residuals",
        na.action = na.remove)
    pacf(data, main = paste("PACF of", x$series))
    pacf(x$residuals, main = "PACF of Residuals",
         na.action = na.remove)
    par(ask = op$ask, mfrow = op$mfrow)
    invisible(x)
}
## Copyright (C) 1997-2003  Adrian Trapletti
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2, or (at your option)
## any later version.
##
## This program 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
## General Public License for more details.
##
## A copy of the GNU General Public License is available via WWW at
## http://www.gnu.org/copyleft/gpl.html.  You can also obtain it by
## writing to the Free Software Foundation, Inc., 59 Temple Place,
## Suite 330, Boston, MA  02111-1307  USA.

##
## Financial time series analysis
##

portfolio.optim <- function (x, ...) UseMethod ("portfolio.optim")

portfolio.optim.ts <-
function (x, ...)
{
    if(!is.ts(x))
        stop("method is only for time series")
    if(NCOL(x) == 1)
        stop("x is not a multivariate time series")
    res <- portfolio.optim.default(as.matrix(x), ...)
    res$px <- ts(res$px, start = start(x), frequency = frequency(x))
    return(res)
}

portfolio.optim.default <-
function(x, pm = mean(x), riskless = FALSE, shorts = FALSE,
         rf = 0.0, reslow = NULL, reshigh = NULL, covmat = cov(x), ...) 
{
    if(!require("quadprog", quietly=TRUE))
        stop("package", sQuote("quadprog"), "is needed.  Stopping")
    if(NCOL(x) == 1)
        stop("x is not a matrix")
    if(any(is.na(x)))
        stop("NAs in x")
    k <- dim(x)[2]
    if(!is.matrix(covmat)) {
        stop("covmat is not a matrix")
    }
    if((dim(covmat)[1] !=k) | (dim(covmat)[2] !=k)) {
      stop("covmat has not the right dimension")
    }
    Dmat <- covmat
    dvec <- rep(0, k)
    big <- 1e+100
    if(!is.null(reslow) & is.null(reshigh)) {
        reshigh <- rep(big, k)
    }
    if(is.null(reslow) & !is.null(reshigh)) {
        reslow <- -rep(big, k)
    }
    if(!is.null(reslow)) {
        if(!is.vector(reslow)) {
            stop("reslow is not a vector")
        }
        if(length(reslow) != k) {
            stop("reslow has not the right dimension")
        }
    }
    if(!is.null(reshigh)) {
        if(!is.vector(reshigh)) {
            stop("reshigh is not a vector")
        }
        if(length(reshigh) != k) {
            stop("reshigh has not the right dimension")
        }
    }
    if(riskless) {
        a1 <- apply(x, 2, mean)-rf
        if(shorts) {
            a2 <- NULL
            b2 <- NULL
        }
        else {
            a2 <- matrix(0, k, k)
            diag(a2) <- 1
            b2 <- rep(0, k)
        }
        if(!is.null(reslow) & !is.null(reshigh)) {
            a3 <- matrix(0, k, k)
            diag(a3) <- 1
            Amat <- t(rbind(a1, a2, a3, -a3))
            b0 <- c(pm-rf, b2, reslow, -reshigh)
        }
        else {
            Amat <- t(rbind(a1, a2))
            b0 <- c(pm-rf, b2)
        }
        res <- solve.QP(Dmat, dvec, Amat, bvec=b0, meq=1)
    }
    else {
        a1 <- rep(1, k)
        a2 <- apply(x, 2, mean)
        if(shorts) {
            if(!is.null(reslow) & !is.null(reshigh)) {
                a3 <- matrix(0, k, k)
                diag(a3) <- 1
                Amat <- t(rbind(a1, a2, a3, -a3))
                b0 <- c(1, pm, reslow, -reshigh)
            }
            else {
                Amat <- t(rbind(a1, a2))
                b0 <- c(1, pm)
            }              
        }
        else {
            a3 <- matrix(0, k, k)
            diag(a3) <- 1
            b3 <- rep(0, k)
            if(!is.null(reslow) & !is.null(reshigh)) {
                Amat <- t(rbind(a1, a2, a3, a3, -a3))
                b0 <- c(1, pm, b3, reslow, -reshigh)
            }
            else {
                Amat <- t(rbind(a1, a2, a3))
                b0 <- c(1, pm, b3)
            }
        }
        res <- solve.QP(Dmat, dvec, Amat, bvec=b0, meq=2)
    }
    y <- t(res$solution%*%t(x))
    ans <- list(pw = res$solution, px = y, pm = mean(y), ps = sd(y))
    return(ans)
}

get.hist.quote <-
function (instrument = "^gdax", start, end,
          quote = c("Open", "High", "Low", "Close"),
          provider = "yahoo", method = "auto", origin = "1899-12-30") 
{
    if(missing(start)) start <- "1991-01-02"
    if(missing(end)) end <- format(Sys.time() - 86400, "%Y-%m-%d")
  
    provider <- match.arg(provider)

    start <- as.POSIXct(start, tz = "GMT")
    end <- as.POSIXct(end, tz = "GMT")

    if(provider == "yahoo") {
        url <-
            paste("http://chart.yahoo.com/table.csv?s=",
                  instrument, 
                  format(start,
                         paste("&a=",
                               as.character(as.numeric(format(start, "%m"))-1),
                               "&b=%d&c=%Y",
                               sep = "")),
                  format(end,
                         paste("&d=",
                               as.character(as.numeric(format(end, "%m"))-1),
                               "&e=%d&f=%Y",
                               sep = "")), 
                  "&g=d&q=q&y=0&z=",
                  instrument,
                  "&x=.csv",
                  sep = "")
        destfile <- tempfile()
        status <- download.file(url, destfile, method = method)
        if(status != 0) {
            unlink(destfile)
            stop(paste("download error, status", status))
        }
        nlines <- length(count.fields(destfile, sep = "\n"))
        if(nlines == 1) {
            unlink(destfile)
            stop(paste("No data available for", instrument))
        }
        
        ## Yahoo includes rows concerning dividends,
        ## hence need fill = TRUE and na.omit
        x <- read.table(destfile, header = TRUE, sep = ",", as.is = TRUE, fill = TRUE)
        x <- na.omit(x)
        
        ## Debug
        ## cat("read.table: start =", x[NROW(x),"Date"], "\n")
        ## cat("read.table: end   =", x[1,"Date"], "\n")
        
        unlink(destfile)

        nser <- pmatch(quote, names(x)[-1]) + 1
        if(any(is.na(nser)))
            stop("This quote is not available")
        n <- nrow(x)

        ## Yahoo currently formats dates as '26-Jun-01', hence need C
        ## LC_TIME locale for getting the month right.
        lct <- Sys.getlocale("LC_TIME")
        Sys.setlocale("LC_TIME", "C")
        on.exit(Sys.setlocale("LC_TIME", lct))

        dat <- gsub(" ", "0", as.character(x[, 1])) # Need the gsub?
        dat <- as.POSIXct(strptime(dat, "%d-%b-%y"), tz = "GMT")
        if(dat[n] != start)
            cat(format(dat[n], "time series starts %Y-%m-%d\n"))
        if(dat[1] != end)
            cat(format(dat[1], "time series ends   %Y-%m-%d\n"))
        jdat <-
            unclass(julian(dat, origin = as.POSIXct(origin, tz = "GMT")))
        ## We need unclass() because 1.7.0 does not allow adding a number
        ## to a "difftime" object. 
        ind <- jdat - jdat[n] + 1
        y <- matrix(NA, nr = max(ind), nc = length(nser))
        y[ind, ] <- as.matrix(x[, nser, drop = FALSE])
        colnames(y) <- names(x)[nser]
        return(ts(y, start = jdat[n], end = jdat[1]))
    }
    else stop("provider not implemented")
}

maxdrawdown <-
function(x)
{
    if(NCOL(x) > 1)
        stop("x is not a vector or univariate time series")
    if(any(is.na(x)))
        stop("NAs in x")
    cmaxx <- cummax(x)-x
    mdd <- max(cmaxx)
    to <- which(mdd == cmaxx)
    from <- double(NROW(to))
    for (i in 1:NROW(to))
        from[i] <- max(which(cmaxx[1:to[i]] == 0))
    return(list(maxdrawdown = mdd, from = from, to = to))
}

sharpe <-
function(x, r = 0, scale = sqrt(250))
{
    if(NCOL(x) > 1)
        stop("x is not a vector or univariate time series")
    if(any(is.na(x)))
        stop("NAs in x")
    if(NROW(x) == 1)
        return(NA)
    else {
        y <- diff(x)
        return(scale * (mean(y)-r)/sd(y))
    }
}

sterling <-
function(x)
{
    if(NCOL(x) > 1)
        stop("x is not a vector or univariate time series")
    if(any(is.na(x)))
        stop("NAs in x")
    if(NROW(x) == 1)
        return(NA)
    else {
        return((x[NROW(x)]-x[1]) / maxdrawdown(x)$maxdrawdown)
    }
}

plotOHLC <-
function(x, xlim = NULL, ylim = NULL, xlab = "Time", ylab,
         col = par("col"), bg = par("bg"), axes = TRUE,
         frame.plot = axes, ann = par("ann"), main = NULL,
         date = c("calendar", "julian"), format = "%Y-%m-%d",
         origin = "1899-12-30", ...)
{
  if ((!is.mts(x)) ||
      (colnames(x)[1] != "Open") ||
      (colnames(x)[2] != "High") ||
      (colnames(x)[3] != "Low") ||
      (colnames(x)[4] != "Close"))
      stop("x is not a open/high/low/close time series")
  xlabel <- if (!missing(x)) 
      deparse(substitute(x))
  else NULL
  if (missing(ylab)) 
      ylab <- xlabel
  date <- match.arg(date)
  time.x <- time(x)
  dt <- min(lag(time.x)-time.x)/3
  if (is.null(xlim)) 
      xlim <- range(time.x)
  if (is.null(ylim)) 
      ylim <- range(x[is.finite(x)])
  plot.new()
  plot.window(xlim, ylim, ...)
  segments(time.x, x[, "High"], time.x, x[, "Low"], col = col[1], bg = bg)
  segments(time.x - dt, x[, "Open"], time.x, x[, "Open"],
           col = col[1], bg = bg)
  segments(time.x, x[, "Close"], time.x + dt, x[, "Close"],
           col = col[1], bg = bg)
  if (ann) 
      title(main = main, xlab = xlab, ylab = ylab, ...)  
  if (axes) {
      if (date == "julian") {
          axis(1, ...)
          axis(2, ...)
      }
      else {
          n <- NROW(x)
          lab.ind <- round(seq(1, n, length=5))
          D <- as.vector(time.x[lab.ind]*86400) + as.POSIXct(origin, tz = "GMT")
          DD <- format.POSIXct(D, format = format, tz ="GMT")
          axis(1, at=time.x[lab.ind], lab=DD, ...)
          axis(2, ...)
      }
  }
  if (frame.plot) 
      box(...)
}
## Copyright (C) 1997-1999  Adrian Trapletti
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2, or (at your option)
## any later version.
##
## This program 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
## General Public License for more details.
##
## A copy of the GNU General Public License is available via WWW at
## http://www.gnu.org/copyleft/gpl.html.  You can also obtain it by
## writing to the Free Software Foundation, Inc., 59 Temple Place,
## Suite 330, Boston, MA  02111-1307  USA.

##
## GARCH class
##

garch <-
function (x, order = c(1, 1), coef = NULL, itmax = 200, eps = NULL,
          grad = c("analytical","numerical"), series = NULL,
          trace = TRUE, ...)
{
    if(NCOL(x) > 1)
        stop("x is not a vector or univariate time series")
    if(!is.vector(order)) stop("order is not a vector")
    grad <- match.arg(grad)
    switch(grad,
           analytical = (agrad <- TRUE),
           numerical = (agrad <- FALSE))
    if(is.null(series)) series <- deparse(substitute(x))
    ists <- is.ts(x)
    x <- 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 <- nrow(x)
    e <- double(n)
    ncoef <- order[1]+order[2]+1
    hess <- matrix(0.0, ncoef, ncoef)
    small <- 0.05
    if(is.null(coef))
        coef <- c(var(x)*(1.0-small*(ncoef-1)),rep(small,ncoef-1))
    if(!is.vector(coef)) stop("coef is not a vector")
    if(ncoef != length(coef)) stop("incorrect length of coef")
    if(is.null(eps)) eps <- .Machine$double.eps
    nlikeli <- 1.0e+10
    fit <- .C("fit_garch",
              as.vector(x, mode = "double"),
              as.integer(n),
              coef = as.vector(coef, mode = "double"),
              as.integer(order[1]),
              as.integer(order[2]),
              as.integer(itmax),
              as.double(eps),
              nlikeli = as.double(nlikeli),
              as.integer(agrad),
              as.integer(trace),
              PACKAGE="tseries")
    pred <- .C("pred_garch",
               as.vector(x, mode = "double"),
               e = as.vector(e, mode = "double"),
               as.integer(n),
               as.vector(fit$coef, mode = "double"),
               as.integer(order[1]),
               as.integer(order[2]),
               as.integer(FALSE),
               PACKAGE = "tseries")
    com.hess <- .C("ophess_garch",
                   as.vector(x, mode = "double"),
                   as.integer(n),
                   as.vector(fit$coef, mode = "double"),
                   hess = as.matrix(hess),
                   as.integer(order[1]),
                   as.integer(order[2]),
                   PACKAGE="tseries")
    rank <- qr(com.hess$hess, ...)$rank
    if(rank != ncoef) {
        se.garch <- rep(NA, ncoef)
        cat("Warning: singular information\n")
    }
    else
        se.garch <- sqrt(diag(solve(com.hess$hess)))
    sigt <- sqrt(pred$e)
    sigt[1:max(order[1],order[2])] <- rep(NA, max(order[1],order[2]))
    f <- cbind(sigt,-sigt)
    colnames(f) <- c("sigt","-sigt")
    e <- as.vector(x)/sigt  
    if(ists) {
        attr(e, "tsp") <-  attr(f, "tsp") <- xtsp
        attr(e, "class") <- attr(f, "class") <- "ts"
    }
    names(order) <- c("p","q")
    coef <- fit$coef
    nam.coef <- "a0"
    if(order[2] > 0)
        nam.coef <- c(nam.coef, paste("a", seq(order[2]), sep = ""))
    if(order[1] > 0)
        nam.coef <- c(nam.coef, paste("b", seq(order[1]), sep = ""))
    names(coef) <- nam.coef
    names(se.garch) <- nam.coef
    garch <- list(order = order,
                  coef = coef,
                  n.likeli = fit$nlikeli,
                  n.used = n,
                  residuals = e,
                  fitted.values = f,
                  series = series, 
                  frequency = xfreq,
                  call = match.call(),
                  asy.se.coef = se.garch)
    class(garch) <- "garch"
    return(garch)
}

coef.garch <-
function(object, ...)
{
    if(!inherits(object, "garch"))
        stop("method is only for garch objects")
    return(object$coef)
}

residuals.garch <-
function(object, ...)
{
    if(!inherits(object, "garch"))
        stop("method is only for garch objects")
    return(object$residuals)
}

fitted.garch <-
function(object, ...)
{
    if(!inherits(object, "garch"))
        stop("method is only for garch objects")
    return(object$fitted.values)
}

print.garch <-
function(x, digits = max(3, getOption("digits") - 3), ...)
{
    if(!inherits(x, "garch"))
        stop("method is only for garch objects")
    cat("\nCall:\n", deparse(x$call), "\n\n", sep = "")
    cat("Coefficient(s):\n")
    print.default(format(coef(x), digits = digits), print.gap = 2,
                  quote = FALSE)
    cat("\n")
    invisible(x)
}

summary.garch <-
function(object, ...)
{
    if(!inherits(object, "garch"))
        stop("method is only for garch objects")
    ans <- NULL
    ans$residuals <- na.remove(object$residuals)
    tval <- object$coef / object$asy.se.coef
    ans$coef <- cbind(object$coef, object$asy.se.coef, tval,
                      2*(1-pnorm(abs(tval))))
    dimnames(ans$coef) <-
        list(names(object$coef),
             c(" Estimate"," Std. Error"," t value","Pr(>|t|)"))
    ans$call <- object$call
    ans$order <- object$order
    Residuals <- ans$residuals
    ans$j.b.test <- jarque.bera.test(Residuals)
    Squared.Residuals <- ans$residuals^2
    ans$l.b.test <- Box.test(Squared.Residuals, type = "Ljung-Box")
    class(ans) <- "summary.garch"
    return(ans)
}

plot.garch <- function(x, ask = interactive(), ...)
{
    if(!inherits(x, "garch"))
        stop("method is only for garch objects")
    op <- par()
    par(ask = ask, mfrow = c(2,1))
    data <- eval.parent(parse(text=x$series))
    if(any(is.na(data))) stop(paste("NAs in", x$series))
    plot(data, main = x$series, ylab = "Series")
    plot(x$residuals, main = "Residuals", ylab = "Series")
    hist(data, main = paste("Histogram of", x$series), xlab = "Series")
    hist(x$residuals, main = "Histogram of Residuals", xlab = "Series")
    qqnorm(data, main = paste("Q-Q Plot of", x$series),
           xlab = "Normal Quantiles")
    qqnorm(x$residuals, main = "Q-Q Plot of Residuals",
           xlab = "Normal Quantiles")
    acf(data^2, main = paste("ACF of Squared", x$series))
    acf(x$residuals^2, main = "ACF of Squared Residuals",
        na.action = na.remove)
    par(ask = op$ask, mfrow = op$mfrow)
    invisible(x)
}

print.summary.garch <-
function(x, digits = max(3, getOption("digits") - 3),
         signif.stars = getOption("show.signif.stars"), ...)
{
    if(!inherits(x, "summary.garch"))
        stop("method is only for summary.garch objects")
    cat("\nCall:\n", deparse(x$call), "\n", sep = "")
    cat("\nModel:\nGARCH(", x$order[1], ",", x$order[2], ")", "\n",
        sep = "")
    cat("\nResiduals:\n")
    rq <- structure(quantile(x$residuals),
                    names = c("Min","1Q","Median","3Q","Max"))
    print(rq, digits = digits, ...)
    cat("\nCoefficient(s):\n")
    printCoefmat(x$coef, digits = digits,
                 signif.stars = signif.stars, ...)
    cat("\nDiagnostic Tests:")
    print(x$j.b.test)
    print(x$l.b.test)
    invisible(x)
}

predict.garch <-
function(object, newdata, genuine = FALSE, ...)
{
    if(!inherits(object, "garch"))
        stop("method is only for garch objects")
    if(missing(newdata)) {
        newdata <- eval.parent(parse(text=object$series))
        if(any(is.na(newdata))) stop("NAs in newdata")
    }
    if(NCOL(newdata) > 1)
        stop("newdata is not a vector or univariate time series")
    ists <- is.ts(newdata)
    if(ists) newdata.tsp <- tsp(newdata)
    newdata <- as.matrix(newdata)
    n <- nrow(newdata)
    if(genuine) h <- double(n+1)
    else h <- double(n)
    pred <- .C("pred_garch",
               as.vector(newdata, mode = "double"),
               h = as.vector(h, mode = "double"),
               as.integer(n),
               as.vector(object$coef, mode = "double"),
               as.integer(object$order[1]),
               as.integer(object$order[2]),
               as.integer(genuine),
               PACKAGE="tseries")
    pred$h <- sqrt(pred$h)
    pred$h[1:max(object$order[1],object$order[2])] <-
        rep(NA, max(object$order[1],object$order[2]))
    pred$h <- cbind(pred$h,-pred$h)
    if(ists) {
        attr(pred$h, "tsp") <-
            if(genuine)
                c(newdata.tsp[1],
                  newdata.tsp[2] + 1 / newdata.tsp[3],
                  newdata.tsp[3])
            else
                newdata.tsp
        attr(pred$h, "class") <- "ts"
    }
    return(pred$h)
}

logLik.garch <-
function(object, ...)
{
    if(!inherits(object, "garch"))
        stop("method is only for garch objects")
    n <- length(na.remove(object$residuals))
    val <- (-object$n.likeli) - 0.5*n*log(2*pi)
    attr(val, "df") <- length(object$coef)
    class(val) <- "logLik"
    return(val)
}

## Copyright (C) 1997-2003  Adrian Trapletti
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2, or (at your option)
## any later version.
##
## This program 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
## General Public License for more details.
##
## A copy of the GNU General Public License is available via WWW at
## http://www.gnu.org/copyleft/gpl.html.  You can also obtain it by
## writing to the Free Software Foundation, Inc., 59 Temple Place,
## Suite 330, Boston, MA  02111-1307  USA.

##
## Irregular time series objects
##

value <- function (x, ...) UseMethod ("value")
                    
irts <-
function(time, value)
{
    if(inherits(time, "POSIXct")) {
        time <- as.numeric(time)
    }
    if(!is.vector(time))
        stop("time is not a vector")
    if(!is.vector(value) & !is.matrix(value))
        stop("value is not a vector and not a matrix")
    if(length(time) != NROW(value))
        stop("time and value have not the same number of rows")
    class(time) <- c("POSIXt", "POSIXct")
    irts <- list(time = time, value = value)
    class(irts) <- "irts"
    return(irts)
}

is.irts <-
function(object)
{
    return(inherits(object, "irts"))
}

as.irts <-
function(object)
{
    return(irts(object[,1], object[,-1]))
}

value.irts <-
function(x, ...)
{
    if(!inherits(x, "irts"))
        stop("method is only for irts objects")
    return(x$value)
}

time.irts <-
function(x, ...)
{
    if(!inherits(x, "irts"))
        stop("method is only for irts objects")
    return(x$time)
}

print.irts <-
function(x, format = "%Y-%m-%d %H:%M:%S", tz = "GMT",
         usetz = TRUE, format.value = NULL, ...)
{
    if(!inherits(x, "irts"))
        stop("method is only for irts objects")
    n <- length(x$time)
    for(i in 1:n) {
        cat(format(x$time[i], format = format, tz = tz, usetz = usetz))
        cat(" ")
        if(is.vector(x$value))
            cat(formatC(x$value[i], format = format.value, ...))
        else
            cat(formatC(x$value[i,], format = format.value, ...))
        cat("\n")
    }
    invisible(x)
}

read.irts <-
function(file, format = "%Y-%m-%d %H:%M:%S", tz = "GMT", ...)
{
    seqN <- function(from, to) {
        if((0 == length(from)) | (0 == length(to)))
            NULL
        else if(to-from+1 <= 0) 
            NULL
        else seq(from, to)
    }
    
    data <- read.table(file, as.is = TRUE, ...)
    n <- length(unlist(strsplit(format, split = " ")))
    tmp <- data[,1]
    j <- 2
    while(j <= n) {
        tmp <- paste(tmp, data[,j])
        j <- j+1
    }
    time <- as.numeric(as.POSIXct(strptime(tmp, format = format), tz = tz))
    value <- as.matrix(data[,-seqN(1, n)])
    return(irts(time, value[,,drop = TRUE]))
}

write.irts <-
function(object, file = "", append = FALSE, quote = FALSE, sep = " ", eol = "\n",
         na = "NA", dec = ".", row.names = FALSE, col.names = FALSE, qmethod = "escape",
         format = "%Y-%m-%d %H:%M:%S", tz = "GMT", usetz = FALSE, format.value = NULL, ...)
{
    dataframe <- data.frame(time = format(object$time, format = format, tz = tz, usetz = usetz),
                            value = formatC(object$value, format = format.value, ...))
    write.table(dataframe, file = file, append = append, quote = quote, sep = sep, eol = eol,
                na = na, dec = dec, row.names = row.names, col.names = col.names, qmethod = qmethod)
    invisible(object)
}

weekday <-
function(object, tz = "GMT")
{
    if(!inherits(object, "irts"))
        stop("function is only for irts objects")
    return(as.POSIXlt(object$time, tz = tz)$wday)
}

daysecond <-
function(object, tz = "GMT")
{
    if(!inherits(object, "irts"))
        stop("function is only for irts objects")
    hour <- as.POSIXlt(object$time, tz = tz)$hour
    min <- as.POSIXlt(object$time, tz = tz)$min
    sec <- as.POSIXlt(object$time, tz = tz)$sec
    return(3600*hour+60*min+sec)
}

is.businessday <-
function(object, tz = "GMT")
{
    if(!inherits(object, "irts"))
        stop("function is only for irts objects")
    wday <- as.POSIXlt(object$time, tz = tz)$wday
    return((0 < wday) & (wday < 6))
}

is.weekend <- function(object, tz = "GMT")
{
    if(!inherits(object, "irts"))
        stop("function is only for irts objects")
    wday <- as.POSIXlt(object$time, tz = tz)$wday
    return((0 == wday) | (wday == 6))
}

"[.irts" <-
function(x, i, j, ...)
{
    if(!inherits(x, "irts"))
        stop("method is only for irts objects")
    if(is.vector(x$value)) {
        if(nargs() > 2) {   
            stop("incorrect number of dimensions")
        }
        if(missing(i)) {
            return(x)
        } else {
            return(irts(as.numeric(x$time)[i], x$value[i]))
        }
    } else {
        if(missing(i)) {
            if(missing(j)) {
                return(x)
            } else {
                return(irts(as.numeric(x$time), x$value[,j,drop = FALSE]))
            }
        } else {
            if(missing(j)) {
                return(irts(as.numeric(x$time)[i], x$value[i,,drop = FALSE]))
            } else {
                return(irts(as.numeric(x$time)[i], x$value[i,j,drop = FALSE]))
            }
        }
    }
}

approx.irts <-
function(object, time, ...)
{
    if(!inherits(object, "irts"))
        stop("function is only for irts objects")
    if(!inherits(time, "POSIXct"))
        stop("time is not of class POSIXct")
    ovalue <- as.matrix(object$value)
    otime <- as.numeric(object$time)
    time <- as.numeric(time)
    value <- matrix(0, NROW(time), NCOL(ovalue))
    for(i in 1:NCOL(ovalue)) {
        result <- approx(otime, ovalue[,i,drop = TRUE], time, ...)
        value[,i] <- result$y
    }
    return(irts(time, value[,,drop = TRUE]))
}

plot.irts <-
function(x, type = "l", plot.type = c("multiple", "single"),
         xlab = "Time", ylab = NULL, main = NULL, ylim = NULL,
         oma = c(6, 0, 5, 0), ...)
{
    seqN <- function(from, to) {
        if((0 == length(from)) | (0 == length(to)))
            NULL
        else if(to-from+1 <= 0) 
            NULL
        else seq(from, to)
    }
    
    addmain <- function(main, cex.main = par("cex.main"),
                        font.main = par("font.main"), 
                        col.main = par("col.main"), ...) {
        mtext(main, 3, 3, cex = cex.main, font = font.main, col = col.main, ...)
    }
    
    if(!inherits(x, "irts"))
        stop("method is only for irts objects")
    t <- time(x)
    v <- value(x)
    nser <- NCOL(v)
    if(is.null(main)) 
        main <- deparse(substitute(x))
    if(nser == 1) {
        if(is.null(ylab))
            ylab <- "Series"
        if(is.null(ylim))
            ylim <- range(v[is.finite(v)])
        plot(t, v, type = type, xlab = xlab, ylab = ylab,
             main = main, ylim = ylim, ...)
    } else if(nser <= 10) {
        plot.type <- match.arg(plot.type)
        if(is.null(ylab)) {
            ylab <- colnames(v)
            if(is.null(ylab)) 
                ylab <- paste("Series", 1:nser)
        }
        if(plot.type == "single") {
            if(is.null(ylim))
                ylim <- range(v[is.finite(v)])
            plot.default(t, v[,1], type = type, xlab = xlab, ylab = ylab,
                         main = main, ylim = ylim, xaxt = "n", ...) 
            for(i in seqN(2, nser)) {
                points(t, v[,i], type = type, xaxt = "n") 
            }
            axis.POSIXct(1, t)
        } else if(plot.type == "multiple") {
            oldpar <- par("mar", "oma", "mfcol")
            on.exit(par(oldpar))
            par(mar = c(0, 5.1, 0, 2.1), oma = oma)
            nc <- if(nser > 4) 2 else 1
            nr <- ceiling(nser/nc)
            par(mfcol = c(nr, nc))
            for(i in seqN(1, nser)) {
                plot.default(t, v[,i], type = type, xlab = xlab, ylab = "", xaxt = "n", ...)
                mtext(ylab[i], 2, 3)
                if((i%%nr == 0) || (i == nser))
                    axis.POSIXct(1, t)
            }
            if(!is.null(main)) {
                par(mfcol = c(1, 1))
                addmain(main, ...)
            }
        }
    } else {
        stop("Can't plot more than 10 series")
    }
    invisible(x)
}

lines.irts <-
function(x, type = "l", ...)
{
    if(!inherits(x, "irts"))
        stop("method is only for irts objects")
    t <- time(x)
    v <- value(x)
    nser <- NCOL(v)
    if(nser == 1) {
        lines(t, v, type = type, ...)
    } else {
        stop("Can't plot multivariate irregular time-series object")
    }
    invisible(x)
}

points.irts <-
function(x, type = "p", ...)
{
    if(!inherits(x, "irts"))
        stop("method is only for irts objects")
    t <- time(x)
    v <- value(x)
    nser <- NCOL(v)
    if(nser == 1) {
        points(t, v, type = type, ...)
    } else {
        stop("Can't plot multivariate irregular time-series object")
    }
    invisible(x)
}

## Copyright (C) 1997-2002  Adrian Trapletti
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2, or (at your option)
## any later version.
##
## This program 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
## General Public License for more details.
##
## A copy of the GNU General Public License is available via WWW at
## http://www.gnu.org/copyleft/gpl.html.  You can also obtain it by
## writing to the Free Software Foundation, Inc., 59 Temple Place,
## Suite 330, Boston, MA  02111-1307  USA. 

##
## Mostly time series tests
##

runs.test <-
function (x, alternative = c("two.sided", "less", "greater"))
{
    if(!is.factor(x))
        stop("x is not a factor")
    if(any(is.na(x)))
        stop("NAs in x")
    if(length(levels(x)) != 2)
        stop("x does not contain dichotomous data")
    alternative <- match.arg(alternative)
    DNAME <- deparse(substitute(x))
    n <- length(x)
    R <- 1 + sum(as.numeric(x[-1] != x[-n]))
    n1 <- sum(levels(x)[1] == x)
    n2 <- sum(levels(x)[2] == x)
    m <- 1 + 2*n1*n2 / (n1+n2)
    s <- sqrt(2*n1*n2 * (2*n1*n2 - n1 - n2) / ((n1+n2)^2 * (n1+n2-1)))
    STATISTIC <- (R - m) / s
    METHOD <- "Runs Test"
    if(alternative == "two.sided")
        PVAL <- 2 * pnorm(-abs(STATISTIC))
    else if(alternative == "less")
        PVAL <- pnorm(STATISTIC)
    else if(alternative == "greater")
        PVAL <- pnorm(STATISTIC, lower.tail = FALSE)
    else stop("irregular alternative")
    names(STATISTIC) <- "Standard Normal"
    structure(list(statistic = STATISTIC,
                   alternative = alternative,
                   p.value = PVAL,
                   method = METHOD,
                   data.name = DNAME),
              class = "htest")
}

bds.test <-
function(x, m = 3, eps = seq(0.5*sd(x),2*sd(x),length=4), trace = FALSE)
{
    if(NCOL(x) > 1)
        stop("x is not a vector or univariate time series")
    if(any(is.na(x)))
        stop("NAs in x")
    if(m < 2)
        stop("m is less than 2")
    if(length(eps) == 0)
        stop("invalid eps")
    if(any(eps <= 0))
        stop("invalid eps")
    DNAME <- deparse(substitute(x))
    n <- length(x)
    k <- length(eps)
    cc <- double(m+1)
    cstan <- double(m+1)
    STATISTIC <- matrix(0,m-1,k)
    for(i in (1:k)) {
        res <- .C("bdstest_main",
                  as.integer(n),
                  as.integer(m),
                  as.vector(x, mode="double"),
                  as.vector(cc),
                  cstan = as.vector(cstan),
                  as.double(eps[i]),
                  as.integer(trace),
                  PACKAGE="tseries")
        STATISTIC[,i] <- res$cstan[2:m+1]
    }
    colnames(STATISTIC) <- eps
    rownames(STATISTIC) <- 2:m
    PVAL <- 2 * pnorm(-abs(STATISTIC))
    colnames(PVAL) <- eps
    rownames(PVAL) <- 2:m
    METHOD <- "BDS Test"
    PARAMETER <- list(m = 2:m, eps = eps)
    structure(list(statistic = STATISTIC,
                   p.value = PVAL,
                   method = METHOD,
                   data.name = DNAME,
                   parameter = PARAMETER), 
              class = "bdstest")
}

print.bdstest <-
function(x, digits = 4, ...)
{
    if(!inherits(x, "bdstest"))
        stop("method is only for bdstest objects")
    cat("\n\t", x$method, "\n\n")
    cat("data: ", x$data.name, "\n\n")
    if(!is.null(x$parameter)) {
        cat("Embedding dimension = ",
            format(round(x$parameter$m, digits)), sep = " ", "\n\n")
        cat("Epsilon for close points = ",
            format(round(x$parameter$eps, digits)), sep = " ", "\n\n")
    }
    if(!is.null(x$statistic)) {
        colnames(x$statistic) <-
            round(as.numeric(colnames(x$statistic)), digits)
        colnames(x$statistic) <- paste("[",colnames(x$statistic),"]")
        rownames(x$statistic) <-
            round(as.numeric(rownames(x$statistic)), digits)
        rownames(x$statistic) <- paste("[",rownames(x$statistic),"]")
        cat("Standard Normal = \n")
        print(round(x$statistic, digits))
        cat("\n")
    }
    if(!is.null(x$p.value)) {
        colnames(x$p.value) <-
            round(as.numeric(colnames(x$p.value)), digits)
        colnames(x$p.value) <- paste("[",colnames(x$p.value),"]")
        rownames(x$p.value) <-
            round(as.numeric(rownames(x$p.value)), digits)
        rownames(x$p.value) <- paste("[",rownames(x$p.value),"]")
        cat("p-value = \n")
        print(round(x$p.value, digits))
        cat("\n")
    }
    cat("\n")
    invisible(x)
}

adf.test <-
function(x, alternative = c("stationary", "explosive"),
         k = trunc((length(x)-1)^(1/3)))
{
    if(NCOL(x) > 1)
        stop("x is not a vector or univariate time series")
    if(any(is.na(x)))
        stop("NAs in x")
    if(k < 0)
        stop("k negative")
    alternative <- match.arg(alternative)
    DNAME <- deparse(substitute(x))
    k <- k+1
    y <- diff(x)
    n <- length(y)
    z <- embed(y, k)
    yt <- z[,1]
    xt1 <- x[k:n]
    tt <- k:n
    if(k > 1) {
        yt1 <- z[,2:k] 
        res <- lm(yt ~ xt1 + 1 + tt + yt1)
    }
    else
        res <- lm(yt ~ xt1 + 1 + tt)
    res.sum <- summary(res)
    STAT <- res.sum$coefficients[2,1] / res.sum$coefficients[2,2]
    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
    interpol <- approx(tableipl, tablep, STAT, rule=2)$y
    if(is.na(approx(tableipl, tablep, STAT, rule=1)$y))
        if(interpol == min(tablep))
            warning("p-value smaller than printed p-value")
        else
            warning("p-value greater than printed p-value")
    if(alternative == "stationary")
        PVAL <- interpol
    else if(alternative == "explosive")
        PVAL <- 1 - interpol
    else stop("irregular alternative") 
    PARAMETER <- k-1
    METHOD <- "Augmented Dickey-Fuller Test"
    names(STAT) <- "Dickey-Fuller"
    names(PARAMETER) <- "Lag order"
    structure(list(statistic = STAT,
                   parameter = PARAMETER,
                   alternative = alternative,
                   p.value = PVAL,
                   method = METHOD,
                   data.name = DNAME), 
            class = "htest")
}

white.test <- function(x, ...) UseMethod("white.test")

white.test.default <-
function(x, y, qstar = 2, q = 10, range = 4,
         type = c("Chisq","F"), scale = TRUE, ...)
{
    DNAME <- paste(deparse(substitute(x)),
                   "and",
                   deparse(substitute(y)))
    x <- as.matrix(x)
    y <- as.matrix(y)
    if(any(is.na(x))) stop("NAs in x")
    if(any(is.na(y))) stop("NAs in y")
    nin <- dim(x)[2]
    t <- dim(x)[1]
    if(dim(x)[1] != dim(y)[1]) 
        stop("number of rows of x and y must match")
    if(dim(x)[1] <= 0) 
        stop("no observations in x and y")
    if(dim(y)[2] > 1)
        stop("handles only univariate outputs")
    if(!missing(type) && !is.na(pmatch(type, "chisq"))) {
        warning(paste("value `chisq' for `type' is deprecated,",
                      "use `Chisq' instead"))
        type <- "Chisq"
    }
    else
        type <- match.arg(type)
    if(scale) {
        x <- scale(x)
        y <- scale(y)
    }
    xnam <- paste("x[,", 1:nin, "]", sep="")
    fmla <- as.formula(paste("y~",paste(xnam,collapse= "+")))
    rr <- lm(fmla)
    u <- residuals(rr)
    ssr0 <- sum(u^2)
    max <- range/2
    gamma <- matrix(runif((nin+1)*q,-max,max),nin+1,q)
    phantom <- (1+exp(-(cbind(rep(1,t),x)%*%gamma)))^(-1)
    phantomstar <- as.matrix(prcomp(phantom,scale=TRUE)$x[,2:(qstar+1)])
    xnam2 <- paste("phantomstar[,", 1:qstar, "]", sep="")
    xnam2 <- paste(xnam2,collapse="+")
    fmla <- as.formula(paste("u~",paste(paste(xnam,collapse= "+"),
                                          xnam2,sep="+")))
    rr <- lm(fmla)
    v <- residuals(rr)
    ssr <- sum(v^2)
    if(type == "Chisq") {
        STAT <- t*log(ssr0/ssr)
        PVAL <- 1-pchisq(STAT,qstar)
        PARAMETER <- qstar
        names(STAT) <- "X-squared"
        names(PARAMETER) <- "df"
    }
    else if(type == "F") {
        STAT <- ((ssr0-ssr)/qstar)/(ssr/(t-qstar-nin))
        PVAL <- 1-pf(STAT,qstar,t-qstar-nin)
        PARAMETER <- c(qstar,t-qstar-nin)
        names(STAT) <- "F"
        names(PARAMETER) <- c("df1","df2")
    }
    else
        stop("invalid type")
    ARG <- c(qstar,q,range,scale)
    names(ARG) <- c("qstar","q","range","scale")
    METHOD <- "White Neural Network Test"
    structure(list(statistic = STAT,
                   parameter = PARAMETER,
                   p.value = PVAL,  
                   method = METHOD,
                   data.name = DNAME,
                   arguments = ARG),
              class = "htest")
}

white.test.ts <-
function(x, lag = 1, qstar = 2, q = 10, range = 4,
         type = c("Chisq","F"), scale = TRUE, ...)
{
    if(!is.ts(x))
        stop("method is only for time series")
    if(NCOL(x) > 1)
        stop("x is not a vector or univariate time series")
    if(any(is.na(x)))
        stop("NAs in x")
    if(lag < 1) 
        stop("minimum lag is 1")
    if(!missing(type) && !is.na(pmatch(type, "chisq"))) {
        warning(paste("value `chisq' for `type' is deprecated,",
                      "use `Chisq' instead"))
        type <- "Chisq"
    }
    else
        type <- match.arg(type)
    DNAME <- deparse(substitute(x))
    t <- length(x)
    if(scale) x <- scale(x)
    y <- embed(x, lag+1)
    xnam <- paste("y[,", 2:(lag+1), "]", sep="")
    fmla <- as.formula(paste("y[,1]~",paste(xnam,collapse= "+")))
    rr <- lm(fmla)
    u <- residuals(rr)
    ssr0 <- sum(u^2)
    max <- range/2
    gamma <- matrix(runif((lag+1)*q,-max,max),lag+1,q)
    phantom <- (1+exp(-(cbind(rep(1,t-lag),y[,2:(lag+1)])%*%gamma)))^(-1)
    phantomstar <- as.matrix(prcomp(phantom,scale=TRUE)$x[,2:(qstar+1)])
    xnam2 <- paste("phantomstar[,", 1:qstar, "]", sep="")
    xnam2 <- paste(xnam2, collapse="+")
    fmla <- as.formula(paste("u~",paste(paste(xnam,collapse= "+"),
                                          xnam2,sep="+")))
    rr <- lm(fmla)
    v <- residuals(rr)
    ssr <- sum(v^2)
    if(type == "Chisq") {
        STAT <- t*log(ssr0/ssr)
        PVAL <- 1-pchisq(STAT,qstar)
        PARAMETER <- qstar
        names(STAT) <- "X-squared"
        names(PARAMETER) <- "df"
    } else if(type == "F") {
        STAT <- ((ssr0-ssr)/qstar)/(ssr/(t-lag-qstar))
        PVAL <- 1-pf(STAT,qstar,t-lag-qstar)
        PARAMETER <- c(qstar,t-lag-qstar)
        names(STAT) <- "F"
        names(PARAMETER) <- c("df1","df2")
    }
    else
        stop("invalid type")
    ARG <- c(lag,qstar,q,range,scale)
    names(ARG) <- c("lag","qstar","q","range","scale")
    METHOD <- "White Neural Network Test"
    structure(list(statistic = STAT,
                   parameter = PARAMETER,
                   p.value = PVAL, 
                   method = METHOD,
                   data.name = DNAME,
                   arguments = ARG),
              class = "htest")
}

terasvirta.test <- function(x, ...) UseMethod("terasvirta.test")

terasvirta.test.default <-
function(x, y, type = c("Chisq", "F"), scale = TRUE, ...)
{
    DNAME <- paste(deparse(substitute(x)),
                   "and",
                   deparse(substitute(y)))
    x <- as.matrix(x)
    y <- as.matrix(y)
    if(any(is.na(x))) stop("NAs in x")
    if(any(is.na(y))) stop("NAs in y")
    nin <- dim(x)[2]
    if(nin < 1)
        stop("invalid x")
    t <- dim(x)[1]
    if(dim(x)[1] != dim(y)[1]) 
        stop("number of rows of x and y must match")
    if(dim(x)[1] <= 0) 
        stop("no observations in x and y")
    if(dim(y)[2] > 1)
        stop("handles only univariate outputs")
    if(!missing(type) && !is.na(pmatch(type, "chisq"))) {
        warning(paste("value `chisq' for `type' is deprecated,",
                      "use `Chisq' instead"))
        type <- "Chisq"
    }
    else
        type <- match.arg(type)
    if(scale) {
        x <- scale(x)
        y <- scale(y)
    }
    xnam <- paste("x[,", 1:nin, "]", sep="")
    fmla <- as.formula(paste("y~",paste(xnam,collapse= "+")))
    rr <- lm(fmla)
    u <- residuals(rr)
    ssr0 <- sum(u^2)
    xnam2 <- NULL
    m <- 0
    for(i in (1:nin)) {
        for(j in (i:nin)) {
            xnam2 <- c(xnam2,paste("I(x[,",i,"]*x[,",j,"])",sep=""))
            m <- m+1
        }
    }
    xnam2 <- paste(xnam2,collapse="+")
    xnam3 <- NULL
    for(i in (1:nin)) {
        for(j in (i:nin)) {
            for(k in (j:nin)) {
                xnam3 <- c(xnam3, paste("I(x[,", i, "]*x[,", j, "]*x[,",
                                        k ,"])", sep=""))
                m <- m+1
            }
        }
    }
    xnam3 <- paste(xnam3,collapse="+")
    fmla <- as.formula(paste("u~",paste(paste(xnam,collapse= "+"),
                                          xnam2,xnam3,sep="+")))
    rr <- lm(fmla)
    v <- residuals(rr)
    ssr <- sum(v^2)
    if(type == "Chisq") {
        STAT <- t*log(ssr0/ssr)
        PVAL <- 1-pchisq(STAT,m)
        PARAMETER <- m
        names(STAT) <- "X-squared"
        names(PARAMETER) <- "df"
    } else if(type == "F") {
        STAT <- ((ssr0-ssr)/m)/(ssr/(t-nin-m))
        PVAL <- 1-pf(STAT,m,t-nin-m)
        PARAMETER <- c(m,t-nin-m)
        names(STAT) <- "F"
        names(PARAMETER) <- c("df1","df2")
    }
    else
        stop("invalid type")
    METHOD <- "Teraesvirta Neural Network Test"
    ARG <- scale
    names(ARG) <- "scale"
    structure(list(statistic = STAT,
                   parameter = PARAMETER,
                   p.value = PVAL, 
                   method = METHOD,
                   data.name = DNAME,
                   arguments = ARG),
              class = "htest")
}

terasvirta.test.ts <-
function(x, lag = 1, type = c("Chisq", "F"), scale = TRUE, ...)
{
    if(!is.ts(x))
        stop("method is only for time series")
    if(NCOL(x) > 1)
        stop("x is not a vector or univariate time series")
    if(any(is.na(x)))
        stop("NAs in x")
    if(lag < 1) 
        stop("minimum lag is 1")
    if(!missing(type) && !is.na(pmatch(type, "chisq"))) {
        warning(paste("value `chisq' for `type' is deprecated,",
                      "use `Chisq' instead"))
        type <- "Chisq"
    }
    else
        type <- match.arg(type)
    DNAME <- deparse(substitute(x))
    t <- length(x)
    if(scale) x <- scale(x)
    y <- embed(x, lag+1)
    xnam <- paste("y[,", 2:(lag+1), "]", sep="")
    fmla <- as.formula(paste("y[,1]~",paste(xnam,collapse= "+")))
    rr <- lm(fmla)
    u <- residuals(rr)
    ssr0 <- sum(u^2)
    xnam2 <- NULL
    m <- 0
    for(i in (1:lag)) {
        for(j in (i:lag)) {
            xnam2 <- c(xnam2,paste("I(y[,",i+1,"]*y[,",j+1,"])",sep=""))
            m <- m+1
        }
    }
    xnam2 <- paste(xnam2,collapse="+")
    xnam3 <- NULL
    for(i in (1:lag)) {
        for(j in (i:lag)) {
            for(k in (j:lag)) {
                xnam3 <- c(xnam3, paste("I(y[,", i+1, "]*y[,", j+1,
                                        "]*y[,", k+1, "])", sep=""))
                m <- m+1
            }
        }
    }
    xnam3 <- paste(xnam3,collapse="+")
    fmla <- as.formula(paste("u~",paste(paste(xnam,collapse= "+"),
                                          xnam2,xnam3,sep="+")))
    rr <- lm(fmla)
    v <- residuals(rr)
    ssr <- sum(v^2)
    if(type == "Chisq") {
        STAT <- t*log(ssr0/ssr)
        PVAL <- 1-pchisq(STAT,m)
        PARAMETER <- m
        names(STAT) <- "X-squared"
        names(PARAMETER) <- "df"
    }
    else if(type == "F") {
        STAT <- ((ssr0-ssr)/m)/(ssr/(t-lag-m))
        PVAL <- 1-pf(STAT,m,t-lag-m)
        PARAMETER <- c(m,t-lag-m)
        names(STAT) <- "F"
        names(PARAMETER) <- c("df1","df2")
    }
    else
        stop("invalid type")
    METHOD <- "Teraesvirta Neural Network Test"
    ARG <- c(lag,scale)
    names(ARG) <- c("lag","scale")
    structure(list(statistic = STAT,
                   parameter = PARAMETER,
                   p.value = PVAL, 
                   method = METHOD,
                   data.name = DNAME,
                   arguments = ARG),
              class = "htest")
}

jarque.bera.test <-
function(x)
{
    if(NCOL(x) > 1)
        stop("x is not a vector or univariate time series")
    if(any(is.na(x)))
        stop("NAs in x")
    DNAME <- deparse(substitute(x))
    n <- length(x)
    m1 <- sum(x)/n
    m2 <- sum((x-m1)^2)/n
    m3 <- sum((x-m1)^3)/n
    m4 <- sum((x-m1)^4)/n
    b1 <- (m3/m2^(3/2))^2
    b2 <- (m4/m2^2)
    STATISTIC <- n*b1/6+n*(b2-3)^2/24
    names(STATISTIC) <- "X-squared"
    PARAMETER <- 2
    names(PARAMETER) <- "df"
    PVAL <- 1 - pchisq(STATISTIC,df = 2)
    METHOD <- "Jarque Bera Test"
    structure(list(statistic = STATISTIC,
                   parameter = PARAMETER,
                   p.value = PVAL,
                   method = METHOD,
                   data.name = DNAME),
              class = "htest")
}

pp.test <-
function(x, alternative = c("stationary", "explosive"),
         type = c("Z(alpha)", "Z(t_alpha)"), lshort = TRUE)
{
    if(NCOL(x) > 1)
        stop("x is not a vector or univariate time series")
    type <- match.arg(type)
    alternative <- match.arg(alternative)
    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)
    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),
                 ssqrtl=as.double(ssqru),
                 PACKAGE="tseries")$ssqrtl
    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
    if(type == "Z(alpha)") {
        alpha <- res.sum$coefficients[3,1]
        STAT <- n*(alpha-1)-(n^6)/(24*Dx)*(ssqrtl-ssqru)
        table <- cbind(c(22.5, 25.7, 27.4, 28.4, 28.9, 29.5),
                       c(19.9, 22.4, 23.6, 24.4, 24.8, 25.1),
                       c(17.9, 19.8, 20.7, 21.3, 21.5, 21.8),
                       c(15.6, 16.8, 17.5, 18.0, 18.1, 18.3),
                       c(3.66, 3.71, 3.74, 3.75, 3.76, 3.77),
                       c(2.51, 2.60, 2.62, 2.64, 2.65, 2.66),
                       c(1.53, 1.66, 1.73, 1.78, 1.78, 1.79),
                       c(0.43, 0.65, 0.75, 0.82, 0.84, 0.87))
    }
    else if(type == "Z(t_alpha)") {
        tstat <-
            (res.sum$coefficients[3,1]-1)/res.sum$coefficients[3,2]
        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))
    }
    else
        stop("irregular type")
    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
    interpol <- approx(tableipl, tablep, STAT, rule=2)$y
    if(is.na(approx(tableipl, tablep, STAT, rule=1)$y))
        if(interpol == min(tablep))
            warning("p-value smaller than printed p-value")
        else
            warning("p-value greater than printed p-value") 
    if(alternative == "stationary")
        PVAL <- interpol
    else if(alternative == "explosive")
        PVAL <- 1 - interpol
    else stop("irregular alternative")
    PARAMETER <- l
    METHOD <- "Phillips-Perron Unit Root Test"
    names(STAT) <- paste("Dickey-Fuller", type)
    names(PARAMETER) <- "Truncation lag parameter"
    structure(list(statistic = STAT,
                   parameter = PARAMETER,
                   alternative = alternative,
                   p.value = PVAL,
                   method = METHOD,
                   data.name = DNAME),
              class = "htest")
}

po.test <-
function(x, demean = TRUE, lshort = TRUE)
{
    if(NCOL(x) <= 1)
        stop("x is not a matrix or multivariate time series")
    DNAME <- deparse(substitute(x))
    x <- as.matrix(x)
    dimx <- ncol(x)
    if(dimx > 6) stop("no critical values for this dimension")
    if(demean)
        res <- lm(x[,1]~x[,-1])
    else
        res <- lm(x[,1]~x[,-1]-1)
    z <- embed(residuals(res), 2)
    ut <- z[,1]
    ut1 <- z[,2]
    n <- length(ut)
    res <- lm(ut ~ ut1 - 1)
    if(res$rank < 1)
        stop("Singularities in regression")
    res.sum <- summary(res)
    k <- residuals(res)
    ssqrk <- sum(k^2)/n
    if(lshort)
        l <- trunc(n/100)
    else
        l <- trunc(n/30)
    ssqrtl <- .C("R_pp_sum",
                 as.vector(k, mode="double"),
                 as.integer(n),
                 as.integer(l),
                 ssqrtl=as.double(ssqrk),
                 PACKAGE="tseries")$ssqrtl
    alpha <- res.sum$coefficients[1,1]
    STAT <- n*(alpha-1)-0.5*n^2*(ssqrtl-ssqrk)/(sum(ut1^2))
    if(demean) {
        table <- cbind(c(28.32, 34.17, 41.13, 47.51, 52.17),
                       c(23.81, 29.74, 35.71, 41.64, 46.53),
                       c(20.49, 26.09, 32.06, 37.15, 41.94),
                       c(18.48, 23.87, 29.51, 34.71, 39.11),
                       c(17.04, 22.19, 27.58, 32.74, 37.01),
                       c(15.93, 21.04, 26.23, 31.15, 35.48),
                       c(14.91, 19.95, 25.05, 29.88, 34.20))
    }
    else {
        table <- cbind(c(22.83, 29.27, 36.16, 42.87, 48.52),
                       c(18.89, 25.21, 31.54, 37.48, 42.55),
                       c(15.64, 21.48, 27.85, 33.48, 38.09),
                       c(13.81, 19.61, 25.52, 30.93, 35.51),
                       c(12.54, 18.18, 23.92, 28.85, 33.80),
                       c(11.57, 17.01, 22.62, 27.40, 32.27),
                       c(10.74, 16.02, 21.53, 26.17, 30.90))
    }
    table <- -table
    tablep <- c(0.01, 0.025, 0.05, 0.075, 0.10, 0.125, 0.15)
    PVAL <- approx(table[dimx-1,], tablep, STAT, rule=2)$y   
    if(is.na(approx(table[dimx-1, ], tablep, STAT, rule=1)$y))
        if(PVAL == min(tablep))
            warning("p-value smaller than printed p-value")
        else
            warning("p-value greater than printed p-value") 
    PARAMETER <- l
    METHOD <- "Phillips-Ouliaris Cointegration Test"
    if(demean)
        names(STAT) <- "Phillips-Ouliaris demeaned"
    else
        names(STAT) <- "Phillips-Ouliaris standard"
    names(PARAMETER) <- "Truncation lag parameter"
    structure(list(statistic = STAT,
                   parameter = PARAMETER,
                   p.value = PVAL,
                   method = METHOD,
                   data.name = DNAME),
              class = "htest")
}

kpss.test <-
function(x, null = c("Level", "Trend"), lshort = TRUE)
{
    if(NCOL(x) > 1)
        stop("x is not a vector or univariate time series")
    DNAME <- deparse(substitute(x))
    null <- match.arg(null)
    n <- length(x)
    if(null == "Trend") {
        t <- 1:n
        e <- residuals(lm(x ~ t))
        table <- c(0.216, 0.176, 0.146, 0.119)
    }
    else if(null == "Level") {
        e <- residuals(lm(x ~ 1))
        table <- c(0.739, 0.574, 0.463, 0.347)
    }
    tablep <- c(0.01, 0.025, 0.05, 0.10)
    s <- cumsum(e)
    eta <- sum(s^2)/(n^2)
    s2 <- sum(e^2)/n
    if(lshort)
        l <- trunc(3*sqrt(n)/13)
    else
        l <- trunc(10*sqrt(n)/14)
    s2 <- .C("R_pp_sum",
             as.vector(e, mode="double"),
             as.integer(n),
             as.integer(l),
             s2=as.double(s2),
             PACKAGE="tseries")$s2
    STAT <- eta/s2
    PVAL <- approx(table, tablep, STAT, rule=2)$y
    if(is.na(approx(table, tablep, STAT, rule=1)$y))
        if(PVAL == min(tablep))
            warning("p-value smaller than printed p-value")
        else
            warning("p-value greater than printed p-value") 
    PARAMETER <- l
    METHOD <- paste("KPSS Test for", null, "Stationarity")
    names(STAT) <- paste("KPSS", null)
    names(PARAMETER) <- "Truncation lag parameter"
    structure(list(statistic = STAT,
                   parameter = PARAMETER,
                   p.value = PVAL,
                   method = METHOD,
                   data.name = DNAME),
              class = "htest")
}
## Copyright (C) 1997-2001 Adrian Trapletti
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2, or (at your option)
## any later version.
##
## This program 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
## General Public License for more details.
##
## A copy of the GNU General Public License is available via WWW at
## http://www.gnu.org/copyleft/gpl.html.  You can also obtain it by
## writing to the Free Software Foundation, Inc., 59 Temple Place,
## Suite 330, Boston, MA  02111-1307  USA.

##
## Various time series related routines
##

read.ts <-
function(file, header = FALSE, sep = "", skip = 0, ...)
{
    x <- read.matrix(file, header = header, sep = sep, skip = skip)
    x <- ts(x, ...)
    return(x)
}

fftsurr <-
function(x)
{
    ## This is algorithm 1, p. 183 from "Theiler et al. (1992): Using
    ## Surrogate Data to Detect Nonlinearity in Time Series, in
    ## Nonlinear Modelling and Forecasting, Editors Casdagli & Eubank,
    ## Santa Fe Institute, Addison Wesley". Note that Step 7. and 8. are
    ## only for t = 2,...,N.
    z <- fft(x)
    zz <- z*exp(1i*runif(z, max=2*pi))
    re <- Re(zz[2:length(zz)]+zz[length(zz):2])/2
    im <- Im(zz[2:length(zz)]-zz[length(zz):2])/2
    zzz1 <- Re(zz[1]+zz[1])/2+1i*Im(zz[1]-zz[1])/2 
    zzz <- c(zzz1,re+1i*im)
    return(Re(fft(zzz, inverse=TRUE)))
}

ampsurr <-
function(x)
{
    ## This is algorithm 2, pp. 183, 184 from "Theiler et al. (1992):
    ## Using Surrogate Data to Detect Nonlinearity in Time Series, in
    ## Nonlinear Modelling and Forecasting, Editors Casdagli & Eubank,
    ## Santa Fe Institute, Addison Wesley".
    sx <- sort(x)
    rx <- rank(x)
    g <- rnorm(x)
    sg <- sort(g)
    y <- sg[rx]
    yy <- fftsurr(y)
    ryy <- rank(yy)
    return(sx[ryy])
}

surrogate <-
function(x, ns = 1, fft = FALSE, amplitude = FALSE, statistic = NULL, ...)
{
    call <- match.call()
    if(NCOL(x) > 1)
        stop("x is not a vector or univariate time series")
    if(any(is.na(x)))
        stop("NAs in x")
    if(ns < 1)
        stop("ns is not positive")
    n <- length(x)
    if(is.null(statistic)) {
        ists <- is.ts(x)
        if(ists) xtsp <- tsp(x)
        surrogate <- matrix(x, nrow=n, ncol=ns)
        if(fft) {
            if(amplitude)
                surrogate <- apply(surrogate, 2, ampsurr)
            else
                surrogate <- apply(surrogate, 2, fftsurr)
        }
        else
            surrogate <- apply(surrogate, 2, sample, replace = FALSE)
        if(ists) {
            attr(surrogate, "tsp") <- xtsp
            attr(surrogate, "class") <- "ts"
        }
        return(drop(surrogate))
    }
    else {
        orig.statistic <- statistic(x, ...)
        l.stat <- length(orig.statistic)
        names(orig.statistic) <- paste("t", 1:l.stat, sep="")
        stat <- matrix(0, ns, l.stat)
        if(fft) {
            if(amplitude)
                for(i in 1:ns)
                    stat[i,] <- statistic(ampsurr(x), ...)
            else
                for(i in 1:ns)
                    stat[i,] <- statistic(fftsurr(x), ...)
        }
        else
            for(i in 1:ns)
                stat[i,] <- statistic(sample(x, replace=FALSE), ...)
        colnames(stat) <- names(orig.statistic)
        bias <- apply(stat, 2, mean) - orig.statistic
        se <- apply(stat, 2, sd)
        res <- list(statistic = drop(stat),
                    orig.statistic = drop(orig.statistic),
                    bias = drop(bias),
                    se = drop(se),
                    call = call)
        attr(res, "class") <- "resample.statistic"
        return(res)
    }
}

quadmap <-
function(xi = 0.2, a = 4.0, n = 1000)
{
    if(n < 1) stop("n is not positive")
    if((xi < 0) | (xi > 1)) stop("xi is not in [0,1]")
    if((a < 0) | (xi > 4)) stop("a is not in [0,4]")
    x <- double(n)
    res <- .C("R_quad_map",
              x = as.vector(x),
              as.double(xi),
              as.double(a),
              as.integer(n),
              PACKAGE="tseries")
    return(ts(res$x))
}

read.matrix <-
function(file, header = FALSE, sep = "", skip = 0)
{
    row.lens <- count.fields(file, sep = sep, skip = skip)
    if(any(row.lens != row.lens[1])) 
        stop("number of columns is not constant")
    if(header) {
        nrows <- length(row.lens) - 1
        ncols <- row.lens[2]
        col.names <- scan(file, what = "", sep = sep, nlines = 1,
                          quiet = TRUE, skip = skip)
        x <- scan(file, sep = sep, skip = skip + 1, quiet = TRUE)
    }
    else {
        nrows <- length(row.lens)
        ncols <- row.lens[1]
        x <- scan(file, sep = sep, skip = skip, quiet = TRUE)
        col.names <- NULL
    }
    x <- as.double(x)
    if(ncols > 1) {
        dim(x) <- c(ncols,nrows)
        x <- t(x)
        colnames(x) <- col.names
    }
    else if(ncols == 1)
        x <- as.vector(x)
    else
        stop("wrong number of columns")
    return(x)
}

na.remove <- function(object, ...) UseMethod("na.remove")

na.remove.ts <-
function(object, ...)
{
    x <- object                         # generic/method
    if(!is.ts(x)) stop("method is only for time series")
    if(any(is.na(x))) {
        y <- na.remove.default(x)
        ok <- seq(1,NROW(x))[-attr(y,"na.removed")]
        xfreq <- frequency(x)
        start <- tsp(x)[1]+(ok[1]-1)/xfreq
        end <- tsp(x)[1]+(ok[length(ok)]-1)/xfreq
        yfreq <- (NROW(y)-1)/(end-start)
        attr(y, "tsp") <- c(start,end,yfreq)
        attr(y, "class") <- attr(x, "class")
        return(y)
    }
    else return(x)
}

na.remove.default <-
function(object, ...)
{
    x <- object                         # generic/method
    if(any(is.na(x))) {
        if(is.matrix(x)) {
            nas <- apply(is.na(x),1,any)
            y <- matrix(as.vector(x)[rep(!nas,ncol(x))],ncol=ncol(x))
            dimnames(y) <- dimnames(x)
            nas <- which(nas)
        }
        else {
            nas <- which(is.na(x))
            y <- x[-nas]
        }
        attr(y, "na.removed") <- nas
        return(y)
    }
    else return(x)
}

seqplot.ts <-
function(x, y, colx = "black", coly = "red", typex = "l",
         typey = "l", pchx = 1, pchy = 1, ltyx = "solid",
         ltyy = "solid", oma = c(6, 0, 5, 0), ann = par("ann"),
         xlab = "Time", ylab = deparse(substitute(x)), main = NULL)
{
    if(!is.ts(x) || !is.ts(y))
        stop("x or y is not a time series")
    if(abs(frequency(x)-frequency(y)) > getOption("ts.eps"))
        stop("x and y don't have the same frequency")
    nser <- NCOL(x)
    nsery <- NCOL(y)
    if(nser != nsery) stop("x and y don't have consistent dimensions")
    if(nser == 1) {
        xlim <- range(time(x), time(y))
        ylim <- range(x[is.finite(x)], y[is.finite(y)])
        plot(x, xlim = xlim, ylim = ylim, col = colx, type = typex, pch
             = pchx, lty = ltyx, xlab = "", ylab = ylab)
        points(y, col = coly, type = typey, pch = pchy, lty = ltyy)
        if(ann) {
            mtext(xlab, 1, 3)
            if(!is.null(main)) title(main)
        }
    }
    else {
        if(nser > 10) stop("Can't plot more than 10 series")
        if(is.null(main)) main <- deparse(substitute(x))
        nm <- colnames(x)
        if(is.null(nm)) nm <- paste("Series", 1:nser)
        nc <- if(nser >  4) 2 else 1
        oldpar <- par("mar", "oma", "mfcol")
        on.exit(par(oldpar))
        par(mar = c(0, 5.1, 0, 2.1), oma = oma)
        nr <- ceiling(nser %/% nc)
        par(mfcol = c(nr, nc))
        for(i in 1:nser) {
            xlim <- range(time(x[,i]), time(y[,i]))
            ylim <- range(x[is.finite(x[,i]),i], y[is.finite(y[,i]),i])
            plot(x[,i], xlim = xlim, ylim = ylim, col = colx, type =
                 typex, pch = pchx, lty = ltyx, axes = FALSE, xlab = "",
                 ylab = "")
            points(y[,i], col = coly, type = typey, pch = pchy, lty =
                   ltyy)
            box()
            axis(2, xpd = NA)
            mtext(nm[i], 2, 3)
            if((i%%nr==0) || (i==nser))
                axis(1, xpd = NA)
        }
        if(ann) {
            mtext(xlab, 1, 3)
            if(!is.null(main)) {
                par(mfcol = c(1,1))
                mtext(main, 3, 3, cex=par("cex.main"),
                      font=par("font.main"), col=par("col.main"))
            }
        }
    }
    invisible()
}

boot.sample <-
function(x, b, type)
{
    return(.C("boot",
              as.vector(x, mode = "double"),
              x = as.vector(x, mode = "double"),
              as.integer(length(x)),
              as.double(b),
              as.integer(type),
              PACKAGE = "tseries")$x)
}

tsbootstrap <- function(x, nb = 1, statistic = NULL, m = 1, b = NULL,
                        type = c("stationary","block"), ...)
{
    call <- match.call()
    type <- match.arg(type)
    if(NCOL(x) > 1)
        stop("x is not a vector or univariate time series")
    if(any(is.na(x)))
        stop("NAs in x")
    if(nb < 1)
        stop("nb is not positive")
    n <- NROW(x)
    if (n <= m)
        stop("x should contain more than m observations")
    const <- 3.15
    if(type == "stationary") {
        type <- 0
        if(is.null(b)) b <- const*n^(1/3)
        b <- 1/b
        if((b <= 1/n) || (b >= 1))
            stop(paste("b should be in (1,length(x))",
                       "for the stationary bootstrap"))
    }
    else {
        type <- 1
        if(is.null(b)) b <- const*n^(1/3)
        if((b < 1) || (b > n))
            stop(paste("b should be in [1,length(x)]",
                       "for the blockwise bootstrap"))
    }
    if(is.null(statistic)) {
        if (m > 1)
            stop("Can only return bootstrap data for m = 1")
        ists <- is.ts(x)
        if(ists) xtsp <- tsp(x)
        boot <- matrix(x, nrow=n, ncol=nb)
        boot <- apply(boot, 2, boot.sample, b, type)    
        if(ists) {
            attr(boot, "tsp") <- xtsp
            attr(boot, "class") <- "ts"
        }
        return(drop(boot))
    }
    else {
        y <- embed(x, m)
        yi <- 1:NROW(y)
        orig.statistic <- statistic(drop(y), ...)
        l.stat <- length(orig.statistic)
        names(orig.statistic) <- paste("t", 1:l.stat, sep="")
        stat <- matrix(0, nb, l.stat)
        for(i in 1:nb)
            stat[i,] <- statistic(y[boot.sample(yi, b, type), , drop=TRUE], ...)
        colnames(stat) <- names(orig.statistic)
        bias <- apply(stat, 2, mean)-orig.statistic
        se <- apply(stat, 2, sd)
        res <- list(statistic = drop(stat),
                    orig.statistic = drop(orig.statistic),
                    bias = drop(bias),
                    se = drop(se),
                    call = call)
        attr(res, "class") <- "resample.statistic"
        return(res)
    }
}

print.resample.statistic <-
function(x, digits = max(3, getOption("digits") - 3), ...)
{
    cat("\nCall:", deparse(x$call), "", sep = "\n")
    nam <- c("original", "bias", "std. error")
    stat <- cbind(x$orig.statistic, x$bias, x$se)
    colnames(stat) <- nam
    cat("Resampled Statistic(s):\n")
    print(drop(stat), digits = digits, ...)
    cat("\n")
    invisible(x)
}
.First.lib <-
function(lib, pkg)
{
    library.dynam("tseries", pkg, lib)
    require("stats", quietly = TRUE)
    mylib <- dirname(system.file(package = "tseries"))
    ver <- packageDescription("tseries", lib = mylib)["Version"]
    txt <- c("\n",
             paste(sQuote("tseries"), "version:", ver),
             "\n",
             paste(sQuote("tseries"),
                   "is a package for time series analysis",
                   "and computational finance."),
             "\n",
             paste("See",
                   sQuote("library(help=\"tseries\")"),
                   "for details."),
             "\n")
    if(interactive() || getOption("verbose"))
        writeLines(strwrap(txt, indent = 4, exdent = 4))
}
