Costationarity of locally stationary time series

This article describes the R package costat. This package enables a user to (i) perform a test for time series stationarity; (ii) compute and plot time-localized autocovariances, and (iii) to determine and explore any costationary relationship between two locally stationary time series. Two locally stationary time series are said to be costationary if there exists two time-varying combination functions such that the linear combination of the two series with the functions produces another time series which is stationary. Costationarity existing between two time series indicates a relationship between the series that might be usefully exploited in a number of ways. Sometimes the relationship itself is of interest, sometimes the derived stationary series is of interest and useful as a substitute for either of the original stationary series in some applications.


Summary
The costat package is designed for the analysis of locally stationary time series, particularly locally stationary wavelet time series.The package includes functionality for • computing tests of stationarity, BootTOS; • computing localized autocovariances, lacv; • discovering costationarity between two time series, findstysols Several method functions exist that print, summarize or plot the various outputs of these functions.

Introduction: Locally stationary time series
This article concerns the analysis of discrete time series, X t .That is, a sequence of observations taken through time where t could be an integer or some other regular spacing.Many are familiar with the concept of a stationary time series: that is, loosely speaking, a series whose statistical properties do not change over time.In much of classical time series analysis stationary means second-order stationary where the mean and variance of a series, X t , do not depend on time and the autocovariance: only depends on the lag τ between two different variates X t and X t+τ , but not the time point t = 1, . . ., T .There are several excellent books dealing with (stationary) time series analysis, for example, Brockwell and Davis (1991), Chatfield (2003), Hamilton (1994), Hannan (1960), Priestley (1983) or Shumway and Stoffer (2006).
Mathematical theory demands that a stationary time series X t possesses the following representation: where A(ω) is an amplitude function, {e itω } is the usual system of harmonic complex exponentials and dz(ω) is an orthonormal increments process.The amplitude function controls the second-order properties of the time series.The usual spectrum f (ω) = |A(ω)| 2 and the spectrum and autocovariance are a Fourier transform pair.The amplitude, spectrum and autocovariance here reflect the time-invariant nature of the time series X t : none of them depend on time, t.
However, many actual time series, arising in many disciplines, are not stationary.Sometimes the nonstationarities can be detected via a simple time series plot.Another way of detecting nonstationarity is to compute a simple statistical time series statistic on different sections of the time series and look for differences.For example, using the sample autocovariance function acf in R, R Development Core Team (2009).Such exploratory methods for detecting nonstationarities are undoubtedly useful, but somewhat ad hoc, and more statistically reliable methods are necessary.
Over the years several extensions to the basic stationary model (2) have been proposed to deal with nonstationary time series.Important general classes of models are the locally stationary time series models.One locally stationary extension of (2) replaces the time-constant A(ω) amplitude function by one which explicitly depends on time, for example A t (ω).For example, the oscillatory processes of Priestley (1983) or the locally stationary Fourier (LSF) processes introduced by Dahlhaus (1997), and more recently Dahlhaus and Polonik (2006) and Dahlhaus and Polonik (2009).A comprehensive and recent review of locally stationary processes can be found in Dahlhaus (2012).
Another locally stationary extension, developed by Nason, von Sachs, and Kroisandt (2000), replaced the exponential harmonics by nondecimated discrete wavelets as well as introducing time-locality into the expansion.Their locally stationary wavelet (LSW) model was given by where w j,k is the amplitude sequence, {ψ j, } are the nondecimated discrete wavelets and ξ j,k is a sequence of random variables satisfying COV(ξ j,k , ξ ,m ) = δ j, δ k,m , where δ k,m is the Kronecker delta function.The ξ j,k also satisfy E(ξ j,k ) = 0 which means that EX t = 0: any real time series that does not have zero mean can be detrended, see Chatfield (2003), for example.Additionally, deterministic trend needs to be removed before modelling using LSW processes.Additional details on the LSW processes and examples of their use can be found in Nason (2008).As with stationary and LSF processes the amplitude sequence w j,k is linked to a spectral quantity by w 2 j,k ≈ S j (z), where z = k/T is rescaled time, z ∈ (0, 1) and k = 1, . . ., T .The quantity {S j (z)} ∞ j=1 is known as the evolutionary wavelet spectrum.Assuming that a time series X t is appropriately modeled by a LSW process, a key task is to estimate its spectrum using an estimator { Ŝj (z)} J j=1 } where 2 J = T .For a finite length T time series the spectrum can only be computed for a finite set of scales 1, . . ., J. The wavethresh package, Nason (2012a), contains a function ewspec which estimates the evolutionary wavelet spectrum of a time series and there are various other functions to both plot the time-varying spectrum and extract its values for further processing.Again, many more details on these kinds of practical analyses can be found in Nason (2008).
Much recent research has considered the case of nonstationary models of one kind or another, and methods of estimating the associated spectral quantities.Fewer articles consider the tricky problem of given an actual time series how does one decide whether it is a stationary or not?Some examples of stationarity tests are the early paper by Priestley and Subba Rao (1969), tests that use wavelets by von Sachs and Neumann (2000) and Nason (2013), Ahamada and Boutahar (2002) in economics, Starica and Granger (2005) and Paparoditis (2009) both which measure the difference between a periodogam and a model spectrum on intervals and Dwivedi and Subba Rao (2009).
We should mention at this point that we are considering the testing of second-order stationarity of locally stationary time series.That is, we assume that the mean is constant and there is no trend (or the series has been detrended and deseasonalised).There is a substantial literature on the topic of stationarity testing for time series exhibiting trends, see, for example, the KPSS test, Kwiatowski et al. (1992), and the related problem of unit-root testing, see, e.g.Dickey and Fuller (1979), Elliott et al. (1996) or Phillips and Perron (1988).
In R the package fractal includes an implementation of the Priestley-Subba Rao test, called stationarity, modified with improved spectral density function estimators based on averaging multitaper estimators, see Constantine and Percival (2007).This article describes the (new) test of stationarity introduced in Cardinali and Nason (2010) and presented in the costat package: this is described in Section 3.
The "flip" side to spectral quantities are covariances.In the stationary theory the autocovariance function is the inverse Fourier transform of the spectrum.For locally stationary processes attention has focussed on spectral quantities such as f (z, ω) or S j (z) although some works have concentrated on the covariance side notably Sanderson, Fryzlewicz, and Jones (2010), parts of Nason et al. (2000) and Nason (2013).Section 4 provides details on how to compute both localized autocovariance (LACV) from Nason et al. (2000) and how to present the results in attractive and useful formats.
Finally, in contrast to the stationary or nonlinear time series world, little work has been conducted on the case of where there is more than one locally stationary time series to be considered, in other words, multivariate time series.Again, Sanderson et al. (2010) is a notable exception.This article considers the implementation of the costationary methods introduced by Cardinali and Nason (2010) in Section 5. Given two locally stationary time series the goal of costationarity determination is to investigate whether there are two sequences α t , β t that can be found such that Z t is a classically (second-order) stationary time series where The concept of costationarity is inspired by the cointegration of Engle and Granger (1987) except that order-1 nonstationarity is replaced by order-2, and the costationary vectors, (α t , β t ) are time-varying (necessarily for a non-trivial theory).In cointegration the X t , Y t processes are assumed to be integrated (differencing required to achieve stationarity) and a linear combination (α, β) is sought such that Z t = αX t +βY t is a stationary series.In our work X t , Y t have time-varying second-order statistics and the linear combination in ( 4) is sought to produce a second-order stationary (time invariant second-order statistics) series.
The rationale for costationary detection/determination are (i) learning of any costationary relationship itself, as this could be interesting, (ii) estimating the strength of such a relationship, (iii) using the derived stationary series Z t , in some applications in preference to either of the original series and (iv) using the relationship to learn about X t from data on Y t or vice versa.Cardinali and Nason (2010) give two examples where a costationary relation is discovered.One example, from financial portfolio planning, obtains asset portfolios that are shown to be preferable to classically defined ones or others defined using classical methods augmented using conditional variance estimation techniques (essentially VECH GARCH).Their other example, investigates costationary combinations of volatile wind power series to obtain a less volatile series.Volatility is a real problem for wind power as generators prefer stable and predictable sources and wind is anything but.Usable power from wind can be obtained by aggregating wind power over geographical areas and costationarity can help discover useful combinations to enhance stability.

The Test
Given a time series one might wish to test whether it is stationary or not.The costat package incorporates a function BootTOS which performs a simple bootstrap based test of stationarity.Given a time series X t , assumed to be modeled by a LSW process with spectrum {S j (z)} J j=1 , the null hypothesis of the costat test of stationarity is: versus the alternative H A : S j (z) is not constant for some j.Details of the test appear in Cardinali and Nason (2010) but the test is based on the metric: where Sj = 1 0 S j (z)dz.Clearly, if S j (z) is constant, for all j, then T ({S j (z)}) = 0 and T is nonzero for nonconstant spectra.The statistical significance of the test statistic T computed on the time series realization is assessed via a parametric bootstrap calculation.Multiple comparator simulations are produced under the null hypothesis assuming constant spectral quantities estimated from the data and then the test statistic evaluated on those simulations.
To demonstrate BootTOS we use the SP500FTSElr data distributed with costat which contains log returns of both the well-known SP500 and FTSE financial index series between 21st June 1995 and 2nd October 2002.For the purposes of this article we will examine a smaller segment of this data called sret and fret which are the SP500 and FTSE log returns extracted from the SP500FTSElr object from index 256 to 767.The methods we describe below do work on the whole series contained in SP500FTSElr but can take quite a long time on single-processor machines especially when the main purpose here is pedagogical.
Should the reader wish, it is quite easy to convert the SP500FTSElr object into a multivariate time series objects present in R. For example, a timeSeries class object by loading the timeSeries package and typing and then this new times series, newts, can be plotted or analyzed using the functions present in the timeSeries package.
The log returns of a sequence X t are given by log(X t /X t−1 ) and is a common quantity of analysis for economic and financial data.The reader will note that the extracted series are of length 512 which is a power of two.Like many wavelet-based codes the costat software works on data sets that are of dyadic length (i.e. 2 J for some natural number J).It should be made clear that this is not a limitation of wavelets per se, but of the computationally efficient algorithms used to compute the intended quantities.Data sets of other lengths can be handled by zero-padding or truncation.The data can be accessed by typing: R> library("costat") R> data("sret", "fret") As with many data analyses it is useful to first produce a plot of the series.We can do this with the following code.

R> ts.plot(sret) R> lines(fret, col=2)
This code produces the plot shown in Figure 1 which indicates that both series are probably not stationary as there appears to be clear variance changes in the series.The following code performs the bootstrap test of stationarity mentioned above: The number of bootstrap simulations is controlled by the Bsims argument which is set to 100 by default.The BootTOS test is computationally intensive.However, it can use the mclapply function of the multicore package, Urbanek (2011), to execute the bootstrap simulations in parallel.On a dual core (Intel Core 2 Duo, 2.8Ghz) iMac the bootstrap test took approximately 10.4 seconds, on a twelve core (Intel Xeon, 3.07Ghz) Linux machine the test took 1.3 seconds.
The BootTOS function returns an object of class c("BootTOS", "htest").One can see the results of the significance test by typing the object name as follows: The function output indicates that the p-value of the test is less than 0.01 and hence the sret series is nonstationary as its spectrum is much more variable that could be obtained by a stationary equivalent (i.e. with the same marginal S j values). As

Size and Power: Simulation Results
With any statistical test it is important to have some understanding of its size and power characteristics which we do here via simulation.
Our size and power results involve three tests of second-order stationarity.The Priestley and Subba Rao (1969) test as implemented in the package fractal, the Haar wavelet test of stationarity from Nason (2013) and BootTOS as described above.We shall abbreviate the names of these tests in our tables as PSR, HWTOS and BootTOS.
To assess size, we simulate data from a number of stationary models using Gaussian innovations and assess how often our tests of stationarity reject the null hypothesis.The models are: S1 iid standard normal; S2 AR(1) model with AR parameter of 0.9 with standard normal innovations; S3 As S2 but with AR parameter of −0.9; S4 MA(1) model with parameter of 0.8; S5 As S4 but with parameter of −0.8.
Models S1-S7 are the same as in Nason (2013) to aid comparisons with another test of stationarity.
The empirical size values are shown in Table 1.Note: the PSR column is a direct copy from Nason (2013) whereas those for HWTOS and BootTOS are new for this paper (although column HWTOS shows similar values to that in Nason (2013)).Overall, the size characteristics of BootTOS are extremely good and conservative apart from S3 (where actually none of the tests meet the 5% nominal size, although BootTOS performs much the worst) and S7 where PSR Power.To explore statistical power we create processes that are nonstationary and then count the number of times each test reckons a realisation is not stationary over multiple simulations.Again, we use models from Nason (2013) as follows.
P1 Time-varying AR model X t = α t X t−1 + t with iid standard normal innovations and the AR parameter evolving linearly from 0.9 to -0.9 over the 512 observations.
P3 A LSW process based on Haar wavelets with spectrum S j (z) = 0 for j > 2 and S 1 (z) as for P2 and S 2 (z) = S 1 (z + 1 2 ) using periodic boundaries (for the construction of the spectrum only).
The spectra and single realisations for these processes are displayed in Nason (2012b).All tests have excellent power for model P4.BootTOS has excellent power on P1 and P2 where one of PSR (for P1) or HWTOS (for P2) do not perform at all well.No test performs particularly well for P3 although PSR does best and BootTOS is second both detecting some nonstationarity.HWTOS for P3 performs particularly poorly.
Overall, although no test can be declared the 'winner' in these experiments the empirical performance of BootTOS is creditable.

Local autocovariance
This section introduces functions that compute and plot localized autocovariances (LACV).

LACV on a stationary series
Before analyzing a nonstationary example, we use the localized autocovariance on a stationary series, in order to assess the correctness of the answer given by comparing the result to that produced by the standard R function acf.
Even so, we should not expect the LACV to be as accurate as the regular acf function on a stationary series.This is because the acf function will produce its estimate by averaging over the full time period of quantities whose expectation is constant over that time period.The localized function will use local smoothing and not take advantage of our prior knowledge of the constancy of the autocorrelations in this specific case.
Computing the regular acf.For our stationary series we choose to generate the AR time series X t = 0.8X t−1 + Z t , and also print the autocovariance.Note that the lag one autocorrelation of 0.822 estimates the AR parameter of 0.8.

R> vsim
Computing the localized autocovariance.We can use costat's function, lacv, for computing the localized autocovariance, c(z, τ ), which is an implementation of the formula: from Nason et al. (2000), where Ψ j (τ ) = k ψ j,k ψ j,k−τ is the autocorrelation wavelet of the discrete nondecimated wavelet ψ j,k .
The quantity c(z, τ ) in ( 7) is the autocovariance of X t at lag τ and at rescaled time z ∈ (0, 1).Rescaled time is related to real time by the formula z = t/T for time points t = 1, . . ., T where T is the length of the time series.The lacv function produces an object of class lacv which contains the localized autocovariance and autocorrelation information.The second command, where we typed the name of the object vsim.lacv,invokes the print.lacvfunction which informs the user of the type of object, subsequently calls summary.lacvwhich gives information on when the object was produced and its size.
The vsim data series contained T = 1024 data points.Without specifying the lag.max=100 argument the vsim.lacvwould have contained large number of lags: 7162.We have also chosen to compute the localized autocovariance with the filter.number=4argument which chooses a particular wavelet.Using different wavelets will result in different results for the autocovariance estimators but the differences will typically be small.The wavelet with four vanishing moments can be produced using the draw.defaultfunction in wavethresh and has intermediate smoothness between the piecewise constant Haar wavelet and the very smooth wavelet with 10 vanishing moments (selected by filter.select=1and filter.select=10arguments respectively).More information on the different types and smoothnesses of wavelets can be found in the wavethresh help function on filter.select.
The call to lacv to produce vsim.lacvtook 3.8 seconds on the iMac and 1.9 seconds on the Linux machine.
Plotting the localized autocovariance.The LACV is a function of two arguments and hence there are various ways it can be plotted to elicit information about its form.
Since the LACV is potentially a quantity that can vary over time our first plot will draw the ĉ(z, τ ) estimates over all time z ∈ (0, 1) for a fixed number of lags ranging from lag τ = 0 to τ = 7. Remember z is rescaled time and can be mapped back to the actual time points by z = t/T .

Execute the following commands:
R> plot(vsim.lacv,lags=0:7, lcol=1:8) R> abline(h=0.8,lty=2) The output from this code is shown in Figure 3.The autocorrelation at each lag τ = 0 to τ = 7 is plotted as a solid line in Figure 3 with the lag associated with each line plotted as integers superimposed on the line.Note, the lag τ = 0 line is completely straight because autocorrelation is always one at τ = 0 for regular and localized versions.
At lag τ = 1 the theoretical autocorrelation value for this process is 0.8.The regular acf estimate computed above by acf was 0.822.Figure 3 shows the localized autocorrelation varies: on the first half of the series it is above 0.8 and for the second half it tends to be below.Although the underlying process is stationary, and all theoretical autocorrelations constant, the lacv function does not assume this and so all the estimated autocorrelations vary (apart from the lag zero).However, the values of ĉ(z, τ ) are not highly variable.
Would we produce a plot such as Figure 3 in practice?Well, normally we would perform a test of stationarity first.Applying BootTOS to vsim results in a p-value of 0.97 which means that we would accept the hypothesis of stationarity for this series and, in these circumstances not compute the localized version.
The plot.lacv method for lacv objects also permits plotting of ĉ(z, τ ) for a given fixed value of z (or t) in the style of the usual plot produced by the R function acf.Type: R> plot(vsim.lacv,type="acf",the.time=512)The plot.lacv method can also be used to produce perspective plots of the autocorrelation function using the type="persp" argument.

LACV on a simulated locally stationary series
We first construct a function, tvar1sim, to generate realisations from a time-varying autoregressive process TVAR(1) as follows: n <-512 arvec <-seq(from = 0.9, to = -0.9,length = 512) x <-c(rnorm(1, mean = 0, sd = sd), rep(0, n -1)) for (i in 2 The model behind this process is X t = α t X t−1 + t , where { t } is an iid white noise process with variance of σ 2 .The AR parameter of this process varies from α = 0.9 to α = −0.9 over the length of the process. We simulate a realisation from this TVAR(1) process, and plot it using R> tvar1.sim<-tvar1sim() R> ts.plot(tvar1.sim)The plot produced is shown in Figure 5.The LACV can be computed and plotted using the following commands: R> tvar1.lacv<-lacv(tvar1.sim,filter.number=4,lag.max=50)R> plot(tvar1.lacv,lags=0:2, lcol=1:3) Figure 6 shows the LACV plot where we have chosen to show up to lag two for clarity.The lag one autocorrelation decays nicely from positive 0.81 to negative 0.71 across the length of the series.This decay is not perfectly linear from 0.9 to -0.9 (as specified in the function tvar1sim()) but the estimate is from a finite length realisation and even stationary ACF estimation is not perfect (For example, in the stationary AR(1) α = 0.8 model in Section 4.1 the estimate was 0.822 and that result was achieved by averaging over 1023 pairs all with the same distribution!)

LACV on a locally stationary series
We now compute and plot the localized autocovariance function on the sret SP500 log-returns data set using the following commands R> sret.lacv<-lacv(sret, filter.number=4,lag.max=50)R> plot(sret.lacv,lags=0:5, lcol=1:6) which produces the picture displayed in Figure 7.We have only chosen to display the first five lags so as to provide an uncluttered plot.Figure 7 shows that the autocorrelations are all reasonably small (less than 0.2 in magnitude) apart from a big spike at just prior to t = 350.Why do we see the spike?The answer is to look back at the actual time series and one indeed can see a negative outlier at time t = 332 in Figure 1.Such jump changes are permitted within the locally stationary framework and the jump is effectively localized in that the remaining autocovariances in Figure 7 will be relatively unaffected.
To see this, we remove this outlier, and a few others, using the commands: The LACV plot is shown in Figure 8. Log-returns time series of this kind are often modeled with a GARCH process.One of the features of GARCH is that the empirical autocovariance function is zero for all non-zero lags and this seems to be the case here.Repeated regular applications of acf were applied to different portions of the series and most autocovariances on most occasions were deemed to be not significantly different from zero.
One could argue that the outlier means that the autocorrelation spike contradicts the zero autocorrelation finding.However, GARCH models often admit heavy tails and so the large negative observation would not usually be thought of as an outlier.So, although the outlier only disturbs our LACV locally it would be better, maybe, to construct a LACV estimator which is more robust to this kind of disturbance.Such a concept is desirable but beyond the scope of the current article.

Commands to discover costationarity
Finally, we examine the case where we wish to discover any potential costationary solutions.That is, given time series X t , Y t find α t , β t such that we obtain (4) where Z t is second-order stationary.The main function to discover costationarity is findstysols.This starts with a random choice of the vectors α t , β t and, on any one run, numerically optimizes the vectors to try and find the combination which results in the flattest spectrum of Z t over time.Finally, the 'optimal' Z t is tested using the test of stationarity mentioned above to see whether the resultant Z t can be statistically assessed to be stationary.Since the output depends on the random input it is good practice to restart the algorithm using many random starts.The number of random starts is controlled by the Nsims argument.In general, you should run as many invocations with different random starts as your computational resources will allow.
There is a function, mergexy that can merge the outputs from many runs of the findstysols function into a single output file.This means that parallel computation can be used to great effect by using multiple serial jobs to run from the repeated random starts.In the examples below Nsims is set to be 10, which is really too small.Also, Nsims should grow with the length of the time series.For longer series, there are typically more ways, and more flexibility, in finding stationary solutions.However, there are no hard and fast rules concerning the choice of Nsims.
Let us try and discover costationary solutions between the sret and fret time series.To repeat this discovery process using ten random starts we can issue the command: R> sretfret.fss<-findstysols(Nsims=10, tsx=sret, tsy=fret) Note: this command takes quite some time to execute.On our dual core iMac it took 8 minutes and on the twelve core Linux machine it took 48 seconds with this using the package multicore and additionally using the lapplyfn="mclapply" within the function call above.The summary output contains some basic information on the names and dimensions of the input time series and solutions sought.In this case we asked for ten optimization runs to be executed and indeed ten sets of solutions were produced.However, the summary indicates that not all optimization runs converged.The reason why some did not converge can be gained by looking at the convergence component of the csFSS object.For example,

R> sretfret.fss$convergence
[1] 0 0 1 0 0 0 10 0 1 0 As the text in the summary output indicates these numbers correspond to the status code returned by the optim function that was used to perform the optimizations.The meaning of the status code can be revealed by examining the help page for optim and looking at its convergence output.In particular, a code of zero indicates successful convergence, a code of one indicates that optim ran out of iterations before it deemed convergence had occured.
The single occurence of code ten indicates degeneracy of the Nelder-Mead simplex in that method.For longer, and possibly more complex series, one often sees multiple error codes of one indicating that the optimizers require more iterations.Extra iterations can be given by increasing the my.maxit argument or directly passing control information to the optim calls by the optim.controlargument of findstysols.
In any case the solution above did converge seven times and all of them appear to be stationary (in that all p-values were greater than 0.05.The value that controls this assessment for stationarity, the size argument in the call to summary can be changed, if required).So far, we have performed ten optimizations, seven have converged and look interesting.The commands described next allow us to further investigate the nature of the solutions.

Plots for obtaining an overview of all solutions
Each set of solutions consists of a set of of parameters that encode the α t and β t functions (actually wavelet coefficients of the functions, the number of such coefficients controlled by the Ncoefs argument to findstysols) and this is stored as a multivariate data set with the number of cases (rows) equal to the number of solutions and the number of variables equal to twice Ncoefs as there are coefficients associated with each of α t and β t .
The costat package provides two functions to display and interpret this information.One method uses hierarchical cluster analysis and R's hclust function which clusters similar solutions together.For example, for the current example, the command R> plot(sretfret.fss,ALLplotscale=FALSE) produces Figure 9: a dendrogram showing how solutions 1 and 4 are very similar as well as solutions 5 and 10.Running the plot with the ALLplotclust argument set to FALSE causes a plot of a multidimensional scaling solution applied to a Euclidean distance calculation on the set of solutions (applying cmdscale to dist on the set of solution vectors).This latter plot would show the solution ids on a two-dimensional plot which again gives some indication of how similar different solution sets are.Suppressing both the ALLplotclust and ALLplotscale arguments causes both plots to be produced.
The idea of these 'complete set' plots is to identify which solutions are worth examining further.For example, from the above plot there is little point examining both solution 1 and 4 in detail since they are very similar, so one can examine one of them and that will be representative.Typically, one would compute many more solutions.In the initial call to findstysols the Nsims argument might be a lot higher, e.g.100.In such a case one will find that the solutions tend to group into clusters of like-valued solutions and one need only examine a single representative from each cluster, maybe using the plots described in the next section.

Plotting information about individual solutions
After the solutions are examined collectively one can then choose to investigate specific solutions.Individual solutions can be examined by using again the plot function with the additional solno argument.The solno argument can be a vector of solution numbers and a series of four plots is produced for each solution.It can be convenient to produce all four graphics pertaining to a solution on the same plot.This can be achieved using par(mfrow=c(2,2)) or using layout such as in the following example: R> def.par <-par(no.readonly= TRUE) # For resetting R> nf <-layout(matrix(1:4, 2,2, byrow=TRUE)) R> plot(sretfret.fss,solno=6) R> par(def.par)# Resets plot to single plot, if required produces Figure 10 for solution number 6.This figure shows the costationary vectors (α t , β t ), that were used to produce the combined series Z t , which is also shown along with its associated local spectral estimate.Figure 10 shows how the costationary vectors change over time.The bottom row of the plot demonstrates that the resultant Z t series looks stationary with a fairly Figure 10: Top left: optimal α t costationary vector; Top right: optimal β t costationary vector; Bottom left: the combined series Z t ; Bottom right: the evolutionary wavelet spectral estimate of the series Z t with the computed p-value for the test of stationarity for the Z t series.
flat spectrum.More information can be obtained using the plotstystat=TRUE argument.Using this argument additionally plots the combined time series, the regular and partial autocorrelations of the combined series and a spectral estimate.These plots can be produced as the test of stationarity has deemed that this Z t is stationary (or at least, no evidence that it is not)

Figure 1 :
Figure 1: Time series plot of the sret (black) and fret (red) time series.

Figure 3 :
Figure 3: Localized autocorrelation plot from a realisation of a stationary AR(1) process with α = 0.8.The horizontal dashed line is plotted at height of 0.8

Figure 4 :
Figure 4: Localized autocorrelations plotted on an acf style plot at time t = 512 from a realisation of a stationary AR(1) process with α = 0.8.The stars are the values of the true autocorrelations for this process.

Figure 6 :
Figure 6: The LACV plot of the TVAR realisation shown in Figure 5.

Figure 8 :
Figure 8: Localized autocorrelation values for SP500 log-returns time series with outliers removed.

Table 1 :
Nason (2013)ple, let us run the same test but on a series we know to be stationary.Let us create a new series which samples 512 observations from a standard Gaussian distribution and repeat the test of stationarity.Here one can see that the p-value of the test is 0.51, and hence we can feel comfortable in assessing the series x to be stationary (as we knew it would because we set it up that way).The plot from the plot(x.TOS) command is shown in Figure2.The histogram shows all the test statistic values computed on the data and the bootstrap simulations.The vertical line corresponds to the test statistic computed on the original data.It can be seen that roughly Simulated size estimates (%) for stationary Gaussian models T = 512.Based on 1000 simulations with 1000 bootstrap simulations per main simulation.Note: column PSR corresponds to simulation results fromNason (2013)for 100 simulations.
R> x <-rnorm(512) R> x.TOS <-BootTOS(x) 51% of the bootstrap samples are higher in value.Hence, the sample test statistic is nothing unusual compared to a bootstrap sample on comparable distributional setup and hence we have no reason to reject the null hypothesis.

Table 2 :
Simulated power estimates (%) for models P1-P4 with nominal size of 5% BootTOS better, but still not near to 5%.Note, PSR does not do very well in S2 or S6 where HWTOS and BootTOS do well.
The returned sretfreet.fssobject is of class csFSS and contains a lot of information concerning the solutions returned by the multiple optimizations.Printing out the results displays: