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

#### Documented in periodogram

#' @title Periodogram function
#' @description This function computes the periodogram from a stationary time
#' serie. Returns the periodogram, its graph and the Fourier frequency.
#' @details
#' The tapered periodogram it is given by
#' \deqn{I(\lambda) = \frac{|D_n(\lambda)|^2}{2\pi
#' H_{2,n}(0)}} with \eqn{D(\lambda) = \sum_{s=0}^{n-1} h
#' \left(\frac{s}{N}\right) y_{s+1}\,
#' e^{-i\,\lambda\,s}}, \eqn{H_{k,n} = \sum_{s=0}^{n-1}h
#' \left(\frac{s}{N}\right)^k\,
#' e^{-i\,\lambda\,s}} and \eqn{\lambda} are Fourier frequencies defined as
#' \eqn{2\pi k/n}, with \eqn{k = 1,\,\ldots,\, n}.
#' The data taper used is the cosine bell function,
#' \eqn{h(x) = \frac{1}{2}[1-\cos(2\pi x)]}. If the series has missing data,
#' these are replaced by the average of the data and \eqn{n} it is corrected by
#' $n-N$, where \eqn{N} is the amount of missing values of serie. The plot of
#' the periodogram is \code{periodogram} values vs. \eqn{\lambda}.
#' @param y (type: numeric) data vector
#' @param plot (type: logical) logical argument which allows to plot the
#' periodogram. Defaults to TRUE.
#' @param include.taper (type: logical) logical argument which by default is
#' \code{FALSE}. If \code{include.taper=TRUE} then \code{y} is multiplied by
#' \eqn{0.5(1 - \cos(2\pi(n-1)/n))} (cosine bell).
#' @references
#' \insertRef{brockwell2002introduction}{LSTS}
#' \insertRef{dahlhaus1997fitting}{LSTS}
#' @examples
#' # AR(1) simulated
#' set.seed(1776)
#' ts.sim <- arima.sim(n = 1000, model = list(order = c(1, 0, 0), ar = 0.7))
#' per <- periodogram(ts.sim)
#' per\$plot
#' @return
#' A list with with the periodogram and the lambda values.
#' @importFrom stats fft
#' @importFrom ggplot2 geom_line
#' @importFrom scales math_format
#' @export
periodogram <- function(y, plot = TRUE, include.taper = FALSE) {
series <- y - mean(y, na.rm = TRUE)
N <- sum(is.na(series))
series[is.na(series)] <- 0
n <- length(series)
if (include.taper == TRUE) {
a <- 0.5 * (1 - cos(2 * pi * (0:(n - 1)) / n))
series <- series * a
}
aux <- Mod(fft(series))^2
m <- n / 2
periodogram <- (aux[2:(m + 1)]) / (2 * pi * (n - N))
if (include.taper == TRUE) {
periodogram <- (aux[2:(m + 1)]) / (3 * pi * (n - N) / 4)
}
lambda <- (2 * pi * (1:m)) / n
if (plot == TRUE) {
pi_scales <- math_format(.x * pi, format = function(x) x / pi)
g <- ggplot(data.frame(x = lambda, y = periodogram)) +
geom_line(aes(x = as.numeric(x), y = as.numeric(y)), color = "#1f77b4") +
labs(x = "Frequency", y = "Periodogram") +
scale_x_continuous(labels = pi_scales, breaks = seq(0, pi, pi / 4)) +
theme_minimal()
} else {
g <- NULL
}
return(list(periodogram = periodogram, lambda = lambda, plot = g))
}


## 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.