R/ubgenerator.r

In jlivsey/sigex: Ecce Signum: Multivariate Signal Extraction

#' Compute the product of unit root differencing operators.
#'
#' @param period  Possibly non-integer period of sinusoids
#' @param trunc.len  Gives the maximal index of k, and is set to [period/2]
#'     when trunc.len = NULL
#' @param m  Large integer giving number of cepstral coefficients
#' @param rho Parameter between 0 and 1, can equal 1 for unit root case
#'
#' @return wolds: coefficients of product polynomial.
#' @export
#'

ubgenerator <- function(period,trunc.len,m,rho)
{

  ##########################################################################
  #
  #	ubgenerator
  # 	    Copyright (C) 2020  Tucker McElroy
  #
  #    This program is free software: you can redistribute it and/or modify
  #    it under the terms of the GNU General Public License as published by
  #    the Free Software Foundation, either version 3 of the License, or
  #    (at your option) any later version.
  #
  #    This program is distributed in the hope that it will be useful,
  #    but WITHOUT ANY WARRANTY; without even the implied warranty of
  #    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  #    GNU General Public License for more details.
  #
  #    You should have received a copy of the GNU General Public License
  #    along with this program.  If not, see <https://www.gnu.org/licenses/>.
  #
  ############################################################################

  ################# Documentation #####################################
  #
  #	Purpose: compute the product of differencing operators.
  # Background: the goal is to annihilate periodic components of
  #   frequency 2*pi*k/period for k = 1, 2, ..., [period/2].
  #   Compute product_k  (1 - 2 rho cos(2*pi*k/period) z + rho^2 z^2) by using
  #   cepstral method for higher accuracy.
  #	Inputs:
  #   period: possibly non-integer period of sinusoids
  #   trunc.len: gives the maximal index of k, and is set to [period/2]
  #     when trunc.len = NULL
  #   m: large integer giving number of cepstral coefficients
  #   rho: parameter between 0 and 1, can equal 1 for unit root case
  #	Outputs:
  #	  wolds: coefficients of product polynomial.
  #
  ##########################################################################################

  ceps2wold <- function(ceps,q)
  {
    m <- length(ceps)
    if(q > m) {	ceps <- c(ceps,rep(0,q-m)) }
    wold <- 1
    wolds <- wold
    for(j in 1:q)
    {
      wold <- sum(seq(1,j)*ceps[1:j]*wolds[j:1])/j
      wolds <- c(wolds,wold)
    }
    return(wolds)
  }

  half.len <- floor(period/2)
  if(length(trunc.len)==0) { trunc.len <- half.len }
  ceps <- rep(0,m)

  for(ell in 1:m)
  {
    ceps[ell] <- -2*(rho^ell)*sum(cos(2*pi*ell*seq(1,trunc.len)/period))/ell
  }
  wolds <- ceps2wold(ceps,2*trunc.len)

  return(wolds)
}