Bayesian Lattice Filters for Time-Varying Autoregression and Time-Frequency Analysis

Modeling nonstationary processes is of paramount importance to many scientific disciplines including environmental science, ecology, and finance, among others. Consequently, flexible methodology that provides accurate estimation across a wide range of processes is a subject of ongoing interest. We propose a novel approach to model-based time-frequency estimation using time-varying autoregressive models. In this context, we take a fully Bayesian approach and allow both the autoregressive coefficients and innovation variance to vary over time. Importantly, our estimation method uses the lattice filter and is cast within the partial autocorrelation domain. The marginal posterior distributions are of standard form and, as a convenient by-product of our estimation method, our approach avoids undesirable matrix inversions. As such, estimation is extremely computationally efficient and stable. To illustrate the effectiveness of our approach, we conduct a comprehensive simulation study that compares our method with other competing methods and find that, in most cases, our approach performs superior in terms of average squared error between the estimated and true time-varying spectral density. Lastly, we demonstrate our methodology through three modeling applications; namely, insect communication signals, environmental data (wind components), and macroeconomic data (US gross domestic product (GDP) and consumption).


Introduction
Recent advances in technology have lead to the extensive collection of complex high-frequency nonstationary signals across a wide array of scientific disciplines.In contrast to the timedomain, the time-varying spectrum may provide better insight into important characteristics of the underlying signal (e.g., Holan et al., 2010Holan et al., , 2012;;Rosen et al., 2012;Yang et al., 2013, and the references therein).For example, Holan et al. (2010)  In general, time-frequency analyses can either proceed using a nonparametric or modelbased (parametric) approach.The most common nonparametric approach is the short-time Fourier transform (i.e., windowed Fourier transform) which produces a time-frequency representation characterizing local signal properties (Gröchenig, 2001;Oppenheim and Schafer, 2009).Another path to time-frequency proceeds using smoothing splines (Rosen et al., 2009(Rosen et al., , 2012) ) or by parameterizing the spectral density to estimate the local spectrum via the Whittle likelihood (Everitt et al., 2013).Similarly, time-frequency can be achieved by applying smooth localized complex exponential (SLEX) functions to the observed signal (Ombao et al., 2001).In contrast to window based approaches, the SLEX functions are produced using a projection operator and are, thus, simultaneously orthogonal and localized in both time and frequency.Alternatively, one can use the theory of frames and over-complete bases to produce a time-frequency representation.For example, continuous wavelet transforms (Vidakovic, 1999;Percival and Walden, 2000;Mallat, 2008) or Gabor frames (Wolfe et al., 2004;Feichtinger and Strohmer, 1998;Fitzgerald et al., 2000) could be used.By introducing redundancy into the basis functions, these representations may provide better simultaneous resolution over both time and frequency.
Model-based approaches typically proceed through the time-domain in order to produce a time-frequency representation for a given nonstationary signal.In this setting common approaches include fitting piecewise autoregressive (AR) models as well as time-varying autoregressive (TVAR) models.The former approach assumes that the nonstationary signal is piecewise stationary.Consequently, the estimation procedure attempts to identify the order of the AR models along with the location of each piecewise stationary series.For example, Davis et al. (2006) propose the AutoPARM method using minimum description length (MDL) in conjunction with a genetic algorithm (GA) to automatically locate the break points and AR model order within each segmentation.In addition to providing a time-frequency representation, this approach also locates changepoints.Wood et al. (2011) propose fitting mixtures of AR models within each segment via Markov chain Monte Carlo (MCMC) methods.Their approach selects a common segment length and then divides the signal into these segments prior to implementation of the fitting procedure.Although such approaches may accommodate signals with several piecewise stationary structures, they lack the capability of capturing momentary shocks to the system (i.e., changes to the evolutionary structure that only occur over relatively few time points).
For many processes, TVAR models may provide superior resolution within the timefrequency domain for both large and small scale features through modeling time-varying parameters.To estimate the TVAR model coefficients, Kitagawa and Gersch (1996) and Kitagawa (2010) treat the coefficients as a stochastic process and model them using difference equations under the assumption of a maximum fixed order of the TVAR model.Their estimation procedure for the coefficients is then based on state-space models with smoothness priors.In this context, the innovation variances are treated as constant and estimated using a maximum likelihood approach.Subsequently, West et al. (1999) propose a fully Bayesian TVAR framework that simultaneous models the coefficients and the innovation variances using random walk models.Alternatively, by assuming a constant innovation variance, Prado and Huerta (2002) model the coefficients and order of the TVAR model using random walk models.Further, to make the TVAR models stable, the constraint that the roots of the characteristic polynomial lie within the unit circle could be imposed.However, such an added condition makes estimation more complicated and computationally expensive.
To avoid these issues, we instead work with the partial autocorrelation coefficients (i.e., in the partial correlation (PARCOR) coefficient domain) and then use the Levinson recursion to connect the PARCOR coefficients and TVAR model coefficients (Kitagawa and Gersch, 1996;Godsill et al., 2004).Godsill et al. (2004) model the PARCOR coefficients and innovation standard deviations using a truncated normal first-order autoregression and a log Gaussian first-order autoregression, respectively, with a given constant order.To estimate these values, a sequential Monte Carlo algorithm is used.Alternatively, as previously alluded to, by assuming a constant innovation variance, Kitagawa and Gersch (1996) implement the smoothness prior within a lattice filter to estimate the PARCOR coefficients.After the PAR-COR coefficients have been estimated, a constant innovation variance is estimated using a maximum likelihood approach.However, the former approach is computationally expensive and may suffer from the degeneracy problem (i.e., the collapse of approximations of the marginal distributions) when TVAR model order is large.In addition, certain hyperparameter values (i.e., the TVAR coefficients associated with the two latent models) may be sensitive to starting values and may require prior knowledge or expert supplied subjective information to achieve convergence (Godsill et al., 2004).Although the latter takes advantage of the lattice form to estimate the PARCOR coefficients, estimation of the innovation variance is achieved outside of the lattice structure; that is, the estimation procedure is a two-stage method.Importantly, this approach is designed for a constant innovation variance and cannot deal with time-dependent innovation variances.As such, we propose a novel approach that addresses these issues within a fully Bayesian context.Kitagawa (1988) and Kitagawa and Gersch (1996) constitute the first attempts at utilizing the lattice structure to estimate the coefficients of a TVAR model when the innovations follow a Gaussian distribution with mean zero and constant variance.At each stage of the lattice filter, they assume that the residual at each time between the forward and backward prediction errors follows a Cauchy distribution, and that the PARCOR coefficient is modeled as a Gaussian random walk.This produces a non-Gaussian state space model at each stage and thus, a numeric algorithm is conducted for estimation.Moreover, the assumptions of their approach ignore an implicit connection between the innovation term of the TVAR model and the residual term between the forward and backward prediction errors.That is, the distribution of residual term should be a Gaussian distribution rather than a Cauchy distribution.Consequently, their approach leads to a TVAR model with innovation terms following a Cauchy distribution; hence the innovation variances do not exist.In contrast, our approach assumes the residual term follows a Gaussian distribution.
We propose a fully Bayesian approach to efficiently estimate the TVAR coefficients and innovation variances within the lattice structure.One novel aspect of our approach is that we model both the PARCOR coefficients and the TVAR innovation variances within the lattice structure and then estimate them simultaneously.This is different from the frequentist twostage method of Kitagawa (1988) and Kitagawa and Gersch (1996).Another novel aspect is that we take advantage of dynamic linear model (DLM) theory (West and Harrison, 1997;Prado and West, 2010) to regularize the PARCOR coefficients instead of using truncated distributions.Thus, our method provides marginal posterior distributions with standard forms for both the PARCOR coefficients and innovation variances.Since our approach takes advantage of the lattice structure, the computational efficiency of our approach is not affected by the order of the TVAR model; that is, our approach avoids having to calculate higher dimensional inverse matrices.To select the TVAR model, we provide both a visual and a numerical method.Importantly, the simulation study we provide demonstrates that our approach leads to superior performance in terms of estimating the time-frequency representation of various nonstationary signals, as measured by average squared error.Thus, our approach provides a stable and computationally efficient way to fit TVAR models for time-frequency analysis.
The remainder of this paper is organized as follows.Section 2 briefly introduces the lattice structure and describes our methodology along with prior specification.Section 3 presents a comprehensive simulation study that illustrates the effectiveness of our approach across an expansive array of nonstationary processes.Subsequently, in Section 4, our methodology is demonstrated through three modeling applications; namely, insect communication signals, environmental data (i.e., wind components), and macroeconomic data (i.e., US gross domestic product (GDP) and consumption).Lastly, Section 5 concludes with discussion.
For convenience of exposition, details surrounding the estimation algorithms and additional figures are left to an Appendix.

Time-Varying Coefficient Autoregressive Models
The TVAR model of order P for a nonstationary univariate time series x t , t = 1, . . ., T , can be expressed as where a (P ) t,m and ǫ t are the TVAR coefficients associated with time lag m at time t and the innovation at time t, respectively.Typically, the innovations are assumed to be uncorrelated mean-zero Gaussian random variables (i.e., ǫ t ∼ N(0, σ 2 t ), with time-varying variance σ 2 t ).Therefore, the TVAR model corresponds to a nonstationary AR model with the AR coefficients and variances evolving through time.In such settings, the model is locally stationary but nonstationary globally.As will be illustrated, the assumption of local stationarity is not required for our approach; that is, the forward and backward partial autocorrelations (defined in Section 2.2) need not be equal.Because this model generally allows both slow and rapid changes in the parameters, it can flexibly model the stochastic pattern changes often exhibited by complex nonstationary signals.

Lattice Structures
The Levinson-Durbin algorithm yields a unique correspondence between the PARCOR coefficients and the AR coefficients (Shumway and Stoffer, 2006;Kitagawa, 2010).Therefore, a disciplined approach to fitting AR models can be achieved through estimation of the PAR-COR coefficients.The lattice structure described below provides a direct way of associating the PARCOR coefficients with the observed time series (see Hayes (1996, Page 225) and the Supplementary Appendix for additional discussion).As such, the lattice structure provides an effective path to AR model estimation.In fact, the Levinson-Durbin algorithm for a stationary time series can be derived using the lattice structure (see Kitagawa, 2010, Appendix B).Then, the m-th stage of the lattice filter can be characterized by the pair of input-output relations between the forward and backward predictions, with the initial condition, f on the right hand side of ( 2) and (3) are the prediction errors of the forward and backward AR models, respectively, we consider an example when P = 2.In this case, we first derive that the difference between f (1) t and α (2) 2 b (1) t−2 of ( 2) is equal to the forward prediction errors of an AR(2) model as follows The above derivation shows that the difference between f (1) t and α (2) 2 b (1) t−2 is equal to the prediction errors of the forward AR(2).It can also be shown that (3) is true.Moreover, when the signal is stationary, the forward and backward PARCOR coefficients are equal (i.e., α In such cases, we can change the second term of the third row to (α The PARCOR coefficients α (4) This equation implies that once the PARCOR coefficient α (P ) P is estimated, then all of the other coefficients are immediately determined as well.

The Lattice Structure of the TVAR model
Given the assumption of second-order stationarity, the forward and backward PARCOR coefficients are constant over time (i.e., shift-invariant).However, since most real-world signals are nonstationary, the shift-invariant PARCOR coefficients are typically inappropriate.In such cases we can modify ( 2) and (3) as follows with both the forward and backward PARCOR coefficients α t,m and β t,m may not be identical at each time t.Therefore, our approach will proceed without this constraint.Also, the residual terms, f t , are assumed to follow zero-mean Gaussian distributions, N(0, σ 2 f,m,t ) and N(0, σ 2 b,m,t ), respectively.Importantly, when the true process is TVAR(P ), the variance σ 2 f,P,t is equal to the innovation variance σ 2 t .To verify this statement, we use the fact that PARCOR coefficient of AR(P ) is equal to zero when the lag is larger than P .This property can be applied to TVAR models since TVAR models correspond to AR models at each time t.Using such property, α for m ≥ P .The forward prediction error f (P ) t is identical to ǫ t ; i.e., f (P ) t and ǫ t are identically distributed.Therefore, as mentioned in Section 1, the Gaussian distribution provides a more reasonable assumption for the target model (1) than the Cauchy distributions used by Kitagawa (1988) and Kitagawa and Gersch (1996).
For each stage m of the lattice structure, we construct the following equations to obtain the coefficients, a t,k , of the forward and backward TVAR models (Hayes, 1996;Kitagawa and Gersch, 1996;Haykin, 2002) t,m .Equations ( 7) and ( 8) describe the relationship between the coefficients of the forward and backward TVAR models.In particular, these relations illustrate that the forward coefficients at the current stage are a linear combination of the forward and backward coefficients of the previous stage, with the weights equal to the PAR-COR coefficients.Importantly, such a combination also includes the stationary and locally stationary cases.For the stationary case, since α t,m are constant over time, the general equations ( 7) and ( 8) can be reduced to (4).For locally stationary cases, since t,m at time t, (7) and ( 8) are identical.

Model Specification and Bayesian Inference
Since both the forward and backward PARCOR coefficients of ( 5) and ( 6) as well as the corresponding innovation variances require time-varying structures, we consider random walk models for their evolutions.In such cases, the following two hierarchical components are added to ( 5) and ( 6).The evolution of the forward and backward PARCOR coefficients are modeled, respectively, as follows: where w α,m,t and w β,m,t are time dependent system variances.These system variances are then defined in terms of hyperparameters γ f,m and γ b,m , so called discount factors with range (0, 1), respectively (West and Harrison, 1997); see the Appendix for further details.Usually, we treat γ f,m = γ b,m = γ m at each stage m.These two equations also imply a sequential update form that the PARCOR coefficient at time t + 1 is equal to the sum of the PARCOR coefficient at time t plus a correction.
Similarly, both the evolution innovation variances, σ 2 f,m,t and σ 2 b,m,t , are modeled through multiplicative random walks as follows where δ f,m and δ b,m are hyperparameters (i.e., discount factors on the range (0,1)), and the multiplicative innovations, η f,t,m and η b,t,m follow beta distributions with parameters, (g f,m,t , h f,m,t ) and (g b,m,t , h b,m,t ), respectively (West et al., 1999).These parameters are defined at each time t by the discount factors, δ f,m and δ b,m , as detailed in the Appendix.In many cases, we also assume that δ f,m = δ b,m = δ m at each stage m.The series of stochastic error terms ǫ α,m,t , ǫ β,m,t , η f,m,t , and η b,m,t are mutually independent, and independent of the forward and backward innovations, f (m) t of ( 5) and ( 6).
We specify conjugate initial priors for α (m) 0,m and σ 2 f,m,0 at each stage m as follows p(α where D f,m,0 denotes the information set at the initial time t = 0, G(•, •) is the gamma distribution, µ f,m,0 and c f,m,0 are the mean and variance for a normal distribution, and v f,m,0 /2 and κ f,m,0 /2 are the shape and scale parameters for a gamma distribution.Usually, we treat the starting values µ f,m,0 , c f,m,0 , v f,m,0 , and κ f,m,0 as common constants over stage m.Typically, we choose µ f,m,0 and c f,m,0 to be zero and one, respectively.In addition, to set v f,m,0 and v f,m,0 , we first fix v f,m,0 = 1 and calculate the sample variance of the initial components of the signal.Given these two values, we can obtain κ f,m,0 through the formula for the expectation of the gamma distribution.In such prior settings, the DLM sequential filtering and smoothing algorithms provide the necessary components for the marginal posterior distributions (West and Harrison, 1997).Specifically, for t = 1, . . ., T , with the information set

Model Selection
Selection of the model order and set of discount factors {P, γ m , δ m ; m = 1, . . ., P } is essential for our approach.First, one can assume γ m = γ and δ m = δ, for m = 1, . . ., P .Then, the analysis will proceed using a set of various pre-specified combinations of (P, γ, δ).Since γ is related to the variability of the PARCOR coefficients, it also affects the variability of the TVAR coefficients.Hence, one can model the variance of the time-varying coefficients and the innovation variances of the TVAR models using discount factors γ and δ, respectively (West et al., 1999).Note that, in our context, the discount factors (γ m , δ m ) are a function of m (the lattice filter stage), whereas the West et al. (1999) setting does not make use of the lattice filter and, thus, there is only one set of discount factors (γ, δ) that need to be estimated.However, estimation of (P, γ, δ) using the approach of West et al. (1999) entails repeatedly having to calculate inverse matrices in the sequential filtering process.Our approach allows γ m and δ m to vary by stage.Thus, we first specify a potential maximum value of P and a set of combinations of {γ m , δ m } for each stage m.Given a value of P , we search for the combination of {γ 1 , δ 1 } maximizing the log likelihood of (5) at stage one.
Using the selected γ 1 and δ 1 , we can obtain the corresponding series {f (2) t } and {b (2) t }, for t = 1, . . ., T , as well as the value, L 1 , of log maximum likelihood of (5).We then, repeat the above search procedure for stage two using the output {f (2) t } and {b (2) t } obtained from implementing the selected hyperparameters γ 1 and δ 1 .In turn, this produces a new series of {f (3) t } and {b (3) t }, for t = 1, . . ., T , as well as a value L 2 .We repeat the procedure until the set of {γ m , δ m , L m }, m = 1, . . ., P , has been selected.
Here, we provide both a visual and numerical method to select the order.Similar to the scree plot widely used in multivariate analysis (Rencher, 2002), we can plot L m against the order m.When the observed series follows an AR or TVAR model, the values of L m will stop increasing after a specific lag, this lag can be chosen as the order for the estimated model.
Henceforth, this plot is referred to as "BLF-scree."This type of visual order determination can be directly quantified through the relative change of L m .Specifically, we provide a numerical method of order selection based on calculating the percent change in going from Based on simulation of various TVAR models, we choose τ = 0.5 with m − 1 reflecting the "best" value for the order.That is, we have found that 0.5 provides an effective cutoff for choosing the order.Although this approach provides a good guide to order selection, other model selection methods could be considered (e.g., shrinkage through Bayesian variable selection, reversible jump MCMC, or by minimizing an information criteria).Development of alternative model selection approaches in this setting constitutes an area of future research.

Simulation Studies
In this section, we simulate various nonstationary time series in order to compare the performance of our approach with four other approaches used to estimate the time-frequency representation.The first approach is AdaptSPEC proposed by Rosen et al. (2012).This approach adaptively segments the signal into finite pieces and then estimates the time-frequency representation using smoothing splines to fit local spectra via the Whittle likelihood approximation.The size of a segment and the number of the spline basis functions are two essential parameters for this approach.To reduce any subjectivity in our comparisons, we choose settings for these two parameters similar to those considered in Rosen et al. (2012) (with their t min = 40), as well as the same settings for MCMC iterations and burn-in.However, rather than using 10 spline basis functions we use 15, as this provides slightly better results along with superior computational stability.We note that the approach of Everitt et al. (2013) is not considered here due to the fact that the parameterization of the spectral density through the Wittle likelihood requires subjective knowledge.
The second approach is the AutoPARM method (Davis et al., 2006).Although this approach combines the GA and MDL to automatically search for potential break points along with the AR orders for each segment, four parameters are crucial for the GA: the number of islands, the number of chromosomes in each island, the number of generations for migration, and the number of chromosomes replaced in a migration; see Davis et al. (2006) for a comprehensive discussion.All of these parameters were chosen identical to those used in Davis et al. (2006).
The third approach is the AutoSLEX method (Ombao et al., 2001).Given a fixed value for the complexity penalty parameter of the cost function, AutoSLEX can automatically segment a given signal and choose a smoothing parameter.Following the suggestion of Ombao et al. (2001), we set this parameter equal to one.The last method we consider is the approach of West et al. (1999), referred to as WPK1999.This approach requires specification of three parameters: the TVAR order and two discount factors -one associated with the variance of the time-varying coefficients and the other with the innovation variances.In general, the discount factor values are in the range 0.9−0.999(West et al., 1999).Therefore, for our simulations, we give each discount factor a set of values from 0.8 to 1 (with equal spacing of 0.02) and, further, a set of values for the TVAR order from 1 to 15.Given these values, we choose the combination that achieves the maximum likelihood (West et al., 1999).
Our approach uses the two selection methods discussed in Section 2.5 to search for the TVAR order with appropriate discount factor values.The selected combination of (P, γ, δ) with γ and δ held fixed over all stages of the Bayesian lattice filter is referred to as BLFFix.
The selected combination of (P, γ m , δ m ), for m = 1, . . ., P , is referred to as BLFDyn.Again, the candidate space of parameters for both discount factors are from 0.8 to 1 (with equal spacing of 0.02), along with orders from 1 to 15.
We consider four types of nonstationary signals: 1) TVAR of order 2 with constant innovation variance; 2) TVAR of order 6 with constant innovation variance; 3) a piecewise AR process with constant innovation variance; and 4) simulated signals based on an Enchenopa treehopper communication signal (Holan et al., 2010); see Section 4.1.Each simulation consists of 200 realizations.To evaluate the performance in estimating the various timefrequency representations, we calculate the average squared error (ASE) for each realization as follows (Ombao et al., 2001) where n = 1, . . ., 200, ω l = 0, 0.005, . . ., 0.5, and T denotes the length of the simulated series.Lastly, we denote ASE = (1/200) 200 n=1 ASE n .For each simulation study, Table 1 summarizes the mean values and standard deviations for ASE n .
In the case of the AutoSLEX method the number of frequencies in (17) differs from the other approaches considered.In particular, the AutoSLEX approach dyadically segments the signal up to a given maximum scale J such that 2 J is less than signal length, T .Subsequently, AutoSLEX automatically determines whether a segment at a particular scale will be included in final segmentation.Once this has been completed, the frequency resolution for the AutoSLEX approach is equal to T /2 (j+1) where j is the scale of the largest segment included in the final segmentation.

Time-Varying AR(2) Process
We simulate signals from the same time-varying AR(2) process (TVAR2), used in Davis et al. (2006) and Rosen et al. (2009Rosen et al. ( , 2012)), which is defined as follows where ǫ t iid ∼ N(0, 1) and t = 1, . . ., 1024. Figure 1 shows the BLF-scree plot, suggesting that order two is the appropriate choice for all 200 realizations.Since the time-varying coefficient a t , varies slowly with time, this process naturally exhibits a slowly evolving time-varying spectrum (Figure 2).The box-plots of the ASE values in Figure 2 show that the group of TVAR-based models (i.e., WPK1999, BLFFix, and BLFDyn) perform superior to the group of non-TVAR-based models (i.e., AdaptSPEC, AutoPARM, and AutoSLEX), with BLFDyn performing the best.In particular, there is a significant percent reduction in ASE for BLFDyn relative to the other methods considered (Table 1).
The BLF-scree plot (not shown) suggests that order six is the appropriate choice for all 200 realizations.The TVAR(6) contains three pairs of time-varying conjugate complex roots.
Figure 3 illustrates that the TVAR(6) has a time-varying spectrum with three peaks.Similar to the TVAR(2) analysis (Section 3.1), the TVAR-based models outperform the group of non-TVAR-based models, with BLFDyn performing superior to the others.Again, the percent reduction in ASE is significant relative to the other approaches considered (Table 1).

Piecewise Stationary AR Process
The signals simulated here are based on the same piecewise stationary AR process, used by Davis et al. (2006) and Rosen et al. (2009Rosen et al. ( , 2012) ) and is defined as follows where ǫ t iid ∼ N(0, 1).These generated signals are referred to as PieceAR.Since it is difficult to choose the order for some realizations visually using the BLF-scree plot, we use ( 15), with τ = 0.5, to choose the order.The numerical method suggests order two for some realizations and order three for the others.The true process includes three segments, with each of the segments mutually independent.The piecewise nature of this process is clearly depicted by its time-varying spectrum (Figure 4b).The box-plot (Figure 4c) shows that AutoPARM exhibits superior performance in terms of the smallest median ASE.However, WPK1999 and BLFFix may perform more robustly (i.e., less outlying ASE values).Although we see a 23.78% reduction in ASE for AutoPARM versus BLFFix, we find that BLFFix performs superior to the remainder of the approaches and has the smallest standard deviation across the 200 simulations.Table 1 summarizes the mean values and standard deviations for this simulation.The findings here are not surprising as AutoPARM is ideally suited toward identifying and estimating piecewise AR processes.For our approach, taking (γ, δ) fixed is advantageous for processes that are not slowly-varying.

Simulated Insect Communication Signals
The signals considered in this simulation are formulated such that they exhibit the same properties as an Enchenopa treehopper mating signal; see Section 4.1 for a complete discussion.Specifically, we fit a TVAR(6) model to the signal x t , t = 1, . .  Figure 5b, we see that BLFFix and BLFDyn perform better than the other approaches, in terms of median ASE.Further, we find that BLFFix and BLFDyn are similar, in terms of ASE, although the median of BLFDyn (0.2675) is smaller than that of BLFFix (0.3176).
Table 1 summarizes the mean values and standard deviations for this simulation.From this table we see that the reduction in ASE for BLFDyn is 6.01% over BLFFix and that both exhibit a substantial percent reduction in ASE over the other methods considered.

Animal Communication Signals
Understanding the dynamics of populations is an important component of evolutionary biology.Many organisms exhibit complex characteristics that intricately relate to fitness.For example, the mating signal of the Enchenopa treehopper represents a phenotype of the insect that is used in mate selection.During the mating season, males in competition deliver their vibrational signals through stems of plants to females (see Cocroft and McNett, 2006, and the references therein).The data considered here comes from an experiment that was previously analyzed in Cocroft and McNett (2006) and Holan et al. (2010).The experiment was designed with the goal of reducing potential confounding effects between environmental and phenotypical variation.In this experiment, males signals were recorded one week prior to the start of mating.Figure 6a displays a typical signal of from a successful mater, with length 4,739 downsampled from registered signals of length 37,912.Justification for the appropriateness of downsampling the original signal, in this context, can be found in Holan et al. (2010).Also, as discussed in Holan et al. (2010), this signal shows a series of broadband clicks preceding a frequency-modulated sinusoidal component, followed by a series of pulses.
For this analysis, we used the BLFDyn approach to search for a model having both discount factors in the range of 0.8 to 1 (with equal spacing of 0.02) and an order between 1 and 25. Figure 7a shows an increase in L m along with the order, which is different from the simulated TVAR(6) model in Section 3.4.Hence, we use (15) with the model order chosen by τ = 0.5.This rule yields a TVAR(6) model.In Figure 7b, the PARCOR coefficients of lag larger than two are close to zero following time around 0.3 (where the time axis has been normalized such that t ∈ (0, 1)).Thus, the last four TVAR coefficients after time 0.3 are close to zero (Figure 7c).Such phenomena suggests that the period before time 0.3 has a more complex dependence structure.Figure 7d illustrates that the innovation variance exhibits higher volatility at the beginning signal.These bursts in the innovation variance are related to the series of broadband clicks at the beginning of the signal.Finally, Figure 6b presents the time-frequency representation of the treehopper signal using the TVAR(6) model and corroborates the significance of the broadband clicks at the beginning of the signal.Figure 6c shows the posterior standard deviation of the time-frequency representation using 2,000 MCMC samples from 210,000 MCMC iterations by discarding the first 10,000 iterations as burn-in and keeping every 100-th iteration of the remainder.

Wind Components
We study the time-frequency representations of the east/west and north/south wind components, recorded daily at Chuuk Island in the tropical Pacific during the period of 1964 to 1994 (see Cressie and Wikle, 2011, Sections 3.5.3 and 3.5.4).The data studied are at the level of 70 hPa, which is important scientifically due to the likely presence of westward and eastward propagating tropical waves, and the presence of the quasi-biennial oscillation (QBO) (Wikle et al., 1997).Figure 8 shows the two wind component time series from which we can discern visually that the east/west series clearly exhibits the QBO signal, but no discernible smaller-scale oscillations are present in either series.Our interest is then whether the time-varying spectra for these series suggest the presence of time-varying oscillations, which are theorized to be present.
We consider the same search space for the model order and discount factors as that used for the treehopper communication signal (Section 4.1).The BLF-scree plot for the east/west component shows an increase of L m along with the order.Therefore, we use ( 15) with τ = 0.5 and choose the order equal to three.The PARCOR coefficient of lag one depends on time but the PARCOR coefficient of lag two and three appear to be constant.The innovation variances of TVAR(3) model are time dependent.On the other hand, the BLF-scree plot for the north/south component shows a turning point at order four so that we choose the order equal to four.The PARCOR coefficients are time independent but the innovation variances of the TVAR(4) model are time dependent.The preceding findings are illustrated in Figures 13 and 14 of the Appendix.Figures 9a and 9b show the estimated time-frequency representations of the wind east/west and north/south components.Figures 9c and 9d present the associated posterior standard deviations of the time-frequency representation using 2,000 MCMC samplers from 210,000 MCMC iterations by discarding the first 10,000 iterations as burn-in and keeping every 100-th iteration of the remainder.The east/west component time-varying spectrum does suggest that the QBO intensity varies considerably as evidenced by the power in the low-frequencies.Perhaps more interesting is the suggestion of time-varying equatorial waves in the north/south wind component time-varying spectrum.
In particular, the lower-frequency (Kelvin and Rossby) waves with frequencies between 0.1 and 0.2 show considerable variation in duration of wave activity, as well as intensity.One also sees time-variation in the likely mixed-Rossby gravity waves in the frequency band between 0.2 and 0.35.Interestingly, in some cases these are in phase with the lower-frequency wave activity but more often act in opposition.We also note the almost complete collapse of the equatorial wave activity centered on 1984.

Economic Index
Koopman and Wong ( 2011) study the time-frequency representation of the log difference of the US gross domestic product (GDP), consumption, and investment from the first quarter of 1947 to the first quarter of 2010.For comparison, we obtained the series from the Federal Reserve Bank of St. Louis (http://research.stlouisfed.org/fred2).Since the investment series is unavailable from the website, we only consider the GDP and consumption in our study.
Figure 10a illustrates the trend of the logarithm of the GDP and consumption series.Similar to Koopman and Wong (2011), we analyzed the log difference of these two series (Figures 10b   and 10c).
We consider the same search space for the model order and discount factors as that used for the treehopper communication signal (Section 4.1).By using (15) with τ = 0.5, we choose TVAR(1) for the log difference of GDP and TVAR(2) for the log difference of consumption.

Discussion
This paper develops a computationally efficient method for model-based time-frequency analysis.Specifically, we consider a fully Bayesian lattice filter approach to estimating timevarying autoregressions.By taking advantage of the partial autocorrelation domain, our approach is extremely stable.That is, the PARCOR coefficients and the TVAR innovation variances are specified within the lattice structure and then estimated simultaneously.Notably, the full conditional distributions arising from our approach are all of standard form and, thus, facilitate easy estimation.
The framework we propose extends the current model-based approaches to time-frequency analysis and, in most cases, provides superior performance, as measured by the average squared error between the true and estimated time-varying spectral density.In fact, for slowly-varying processes we have demonstrated significant estimation improvements from using our approach.In contrast, when the true process comes from a piecewise AR model the approach of Davis et al. (2006) performed best, with our approach a close competitor and performing second best.This is not unexpected as the AutoPARM method is a model-based segmented approach and more closely mimics the behavior of a piecewise AR.

Acknowledgments
demonstrated that features in the time-frequency domain of nonstationary Enchenopa treehopper mating signals may describe crucial phenotypes of sexual selection.
error at time t for a forward and backward AR(

m
are the lag m forward and backward PARCOR coefficients, respectively.Equation (2) shows that the forward PAR-COR coefficient of lag m is a regression coefficient of the forward prediction error f (m−1) t regressed on the backward prediction error b the forward AR(m) model.Similarly, (3) shows that the backward PAR-COR coefficient of lag m is a regression coefficient of the backward prediction error b the backward AR(m) model.Using (2) and (3) recursively, we can derive the PARCOR coefficients for a given lag.In the stationary case, the forward and backward PARCOR coefficients are equivalent; i.e., α

P
the last component of the coefficients of the forward AR(P ) model; i.e., α the Levinson-Durbin algorithm, the remainder of the AR coefficients and the innovation variance can be obtained as follows a (P ) m = a (P −1) −m , m = 1, . . ., P − 1.
now time dependent.Note that for notational simplicity, f (m−1) t and b (m−1) t here denote the prediction error at time t of the forward and backward TVAR(m − 1).For locally stationary signals, we may impose the constraint that α at each time t.However, for general nonstationary cases, α

D
f,m,T up to time T , the marginal posterior distributions p(α (m) t,m |D f,m,T ) and p(σ −2 f,m,t |D f,m,T ) are the t-distribution and gamma distribution, respectively.Analogous to α (m) 0,m and σ 2 f,m,0 , the same conjugate initial priors for β (m) 0,m and σ 2 b,m,0 are specified at each stage m.Details of the sequential filtering and smoothing for the PARCOR coefficients and innovation variances for each stage m are discussed in the Appendix.
}, { β (P ) t,P }, { σ 2 f,P,t }, and { σ 2 b,P,t } have been obtained.Then, recursively plugging the estimates of {α . , 4096, to obtain time-varying AR coefficients and innovation variances.Typically, with these type of nonstationary signals, the innovation variances are time dependent, which is markedly different from the previous examples where the innovation variance was constant.The signals generated by these parameters are referred to as SimBugs.As expected, the BLF-scree plot (not shown) suggests that order six may be an appropriate choice for all 200 realizations.

Figure
Figure 5a illustrates one realization of the SimBugs, whereas Figure 5b provides box-plots that characterize the distribution of ASE n over the 200 simulated signals.Specifically, from

Figures
Figures 11a and 11bshow the estimated time-frequency representations and the associated standard deviations of the log difference of GDP and consumption series, respectively.Similar to the results ofKoopman and Wong (2011), both the US macroeconomic indices exhibit This research was partially supported by the U.S. National Science Foundation (NSF) and the U.S. Census Bureau under NSF grant SES-1132031, funded through the NSF-Census Research Network (NCRN) program.The authors would like to thank the Reginald Cocroft lab for use of the insect communication data analyzed in Section 4.1.The authors would like to thank Drs.Richard Davis, Hernando C. Ombao, and Ori Rosen for providing their codes.Finally, we also thank Dr. Mike West for generously providing his code online.

Figure 7 :
Figure 7: (a) shows the BLF-scree plot of the treehopper communication signal.(b) depicts the first six time-varying estimated PARCOR coefficients.(c) and (d) show the estimated time-varying coefficients and innovation variances of the TVAR(6) model.

Figure 8 :
Figure 8: (a) and (b) show daily time series of east/west and north/south components of wind, respectively.Both components are measured in meters per second (m/s).

FrequencyFigure 9 :Figure 10 :Figure 11 :
Figure 9: (a) and (c) display the posterior mean and standard deviation of time-frequency representations of the wind east/west component by fitting a TVAR(3) model.(b) and (d) display the posterior mean and standard deviation of time-frequency representations of the wind north/south component by fitting a TVAR(4) model.

Figure 13 :
Figure 13: (a) shows the BLF-scree plot of the east/west component of the wind signal.(b) depicts the first three time-varying estimated PARCOR coefficients.(c) and (d) show the estimated time-varying coefficients and innovation variances of the TVAR(3) model.

Figure 14 :
Figure 14: (a) shows the BLF-scree plot of the north/south component of the wind signal.(b) depicts the first four time-varying estimated PARCOR coefficients.(c) and (d) show the estimated time-varying coefficients and innovation variances of the TVAR(4).

Figure 15 :Figure 16 :
Figure 15: (a) shows the BLF-scree plot of the difference of logarithm GDP series.(b) depicts the first time-varying estimated PARCOR coefficients.(c) and (d) show the estimated timevarying coefficients and innovation variances of the time-varying AR(1).

Table 1 :
The mean, ASE and standard deviation, sd ASE , of the ASE values for the simulations presented in Section 3. Note that the bold values represent the approach having minimum ASE.