Transfer Function Models with Time-varying Coefficients

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.


Introduction
In this paper, we consider the transfer function model where T is the number of observations and ε t is i.i.d.0, σ 2 .We assume that the error and the input series are independent.The functions ω j u , j 0, 1, 2, . . ., n, and δ j u , j 1, 2, . . ., m, are supported on the interval 0,1 and connected to the underlying series by an appropriate rescaling of time, u t/T .As noted by Nason and Sapatinas 1 , when both time series X t and Y t 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  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 Y t in terms of an explanatory time series X t , both assumed to be nonstationary.Rather than building a model directly between X t and Y t some regression-type model is built between Y t and an NDWPT version of X t .
We consider the problem of estimating ω j u , j 0, 1, 2, . . ., n and δ j u , j 1, 2,. .., m, 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 ω j u and δ j u .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.

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 φ x and the wavelet ψ x , we define infinite collections of translated and scaled versions, φ j,k x 2 j/2 φ 2 j x − k , ψ j,k x 2 j/2 ψ 2 j x − k , j, k ∈ Z {0, ±1, ±2, . ..}.We assume that {φ ,k • } k∈Z ∪ {ψ j,k • } j≥ ;k∈Z forms an orthonormal basis of L 2 R , 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 0, 1 .This will be the case of our functions ω j u and δ j u in 1.1 .So it will be necessary to consider an orthonormal system that spans L 2 0, 1 .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 V j are generated by the ψ j,k .Negative values of j are not necessary, since φ φ 0,0 1 and if j ≤ 0, ψ j,k x 2 −j/2 .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 ψ j,k .
Accordingly, for any function f ∈ L 2 0, 1 , we can expand it in an orthogonal series where we take 0 and I j {k : k 0, . . ., 2 j − 1}.For each j, the set I j brings the values k, such that β jk belongs to scale 2 j .For example, for j 3, we have 8 wavelet coefficients in scale 2 3 , while for j 2 we have half of this number, namely, 4 coefficients in scale 2 2 .The wavelet coefficients are given by α 0,0 f x φ x dx, β j,k f x ψ j,k x dx.

2.3
Often we consider the sum in 2.2 for a maximum level J, in such a way that we approximate f in the V J 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 strategywavelet 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 β jk .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 β jk are assumed to be Gaussian with mean β jk and constant variance, a popular method is to apply the so-called universal threshold λ σ 2 log 2 J , where σ 2 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.

Estimators
The aim is to estimate the functions δ i u , i 1, . . ., m, and ω i u , i 0, 1, . . ., n, appearing in model 1.1 , given T observations of the series X t,T and Y t,T .We assume that the orders m and n are fixed and known.The idea is to expand these functions in wavelet series

3.1
The empirical wavelet coefficients are obtained minimizing with δ i u and ω i u replaced by 3.1 , v max m, n .These empirical wavelet coefficients are then modified using a hard threshold, and finally we build estimators of δ i u and ω i u 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 m 1 and n 0, namely, So we regress Y t,T on X t,T and Y t−1,T , using the expansion 3.1 for δ 1 t/T and ω 0 t/T .Let Δ 2 J − 1.In matrix notation, we have β δ 1 00 . . .where The vector Φ Y is given by The vector Φ X is given by . 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 δ i u and ω i u .Finally, nonlinear smoothed estimators δ i u and ω i u , respectively, are obtained applying some threshold to the empirical detail coefficients β i jk 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.

Properties of Empirical Coefficients
Now we present some properties of the linear empirical wavelet coefficients  The proofs may be obtained from the authors upon request or from the site of the journal.We assume that functions ω i u , i 0, 1, . . ., n and δ i u , i 1, . . ., m belong to some function spaces F i given by Here, i is the smoothness degree of each space, p i and q i 1 ≤ p i , q i ≤ ∞ specify the norm, and C i1 and C i2 are positive constants.For these function spaces, the following result is valid 19 : where s i i 1/2 − 1/ p i , with p i min{p i , 2}.It can be shown, see Donoho et al. 19 , that the loss incurred by truncating at level J is of order T −2 i / 2 i 1 if we choose J such that 2 J−1 ≤ T 1/2 ≤ 2 J .This is achieved if s i > 1.
Concerning the wavelets, we assume that φ and ψ are compactly supported on 0, 1 and have continuous derivatives up to order r > l, with max i .We denote the spectral norm by • 2 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 V J spaces, generated by {φ J,1 , φ J,2 , . . ., φ J,2 J }.With this basis, we can write 3.3 as  with  Analogously to 3.11 , we have the least squares estimates The relationship between where the Γ is a 2 J 1 ×2 J 1 block diagonal matrix.The matrix Γ does the transformation α i 00 , 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 O p 2 J T −1/2 bias .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.If we write γ ε S, we get where the matrix S is an T − 1 × 1 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 0, 1 and have continuous derivatives up to order r > l, with max i .A2 In the estimation procedure, we have used first a linear estimator, truncating the wavelet expansion at scale 2 J .In order to get functions with the appropriate smoothness, J was chosen such that s i > 1.
A3 The matrix Υ in 4.10 satisfies The Daubechies Extremal Phase wavelets D2N, for an appropriate choice of N, 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 −N 1, N .Since β jk is the coefficient of the regressor ψ jk , the squared error of the least squares estimator is β jk − β jk 2 O p T −1 , which implies that the square error δ u − δ u 2 O p 2 J T −1/2 bias , for example.We can improve this rate as in the next result.

Some Simulations
We now present three simulation examples for the linear and nonlinear estimates δ i • , w i • and δ i • , w i • , 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 v max m, n .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 J, how to estimate σ 2 , the noise variance, and the choice of the threshold parameter λ.Concerning J, we have used the rule specified in Section 4, below 4.3 , resulting J 6 for all examples.The variance was estimated using the empirical wavelet coefficients in all scales.Hard thresholding was used with the parameter λ chosen as λ σ 2 log 2 J .
In all cases, the input time series was generated as an AR 2 model, X t,T a 1 t/T X t−1 ,T a 2 t/T X t−2 ,T ς t , where ς t are independent, normally distributed random variables, with zero mean and variance one, for t 1, 2, . . ., T, and the coefficients are given by
As in the Example 5.1, we generated T 2048 data values for X t,T and Y t,T 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 ω 0 in terms of RMSE.This is a well-known fact, see, for example, Percival and Walden 11 .In Figure 8, we use bootstrap with B 300 replications to see that the true transfer function coefficients lie within the 95% confidence regions.
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 δ 1 at some points of the interval 0, 1 .We choose the coefficients ω 0 u and δ 1 u as piecewise constant functions: ω 0 u 2, for 0 < u ≤ 0.5 and ω 0 u −2, for 0.5 < u ≤ 1; δ 1 u 0.6, for 0 < u ≤ 0.25 or 0.50 < u ≤ 0.75, and δ 1 u −0.5, for 0.25 < u ≤ 0.5 or 0.75 < u ≤ 1.We simulated T 2048 data values for X t,T and Y t,T 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.
In Figure 12, we use bootstrap with B 300 replications to see that the true transfer function coefficients lie within the 95% confidence regions.
In Figures 21 and 22 of the appendix, we present histograms for linear and nonlinear estimators of δ 1 at some points of the interval 0, 1 .

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.
In this application, our aim is to model Valley's wind speeds, {Y t,T }, in terms of its lagged values and present and lagged values of wind speeds at Aberporth, {X t,T }.We should remember that both series are not stationary and so classical methods should not be used.We fitted four potential models, namely, 6.1 The Daubechies Extremal Phase D8 wavelet was used in the computations.The model M 4 was marginally superior in terms of RMSE to model M 3 , so it was discarded.The mean residual sums of squares for M 1 , M 2 , and M 3 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 M 3 are shown in Figure 14.
For Model M 3 , Figures 15 and 16 show the 95% confidence regions based on bootstrap, using linear and nonlinear estimates, respectively.
The bootstrap procedure works as follows: a using the real series, we obtain the estimated transfer function coefficients; b we obtain the fitted series Y t,T using these estimated coefficients; c we obtain the residuals from the fitted model; d we resample these residuals and obtain bootstrap samples of the Y t,T series; e we obtain bootstrap estimators of the transfer function coefficients.
We used B 300 bootstrap replications of δ 1 , δ 2 , ω 0 , and δ 1 , δ 2 , ω 0 , to obtain bootstrap standard errors and the 95% confidence regions.We also fitted the model M 3 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.

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.

Journal of Probability and Statistics
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 F i 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 ε t , 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.

Figure 1 :
Figure 1: Example 5.1: a input series and b output series.

Figure 4 :
Figure 4: Example 5.1: Bootstrap 95% confidence regions for linear estimators, a , c , and e and nonlinear estimators, b , d , and f , with true curves inside.

Figure 5 :
Figure 5: Example 5.2: a input series and b output series.

Figure 8 :
Figure 8: Example 5.2: Bootstrap 95% confidence regions for linear estimators, a and c , and nonlinear estimators, b and d , with true curves inside.
. 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 .

Figure 9 :
Figure 9: Example 5.3: a input series and b output series.

Figure 12 :
Figure 12: Example 5.3: Bootstrap 95% confidence regions for linear estimators, a and c , and nonlinear estimators, b and d , with true curves inside.

2 cFigure 14 :
Figure 14: Model M 3 : a δ 1 u full line and δ 1 u dashed line , b δ 2 u full line and δ 2 u dashed line , c ω 0 u full line and ω 0 u dashed line .

Proposition 4 . 1 .Proposition 4 . 2 .− β i jk 2 O T − 1
Under Assumptions (A1)-(A3), we have that ζ − ζ 2 ∞ O p 2 J T −1 log T .This proposition allows us to derive the rate of convergence of the mean square error of the original linear estimators of the wavelet coefficients.If Assumptions (A1)-(A3) hold, then E β i jk holds uniformly in i, k and j < J.

Example 5 . 3 .
Let Let m 1 and n 0; so the transfer function model is as in Example 5.2.
considered a special case of 1.1 , namely,