Abstract
We consider a transfer function model with time-varying coefficients. We propose an estimation procedure, based on the least squares method and wavelet expansions of the time-varying coefficients. We discuss some statistical properties of the estimators and assess the validity of the methodology through a simulation study. We also present an application of the proposed procedure to a real pair of series.
1. Introduction
In this paper, we consider the transfer function model where is the number of observations and is i.i.d. . We assume that the error and the input series are independent. The functions , , and , , are supported on the interval and connected to the underlying series by an appropriate rescaling of time, .
As noted by Nason and Sapatinas [1], when both time series and belong to the class of ARMA models, the transfer function models made popular by Box and Jenkins are appropriate. See Box et al. [2].
Model (1.1) is important because it can be used when either or both time series are not stationary, but it is perhaps more suitable when the series exhibit a locally stationary behaviour, for example, in the sense of Dahlhaus [3] or of Nason et al. [4]. In several areas, such as signal analysis (speech analysis in particular) and economics, empirical and theoretical works suggest that models with time-varying coefficients are useful, reflecting, for example, the evolution of the economy over time.
Chiann and Morettin [5, 6] considered a special case of (1.1), namely, while Dahlhaus et al. [7] entertained the model
These authors used wavelet expansions of the time-varying coefficients in order to obtain estimates.
Related references are Brockwell et al. [8, 9], who used a state-space approach to transfer function modelling, allowing for nonstationary input and output sequences, and Nason and Sapatinas [1] who show how a nondecimated wavelet packet transform (NDWPT) can be used to model a response time series in terms of an explanatory time series , both assumed to be nonstationary. Rather than building a model directly between and some regression-type model is built between and an NDWPT version of .
We consider the problem of estimating , and , , in time domain, using wavelet expansions. We use least squares to obtain the estimators of the wavelet coefficients. Then the empirical detail coefficients are shrunk before the inverse wavelet transform is applied to obtain the final estimates of and . This results in a nonlinear smoothing procedure. See Section 2 for further details.
We will apply the methodology to real data, namely, a pair of hourly wind speeds recorded at two Welsh Meteorological Office stations, Valley and Aberporth. Here the idea is to relate the target wind speeds (Valley) to wind speeds measured at a nearby reference site (Aberporth). See Nason and Sapatinas [1] for a complete description of this data and Section 6.
The plan of the paper is as follows. In Section 2, we introduce some background material on wavelets. In Section 3, we present the estimation procedure, and in Section 4, we give some properties of the empirical wavelet coefficients. In Section 5, we perform some simulations and in Section 6, we apply the methodology to the data described above. The paper ends with some remarks in Section 7.
2. Wavelets
In this section, we discuss some basic ideas on wavelets. For more details, see the books by Vidakovic [10], Percival and Walden [11], and Nason [12]. From two basic functions, the scaling function and the wavelet , we define infinite collections of translated and scaled versions, , , . We assume that forms an orthonormal basis of , for some coarse scale . A key point [13] is that it is possible to construct compactly supported and that generate an orthonormal system and have space-frequency localization, which allows parsimonious representations for wide classes of functions in wavelet series.
In some applications, the functions involved are defined in a compact interval, such as . This will be the case of our functions and in (1.1). So it will be necessary to consider an orthonormal system that spans . Several solutions were proposed, the most satisfactory one being that by Cohen et al. [14]. As in any proposal, close to the boundaries the basis functions are modified to preserve orthogonality. Cohen et al. [14] use a preconditioned step. Periodized wavelets are defined by and these generate a multiresolution level ladder , in which the spaces are generated by the . Negative values of are not necessary, since and if , . See Vidakovic [10] for details and Walter and Cai [15] for a different approach. Other instances are the boundary corrected Meyer wavelets [16]. From now on, we denote the periodized wavelets simply by .
Accordingly, for any function , we can expand it in an orthogonal series where we take and . For each , the set brings the values , such that belongs to scale . For example, for , we have 8 wavelet coefficients in scale 23, while for we have half of this number, namely, 4 coefficients in scale . The wavelet coefficients are given by
Often we consider the sum in (2.2) for a maximum level , in such a way that we approximate in the space.
Wavelet shrinkage provides a simple tool for nonparametric function estimation. It is an active research area where the methodology is based on optimal shrinkage estimators for the location parameters. Some references are Donoho and Johnstone [17, 18] and Donoho et al. [19]. In this paper, we focus on the simplest, yet most important shrinkage strategy-wavelet thresholding. Earlier models considered nonparametric regression with i.i.d. normal errors. These were further extended by Johnstone and Silverman [20] for correlated errors. In this context, the paper by Opsomer at al. [21] is a good survey of methods using kernels, splines, and wavelets.
The thresholding technique consists of reducing the noise included in a signal through the application of a threshold to the estimated wavelets coefficients . Some commonly used forms are the soft and hard thresholds, given by respectively. One crucial issue is the choice of the parameter . If the empirical wavelet coefficients are assumed to be Gaussian with mean and constant variance, a popular method is to apply the so-called universal threshold , where is the noise variance. Other methods include the SURE procedure, cross validation, and Bayesian approaches. See Percival and Walden [11] and Vidakovic [10] for details.
But despite the good properties of these procedures, which are spatially adaptive and nearly optimum over wide classes of function spaces, they are not appropriate when the design is not regular or it is random. Recently, several approaches were proposed dealing with nonregular designs. Delouille [22] introduced the concept of design adapted wavelets, and Kerkyacharian and Picard [23] introduced the warped wavelets. Cai and Brown [24, 25] consider wavelet estimation for models with nonequispaced or random uniform designs, and Sweldens [26] introduced the so-called lifting scheme; now the wavelets are no longer obtained as dilations and translations of some mother wavelet. See also Jansen and Oonincx [27]. Some applications and extensions of these approaches are given in Porto et al. [28, 29] and Morettin and Chiann [30].
In this paper, we will use ordinary wavelets, as in Dahlhaus et al. [7], since these wavelets work well for our purposes. See also a related discussion in Morettin and Chiann [30]. Concerning the choice of the mother wavelet, there are no established rules. In our experience, we agree with Nason and Sapatinas [1], when they say that “how the choice of of wavelets affects the final model and its interpretation is an area of further research.” We will make use in particular of the Daubechies compactly supported wavelets in the simulations and application. Concerning thresholding, we use hard thresholds and the parameter chosen as explained in Sections 4 and 5.
3. Estimators
The aim is to estimate the functions , , and , , appearing in model (1.1), given observations of the series and . We assume that the orders and are fixed and known. The idea is to expand these functions in wavelet series
The empirical wavelet coefficients are obtained minimizing with and replaced by (3.1), . These empirical wavelet coefficients are then modified using a hard threshold, and finally we build estimators of and by applying the inverse wavelet transform to these thresholded coefficients.
For ease of exposition, we restrict our attention from now on to the simple model with and , namely,
So we regress on and , using the expansion (3.1) for and . Let . In matrix notation, we have
or where
The vector is given by and the matrices are given, for , by
The vector is given by and the matrices are given, for , by
Let and .
It follows easily that the least squares estimators of the coefficients are then given by
Having obtained the estimates given by (3.11), we plug them in (3.1), resulting in linear estimates and . Finally, nonlinear smoothed estimators and , respectively, are obtained applying some threshold to the empirical detail coefficients and taking the inverse wavelet transform.
We might think of using other procedures, as the LASSO (least absolute shrinkage and selection operator) procedure of Tibshirani [31] or the nonnegative garotte of Breiman [32]. For the case of orthogonal design, LASSO has a relationship with the soft shrinkage procedure of Donoho and coworkers. But using these procedures would make the derivations of the properties of the estimators above considerably more difficult, or even impossible.
4. Properties of Empirical Coefficients
Now we present some properties of the linear empirical wavelet coefficients. The techniques used to prove the results are quite involved and are based on function space theory. Basically we adapt the results of Dahlhaus et al. [7] for the transfer function model (1.1). The proofs may be obtained from the authors upon request or from the site of the journal.
We assume that functions , and , belong to some function spaces given by where
. Here, is the smoothness degree of each space, and specify the norm, and and are positive constants. For these function spaces, the following result is valid [19]: where , with .
These classes contain Besov, Hölder, and Sobolev spaces, see, for example, Vidakovic [10], Dahlhaus et al. [7], and Triebel [33].
It can be shown, see Donoho et al. [19], that the loss incurred by truncating at level is of order if we choose such that . This is achieved if .
Concerning the wavelets, we assume that and are compactly supported on and have continuous derivatives up to order , with . We denote the spectral norm by and the sup norm by .
In order to analyze the statistical behavior of the estimated coefficients, it is convenient to assume that we take expansions of these coefficients as linear combinations of functions in the spaces, generated by . With this basis, we can write (3.3) as where
Equation (4.4) in matrix form becomes with
In (4.6), and contain the coefficients in the expansion (4.4) and contains lines given by (4.5).
Analogously to (3.11), we have the least squares estimates
The relationship between and is where the is a () block diagonal matrix. The matrix does the transformation ; so each coefficient estimate , in (3.11) can be written as and .
Equation (4.6) can be written in the form with , and the vector as defined above. The least squares estimator of is given by
We notice that the error term is not independent from the regressors , which means that the estimator is biased and its squared error is of order . But this order can be improved, as shown below. To take care of the bias, we might use some robust procedure, as in Martin et al. [34], for example.
From (4.10) and (4.11), we obtain If we write , we get where the matrix is an () matrix, given by
We now state results on the squared error and mean squared error of the linear empirical wavelet coefficients. The following assumptions are needed.(A1) The wavelets and are compactly supported on and have continuous derivatives up to order , with .(A2) In the estimation procedure, we have used first a linear estimator, truncating the wavelet expansion at scale . In order to get functions with the appropriate smoothness, was chosen such that .(A3) The matrix in (4.10) satisfies , for some .
The Daubechies Extremal Phase wavelets , for an appropriate choice of , satisfy (A1). See Härdle et al. [35] for details. The number of continuous derivatives has a connection with the number of null moments of . Also, the support of is contained in the interval .
Since is the coefficient of the regressor , the squared error of the least squares estimator is , which implies that the square error , for example. We can improve this rate as in the next result.
Proposition 4.1. Under Assumptions (A1)–(A3), we have that .
This proposition allows us to derive the rate of convergence of the mean square error of the original linear estimators of the wavelet coefficients.
Proposition 4.2. If Assumptions (A1)–(A3) hold, then holds uniformly in , and .
5. Some Simulations
We now present three simulation examples for the linear and nonlinear estimates , and , , respectively. In the first, the coefficients are very smooth functions, the second presents an intermediary case, and, in the third, the functions are piecewise constants.
The performance of the estimators are assessed via the square root of the mean squared errors (RMSE), namely, the same for the nonlinear estimators. Recall that .
Concerning the choice of the wavelet, we have tried the Daubechies Extremal Phase D8 and D12, the Daubechies Least Asymmetric S8 and S12 and the Haar wavalets. The best results were obtained with the wavelet mentioned in each example. Of course, in a real problem situation, this cannot be done, since we do not know the true function coefficients. In our experience, the Haar and D8 wavelets have proven useful.
Other important issues are the choice of , how to estimate , the noise variance, and the choice of the threshold parameter . Concerning , we have used the rule specified in Section 4, below (4.3), resulting for all examples. The variance was estimated using the empirical wavelet coefficients in all scales. Hard thresholding was used with the parameter chosen as .
In all cases, the input time series was generated as an AR(2) model, , where are independent, normally distributed random variables, with zero mean and variance one, for , and the coefficients are given by
Example 5.1. Let and , so the transfer function model is
We choose , , and as , , and .
For the simulation, we generated data values for and with 5000 replicates. For this example, the Daubechies D12 wavelet is employed for our estimation procedure.
Figure 1 presents the input and output series, and Figure 2 presents the true and estimated coefficients from a typical sample. The boxplots for the 5000 RMSE values are presented in Figure 3. We see that both types of estimators perform well, but clearly the nonlinear estimator perform, better, both visually and in RMSE sense.
(a)
(b)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
In Figure 4, we use bootstrap with replications to see that the true transfer function coefficients lie within the confidence regions. See Section 6 for the description of the bootstrap procedure.
(a)
(b)
(c)
(d)
(e)
(f)
In Figures 17 and 18 of the appendix, we show histograms for linear and nonlinear estimators of the parameter at some points of the interval . Notice that the normal approximation for the estimators is a possibility to be considered under the theoretical point of view.
Example 5.2. Let and ; so the transfer function model is
We choose and as , ; , for and , and , for .
As in the Example 5.1, we generated data values for and with 5000 replicates and the Daubechies Extremal Phase wavelet D8 is employed.
From a typical sample, the input and output series are shown in Figure 5 and the estimated coefficients in Figure 6. The boxplots for the 5000 RMSE values are presented in Figure 7. We see that, visually, both types of estimators perform well, but the thresholded estimator is better for the smoother parameter in terms of RMSE. This is a well-known fact, see, for example, Percival and Walden [11]. In Figure 8, we use bootstrap with replications to see that the true transfer function coefficients lie within the confidence regions.
(a)
(b)
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
In Figures 6(c) and 6(e), we see some unwanted artifacts next to discontinuity points. These do not disappear, even with expansions with many terms. But this Gibbs effect is less severe than in Fourier expansions. See Vidakovic [10] and Percival and Walden [36] for details.
In Figures 19 and 20 of the appendix, we show histograms for linear and nonlinear estimators of the parameter at some points of the interval .
Example 5.3. Let Let and ; so the transfer function model is as in Example 5.2.
We choose the coefficients and as piecewise constant functions: , for and , for ; , for or , and , for or .
We simulated data values for and with 5000 replicates, and, due to the nature of the functions, the Haar wavelet is appropriate in this situation.
Figure 9 presents the input and output series, and Figure 10 presents the true and estimated coefficients from a typical sample. The boxplots for the 5000 RMSE values are presented in Figure 11. We see that both types of estimators perform well, but clearly the nonlinear estimator perform better, both visually and in RMSE sense.
(a)
(b)
(a)
(b)
(c)
(d)
(e)
(f)
In Figure 12, we use bootstrap with replications to see that the true transfer function coefficients lie within the 95% confidence regions.
(a)
(b)
(c)
(d)
In Figures 21 and 22 of the appendix, we present histograms for linear and nonlinear estimators of at some points of the interval .
6. An Application
The data that we analyze are hourly wind speeds recorded from 00:00 on 6th April 1995 at two Welsh Meteorological Office stations: Valley and Aberporth, with 2048 observations. Valley is located approximately 120 km north of Aberporth, and they are mostly separated by Cardigan Bay. Part of these data were analysed by Nason and Sapatinas [1]. Figure 13 shows the data we analyze. We are grateful to N&M Wind Power for providing the data.
(a) Days since 00:00 on 6th April 1995
(b) Days since 00:00 on 6th April 1995
In this application, our aim is to model Valley’s wind speeds, , in terms of its lagged values and present and lagged values of wind speeds at Aberporth, . We should remember that both series are not stationary and so classical methods should not be used. We fitted four potential models, namely,
The Daubechies Extremal Phase D8 wavelet was used in the computations. The model was marginally superior (in terms of RMSE) to model , so it was discarded. The mean residual sums of squares for , , and are 0.02258, 0.02186, and 0.02177, respectively, for linear estimation. For nonlinear estimation, the figures are 0.02391, 0.02359, and 0.02178, respectively. The final coefficient estimates for model are shown in Figure 14.
(a)
(b)
(c)
For Model , Figures 15 and 16 show the confidence regions based on bootstrap, using linear and nonlinear estimates, respectively.
(a)
(b)
(c)
(a)
(b)
(c)
The bootstrap procedure works as follows:(a)using the real series, we obtain the estimated transfer function coefficients;(b)we obtain the fitted series using these estimated coefficients;(c)we obtain the residuals from the fitted model;(d)we resample these residuals and obtain bootstrap samples of the series;(e)we obtain bootstrap estimators of the transfer function coefficients.
We used bootstrap replications of , , , and , , , to obtain bootstrap standard errors and the 95 confidence regions. We also fitted the model with constant coefficients given an MSE value of 0.02354, indicating that the model with varying coefficients is better. We see that is not possible to plot constant lines totally within these regions, showing that, in fact, the coefficients are time varying.
7. Further Remarks
In this paper, we have proposed an estimation procedure for a transfer function model with time-varying coefficients. Basically it is a least squares procedure, with the use of wavelets to expand the function coefficients. Firstly, linear estimators for the time-varying coefficients are obtained, truncating the wavelet expansion at an appropriate scale. Then thresholds are applied to the empirical wavelet coefficients, and inverse wavelet transformation gives the nonlinear smoothed estimators for the function coefficients.
Simulations have shown that this procedure leads to estimators with a good performance. Some statistical properties for the empirical wavelet coefficients in the case of linear estimators were presented.
Further studies are needed on the estimation of the variance and asymptotic normality of the empirical wavelet coefficients in the linear case, on the rate of convergence for the risk of the nonlinear estimator over the smoothness classes and on issues related to the identification, diagnostics, and forecasting for these nonstationary models.
The issue of forecasting was discussed in the literature for some special and related models. Fryzlewicz et al. [37] addressed the problem of how to forecast locally stationary wavelet processes in the sense of Nason et al. [4] by means of nondecimated wavelets. Van Bellegem and von Sachs [38] consider a special modulated locally stationary process, where a stationary process is modulated by a time-varying variance, so the problem reduces to forecast this variance.
Concerning asymptotic distribution of the empirical wavelet estimators, in the linear case, we believe that this can be obtained under further assumptions, namely on the cumulants of the , on the thresholds used and on the roots of the polynomials associated with the model (1.1). This will be the subject of further research.
Appendix
For more details, see Figures 17, 18, 19, 20, 21, and 22.
Acknowledgments
The authors are grateful to two anonymous referees, whose comments improved substantially a previous version of the paper. The first author acknowledges the support of CNPq, and the last three authors acknowledge the suport of Fapesp (Grant 2008/51097-6) and Capes/Procad (Grant 177/2007).