R/kernel_from_eig.R

In tomasrubin/specsimfts: Spectral domain simulation methods for functional time series

#' Takes as its arguments the eigendecomposition of a self-adjoint operator and returns the function of two variables as the kernel of this operator.
#' 
#'
#' @title Produces a kernel function from given eigendecomposition
#' @param eigenvalues The eigendecomposition can be defined either as (i) a list of finite number of eigenvalues/eigenfunctions, or (ii) a function returning the eigenvalue/eigenfunction of given order. In the case (i), the \code{eigenvalues} parameter is simply a vector of eigenvalues. In the case (ii), the \code{eigenvalues} parameter is a function of the variable \code{n} which returns the \code{n}-th eigenvalue.
#' @param eigenfunctions See \code{eigenvalues} before. In the case (i), the \code{eigenfunctions} parameter is a list of functions with the same length as the parameter \code{eigenvalues}. Each function is a function of the variable \code{x} that returns the value of the eigenfunction at point \code{x}. The order of the functions in the list determine the order of eigenfunctions. In the case (ii), the \code{eigenfunctions} parameter is a function of two variables, \code{n} and \code{x}, that returns the value of the \code{n}-th eigenfunction at the point \code{x}.
#' @param n_pc_if_functions In the case (ii) where the eigenvalues and the eigenfunctions are defined as functions representing the infinite rank eigendecomposition, take only \code{n_pc_if_functions} eigenfunctions. In the case (i), this parameter has no effect.
#' @return the function of two variables, \code{x} and \code{y}, returning the value of the kernel at point (\code{x},\code{y})
#' @references Rubin, Panaretos. \emph{Simulation of stationary functional time series with given spectral density}. arXiv, 2020
#' @seealso \code{\link{FARFIMA_simulate}}
#' @examples 
#' # define as infinite rank operator, definition as functions
#' eigenvalues <- function(n) { 1/((n-0.5)*pi)^2 }
#' eigenfunctions <- function(n,x) { sqrt(2)*sin((n-0.5)*pi*x) }
#' 
#' # # Alternative: define as finite rank operator, definition as lists
#' # eigenvalues <- c(1, 0.6, 0.3, 0.1, 0.1, 0.1, 0.05, 0.05, 0.05, 0.05)
#' # eigenfunctions <- list(
#' #   function(x){ sin(2*pi*x) },
#' #   function(x){ cos(2*pi*x) },
#' #   function(x){ sin(4*pi*x) },
#' #   function(x){ cos(4*pi*x) },
#' #   function(x){ sin(6*pi*x) },
#' #   function(x){ cos(6*pi*x) },
#' #   function(x){ sin(8*pi*x) },
#' #   function(x){ cos(8*pi*x) },
#' #   function(x){ sin(10*pi*x) },
#' #   function(x){ cos(10*pi*x) }
#' #  )
#' 
#' 
#' # create the kernel
#' ker <- kernel_from_eig(eigenvalues,eigenfunctions)
#' 
#' # evaluate the kernel at point (x,y) = (0.1, 0.6)
#' ker(0.1, 0.6)
#' @export
kernel_from_eig <- function( eigenvalues, eigenfunctions, n_pc_if_functions=100 ){
  
  if (is.list(eigenfunctions)){
    # if it's lst => low rank specification
    return(function(x,y){
      r <- 0
      for (ii in 1:length(eigenvalues)){
        r <- r + eigenvalues[ii] * eigenfunctions[[ii]](x) * eigenfunctions[[ii]](y)
      }
      return(r)
    })
  } else {
    # if they're functions
    return(function(x,y){
      r<-0
      for (ii in 1:n_pc_if_functions){
        r <- r + eigenvalues(ii) * eigenfunctions(ii,x) * eigenfunctions(ii,y)
      }
      return(r)
    })
  }
}