Quantile-Based Spectral Analysis in an Object-Oriented Framework and a Reference Implementation in R: The quantspec Package

Quantile-based approaches to the spectral analysis of time series have recently attracted a lot of attention. Despite a growing literature that contains various estimation proposals, no systematic methods for computing the new estimators are available to date. This paper contains two main contributions. First, an extensible framework for quantile-based spectral analysis of time series is developed and documented using object-oriented models. A comprehensive, open source, reference implementation of this framework, the R package quantspec, was recently contributed to CRAN by the author of this paper. The second contribution of the present paper is to provide a detailed tutorial, with worked examples, to this R package. A reader who is already familiar with quantile-based spectral analysis and whose primary interest is not the design of the quantspec package, but how to use it, can read the tutorial and worked examples (Sections 3 and 4) independently.

1.A short introduction to quantile-based spectral analysis 1.1.Laplace and copula cumulants and their spectral representation Quantification of serial dependence in a stationary process (X t ) t∈Z is traditionally based on its autocovariance and autocorrelation functions, which measure linear dependencies among observations at different times.Periodicities of a time series are then most commonly analyzed by decomposing the autocovariance function, into a sum of sines and cosines.This approach is referred to as (ordinary) spectral analysis of time series and has been known for decades.As a statistical method, it has been investigated many times and is well understood.In the analysis of centered Gaussian time series this approach is particularly attractive, because the autocovariance function completely characterizes the distribution of the underlying process.If that process is not Gaussian, ordinary spectral analysis suffers from typical weaknesses of L 2 -methods: it is lacking robustness against outliers and heavy tails, and is unable to capture important dynamic features such as changes in the conditional shape (skewness, kurtosis), time-irreversibility, or dependence in the extremes.In addition to this, only time series with an existing second moment can be analyzed at all.All of this was previously realized by many researchers, and various extensions and modifications of the L 2 -periodogram have been quantspec: Quantile-based Spectral Analysis in R proposed to remedy those drawbacks.
Approaches to robustifying the traditional spectral methods against outliers and deviations from the distributional assumptions were taken, among others, by Kleiner, Martin, and Thomson (1979); Klüppelberg and Mikosch (1994); Mikosch (1998); Katkovnik (1998); Hill and Mc-Closkey (2013).To account for more general dynamic features alternative spectral concepts and tools were recently proposed.A first step in that direction was taken by Hong (1999).In order to obtain a complete description of the two-dimensional distributions at lag k, he introduced a generalized spectrum where covariances Cov(X t , X t−k ) are replaced by covariances Cov(e ix 1 Xt , e ix 2 X t−k ) yielding a spectrum closely related to the joint characteristic functions of the pairs (X t , X t−k ).In the quantile-based approach to spectral analysis the objects of interest are the Laplace cross-covariance kernel γ k (q 1 , q 2 ) := Cov I{X t ≤ q 1 }, I{X t−k ≤ q 2 } , q 1 , q 2 ∈ R, k ∈ Z, and the copula cross-covariance kernel where I{A} denotes the indicator function of the event {A} and F the marginal distribution function (the distribution function of any X t , due to the assumed stationarity).Obviously these measures exist without the necessity to make assumptions about moments.Also, when the underlying process is not Gaussian, and the quantile-based measures of serial dependence are considered as functions with arguments q 1 , q 2 , or τ 1 , τ 2 respectively, they provide a much richer picture about the pairwise dependence than would the autocovariances.As in the approach of Hong (1999), a complete description of the joint distributions (or copulas) of the pairs (X t , X t−k ) is available.A particular advantage of the copula cross-covariance kernel is its invariance to monotone transformations.This allows to disentangle the serial features from the marginal features.For a full list of the properties and advantages of those dependence measures the interested reader be refered to Hong (2000); Li (2008Li ( , 2012Li ( , 2013Li ( , 2014)); Hagemann (2011); Lee and Rao (2012); Dette, Hallin, Kley, and Volgushev (2014+); Kley, Volgushev, Dette, and Hallin (2014) and Kley (2014b).
Under sumability conditions on (γ k ) and (γ U k ) the representations of (γ k ) and (γ U k ) in the "frequency domain" take the form of the Laplace spectral density kernel and the copula spectral density kernel where q τ := F −1 (τ ).By the relation e ikω f q 1 ,q 2 (ω)dω, and a similar representation for γ U k (τ 1 , τ 2 ), the representations in the "frequency domain" are seen to be equivalent to the "time domain" quantities.
Sometimes considering the cumulated Laplace or copula spectral density kernels, which can be defined as and is more convenient.
Quantities as γ k and γ U k , and their spectral representations (1)-( 4) naturally come into the picture when the clipped processes (I{X t ≤ q}) t∈Z and (I{F (X t ) ≤ τ }) t∈Z are investigated.Such binary processes have been considered earlier in the literature by e. g.Kedem (1980).Observe that the quantile-based spectral quantities can be interpreted in terms of an orthogonal increment process of a spectral representation of the clipped process which exists for every stictly stationary process; no assumptions about moments are necessary.
Despite the vast amount of theoretical work, a public software solution was so far not available.

Estimators for the quantile-based spectral analysis of time series
In this section of the introduction, various estimators (the so-called quantile periodograms) for the Laplace and copula spectra defined in Section 1.1 are briefly considered.For the upcoming definitions denote by X 0 , . . ., X n−1 an observed time series of the process (X t ) t∈Z , by the empirical distribution function of X 0 , . . ., X n−1 , and by z and z the real and imaginary part of z = z + i z ∈ C, respectively.Further, denotes the so-called check function [cf.Koenker (2005)].
Definition 1 (Quantile-regression based periodograms) , where, for ω = 0 mod π, and τ ∈ (0, 1), If ω = 0 mod π, when the regressor that yields the imaginary part of the estimate vanishes, the estimates need to be adapted as follows: for ω π = 2π(j + 1/2), j ∈ Z, let For ω ∈ 2πZ an adaptation is also possible [cf.Kley (2014b)], but since it is not required for the definition of the smoothed estimates it is omitted for the sake of brevity.Observe that the rank-based periodograms obtained their name due to the fact that n Fn (X t ) is the rank of X t among X 0 , . . ., X n−1 The Laplace periodograms trace back to Katkovnik (1998) who, in the field of signal processing, suggested L p estimators in a harmonic linear model.Li (2008) proved asymptotic normality of the Laplace periodograms for τ 1 = τ 2 = 0.5 and later extended the approach to arbitrary quantiles with 0 < τ 1 = τ 2 < 1 (Li 2012).Dette et al. (2014+) introduced the estimator with distinct quantile levels τ 1 and τ 2 (not necessarily equal), and also considered the rank-based version.
Another estimator is based on the discrete Fourier transform of clipped processes and can be defined as follows: Definition 2 (Periodograms based on clipped time series) For ω ∈ R and q 1 , q 2 ∈ R, the clipped time-series periodogram is defined as For ω ∈ R and τ 1 , τ 2 ∈ [0, 1] the copula rank periodogram (for short CR periodogram) is defined as Note the similarity between all the quantile periodograms and the cross-periodograms in multivariate time series analysis [cf., e. g., Brillinger (1975), p. 235]: each periodogram is a product of two frequency representation objects computed at frequencies (ω and −ω) that sum to zero.
The estimator based on the discrete Fourier transformation of clipped time series was introduced by Hong (2000), who used it for a test of pairwise independence.Hagemann (2011) analyzed the special case of τ 1 = τ 2 in the presence of serial dependence.The case of distinct quantile levels τ 1 and τ 2 (not necessarily equal) was discussed in Kley et al. (2014), where weak convergence to a Gaussian process was established.Lee and Rao (2012) investigated the distributions of Cramér-von Mises type statistics, based on empirical joint distributions.
As in the traditional case the new periodograms are not consistent estimators [cf. the positive variances of the limit distributions in Theorems 3.2 and 3.4 in Dette et al. (2014+) or Proposition 3.4 in Kley et al. (2014)].
A comprehensive description of all estimators and their asymptotic properties is available in Kley (2014b).

An analysis of functional requirements
The quantspec software project was triggered by the development of the quantile-based methods for spectral analysis [cf.Section 1, Dette et al. (2014+) and Kley et al. (2014)].The primary aim has been to make these new methods accessible to a wide range of users.
Before going into the programming-specific details of the project, a conceptual design, nonspecific to any specific programming language was developed.By this procedure, additional insight and a thorough documentation of the computational characteristics of quantile-based spectral analysis could be gained.The conceptual design serves as a blueprint for implementations in (possibly) various environments and can easily be transformed into an implementation plan including the details specific to the respective programming environment.
Aiming for a software system that is most flexible in the ways in which it can be used, that can easily be extended in functionality and also for the ease of its maintenance an objectoriented design was chosen.This type of design also contributes to a structure of the system that can be better understood, both by users and developers.The general structure of the system for performing quantile-based spectral analysis is described using class diagrams of the unified modeling language (UML).In these diagrams, all necessary components and their interrelations are described in a formal manner.
To understand the specification of the framework, recall that in an object-oriented design the components of the system are objects encapsulating both data and behavior of a specific "real-world" entity.The structure of each object can thus be described by a meaningful name (the class name), a collection of data fields (in R these are called slots) and implementations of the behavior (in R these implementations are called methods).In a class diagram each class (i.e., the composite of class name, data files and implemented behavior) is represented as a rectangle subdivided into three blocks.The name of the class is given in the top block, the data fields in the middle block and the methods in the bottom block.Note that in the unified modeling language the data fields and methods are specified in a standardized format.In this format the first symbol is an abbreviation used to specify the visibility of the class member.Here "+" for public and "-" for private members are used, meaning that the member is intended to be seen (and used) from outside the object or from inside the object only, respectively.For a data field the name is then followed by a colon and the type of the field.For a method the parameters are given in parenthesis; optional parameters are denoted by two dots.In the class diagram, relationships between classes are marked as lines connecting them.Currently two different types of relationships are modeled.A line with a triangular shaped tip at one end is used to declare a generalization relationship (sometimes also coined inheritance or "is a" relationship).The class at the end of the line with the triangle is called the superclass or the parent, while the class on the other side is called the subclass or the child.In particular this type of relationship implies that an object that is an instance to the subclass and therefore provides all the subclasses' data fields and methods will also provide the data fields and methods of the superclass (and the superclasses' superclass it there are such, etc.).The second type of relationships used in this framework is that of an aggregation (sometimes called "has a" relationship).A line with a hollow diamond at one end is used to denote this kind of relationship, where objects to the class at the end of the line with the diamond are the ones having objects of the class at the other end of the line as a part of them.
At each end of the line the so-called cardinalities are denoted in the form of two numbers with dots in between them.The left number is the minimum number of objects of that type in the relationship that need to exist, the right number is the maximum number.A star is used to denote an unknown positive integer.If min and max cardinality coincide they are displayed as one number without the dots.
For the class diagrams in this manuscript the classes are arranged in a way such that (whenever possible) generalization relationships are displayed with the superclass on top and the subclasses in the bottom.Aggregations are shown with aggregated classes to the left and/or the right of the class representing "the whole".
A graphical representation of the framework for quantile-based spectral analysis is not given in one holistic diagram, but in two class diagrams that are on display in Figures 1 and 2. The structure of the framework is presented in two, thematically organized class diagrams rather than in one, because the 13 classes and their relations could not be fitted easily onto one page without breaking the above mentioned layout guidelines.On the other hand it was easy to group the classes by two topics.
In the next sections, all classes of the framework and their relations are going to be thoroughly described and motivated.

The base class QSpecQuantity and its successors
Many of the quantities important for the quantile-based spectral analysis of a stationary time series [i.e., the estimators of Definitions 1-4 and the model quantities ( 1)-( 4)] are of the functional form, where To provide a common interface to these objects the abstract class QSpecQuantity was introduced.Its data fields (i.e., the array values, a vector of reals frequencies and a list with two vectors of reals levels), are designed to store the sets and the family Note that the handling of a family of B quantile spectral quantities is necessary when bootstrapping replicates are present.The special case of only one function Q can easily be handled by setting B = 1.
There are four classes inheriting the data structure and the method show1 from the abstract class QSpecQuantity.Two such classes, QuantilePG and SmoothedPG, implement the computation of the various quantile periodograms and smoothed quantile periodograms, respectively.A more detailed description has to include the other related classes and is therefore deferred to a separate section (i.e., Section 2.3).A graphical representation of the relevant parts of the framework can be seen in Figure 1.The other two of the classes generalizing the abstract class QSpecQuantitiy are referred to by the names QuantileSD and IntegrQuantileSD.These two classes implement the model quantities ( 1)-( 4).The graphical representation can be seen in Figure 2. A detailed description is deferred to Section 2.4.

Implementation of the quantile-based (smoothed) periodograms
The components relevant to the implementation of the quantile-based spectral statistics are presented in Figure 1.As alluded to in the previous section the two classes QuantilePG and SmoothedPG will do the job, in conjunction with the superclass QSpecQuantity from which they inherit the data structure to store the computed values.
To better understand the implementation surrounding QuantilePG, observe that the quantilebased periodograms defined in Definitions 1 and 2 share the common structure of an outer  product (scaled with (2πn) −1 ).To compute either one of the four periodograms it suffices to perform the same operation to one of the frequency representation objects bτ respectively.In the framework this fact is incorporated by introduction of the abstract class FreqRep and its two subclasses ClippedFT and QRegEstimator, where the actual computations are implemented via the method initialization().The class FreqRep serves as a common interface to the quantities in ( 6).It provides data fields to store various information, including • the observations Y from which the quantities were computed, • the frequencies and levels for which the computation was performed, • the result of the computation, which is stored in an array values.Further more, a flag isRankBased indicates whether, prior to the main computations, the observations (X t ) were transformed to pseudo-observations ( Fn (X t )).Performing this extra step will yield bτ n,R (ω) instead of bτ n (ω) or d τ n,R (ω) instead of d τ n (ω), respectively.The class BootPos allows for performing a block bootstrap procedure by "shuffling" the observations and repeatedly doing the computations on these bootstrapped observations.Currently only one method, the MovingBlocks bootstrap, is implemented.Now, turning attention to the class SmoothedPG, recall that the various smoothed periodograms are all defined similarly, in the sense that computing the smoothed periodogram for ω j := 2πj/n, j = 1, . . ., n − 1 basically means to do a discrete convolution of the sequence of quantile periodograms computed at ω s with a sequence of appropriately chosen weight functions W n (ω s ).Hence, everything needed for the smoothed periodogram is these two ingredients, which is reflected in the framework by two aggregation relationships involving the class SmoothedPG.The first such relationship links SmoothedPG to the QuantilePGs to be smoothed.The second such relationship links SmoothedPG to a class Weight, which provides a common interface to different weight functions.Currently two implementations are included.Employing weights of type KernelWeight, defined by a kernel W and a scale parameter bw (bandwidth), will yield an estimator for the quantile (i.e, Laplace or copula) spectral density.An alternative is to use weights of type SpecDistrWeight, which yields estimators for the integrated quantile (i.e, Laplace or copula) spectral density.

Implementation of the quantile-based spectral measures
The classes QuantileSD and IntegrQuantileSD were introduced to the framework to make the quantities (1)-( 4) available to the user.
To obtain access to a quantity of the form (1) or ( 2), an instance of QuantileSD can be created.In this case, a number of R independent copies of a time series of length N are obtained by calling the function ts.The function ts is a parameter to specify the model for which to obtain the model quantity and it handles the simulation process.Then, for each of these time series a QuantilePG object is created and their values are averaged: first across the R independent copies, saving the result to meanPG and an estimation of the standard error to stdError.After that the averages are averaged again for each combination of levels across frequencies by smoothing.The result is then made available via the data field values of the superclass.By a call to the method increasePrecision the number of independent copies can be increased at any time to yield a less fluctuating average.
To obtain access to quantities of the form (3) or ( 4), an instance of IntegrQuantileSD can be created.For the computation an object of type QuantileSD is created and subsequently the integral is approximated via a Riemann sum.
3. Reference implementation: The R package quantspec

Overview
The quantspec package (Kley 2014a) is intended to be used both by theoretically oriented statisticians and also by data analysts, who work on a more applied basis.In order to address this broad group of potential users the R system for statistical computing (R Core Team 2012) was chosen as a platform.The R system is particularly well suited for the realization of this project, because it is accessible from many operating systems, without charge, already available to the targeted audience and, in particular, allows to integrate the package's functionality among many other, well developed packages.An important example is that the function rq of the quantreg package (Koenker 2013) could be used.
Both R and the quantspec package are available from the comprehensive R archive network (CRAN, http://cran.r-project.org).The package's development is actively continued with the source code available from a GitHub repository (https://github.com/tobiaskley/quantspec).Besides the source code of the releases, which are also available from the CRAN servers, the GitHub repository additionally contains a detailed history of all changes, including comments, that were applied to the source code since April/10/2014, when the GitHub repository was created.The repository is organized into several branches: the master branch, the develop branch and possibly several topic branches.Since version 1.0-0 the master branch contains the source code of all release candidates and official releases, the develop branch contains the most recent updates and bug fixes that were not yet released.The topic branches contain the source code in which extensions to the package are developed.
To install a package from the source code straight from the repository use the install_github function of the devtools package (Wickham and Chang 2013).More precisely, install the devtools package, if it's not already installed, and call R> devtools::install_github("tobiaskley/quantspec", ref = "master") Instead of using "master", another branch (e. g., develop to pull the most recent updates) or a tag that is pointing to one of the releases (e. g., v1.0-0-rc1, to install the 1st release candidate to version 1.0-0) can be used as ref.
Note that the code in the develop branch is merged into the master branch only for a release (candidate) and after being thoroughly tested, so if you're installing from the develop branch you will potentially be using code that has not been fully tested.
Use the option build_vignettes = FALSE if you don't have L A T E X installed on your system.

R code intended for the user and its documentation
The classes of the quantspec package, their methods, slots, dependencies and inheritance properties are implemented as conceptually designed [cf.Section 2].Recall that the design was presented in form of class diagrams, on display in Figures 1 and 2, and that no specific programming language was assumed.All classes that are intended for the end-user possess a constructor method with the same name as the class itself but beginning with a lower case letter.The classes intended for the end-user and their constructors are listed in Table 1.
For a more detailed description of constructors and classes, documentation within the online help system of R is available.After loading the package, which is done by calling > library("quantspec") the help file of the package, which provides an overview on the design, can be called by executing > help("quantspec") on the R command line.Note that an index of all available functions can be accessed at the very bottom of the page.If for example more information on the constructor of QRegEstimator and on the class itself is desired, then > help("qRegEstimator") > help("QRegEstimator") should be invoked.Using this class to determine the frequency representation b τ n,R (ω), for τ ∈ {0.25, 0.5, 0.75} would look as follows.In a toy example, where eight independent random variables X 0 , . . ., X 7 ∼ N (0, 1) are generated and used to compute b τ n,R (ω), call > Y <-rnorm(8) > bn <-qRegEstimator(Y, levels = c(0.25,0.5, 0.75)) By default the computation is done for all Fourier frequencies ω = 2πj/n ∈ [0, π], n = 8, j = 0, . . ., n/2 .The computed information can then be viewed by typing the name of the variable (i.e., bn) to the R console:

> getParallel(bn) [1] FALSE
To invoke a method f with parameters p1, ..., pk of an object obj the call is f(obj, p1, ..., pk).An example is to invoke the accessor function getValues, which is equipped with parameters to get the values associated with certain frequencies or levels.An exemplary call looks like this: > getValues(bn, levels = c(0.25,0.5)) , , 1 Note that the result is returned as an array of dimension c(J, K, B + 1), where in the present case B = 0 bootstrap replications were performed.For a detailed description on how to use the function getValues in the above mentioned case, access the online help via the command > help("getValues-FreqRep") Note the format method_name-class_name to access the help page of a method and that the attribute values is part of the abstract class FreqRep [cf. Figure 2].
A graphical representation of the data can easily be created by applying the plot command.For example, to compute and plot the frequency representations d τ 32,R (ω), from 32 simulated, standard normally distributed random variables execute the following lines on the R shell: > dn <-clippedFT(rnorm(32), levels = seq(0.05,0.95, 0.05)) > plot(dn, frequencies = 2 * pi * (0:64) / 32, levels = c(0.25,0.5)) The above script will yield the diagrams that are on display in Figure 3.Note that the d τ 32,R (ω) were determined for τ ∈ {0.05, 0.1, . . ., 0.9, 0.95} and, by the default setting, for all Fourier frequencies from [0, π].The plot, however, was parameterized to show only τ ∈ {0.25, 0.5}, but all Fourier frequencies from [0, 4π]; by default all available levels and frequencies would be used.In this example two of the 19 frequencies were selected to yield a plot of a size that is apropriate to fit onto the page.Further more, the plot was parameterized to show d n,R (ω) for all Fourier frequencies from [0, 4π] to illustrate characteristic redundancies in the frequency representation objects, and to point out that the default values are always sufficient.The two relations hold for any ω ∈ R and j ∈ Z, d τ n,R (ω).Therefore, without additional calculations, the plot of d τ n,R (ω) can be determined for any ω ∈ 2πj/n, j ∈ Z, as long as d τ n,R (ω) is known for ω ∈ 2πj/n, j = 0, . . ., n/2 , which is what is determined by the default setting.Note that all of this happens transparently for the user, as the method getValues takes care of it.Another fact that one can presume by inspecting Figure 3 is that d τ n,R (ω) appears to be uncorrelated and centered (for ω = 0 mod 2π).

Additional elements of the package
The quantspec package includes three demos that can be accessed via > demo("sp500") > demo("wheatprices") > demo ("qar-simulation") Several examples explaining how to use the various functions of the package can be found in the online help files or the folder examples in the directory where the package is installed.The package comes with two data sets sp500 and wheatprices that are used in the demos and in the examples.A package vignette amends the online help files.It contains the text of this paper.Unit tests covering all main functions were implemented using the testthat framework (Wickham 2011).

Analysis of the S&P 500 stock index, 2007-2010
In this section the use of the quantspec package from the perspective of a data analysts is explained.To this end an analysis of the returns of the S&P 500 stock index is performed.Note that a similar analysis and the data set used are available in the package.Calling demo("sp500") will start the computations and by sp500 the data set can be referenced to do additional analysis.
For the example the years 2007 through to 2010 were selected to have a time series that, at least to some degree, can be considered stationary.Aside from this more technical consideration, employing the new statistical toolbox will reveal interesting features in the returns collected in the financial crisis that completely escape the analysis with the traditional tools blindly applied.
For a start, use the following R script to plot the data, the autocovariances of the returns and the autocovariances of the squared returns.
> library("zoo") > plot(sp500, xlab = "time t", ylab = "", main = "") > acf(coredata(sp500), xlab = "lag k", ylab = "", main = "") > acf(coredata(sp500)^2, xlab = "lag k", ylab = "", main = "") The three plots are displayed in Figure 4. Inspecting them, it is important to observe that the returns themselves appear to be almost uncorrelated.Therefore, not much insight into the serial dependency structure of the data can be expected from traditional spectral analysis.The squared returns on the other hand are significantly correlated.This observation, typically taken as an argument to fit an ARCH or GARCH model, clearly proves that serial dependency exists.In what follows the copula spectral density will be estimated from the data, using quantile periodograms and smoothing them.It will be seen that using the quantspec package this can be done with only a few lines of code necessary.
First, take a look at the CR periodogram I τ 1 ,τ 2 n,R (ω).In the quantspec package it is represented as a QuantilePG object and can be computed calling the constructor function quantilePG with the parameter type = "clipped".To do the calculation for τ 1 , τ 2 ∈ {0.05, 0.5, 0.95}, all Fourier frequencies ω and with 250 bootstrap replications determined from a moving blocks bootstrap with block length = 32, it suffices to execute the first command of the following script: > CR <-quantilePG(sp500, levels.1 = c(0.05,0.5, 0.95), + type = "clipped", type.boot= "mbb", B = 250, l = 32) > freq <-getFrequencies(CR) > plot(CR, levels = c(0.05,0.5, 0.95), Using the second command it is possible to learn for which frequencies the values of the CR periodogram are available.As pointed out in the previous paragraph it was computed for all Fourier frequencies from the interval [0, 2π), which is the default setting for quantilePG and smoothedPG.With the third command the graphical representation of the CR periodogram, which can be seen in Figure 5, is plotted.The plot seen here is a typical plot of any QSpecQuantity: in a configuration with K levels the plot has the form of a K × K matrix, where the subplots on and below the diagonal display the real part of the CR periodogram I τ 1 ,τ 2 n,R (•), with the levels τ 1 and τ 2 denoted on the left and bottom margins of the plot.Above the diagonal the imaginary parts are shown.
To observe the larger values in the neighborhood of ω = 0 and in the extreme quantile levels more closely a plot showing the CR periodogram only for frequencies ω ∈ [0, π/5] can be generated using the following script: > plot(CR, levels = c(0.05,0.5, 0.95), The plot is shown in Figure 6.
In the next step the computed quantile periodogram CR can be used as the basis to determine a smoothed CR periodogram sCR.In the form of a SmoothedPG object it can be generated by the constructor smoothedPG of that class.Besides the QuantilePG object CR, a KernelWeight object is required, which is easily generated using the constructor kernelWeight.As parameters the constructor kernelWeight requires a kernel W and a bandwidth bw.The quantspec package comes with several kernels already implemented.The Epanechnikov kernel for example can be refered to by the name W1.For a complete list of the available kernels call

> help("kernels")
To compute the smoothed CR periodogram from CR using the Epanechnikov kernel and bandwidth bw = 0.07 the first of the following two commands need to be executed.
> sPG <-smoothedPG(CR, weight = kernelWeight(W = W1, bw = 0.07)) > plot(sPG, levels = c(0.05,0.5, 0.95), type.scaling= "individual", Of course, the second line initiates plotting the smoothed CR periodogram, which is on display in Figure 7.Note that the option type.scalingcan be set to yield a plot with certain subplots possessing the same scale.In Figure 7 pointwise confidence intervals are shown.By default these are determined using a normal approximation to the distribution of the estimator as is suggested by the limit theorem in Kley et al. (2014).An alternative is to use the quantiles of estimates computed from the block bootstrap replicates.These pointwise confidence intervals can be plotted using the option type.CIs = "boot.full",as is shown in the following script: > plot(sPG, levels = c(0.05,0.5, 0.95), type.scaling= "real-imaginary", + ptw.CIs = 0.1, type.CIs = "boot.full", For illustrative purposes a different type of scaling was used for the second plot.A complete description of the options is available in the online help, which can be accessed by calling > help("plot-SmoothedPG") Inspecting the plots in Figures 7 and 8 reveals that serial dependency in the events {Y t ≤ q 0.05 } and {Y t ≤ q 0.95 } is present in the data.This concludes the introduction of the quantspec package for data analysts and we can continue with the presentation of how it can also make the work of a probability theorist easier.

A simulation study: Analysing a quantile autoregressive process
In this section using the quantspec package from the perspective of a probability theorist is explained.The aim is twofold.On the one hand, further insight into a stochastic process shall be gained.Any process for which a function to simulate finite stretches of is available can be studied.On the other hand, the finite sample performance of the new spectral methods are to be evaluated.Note that the example discussed in this section and the functions to simulate QAR(1) processes are available inside the package, by calling demo("qar-simulation") and by referring to the function ts1, which implements the QAR(1) process that was discussed in Dette et al. (2014+) and Kley et al. (2014).Recall that a QAR(1) process is a sequence (X t ) of random variables that fulfills where U t is independent white noise with U t ∼ U[0, 1], and θ 1 , θ 0 : [0, 1] → R are model parameters (Koenker and Xiao 2006).The function ts1 implements the model, where θ 1 (u) = 1.9(u − 0.5), u ∈ [0, 1] and θ 0 = Φ −1 , which was discussed in Dette et al. (2014+) and Kley et al. (2014).A complete list of models included in the package can be seen in the online documentation of the package by calling > help("ts-models") The following, two-line script can be used to generate the graphical representation of the copula spectral density that is on display in Figure 9 > csd <-quantileSD(N = 2^9, seed.init= 2581, type = "copula", + ts = ts1, levels.1 = c(0.25,0.5, 0.75), R = 100, quiet = TRUE) > plot (csd, ylab = expression(f[list(q[tau[1]], q[tau[2]])](omega))) When analysing a time series model the recommended is to compute the quantile spectral density once with high precision, store it to the hard drive, and load it later whenever it is needed.The following two lines of code can be used to do this: > csd <-quantileSD(N=2^12, seed.init= 2581, type = "copula", + ts = ts1, levels.1 = c(0.25,0.5, 0.75), R = 50000) > save(csd, file="csd-qar1.rdata") With the first configuration (N = 2 9 and R = 100) the computation time was only around three seconds.To compute the second csd object (with N = 2 12 and R = 50000) the same machine needed roughly 2.5 hours.Storing the object in a file takes about 1MB of hard disk space.Not only values and stdErrors are stored within the QuantileSD object; also the final state of the pseudo random number generator is stored and the method increasePrecision can be used to add more simulation runs at any time to yield a better approximation to the true quantile spectrum.More information on this method can be found in the online help, which is accessible via > help("increasePrecision-QuantileSD") Once the computation is finished the diagram on display in Figure 10 can be created using the following two lines of code: The parameter frequencies was used when plotting the copula spectral density to create a plot that can be compared to the one in Figure 9.Note that by default the plot would have been created using all available frequencies which yields a grid of 8 times as many points on the x-axis (N = 2 12 vs.N = 2 9 ).Now, to get a first idea of how well the estimator performs, plot the smoothed CR periodogram computed from one simulated QAR(1) time series of length 512: The generated plot is on display in Figure 11.It is worth pointing out that in this example (N = 512) the estimator performs already quite well.Note that a different version of the constructor smoothedPG was used here than in Section 4.1.When computing a smoothed quantile periodogram straight from a time series, the syntax is the same as for quantilePG, but with the additional parameter weight.
Finally, for the simulation study, a number of R = 5000 independent QAR(1) time series are generated.Before the actual simulations, some variables that determine what is to be simulated are defined: > set.seed(2581)> ts <-ts1 > N <-128 > R <-5000 > freq <-2 * pi * (1:16) / 32 > levels <-c(0.25,0.5, 0.75) > J <-length(freq) > K <-length(levels) > sims <-array(0, dim=c(4, R, J, K, K)) > weight <-kernelWeight(W = W1, bw = 0.3) Setting the seed in the very beginning allows for reproducible results.Recall that ts1 is a function to simulate from the QAR(1) model to be studied.N is the length of the time series and also the number of Fourier frequencies for which the quantile periodograms will be computed.By the parameter freq a subset of these Fourier frequencies is specified to be stored; a subset is used to save storage space.The estimates at these frequencies freq and at the specified levels are then stored to the array sims.In this example, the smoothed periodograms are computed using the Epanechnikov kernel and the (rather large) bandwidth of b n = 0.3.Note that the flexible accessor method getValues is used to access the relevant subset of values for the frequencies specified (i.e., freq).Once the array sims is available, many interesting properties of the estimator can be analyzed.Examples include the bias, variance, etc.Here, using the function getValues again, the true copula spectral density is copied to an array trueV.Using the arrays sims and trueV the root integrated mean squared errors are computed as follows: These numbers could now be inspected to observe, for example, that the smoothed quantile periodogram possess smaller root integrated mean squared errors than the quantile periodograms (without smoothing).Further discussion of the numbers is omitted, because the crux of this chapter was to explain how to perform the simulation study, not to actually do it.

Roadmap to future developments and concluding remarks
As the new methodology evolves additional features will be added to the quantspec package.For each new feature an entry to the issue tracker available on the GitHub will be made.Then the new feature will be implemented on a topic branch of the repository.For example, a procedure for data-driven selection of the bandwidth is currently being developed.
Other, more complex extensions to the software include the implementation of functions to perform quantile spectral analysis for locally stationary processes, the computation and smoothing of higher order quantile periodograms for the estimation of quantile polyspectra.Procedures for graphical representations of these objects, possibly animated ones, will amend these planned parts of the package.An overview on the planned extensions will be made available in the issue tracker on GitHub.
Summing up, it can be said that the quantspec package provides a comprehensive and conclusive toolbox to perform quantile-based spectral analysis.Due to the great interest in and active development of the statistical procedures that are quantile-based spectral analysis it was deliberately designed in an object-oriented and extensible fashion.Thus it is well prepared for the many extensions that are sure to come in the near future.The source code is open and extensive documentation of the system freely available.Comments on and contribution to the project is, of course, very much welcome.

Figure 4 :
Figure 4: Returns (Y t ) of the S&P 500 returns example data (left), autocovariances Cov(Y t+k , Y t ) of the returns (middle), and autocovariances Cov(Y 2 t+k , Y 2 t ) of the squared returns (right).