R/plot_spec.R

In boiwsa: Seasonal Adjustment of Weekly Data

#' Compare AR spectra of original and seasonally adjusted series
#'
#' Computes and plots autoregressive (AR) spectral density estimates for the
#' detrended original series and its seasonally adjusted counterpart. Spectra are
#' estimated using \code{\link[stats]{spec.ar}} with AR order set to 60. The plot
#' highlights frequencies corresponding to intramonthly and intrayearly cycles.
#'
#' @importFrom graphics abline legend lines
#' @importFrom stats spec.ar
#' @import ggplot2
#'
#' @param x A \code{boiwsa()} result object. Must contain components \code{x},
#'   \code{sa}, and \code{trend}.
#'
#' @return A \code{ggplot2} object showing the AR spectral density estimates for
#' the detrended original and seasonally adjusted series.
#' @export
#'
#' @examples
#' \donttest{
#' # Not run
#' # Seasonal adjustment of weekly US gasoline production
#' res <- boiwsa(x=gasoline.data$y,dates=gasoline.data$date)
#' plot_spec(res)}
#'
plot_spec=function(x){


  spec0=spec.ar((x$x-x$trend),order = 60,plot = F,n.freq=500)
  spec1=spec.ar((x$sa-x$trend),order = 60,plot = F)

  ggplot2::ggplot()+
    ggplot2::geom_line(ggplot2::aes(x=spec0$freq,y=spec0$spec,color="orig"))+
    ggplot2::geom_line(ggplot2::aes(x=spec0$freq,y=spec1$spec,color="sa"))+
    ggplot2::geom_vline(xintercept = 1:2/4.34,linetype="dashed")+
    ggplot2::geom_text(ggplot2::aes(x=1:2/4.34, label="\n Intra-monthly cycle peaks", y=0.5*max(spec0$spec)), colour="black", angle=90)+
    ggplot2::geom_vline(xintercept = (1:3)/52.1775,linetype="dashed")+
    ggplot2::geom_text(ggplot2::aes(x=3/52.1775, label="\n First three intra-yearly cycle peaks", y=0.5*max(spec0$spec)), colour="black", angle=90)+
    ggplot2::scale_color_manual(name = "",
                                values = c( "orig"='#31a354', "sa"="#3182bd"),
                                labels = c("Original", "Seasonally adjusted"))+
    ggplot2::theme_bw()+
    ggplot2::theme(legend.position="bottom")+
    ggplot2::theme(legend.text = ggplot2::element_text(size = 11))+
    ggplot2::ylab(" ")+
    ggplot2::xlab("Frequency")


}