R/ts.R

In theGreatWhiteShark/climex: Estimating Extreme Events in Climatic Time Series (using Data of the DWD)

### ts.R - Helper functions to handle and preprocess time series in
##    the analysis
##' @title Akaike information criterion
##' @description Calculates the Akaike information criterion of a
##'   \emph{climex.fit.gev} or \emph{climex.fit.gpd} class object or a
##'   list of those.
##' @details The \emph{climex.fit.gev} object is of identical
##'   structure as the output of the \code{\link[stats]{optim}}
##    function.
##'  
##' @param x Either an object returned by the \code{\link{fit.gev}} or
##'   \code{\link{fit.gpd}} function or a list of such elements.
##'
##' @seealso \code{\link{bic}}, \code{\link{fit.gev}},
##'   \code{\link{fit.gpd}}
##' @family ts
##'
##' @export
##' 
##' @return Numeric value, if the input is a single fitting object, or
##'   a numerical vector, if \strong{x} is a list of such objects.
##' @author Philipp Mueller
aic <- function( x ){
  UseMethod( "aic" )
}
##' @title Akaike information criterion
##' @description Calculates the Akaike information criterion of a
##'   list of \emph{climex.fit.gev} or \emph{climex.fit.gpd} class
##'   objects.
##' @details The \emph{climex.fit.gev} object is of identical
##'   structure as the output of the \code{\link[stats]{optim}}
##'   function.
##'  
##' @param x Either a list of objects returned by the
##'   \code{\link{fit.gev}} or \code{\link{fit.gpd}} function.
##'
##' @seealso \code{\link{bic}}, \code{\link{fit.gev}},
##'   \code{\link{fit.gpd}}
##' @family ts
##'
##' @export
##' 
##' @return Numerical vector
##' @author Philipp Mueller
aic.list <- function( x ){
  x.result <- Reduce( c, lapply( x, aic ) )
  return( x.result )
}

##' @title Akaike information criterion
##' @description Calculates the Akaike information criterion of a
##'   \emph{climex.fit.gev} class object.
##' @details The \emph{climex.fit.gev} object is of identical
##'   structure as the output of the \code{\link[stats]{optim}}
##'   function.
##'  
##' @param x An object returned by the \code{\link{fit.gev}} function.
##'
##' @seealso \code{\link{bic}}, \code{\link{fit.gev}}
##' @family ts
##'
##' @export
##' 
##' @return Numeric value
##' @author Philipp Mueller
aic.climex.fit.gev <- function( x ){
  return( 2* x$value + 2* length( x$par ) )
}

##' @title Akaike information criterion
##' @description Calculates the Akaike information criterion of a
##'   \emph{climex.fit.gpd} class object.
##' @details The \emph{climex.fit.gpd} object is of identical
##'   structure as the output of the \code{\link[stats]{optim}}
##'   function.
##'  
##' @param x An object returned by the \code{\link{fit.gpd}} function.
##'
##' @seealso \code{\link{bic}}, \code{\link{fit.gpd}}
##' @family ts
##'
##' @export
##' 
##' @return Numeric value
##' @author Philipp Mueller
aic.climex.fit.gpd <- function( x ){
  return( 2* x$value + 2* length( x$par ) )
}

##' @title Bayesian information criterion
##' @description Calculates the Bayesian information criterion of a
##'   \emph{climex.fit.gev} or \emph{climex.fit.gpd} class object or a
##'   list of those.
##' @details The \emph{climex.fit.gev} object is of identical
##'   structure as the output of the \code{\link[stats]{optim}}
##'   function.
##'  
##' @param x Either an object returned by the \code{\link{fit.gev}} or
##'   \code{\link{fit.gpd}} function or a list of such elements.
##'
##' @seealso \code{\link{aic}}, \code{\link{fit.gev}},
##'   \code{\link{fit.gpd}}
##' @family ts
##'
##' @export
##' 
##' @return Numeric value, if the input is a single fitting object, or
##'   a numerical vector, if \strong{x} is a list of such objects.
##' @author Philipp Mueller
bic <- function( x ){
  UseMethod( "bic" )
}
##' @title Bayesian information criterion
##' @description Calculates the Bayesian information criterion of a
##'   list of \emph{climex.fit.gev} or \emph{climex.fit.gpd} class
##'   objects.
##' @details The \emph{climex.fit.gev} object is of identical
##'   structure as the output of the \code{\link[stats]{optim}}
##'   function.
##'  
##' @param x Either a list of objects returned by the
##'   \code{\link{fit.gev}} or \code{\link{fit.gpd}} function.
##'
##' @seealso \code{\link{aic}}, \code{\link{fit.gev}},
##'   \code{\link{fit.gpd}}
##' @family ts
##'
##' @export
##' 
##' @return Numerical vector
##' @author Philipp Mueller
bic.list <- function( x ){
    x.result <- Reduce( c, lapply( x, bic ) )
  return( x.result )
}

##' @title Bayesian information criterion
##' @description Calculates the Bayesian information criterion of a
##'   \emph{climex.fit.gev} class object.
##' @details The \emph{climex.fit.gev} object is of identical
##'   structure as the output of the \code{\link[stats]{optim}}
##'   function.
##'  
##' @param x An object returned by the \code{\link{fit.gev}} function.
##'
##' @seealso \code{\link{aic}}, \code{\link{fit.gev}}
##' 
##' @family ts
##'
##' @export
##' 
##' @return Numeric value
##' @author Philipp Mueller
bic.climex.fit.gev <- function( x ){
  return( 2* x$value + length( x$par )* log( length( x$x ) ) )
}

##' @title Bayesian information criterion
##' @description Calculates the Bayesian information criterion of a
##'   \emph{climex.fit.gpd} class object.
##' @details The \emph{climex.fit.gpd} object is of identical
##'   structure as the output of the \code{\link[stats]{optim}}
##'   function.
##'  
##' @param x An object returned by the \code{\link{fit.gpd}} function.
##'
##' @seealso \code{\link{aic}}, \code{\link{fit.gpd}}
##' 
##' @family ts
##'
##' @export
##' 
##' @return Numeric value
##' @author Philipp Mueller
bic.climex.fit.gpd <- function( x ){
  return( 2* x$value + length( x$par )* log( length( x$x ) ) )
}

##' @title Remove incomplete years
##' @description Removes all years, which contain either a NA or are
##'   incomplete (derived from the amount of time stamps within a
##'   year). 
##'
##' @details Since incomplete years are detected via the
##'   difference in the time stamps of neighbouring points, the user
##'   has to provide the basic time unit. (Numerical value of the
##'   result of the \code{\link[base]{diff}} function applied to two
##'   consecutive points.) For daily data the value is 1.
##' 
##' @param x Either a \pkg{xts} class object or a list of those.
##' @param time.unit Basic time unit in days. Default = 1.
##' @param mc.cores A numerical input specifying the number of cores
##'   to use for the multi core application of the function (see
##'   \code{\link[parallel]{detectCores}}). This functionality is only
##'   available if the input is a list of different objects. If NULL,
##'   the function will be calculated classically with only one core.
##'   Default = NULL.
##'
##' @seealso \code{\link{remove.seasonality}}
##'
##' @family ts
##' @export
##' @importFrom xts xts
##' @importFrom parallel mclapply
##' @import lubridate
##' 
##' @return Object of the same class as the input \strong{x}.
##' @author Philipp Mueller
remove.incomplete.years <- function( x, time.unit = 1,
                                    mc.cores = NULL ){
  UseMethod( "remove.incomplete.years" )
}
##' @title Remove incomplete years
##' @description Removes all years, which contain either a NA or are
##'   incomplete (derived from the amount of time stamps within a
##'   year). 
##'
##' @details Since incomplete years are detected via the
##'   difference in the time stamps of neighbouring points, the user
##'   has to provide the basic time unit. (Numerical value of the
##'   result of the \code{\link[base]{diff}} function applied to two
##'   consecutive points.) For daily data the value is 1.
##' 
##' @param x A list of \pkg{xts} class objects.
##' @param time.unit Basic time unit in days. Default = 1.
##' @param mc.cores A numerical input specifying the number of cores
##'   to use for the multi core application of the function (see
##'   \code{\link[parallel]{detectCores}}). This functionality is only
##'   available if the input is a list of different objects. If NULL,
##'   the function will be calculated classically with only one core.
##'   Default = NULL.
##'
##' @seealso \code{\link{remove.seasonality}}
##'
##' @family ts
##' @export
##' @importFrom xts xts
##' @importFrom parallel mclapply
##' @import lubridate
##' 
##' @return Object of the same class as the input \strong{x}.
##' @author Philipp Mueller
remove.incomplete.years.list <- function( x, time.unit = 1,
                                         mc.cores = NULL ){
  if ( !is.null( mc.cores ) ){
    return( mclapply( x, remove.incomplete.years,
                     time.unit = time.unit, mc.cores = mc.cores ) )
  } else {
    return( lapply( x, remove.incomplete.years,
                   time.unit = time.unit, mc.cores = mc.cores ) )
  }
}
##' @title Remove incomplete years
##' @description Removes all years, which contain either a NA or are
##'   incomplete (derived from the amount of time stamps within a
##'   year). 
##'
##' @details Since incomplete years are detected via the
##'   difference in the time stamps of neighbouring points, the user
##'   has to provide the basic time unit. (Numerical value of the
##'   result of the \code{\link[base]{diff}} function applied to two
##'   consecutive points.) For daily data the value is 1.
##' 
##' @param x A \pkg{xts} class object.
##' @param time.unit Basic time unit in days. Default = 1.
##' @param mc.cores A numerical input specifying the number of cores
##'   to use for the multi core application of the function (see
##'   \code{\link[parallel]{detectCores}}). This functionality is only
##'   available if the input is a list of different objects. If NULL,
##'   the function will be calculated classically with only one core.
##'   Default = NULL.
##'
##' @seealso \code{\link{remove.seasonality}}
##'
##' @family ts
##' @export
##' @importFrom xts xts
##' @importFrom parallel mclapply
##' @import lubridate
##' 
##' @return Object of the same class as the input \strong{x}.
##' @author Philipp Mueller
remove.incomplete.years.xts <- remove.incomplete.years.default <-
  function( x, time.unit = 1, mc.cores = NULL ){
    if ( !any( class( x ) == "xts" ) )
      stop(
          "Provided element of the wrong class in remove.incomplete.years. Provide object of class xts instead!" )
    ## Removes all values containing a NA, R's general type for
    ## missing data. 
    x <- x[ !is.na( x ) ]

    ## There is data available.
    if ( length( x ) > 0 ){
      ## The missing points are detected by looking for differences in 
      ## the time stamps of neighbouring points. If their are more
      ## than time.unit apart, there are some missing values in
      ## between,
      x.index <- index( x )
      x.index.diff <- diff( x.index )

      ## This numerical vector will hold all years, which will be
      ## removed in the end
      all.incomplete.years <- numeric()
      ## There might be as well several gaps in the data. Therefore we
      ## will perform a loop over all gaps detected via the `diff`
      ## function. 
      for ( ii in seq( 1, ( length( x.index ) - 1 ) )[
                              x.index.diff > time.unit ] ){

        ## Since the gap itself won't tell us which years are missing
        ## in the data, we have to derive them from the entry right
        ## before (ii) and after (ii+1) the gap.
        ## But it can happen that the entry
        ## before the gap is the last value within one year or the one
        ## after the gap the first one within one year. These cases
        ## have to be treated with a little bit more care. For all
        ## others a sequences of years from the one before and after
        ## the gap will be removed.
        if ( year( x[ ii ] ) !=
             year( index( x[ ii ] ) + ddays( time.unit ) ) ){
          ## The entry before the gap is at the end of the year.
          year.start <- year( index( x[ ii ] ) + ddays( time.unit ) )
        } else {
          year.start <- year( x[ ii ] )
        }
        if ( year( x[ ii + 1 ] ) !=
             year( index( x[ ii + 1 ] ) - ddays( time.unit ) ) ){
          ## Entry after the gap is at the very beginning of a year.
          year.end <- year( index( x[ ii + 1 ] ) - ddays( time.unit ) )
        } else {
          year.end <- year( x[ ii + 1 ] )
        }

        all.incomplete.years <- c( all.incomplete.years,
                                  seq( year.start, year.end, 1 ) )
      }
      ## Check whether the first entry is at the beginning of the year
      ## and the last one is at the end of the year. If not, those
      ## years will be considered incomplete and removed too.
      if ( year( x[ 1 ] ) ==
           year( index( x[ 1 ] ) - ddays( time.unit ) ) ){
        ## The first entry is not at the beginning of a year.
        all.incomplete.years <- c( all.incomplete.years,
                                  year( x[ 1 ] ) )
      }
      if ( year( x[ length( x ) ] ) ==
           year( index( x[ length( x ) ] ) + ddays( time.unit ) ) ){
        ## The last entry is not at the end of a year.
        all.incomplete.years <- c( all.incomplete.years,
                                  year( x[ length( x ) ] ) )
      }
      return( x[ !( year( x ) %in% all.incomplete.years ) ] )
    } else {
      ## There is no data left
      return( x )
    }
  }

##' @title Checking for N complete years of data
##' 
##' @description This function tests whether the series do has at
##'   least \strong{number.of.years} complete years of data.
##' 
##' @details It uses the \code{\link{remove.incomplete.years}}
##'   to check for the completeness of the individual years.
##'
##' @param x Either a time series of class \pkg{xts} or a list of
##'   them.
##' @param number.of.years The minimum number of complete years in the
##'   series. Default = 30.
##' @param mc.cores A numerical input specifying the number of cores
##'   to use for the multi core application of the function (see
##'   \code{\link[parallel]{detectCores}}). This functionality is only
##'   available if the input is a list of different objects. If NULL,
##'   the function will be calculated classically with only one core.
##'   Default = NULL.
##'
##' @export
##' 
##' @importFrom parallel mclapply
##'
##' @family ts
##'
##' @return If the input is of type \pkg{xts}, the function returns
##'   TRUE or FALSE. If the input, on the other hand, is a list of
##'   objects of type \pkg{xts}, it returns a trimmed list
##'   containing only those elements, which indeed have more than N
##'   years of complete data.
##' @author Philipp Mueller 
check.completeness <- function( x, number.of.years = 30,
                               mc.cores = NULL){
  UseMethod( "check.completeness" )
}
##' @title Checking for N complete years of data
##' @description Removes all class \pkg{xts} time series from a
##'   list, which are not of at least \strong{number.of.years}
##'   complete years of data.
##'
##' @param x A list of \pkg{xts} class objects.
##' @param number.of.years The minimum number of complete years in the
##'   series. Default = 30.
##' @param mc.cores A numerical input specifying the number of cores
##'   to use for the multi core application of the function (see
##'   \code{\link[parallel]{detectCores}}). This functionality is only
##'   available if the input is a list of different objects. If NULL,
##'   the function will be calculated classically with only one core.
##'   Default = NULL.
##'
##' @export
##'
##' @importFrom parallel mclapply
##'
##' @return Trimmed list containing only those elements, which indeed
##'   have more than N years of complete data.
##' @author Philipp Mueller
check.completeness.list <- function( x, number.of.years = 30,
                                    mc.cores = NULL ){
  if ( !is.null( mc.cores ) ){
    ## Apply the checks to all elements of the list.
    x.list.check <-
      Reduce( c, mclapply( x, check.completeness.xts,
                          number.of.years = number.of.years,
                          mc.cores = mc.cores ) )
  } else {
    ## Apply the checks to all elements of the list.
    x.list.check <-
      Reduce( c, lapply( x, check.completeness.xts,
                        number.of.years = number.of.years,
                        mc.cores = mc.cores ) )
  }
  ## Return only those element matching the conditions
  return( x[ x.list.check ] )
}
##' @title Checking for N complete years of data
##' 
##' @description This function tests whether a \pkg{xts} class object
##'   has at least \strong{number.of.years} complete years of data.
##' 
##' @details It uses the \code{\link{remove.incomplete.years}}
##'   to check for the completeness of the individual years.
##'
##' @param x A time series of class \pkg{xts}.
##' @param number.of.years The minimum number of complete years in the
##'   series. Default = 30.
##' @param mc.cores A numerical input specifying the number of cores
##'   to use for the multi core application of the function (see
##'   \code{\link[parallel]{detectCores}}). This functionality is only
##'   available if the input is a list of different objects. If NULL,
##'   the function will be calculated classically with only one core.
##'   Default = NULL.
##'
##' @importFrom xts is.xts
##' @importFrom parallel mclapply
##'
##' @export
##' 
##' @family data.cleaning
##'
##' @return Logical
##' @author Philipp Mueller 
check.completeness.xts <- check.completeness.default <-
  function( x, number.of.years = 30, mc.cores = NULL ){
    ## Only accept objects of class 'xts'
    if ( !is.xts( x ) ){
      stop( "Only objects of class 'xts' accepted!" )
    }
    if ( !is.numeric( number.of.years ) ){
      stop( "The number of years has to be a 'numerical' value!" )
    }

    ## Remove all incomplete years
    x.complete.years <- remove.incomplete.years( x )

    ## Check whether the remaining number of years are at least of the
    ## specified number.
    x.check <- length( unique(
        lubridate::year( x.complete.years ) ) ) >=
      number.of.years
    return( x.check )
  }

##' @title Remove seasonality
##' @description Calculates the seasonal component of a time series
##'   and subtracts it from the original. 
##'
##' @details The function \code{\link[stats]{stl}} with the argument
##'   \code{s.window = 12} and a 
##'   conversion of the input into a \emph{ts} class object of daily
##'   data is used to calculate the seasonal component. This should be
##'   replaced by a more sophisticated solution as soon I digged
##'   deeper into the field of deseasonalization.
##'   \code{\link{remove.incomplete.years}} is used
##'   to remove incomplete years from the data set. This ensures a
##'   better calculation of the seasonal component but also requires
##'   to forecast it to the length of the original data set and align
##'   it at the right place for subtraction.
##'
##'   This function can also be applied to a list of \pkg{xts} class
##'   objects.
##'
##' @param x Either an object of class \pkg{xts} or a list of those.
##' @param mc.cores A numerical input specifying the number of cores
##'   to use for the multi core application of the function (see
##'   \code{\link[parallel]{detectCores}}). This functionality is only
##'   available if the input is a list of different objects. If NULL,
##'   the function will be calculated classically with only one core.
##'   Default = NULL.
##'
##' @family ts
##'
##' @importFrom xts xts
##' @importFrom parallel mclapply
##' @import lubridate
##' @return Some class as input.
##' @author Philipp Mueller
remove.seasonality <- function( x, mc.cores = NULL ){
  UseMethod( "remove.seasonality" )
}
##' @title Remove seasonality
##' @description Calculates the seasonal component of a time series
##'   and subtracts it from the original. 
##'
##' @details The function \code{\link[stats]{stl}} with the argument
##'   \code{s.window = 12} and a 
##'   conversion of the input into a \emph{ts} class object of daily
##'   data is used to calculate the seasonal component. This should be
##'   replaced by a more sophisticated solution as soon I digged
##'   deeper into the field of deseasonalization.
##'   \code{\link{remove.incomplete.years}} is used
##'   to remove incomplete years from the data set. This ensures a
##'   better calculation of the seasonal component but also requires
##'   to forecast it to the length of the original data set and align
##'   it at the right place for subtraction.
##'
##'   This function can also be applied to a list of \pkg{xts} class
##'   objects.
##'
##' @param x A list of objects of class \pkg{xts}.
##' @param mc.cores A numerical input specifying the number of cores
##'   to use for the multi core application of the function (see
##'   \code{\link[parallel]{detectCores}}). This functionality is only
##'   available if the input is a list of different objects. If NULL,
##'   the function will be calculated classically with only one core.
##'   Default = NULL.
##'
##' @family ts
##'
##' @importFrom xts xts
##' @importFrom parallel mclapply
##' @import lubridate
##' @return Some class as input.
##' @author Philipp Mueller
remove.seasonality.list <- function( x, mc.cores = NULL ){
  if ( !is.null( mc.cores ) ){
    return( mclapply( x, remove.seasonality, mc.cores = mc.cores ) )
  } else {
    return( lapply( x, remove.seasonality, mc.cores = mc.cores ) )
  }
}
##' @title Remove seasonality
##' @description Calculates the seasonal component of a time series
##'   and subtracts it from the original. 
##'
##' @details The function \code{\link[stats]{stl}} with the argument
##'   \code{s.window = 12} and a 
##'   conversion of the input into a \emph{ts} class object of daily
##'   data is used to calculate the seasonal component. This should be
##'   replaced by a more sophisticated solution as soon I digged
##'   deeper into the field of deseasonalization.
##'   \code{\link{remove.incomplete.years}} is used
##'   to remove incomplete years from the data set. This ensures a
##'   better calculation of the seasonal component but also requires
##'   to forecast it to the length of the original data set and align
##'   it at the right place for subtraction.
##'
##'   This function can also be applied to a list of \pkg{xts} class
##'   objects.
##'
##' @param x An object of class \pkg{xts}.
##' @param mc.cores A numerical input specifying the number of cores
##'   to use for the multi core application of the function (see
##'   \code{\link[parallel]{detectCores}}). This functionality is only
##'   available if the input is a list of different objects. If NULL,
##'   the function will be calculated classically with only one core.
##'   Default = NULL.
##'
##' @family ts
##'
##' @importFrom xts xts
##' @importFrom parallel mclapply
##' @import lubridate
##' @return Some class as input.
##' @author Philipp Mueller
remove.seasonality.xts <- remove.seasonality.default <-
  function( x, mc.cores = NULL ){
    if ( !any( class( x ) == "xts" ) )
      stop( "Provided element of the wrong class in remove.seasonality. Provide object of class xts instead!" )  
    ## To ensure reasonal results 
    if ( min( x, na.rm = TRUE ) == -999 )
      x[ x == -999 ] <- NA
    x.clean <- remove.incomplete.years( x )
    
    ## When e.g. no maximal temperatures are provided for a station
    ## the function will exist
    if ( length( x.clean ) == 0 )
      return( list( x = NA, modified = FALSE ) )
    
    ## Was time series modified?
    ifelse( length( x ) == length( x.clean ),
           MODIFIED <- FALSE, MODIFIED <- TRUE )
    
    ## maybe to much data was remove
    if ( length( x ) < 365*5 ){
      warning( "After the removal of the incomplete years less than five years are remaining in the data set in remove.seasonality" )
      return( list( x = NA, modified = MODIFIED ) )
    }
    
    ## Seasonal component of the time series
    x.seas <- stats::stl( stats::ts( as.numeric( x.clean ),
                                    frequency = 365 ),
                         s.window = 12 )$time.series[ , 1 ]
    
    ## Usually the seasonal term is more than just a constant
    ## sinusoidal function but exhibits some kind of trend
    ## itself. That's why it's only logical to use the same or the
    ## nearest year in the seasonal component for subtracting of
    ## original series. Subtracting all years which are also present
    ## in the seasonal component
    x.clean.years <- unique( year( x.clean ) )
    x.years <- unique( year( x ) )
    for ( yy in x.clean.years )
      x[ year( x ) == yy ] <- x[ year( x ) == yy ] -
        x.seas[ which( year( x.clean ) == yy ) ]
    
    ## Looping over all years don't accounted yet and substracting the
    ## year nearst in x.seas to the original one. The alignment of the
    ## to data sets has to be handled with care!
    for ( yy in x.years[ !x.years %in% x.clean.years ] ){
      ## closest seasonal component; Caution, can be a vector!
      ## if the year in x is a leap year one has to search for the
      ## closest leap year.
      ## If the closest year has not 366 days its distance is set to
      ## 999
      if ( length( x[ year( x ) == yy ] ) == 366 ){
        closest.year.vec <-  abs( x.clean.years - yy )
        for ( cc in 1 : length( closest.year.vec ) ){
          if ( length( x.clean[
              year( x.clean ) == x.clean.years[
                                     which( closest.year.vec ==
                                            min( closest.year.vec )
                                           ) ] ] )
              == 366 ){
            closest.year.pos <- which( closest.year.vec ==
                                       min( closest.year.vec ) )
            break
          } else
            closest.year.vec[
                which( closest.year.vec ==
                       min( closest.year.vec ) ) ] <- 999
        }
      } else {
        closest.year.pos <- which( abs( x.clean.years - yy ) ==
                                   min( abs( x.clean.years - yy ) ) )
      }
      closest.year <- x.clean.years[ closest.year.pos[ 1 ] ] 
      
      ## subtracts only those days of the closest year in x.seas which
      ## are also present in the x time series. ydays give the
      ## calender date (1--366)
      days.clean <- yday( x.clean[ year( x.clean ) == closest.year ] )
      days.x <- yday( x[ year( x ) == yy ] )
      ## Even if the length of x is not 366 it can still be a leap
      ## year and the 31.12. will be of the value 366. This way the
      ## seasonal part might not has the some number of entries (the
      ## closest year is not a leap year)
      if ( min( days.x ) == 1 && max( days.x ) == 366 &&
           length( days.x ) < 366 ){
        ## Somewhere at least one day is missing. This "step" has to
        ## be detected and the all days after the step have to be
        ## reduced by one
        step <- which( days.x - c( 1 : length( days.x ) ) > 0 )[ 1 ]
        days.x[ step : length( days.x ) ] <- days.x[
            step : length( days.x ) ] - 1
      } else if ( max( days.x ) == 366 && max( days.clean ) != 366 )
        days.x <- days.x - 1
      
      if ( min( days.x ) == 0 )
        stop( "Wait. Something went wrong in remove.seasonality. Check for the correct alignment of the seasonal part and the original time series!" )
      x[ year( x ) == yy ] <- x[ year( x ) == yy ] - 
        as.numeric( x.seas[ year( x.clean ) == closest.year ][
            days.clean %in% days.x ] )
    }
    return( list( x = x, modified = MODIFIED ) )
  }

##' @title Calculates the mode of a PDF
##' @description Calculates the mode of a PDF
##'
##' @param x Numerical input.
##'
##' @return Numerical
##' @author Philipp Mueller 
mode <- function( x ){
  ## Approximation of the mode (most probable/frequent element) of
  ## a time series 
  ## Working with unique() does not yield good results for a sample of
  ## a continuous distribution, therefore using the Kernel density
  ## approximation and determining the peak of the PDF
  x.pdf <- stats::density( x )
  return( x.pdf$x[ which.max( x.pdf$y ) ] )
}

##' @title Anomalies of a time series
##' @description Calculates the anomalies of an object of class
##'   \pkg{xts} or a list of such objects.
##'
##' @details Construction via the subtraction of the mean
##'   value of the specific date.
##'
##'   Uses the \code{\link[lubridate]{yday}}.
##' 
##' @param x Either a time series of class \pkg{xts} or a list of
##'   them.
##' @param mc.cores A numerical input specifying the number of cores
##'   to use for the multi core application of the function (see
##'   \code{\link[parallel]{detectCores}}). This functionality is only
##'   available if the input is a list of different objects. If NULL,
##'   the function will be calculated classically with only one core.
##'   Default = NULL.
##'
##' @family ts
##' @export
##' @importFrom xts xts
##' @importFrom lubridate yday
##' @importFrom parallel mclapply
##' 
##' @return Same class as the input
##' @author Philipp Mueller
##' 
anomalies <- function( x, mc.cores = NULL ){
  UseMethod( "anomalies" )
}
##' @title Anomalies of a time series
##' @description Calculates the anomalies of an object of class
##'   \pkg{xts} or a list of such objects.
##'
##' @details Construction via the subtraction of the mean
##'   value of the specific date.
##'
##'   Uses the \code{\link[lubridate]{yday}}.
##' 
##' @param x List of time series of class \pkg{xts}
##' @param mc.cores A numerical input specifying the number of cores
##'   to use for the multi core application of the function (see
##'   \code{\link[parallel]{detectCores}}). This functionality is only
##'   available if the input is a list of different objects. If NULL,
##'   the function will be calculated classically with only one core.
##'   Default = NULL.
##'
##' @family ts
##' @export
##' @importFrom xts xts
##' @importFrom lubridate yday
##' @importFrom parallel mclapply
##' 
##' @return Same class as the input
##' @author Philipp Mueller
##' 
anomalies.list <- function( x, mc.cores = NULL ){
  if ( !is.null( mc.cores ) ){
    return( mclapply( x, anomalies, mc.cores = mc.cores ) )
  } else {
    return( lapply( x, anomalies, mc.cores = mc.cores ) )
  }
}
##' @title Anomalies of a time series
##' @description Calculates the anomalies of an object of class
##'   \pkg{xts} or a list of such objects.
##'
##' @details Construction via the subtraction of the mean
##'   value of the specific date.
##'
##'   Uses the \code{\link[lubridate]{yday}}.
##' 
##' @param x A time series of class \pkg{xts}
##' @param mc.cores A numerical input specifying the number of cores
##'   to use for the multi core application of the function (see
##'   \code{\link[parallel]{detectCores}}). This functionality is only
##'   available if the input is a list of different objects. If NULL,
##'   the function will be calculated classically with only one core.
##'   Default = NULL.
##'
##' @family ts
##' @export
##' @importFrom xts xts
##' @importFrom lubridate yday
##' @importFrom parallel mclapply
##' 
##' @return Same class as the input
##' @author Philipp Mueller
anomalies.xts <- anomalies.default <- function( x, mc.cores = NULL ){
  if ( !any( class( x ) == "xts" ) )
    stop( "Only for the class 'xts' the anomalies can be calculated" )
  ## base::ave seems to calculate takes a vector and a factor vector
  ## of the same length. Then it calculates the mean for all values
  ## sharing a factor and is placing the mean value in the element
  ## corresponding to the position of the original one.
  x.anomalies <- x - stats::ave( x, as.factor( yday( x ) ),
                                FUN = function( xx )
                                  mean( xx, na.rm = TRUE ) )
  return( x.anomalies )
}
## End of ts.R