R/isat.R

Defines functions isat

Documented in isat

isat <-
function(y, mc=TRUE, ar=NULL, ewma=NULL, mxreg=NULL,
  iis=FALSE, sis=TRUE, tis=FALSE, uis=FALSE, blocks=NULL,
  ratio.threshold=0.8, max.block.size=30, t.pval=0.001, wald.pval=t.pval,
  vcov.type=c("ordinary", "white", "newey-west"), do.pet=FALSE,
  ar.LjungB=NULL, arch.LjungB=NULL, normality.JarqueB=NULL,
  user.diagnostics=NULL, info.method=c("sc", "aic", "hq"),
  include.gum=FALSE, include.empty=FALSE, tol=1e-07,
  LAPACK=FALSE, max.regs=NULL, print.searchinfo=TRUE, plot=NULL,
  alarm=FALSE)
{

  ##arguments:
  isat.call <- sys.call()
  vcov.type <- match.arg(vcov.type)
  info.method <- match.arg(info.method)
  mod <- arx(y, mc=mc, ar=ar, ewma=ewma, mxreg=mxreg,
    vcov.type=vcov.type, qstat.options=NULL, user.diagnostics=user.diagnostics,
    tol=tol, LAPACK=LAPACK, plot=FALSE)
  y <- mod$aux$y
  y.n <- mod$aux$y.n
  y.index <- mod$aux$y.index
  y.index.as.char <- as.character(y.index)
  y.name <- mod$aux$y.name
  mX <- mod$aux$mX #NULL if is.null(mX)
  mXnames <- mod$aux$mXnames #NULL if is.null(mX)
  colnames(mX) <- mXnames
  mXncol <- mod$aux$mXncol
  vcov.type <- mod$aux$vcov.type
  qstat.options <- mod$aux$qstat.options
  if(is.null(mX)){ mxkeep <- NULL }else{ mxkeep <- 1:mXncol }

  ##indicator saturation matrices:
  ISmatrices <- list()

  if(iis){ #impulse indicators
    mIIS <- matrix(0,y.n,y.n)
    diag(mIIS) <- 1
    colnames(mIIS) <- paste("iis", y.index.as.char, sep="")
    ISmatrices <- c(ISmatrices,list(IIS=mIIS))
  }

  if(sis){ #step-shift indicators
    mSIS <-matrix(0,y.n,y.n)
    loop.indx <- 1:y.n
    tmp <- function(i){ mSIS[i,1:i] <<- 1 }
    tmp <- sapply(loop.indx,tmp)
    colnames(mSIS) <- paste("sis", y.index.as.char, sep="")
    mSIS <- mSIS[,-1]
    ISmatrices <- c(ISmatrices,list(SIS=mSIS))
  }

  if(tis){ #trend indicators
    mTIS <- matrix(0,y.n,y.n)
    v1n <- seq(1,y.n)
    loop.indx <- 1:y.n
    tmp <- function(i){
      mTIS[c(i:y.n),i] <<- v1n[1:c(y.n-i+1)]
    }
    tmp <- sapply(loop.indx,tmp)
    colnames(mTIS) <- paste("tis", y.index.as.char, sep="")
    mTIS <- mTIS[,-1]
    ISmatrices <- c(ISmatrices,list(TIS=mTIS))
  }

  ##user-defined indicators/variables:
  #if uis is a matrix:
  if(!is.list(uis) && !identical(as.numeric(uis),0)){

    ##handle colnames:
    uis <- as.zoo(cbind(uis))
    uis.names <- colnames(uis)
    if(is.null(uis.names)){
      uis.names <- paste("uisxreg", 1:NCOL(uis), sep="")
    }
    if(any(uis.names == "")){
      missing.colnames <- which(uis.names == "")
      for(i in 1:length(missing.colnames)){
       uis.names[i] <- paste("uisxreg", missing.colnames[i], sep="")
      }
    }
    #for the future?: uis.names <- make.names(uis.names)

    ##select sample:
    uis <- na.trim(uis, sides="both", is.na="any")
    uis.index.as.char <- as.character(index(uis))
    t1 <- which(uis.index.as.char==y.index.as.char[1])
    t2 <- which(uis.index.as.char
      ==y.index.as.char[length(y.index.as.char)])
    uis <- coredata(uis)
    uis <- window(uis, start=t1, end=t2)
    uis <- cbind(coredata(as.zoo(uis)))
    colnames(uis) <- uis.names

    #check nrow(uis):
    if(nrow(uis) != y.n) stop("nrow(uis) is unequal to no. of observations")
    ISmatrices <- c(ISmatrices,list(UIS=uis))
  }
  #if uis is a list of matrices:
  if(is.list(uis)){
    #check nrow(uis[[i]]):
    for(i in 1:length(uis)){
      uis[[i]] <- as.matrix(coredata(as.zoo(uis[[i]])))
      if(nrow(uis[[i]]) != y.n){
        stop(paste("nrow(uis[[",i,"]]) is unequal to no. of observations",
          sep=""))
      }
    }
    uis.names <- paste("UIS", 1:length(uis), sep="")
    if(is.null(names(uis))){
      names(uis) <- uis.names
    }else{
      for(i in 1:length(uis)){
        if(names(uis)[i]==""){
          names(uis)[i] <- uis.names[i]
        }else{
          names(uis)[i] <- paste(uis.names[i], ".", names(uis)[i],
            sep="")
        } #close if..else
      } #close for..loop
    }
    ISmatrices <- c(ISmatrices,uis)
  }

  ##check blocks:
  if(is.list(blocks)){
    if(length(ISmatrices)!=length(blocks)){
      stop("No. of IS matrices is unequal to length(blocks)")
    }
    blocks.is.list <- TRUE
    ISblocks <- blocks
  }else{
    blocks.is.list <- FALSE
    ISblocks <- list()
  }

  ##loop on ISmatrices:
  ISfinalmodels <- list()
  for(i in 1:length(ISmatrices)){

    ##blocks:
    if(!blocks.is.list){

      ncol.adj <- NCOL(ISmatrices[[i]])

      if(is.null(blocks)){
        blockratio.value <- ncol.adj/(ratio.threshold*ncol.adj - NCOL(mX))
        blocksize.value <- ncol.adj/max.block.size
        no.of.blocks <- max(2,blockratio.value,blocksize.value)
        no.of.blocks <- ceiling(no.of.blocks)
        no.of.blocks <- min(ncol.adj, no.of.blocks) #ensure blocks < NCOL
#        if(ncol.adj <= 2){ no.of.blocks <- 1 }
      }else{
        no.of.blocks <- blocks
      }

      blocksize <- ceiling(ncol.adj/no.of.blocks)
      partitions.t2 <- blocksize
      for(j in 1:no.of.blocks){
        if( blocksize*j <= ncol.adj ){
          partitions.t2[j] <- blocksize*j
        }
      }
      #check if last block contains last indicator:
      if(partitions.t2[length(partitions.t2)] < ncol.adj){
        partitions.t2 <- c(partitions.t2, ncol.adj)
      }
      blocksadj <- length(partitions.t2)
      partitions.t1 <- partitions.t2 + 1
      partitions.t1 <- c(1,partitions.t1[-blocksadj])

      tmp <- list()
      for(j in 1:blocksadj){
        tmp[[j]] <- partitions.t1[j]:partitions.t2[j]
      }
      ISblocks[[i]] <- tmp

    } #end if(!blocks.is.list)

    ##gets on each block:
    ISspecific.models <- list()
    #for the future?: ISgums <- list(); ISpaths <- list(); ISterminals.results <- list()
    for( j in 1:length(ISblocks[[i]]) ){

      ##print info:
      if(print.searchinfo){
        message("\n", appendLF=FALSE)
        message(names(ISmatrices)[i],
          " block ", j, " of ", length(ISblocks[[i]]), ":",
          appendLF=TRUE)
        message("\n", appendLF=FALSE)
      }

      ##check if block contains 1 regressor:
      if( length(ISblocks[[i]][[j]])==1 ){
        tmp <- colnames(ISmatrices[[i]])[ ISblocks[[i]][[j]] ]
        mXis <- cbind(ISmatrices[[i]][, ISblocks[[i]][[j]] ])
        colnames(mXis) <- tmp
        mXis <- cbind(mX, mXis)
      }else{
        mXis <- cbind(mX,ISmatrices[[i]][, ISblocks[[i]][[j]] ])
      }

      ##if(UIS), then apply dropvar:
#      if( substr(names(ISmatrices)[i],1,3)=="UIS" ){
        mXis <- dropvar(mXis, tol=tol, LAPACK=LAPACK,
          silent=print.searchinfo)
#      }

      ##gum:
      mod <- arx(y, mxreg=mXis, vcov.type=vcov.type,
        qstat.options=qstat.options, user.diagnostics=user.diagnostics,
        tol=tol, LAPACK=LAPACK, plot=FALSE)
#future?: keep=NULL instead of mxkeep?
#      getsis <- getsm(mod, keep=NULL, t.pval=t.pval,
      getsis <- getsm(mod, keep=mxkeep, t.pval=t.pval,
        wald.pval=wald.pval, do.pet=do.pet, ar.LjungB=ar.LjungB,
        arch.LjungB=arch.LjungB, normality.JarqueB=normality.JarqueB,
        user.diagnostics=user.diagnostics, info.method=info.method,
        include.empty=include.empty, max.regs=max.regs,
        estimate.specific=FALSE, print.searchinfo=print.searchinfo,
        plot=FALSE)

      if(is.null(getsis$specific.spec)){
        ISspecific.models[[j]] <- NULL
      }else{
        ISspecific.models[[j]] <- names(getsis$specific.spec)
#For the future?:
#        ISgums[[j]] <- getsis$gum.mean
#        ISpaths[[j]] <- getsis$paths
#        ISterminals.results[[j]] <- getsis$terminals.results
      }

    } #end for(j in 1:length(ISblocks[[i]]))

    ##print info:
    if(print.searchinfo){
      message("\n", appendLF=FALSE)
      message("GETS of union of retained ",
        names(ISmatrices)[i], " variables... ",
        appendLF=TRUE)
    }

    ##if no indicators retained from the blocks:
    if(length(ISspecific.models)==0){
      isNames <- NULL
      ISfinalmodels[[i]] <- NULL
    }

    ##when indicators/variables(uis) retained from the blocks:
    if(length(ISspecific.models)>0){

      isNames <- NULL

      #which indicators/variables(uis) retained?:
      for(j in 1:length(ISspecific.models)){
        #check if mean is non-empty:
        if(!is.null(ISspecific.models[[j]])){
          isNames <- union(isNames, ISspecific.models[[j]])
        }
      } #end for(j) loop
      isNames <- setdiff(isNames, mXnames)
      #isNamesAll[[i]] <- isNames

      #redo gets with union of retained indicators:
      mXisNames <- c(mXnames,isNames)
      mXis <- cbind(mX,ISmatrices[[i]][,isNames])
      colnames(mXis) <- mXisNames
      mXis <- dropvar(mXis, tol=tol, LAPACK=LAPACK,
        silent=print.searchinfo)
      mod <- arx(y, mxreg=mXis, vcov.type=vcov.type,
        qstat.options=NULL, tol=tol, LAPACK=LAPACK,
        plot=FALSE)
      getsis <- getsm(mod, keep=mxkeep, t.pval=t.pval,
        do.pet=do.pet, wald.pval=wald.pval, ar.LjungB=ar.LjungB,
        arch.LjungB=arch.LjungB, normality.JarqueB=normality.JarqueB,
        user.diagnostics=user.diagnostics, info.method=info.method,
        include.gum=include.gum, include.empty=include.empty,
        max.regs=max.regs, print.searchinfo=print.searchinfo,
        estimate.specific=FALSE, plot=FALSE)
      ISfinalmodels[[i]] <- names(getsis$specific.spec)

    } #end if(length(ISspecific.models > 0)

  } #end for(i) loop (on ISmatrices)

  ##add names to ISblocks:
  names(ISblocks) <- names(ISmatrices)

  ##gets of union of retained impulses:
  if(print.searchinfo){
    message("\n", appendLF=FALSE)
    message("GETS of union of ALL retained variables...",
      appendLF=TRUE)
    message("\n", appendLF=FALSE)
  }

  ##no final models estimated:
  if(length(ISfinalmodels)==0){
    ISfinalmodels <- NULL
    if(is.null(mX)){ mXis <- NULL }else{
      mXis <- zoo(cbind(mX), order.by=y.index)
      colnames(mXis) <- mXnames
    }
  }

  ##final models estimated:
  if(length(ISfinalmodels)>0){

    mIS <- NULL #matrix

    #which indicators were retained?
    for(i in 1:length(ISfinalmodels)){
      isNames <- NULL
      #check if non-empty:
      if(!is.null(ISfinalmodels[[i]])){
        isNames <- setdiff(ISfinalmodels[[i]], mXnames)
      }
      if(length(isNames)>0){
        tmp <- cbind(ISmatrices[[i]][, isNames ])
        colnames(tmp) <- isNames
        mIS <- cbind(mIS, tmp)
      }
    } #end for loop

    mXis <- dropvar(cbind(mX,mIS), tol=tol, LAPACK=LAPACK,
      silent=print.searchinfo)
    mXis <- zoo(mXis, order.by=y.index)
  } #end if(length(ISfinalmodels)>0)

  ##gum and gets:
  y <- zoo(y, order.by=y.index)
  mod <- arx(y, mxreg=mXis, vcov.type=vcov.type,
    qstat.options=qstat.options, user.diagnostics=user.diagnostics,
    tol=tol, LAPACK=LAPACK, plot=FALSE)
  getsis <- getsm(mod, keep=mxkeep, t.pval=t.pval,
    do.pet=do.pet, wald.pval=wald.pval, ar.LjungB=ar.LjungB,
    arch.LjungB=arch.LjungB, normality.JarqueB=normality.JarqueB,
    user.diagnostics=user.diagnostics, info.method=info.method,
    include.empty=include.empty, max.regs=max.regs,
    print.searchinfo=print.searchinfo, plot=FALSE)

  ##names of retained impulses, mX colnames:
  ISnames <- setdiff(getsis$aux$mXnames, mXnames)
  if(length(ISnames)==0){ ISnames <- NULL }
  colnames(getsis$aux$mX) <- getsis$aux$mXnames

  ##return:
  getsis$gets.type <- "isat"
  getsis$call <- isat.call
  getsis <- c(list(ISfinalmodels=ISfinalmodels,
    ISnames=ISnames), getsis)
  getsis$aux$t.pval <- t.pval #needed for biascorr
  class(getsis) <- "isat"
  if(alarm){ alarm() }
  if( is.null(plot) ){
    plot <- getOption("plot")
    if( is.null(plot) ){ plot <- FALSE }
  }
  if(plot){ plot.isat(getsis, coef.path=TRUE) }
  return(getsis)
}

Try the gets package in your browser

Any scripts or data that you put into this service are public.

gets documentation built on May 30, 2017, 4:09 a.m.