# Copyright (C) 1997-2000  Adrian Trapletti
#
# 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

#
# 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"
  }
  if (!is.null(lag$ar)) nam.ar <- paste("ar", lag$ar, sep = "")
  else nam.ar <- NULL
  if (!is.null(lag$ma)) nam.ma <- paste("ma", lag$ma, sep = "")
  else nam.ma <- NULL
  if (include.intercept) nam.int <- "intercept"
  else nam.int <- 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 (object, digits = max(3,.Options$digits-3))
{
  if (!inherits(object, "arma")) stop ("method is only for arma objects")
  cat ("\nCall:\n", deparse(object$call), "\n\n", sep = "")
  cat ("Coefficient(s):\n")
  print.default (format(coef(object), digits = digits), print.gap = 2, quote = FALSE)
  cat ("\n")
  invisible (object)
}

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

plot.arma <- function (object, ask = interactive())
{
  if (!inherits(object, "arma")) stop ("method is only for arma objects")
  op <- par()
  par (ask = ask, mfrow=c(2,1))
  x <- eval (parse(text=object$series))
  if (any(is.na(x))) stop ("NAs in x")
  plot(x, main = object$series, ylab = "Series")
  plot(object$residuals, main = "Residuals", ylab = "Series")
  acf (x, main = paste("ACF of",object$series))
  acf (object$residuals, main = "ACF of Residuals", na.action=na.remove)
  pacf (x, main = paste("PACF of",object$series))
  pacf (object$residuals, main = "PACF of Residuals", na.action=na.remove)
  par (ask = op$ask, mfrow = op$mfrow)
  invisible (object)
}

# Copyright (C) 1997-1999  Adrian Trapletti
#
# 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

#
# Descriptive time series analysis in the time domain
#


amif <- function (x, lag.max = NULL, maxbit = 20, confidence = 0.2, ci = 0.95, nsurr = 20,
                  fft = FALSE, amplitude = FALSE, normalized = TRUE, trace = FALSE,
                  plot = TRUE, ...)
{
  if (NCOL(x) > 1) stop ("x is not a vector or univariate time series")
  if ((maxbit < 1) | (maxbit > 25)) stop ("maxbit out of range")
  if ((confidence < 0.01) | (confidence > 0.99)) stop ("confidence out of range")
  if (nsurr < 1) stop ("nsurr is not positive")
  if (ci >= 1) stop ("ci out of range")
  series <- deparse (substitute(x))
  x <- as.ts(x)
  x.freq <- frequency(x)
  x <- as.matrix(x)
  if (any(is.na(x))) stop ("NAs in x")
  sampleT <- nrow(x)
  if (is.null(lag.max))
    lag.max <- floor (10*(log10(sampleT)-log10(1)))
  lag.max <- min (lag.max, sampleT-1)
  if (lag.max < 1) 
    stop ("lag.max must be at least 1")
  lag <- matrix (1, 1, 1)
  lag[lower.tri(lag)] <- -1
  inf <- double (lag.max+1)
  cor <- array(.C ("R_amif", as.vector(x,mode="double"), as.integer(sampleT), inf=as.vector(inf),
                   as.integer(lag.max), as.integer(maxbit), as.double(confidence),
                   as.integer(normalized), as.integer(trace), PACKAGE="tseries")$inf,
               c(lag.max + 1, 1, 1))
  if (ci > 0)
  {
    surrsam <- surrogate (x, ns=nsurr, fft=fft, amplitude=amplitude)
    surrwb <- matrix (0, nrow = nsurr, ncol = lag.max+1)
    for (i in 1:nsurr)
      surrwb[i,] <- .C ("R_amif", as.vector(surrsam[,i],mode="double"), as.integer(sampleT),
                        inf=as.vector(inf), as.integer(lag.max), as.integer(maxbit),
                        as.double(confidence), as.integer(normalized),
                        as.integer(trace), PACKAGE="tseries")$inf
    wb <- apply (surrwb, 2, quantile, ci)
  }
  else
    wb <- NULL
  lag <- outer(0:lag.max, lag/x.freq)
  amif <- structure(.Data = list(acf = cor, type = "covariance", n.used = sampleT, lag = lag,
                      series = series, snames = colnames(x), clim = wb, normalized = normalized),
                    class = c("amif","acf"))
  if (plot)
  {
    plot.amif (amif, ...)
    return (invisible(amif))
  }
  else return (amif)  
}

plot.amif <- function (object, ci.col = "blue", ylab = "AMIF", ...)
{
  if (!inherits(object, "amif")) stop ("method is only for amif objects")
  plot.acf (x = object, ylab = ylab, ...)
  if (!is.null(object$clim))
    lines(object$lag[, 1, 1], object$clim, col = ci.col, lty = 2)
}
# Copyright (C) 1997-1999  Adrian Trapletti
#
# 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

#
# Financial time series analysis
#


portfolio.optim <- function (object, ...) { 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)
{
  if (!require (quadprog, quietly=TRUE))
    stop ("Package 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]
  Dmat <- cov(x)
  dvec <- rep(0,k)
  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)
    }
    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)
    {
      a3 <- NULL
      b3 <- NULL
    }
    else
    {
      a3 <- matrix(0,k,k)
      diag(a3) <- 1
      b3 <- rep(0,k)
    }  
    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", "Volume"),
                            provider = "yahoo", method = "auto")
{
  
  strsplit2 <- function (str)
  {
    unlist(lapply(strsplit(str," "),function(x)x[x!=""]))
  }

  convert.to.julian <- function (dates)
  {
    mn <- c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec")
    ld <- length(dates)
    dm <- matrix(unlist(lapply(dates,strsplit,"-")),ld,3,byrow=T)
    dd <- as.numeric(dm[,1])
    mm <- unlist(lapply(dm[,2],pmatch,mn))
    yy <- as.numeric(dm[,3])
    yy <- (yy<50)*(2000+yy)+(yy>=50)*(1900+yy)
    return (julian(mm,dd,yy))
  }

  if (!require (chron, quietly=TRUE))
    stop ("Package chron is needed. Stopping")
  quote <- match.arg(quote)
  provider <- match.arg(provider)
  mm <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")
  if (missing(start)) start <- "1 2 1991"
  if (missing(end)) end <- paste (pmatch(strsplit2(date())[2],mm),
                                  as.numeric(strsplit2(date())[3]),
                                  strsplit2(date())[5])
  start <- c(as.numeric(strsplit2(start)[1]),
             as.numeric(strsplit2(start)[2]),
             as.numeric(strsplit2(start)[3]))
  end <- c(as.numeric(strsplit2(end)[1]),
           as.numeric(strsplit2(end)[2]),
           as.numeric(strsplit2(end)[3]))
  end <- month.day.year(julian(end[1],end[2],end[3])-1)
  end <- c(end$month,end$day,end$year)
  if (provider == "yahoo")
  {
    url <- paste ("http://chart.yahoo.com/table.csv?s=", instrument, sep="")
    url <- paste (url, "&a=", start[1], "&b=", start[2], "&c=", start[3], sep="")
    url <- paste (url, "&d=", end[1], "&e=", end[2], "&f=", end[3], sep="")
    url <- paste (url, "&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))
    }
    status <- scan (destfile, "", n=1, sep="\n", quiet=TRUE)
    if (substring(status,1,2) == "No")
    {
      unlink (destfile)
      stop (paste("No data available for", instrument))
    }
    x <- read.table (destfile, header=T, sep=",")
    unlink (destfile)
    nser <- pmatch (quote, c("Open", "High", "Low", "Close", "Volume")) + 1
    if (nser > ncol(x)) stop ("This quote is not available")
    n <- nrow(x)
    dat <- gsub(" ", "0", as.character(x[n:1,1]))
    dat <- convert.to.julian (dat)
    ser <- as.vector(x[n:1,nser]) 
    seqdat <- dat[1]:dat[n]
    idx <- match(dat,seqdat)
    newser <- rep(NA,length(seqdat))
    newser[idx] <- ser
    if (dat[1] != julian(start[1],start[2],start[3]))
    {
      st <- month.day.year(dat[1])
      cat (paste("time series starts ",
                 paste(st$month,st$day,st$year,sep= " "), "\n", sep=""))
    }
    if (dat[n] != julian(end[1],end[2],end[3]))
    {
      en <- month.day.year(dat[n])
      cat (paste("time series ends   ",
                 paste(en$month,en$day,en$year,sep= " "), "\n", sep=""))
    }
    return (ts(newser, start = as.numeric(dat)[1], end = as.numeric(dat)[n]))
  }
  else stop ("provider not implemented")
}

# Copyright (C) 1997-1999  Adrian Trapletti
#
# 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

#
# 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]), 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 (object, digits = max(3,.Options$digits-3))
{
  if (!inherits(object, "garch")) stop ("method is only for garch objects")
  cat ("\nCall:\n", deparse(object$call), "\n\n", sep = "")
  cat ("Coefficient(s):\n")
  print.default (format(coef(object), digits = digits), print.gap = 2, quote = FALSE)
  cat ("\n")
  invisible (object)
}

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

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

predict.garch <- function (object, newdata, genuine = FALSE)
{
  if (!inherits(object, "garch")) stop ("method is only for garch objects")
  if (missing(newdata))
  {
    newdata <- eval (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]), 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)
  {
    if (genuine) attr(pred$h, "tsp") <- c(newdata.tsp[1],
                                          newdata.tsp[2]+1/newdata.tsp[3],
                                          newdata.tsp[3])
    else attr(pred$h, "tsp") <- newdata.tsp
    attr(pred$h, "class") <- "ts"
  }
  return (pred$h)
}

# Copyright (C) 1997-2000  Adrian Trapletti
#
# 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

#
# Mostly time series tests
#


runs.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))
  if (any (x == 0.0))
  {
    cat ("Removed", length (x[x==0.0]), "zero(es)\n")
    x <- x[x!=0.0]
  }
  d <- diff (sign(x))
  f <- factor (d)
  sp <- split (d, f)
  resS <- sapply (sp, length)
  resL <- lapply (sp, length)
  n <- length (x)
  sum2 <- sum (resS^2)
  sum3 <- sum (resS^3)
  m <- (n*(n+1)-sum2)/n
  s <- (sum2*(sum2+n*(n+1))-2*n*sum3-n^3)/(n^2*(n-1))
  R <- 1
  if (!is.null(resL$"-2"))
    R <- R+resL$"-2"
  if (!is.null(resL$"2"))
    R <- R+resL$"2"
  STATISTIC <- ((R+0.5)-m)/s
  METHOD <- "Runs Test"
  PVAL <- 2 * pnorm (-abs(STATISTIC))
  names(STATISTIC) <- "Standard Normal"
  structure(list(statistic = STATISTIC,
		 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 (object, digits = 4)
{
  if (!inherits(object, "bdstest")) stop ("method is only for bdstest objects")
  cat("\n\t", object$method, "\n\n")
  cat("data: ", object$data.name, "\n\n")
  if (!is.null(object$parameter))
  {
    cat("Embedding dimension = ", format(round(object$parameter$m, digits)), sep = " ", "\n\n")
    cat("Epsilon for close points = ", format(round(object$parameter$eps, digits)),
        sep = " ", "\n\n")
  }
  if (!is.null(object$statistic))
  {
    colnames(object$statistic) <- round (as.numeric(colnames(object$statistic)), digits)
    colnames(object$statistic) <- paste("[",colnames(object$statistic),"]")
    rownames(object$statistic) <- round (as.numeric(rownames(object$statistic)), digits)
    rownames(object$statistic) <- paste("[",rownames(object$statistic),"]")
    cat("Standard Normal = \n")
    print (round(object$statistic, digits))
    cat("\n")
  }
  if (!is.null(object$p.value))
  {
    colnames(object$p.value) <- round (as.numeric(colnames(object$p.value)), digits)
    colnames(object$p.value) <- paste("[",colnames(object$p.value),"]")
    rownames(object$p.value) <- round (as.numeric(rownames(object$p.value)), digits)
    rownames(object$p.value) <- paste("[",rownames(object$p.value),"]")
    cat("p-value = \n")
    print (round(object$p.value, digits))
    cat("\n")
  }
  cat("\n")
  invisible(object)
}

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 (object, ...) { 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 (!require (mva, quietly=TRUE))
    stop ("Package mva is needed. Stopping")
  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 (!require (mva, quietly=TRUE))
    stop ("Package mva is needed. Stopping")
  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 (object, ...) { 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")
  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")
  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-1999  Adrian Trapletti
#
# 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

#
# 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 (x)
{
  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 (x)
{
  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)) > .Options$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 = F, 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), x=as.vector(x), as.integer(length(x)),
              as.double(b), as.integer(type))$x)
}

bootstrap <- function (x, nb = 1, statistic = NULL, 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 <- length(x)
  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 ("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 ("b should be in [1,length(x)] for the block bootstrap")
  }
  if (is.null(statistic))
  {
    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
  {
    orig.statistic <- statistic (x, ...)
    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 (boot.sample(x, b, type), ...)
    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 (object, digits = max(3,.Options$digits-3), ...)
{
  cat("\nCall:", deparse(object$call), "", sep = "\n")
  nam <- c("original", "bias", "std. error")
  stat <- cbind(object$orig.statistic,object$bias,object$se)
  colnames(stat) <- nam
  cat ("Resampled Statistic(s):\n")
  print (drop(stat), digits, ...)
  cat("\n")
  invisible(object)
}


.First.lib <- function (lib, pkg)
{
  library.dynam("tseries", pkg, lib)
  if (!require (ts, quietly=TRUE))
    stop ("Package ts is needed. Stopping")
  mylib <- .path.package("tseries")
  mylib <- substr(mylib, 1, nchar(mylib)-7)
  ver <- package.description("tseries", lib=mylib)$Version
  vertxt <- paste ("\n      `tseries' version:", ver, "\n")
  introtxt <- paste ("\n      `tseries' is a package for time series analysis with emphasize\n",
                     "       on non-linear modelling.\n",
                     "       See `library (help=tseries)' for details.\n\n", sep="")
  if (interactive() || .Options$verbose)
    cat (paste (vertxt, introtxt))
}
  
