# R/ls_whittle_loglik.R In LSTS: Locally Stationary Time Series

#### Documented in LS.whittle.loglik

#' @title Locally Stationary Whittle log-likelihood Function
#' @description This function computes Whittle estimator for LS-ARMA and
#' LS-ARFIMA models, in data with mean zero. If mean is not zero, then it is
#' subtracted to data.
#' @details
#' The estimation of the time-varying parameters can be carried out by means of
#' the Whittle log-likelihood function proposed by Dahlhaus (1997),
#' \deqn{L_n(\theta) = \frac{1}{4\pi}\frac{1}{M} \int_{-\pi}^{\pi}
#' \bigg\{log f_{\theta}(u_j,\lambda) +
#' \frac{I_N(u_j, \lambda)}{f_{\theta}(u_j,\lambda)}\bigg\}\,d\lambda}
#' where \eqn{M} is the number of blocks, \eqn{N} the length of the series per
#' block, \eqn{n =S(M-1)+N}, \eqn{S} is the shift from block to block,
#' \eqn{u_j =t_j/n}, \eqn{t_j =S(j-1)+N/2}, \eqn{j =1,\ldots,M} and
#' \eqn{\lambda} the Fourier frequencies in the block
#' (\eqn{2\,\pi\,k/N}, \eqn{k = 1,\ldots, N}).
#' @param x (type: numeric) parameter vector.
#' @param series (type: numeric) univariate time series.
#' @param order (type: numeric) vector corresponding to \code{ARMA} model
#' entered.
#' @param ar.order (type: numeric) AR polimonial order.
#' @param ma.order (type: numeric) MA polimonial order.
#' @param sd.order (type: numeric) polinomial order noise scale factor.
#' @param d.order (type: numeric) \code{d} polinomial order, where \code{d} is
#' the \code{ARFIMA} parameter.
#' @param include.d (type: numeric) logical argument for \code{ARFIMA} models.
#' If \code{include.d=FALSE} then the model is an ARMA process.
#' @param N (type: numeric) value corresponding to the length of the window to
#' compute periodogram. If \code{N=NULL} then the function will use
#' \eqn{N = \textrm{trunc}(n^{0.8})}, see Dahlhaus (1998) where \eqn{n} is the
#' length of the \code{y} vector.
#' @param S (type: numeric) value corresponding to the lag with which will go
#' taking the blocks or windows.
#' @param include.taper (type: logical) logical argument that by default is
#' @references
#' \insertRef{brockwell2002introduction}{LSTS}
#' \insertRef{palma2010efficient}{LSTS}
#' @importFrom stats na.omit
#' @export
LS.whittle.loglik <- function(x, series, order = c(p = 0, q = 0),
ar.order = NULL, ma.order = NULL, sd.order = NULL,
d.order = NULL, include.d = FALSE, N = NULL,
S = NULL, include.taper = TRUE) {
y <- series
T. <- length(y)

if (is.null(N)) {
N <- trunc(T.^0.8)
}
if (is.null(S)) {
S <- trunc(0.2 * N)
}

M <- trunc((T. - N) / S + 1)

if (is.null(ar.order)) {
ar.order <- rep(0, order[1])
}
if (is.null(ma.order)) {
ma.order <- rep(0, order[2])
}
if (is.null(sd.order)) {
sd.order <- 0
}
if (is.null(d.order)) {
d.order <- 0
}

p <- na.omit(c(ar.order, ma.order, sd.order))
if (include.d == TRUE) {
p <- na.omit(c(ar.order, ma.order, d.order, sd.order))
}

if (length(x) != sum(p + 1)) {
stop("error in the number of parameters")
} else {
lik <- 0
for (j in 1:M) {
u <- (N / 2 + S * (j - 1)) / T.
aux <- periodogram(y[(1 + S * (j - 1)):(N + S * (j - 1))], include.taper = TRUE, plot = FALSE)
I <- aux$periodogram X <- numeric() k <- 1 for (i in 1:length(p)) { X[i] <- sum(x[k:(k + p[i])] * u^(0:p[i])) k <- k + p[i] + 1 } phi <- numeric() k <- 1 if (order[1] > 0) { phi[is.na(ar.order) == 1] <- 0 phi[is.na(ar.order) == 0] <- X[k:(length(na.omit(ar.order)))] k <- length(na.omit(ar.order)) + 1 } theta <- numeric() if (order[2] > 0) { theta[is.na(ma.order) == 1] <- 0 theta[is.na(ma.order) == 0] <- X[k:(length(na.omit(ma.order)) + k - 1)] k <- length(na.omit(ma.order)) + k } d <- 0 if (include.d == TRUE) { d <- X[k] k <- k + 1 } sigma <- X[k] f <- spectral.density(ar = phi, ma = theta, d = d, sd = sigma, lambda = aux$lambda)

lik <- sum(log(f) + I / f) / N + lik
}

lik <- lik / M

if (is.na(sigma) | sigma <= 0) {
lik <- Inf
}

lik
}
}


## Try the LSTS package in your browser

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

LSTS documentation built on July 29, 2021, 5:07 p.m.