The Numerical Computation of Rational Structures and Asymptotic Standard Deviations in Causal Time Series Data

Problem statement: The specific properties of data series are of primary importance in several sciences. In the field of time series analysis, several researchers have considered the rational approximation theory, particularly the Padé Approximation and Orthogonal Polynomials. Approach: In this study, an approach for the statistical significance of two numerical methods, the r-s and q-d algorithm, was proposed which made possible to identify and compute certain rational structures associated with chronological data. Consideration was given to both univariate and multivariate cases. Results: Both algorithms were illustrated empirically through the use of simulated ARMA and TF models and economic data, some of which were taken from previous studies to compare results. Conclusion: This study highlighted the usefulness of several numerical methods (which were all closely related to the PA and OP) in identifying the rational structures associated with data series.


INTRODUCTION
Over the last two decades, several research programs have developed new procedures and characterization techniques for studying dynamic relations associated with time series. In particular, several authors have considered the use of rational approximation techniques in econometric modeling. Many other research areas also share this concern, such as numerical analysis, control theory, statistics and operations research.
In the univariate case, Autoregressive Moving Average (ARMA) models have been extensively studied over the last two decades [1][2][3] . As for the multivariate case, some research has been devoted to identifying the most appropriate vector ARMA (VARMA) models [4][5][6] and Transfer Function (TF) models [2,7,8] for a given dataset. The latter can be considered as a particular case of the VARMA models and are often used to represent dynamic stochastic systems as a set of rational polynomial expressions in an input-output context.
This study shows the statistical significance of several numerical methods which are closely related to the Padé Approximation (PA) and Orthogonal Polynomials (OP) [9,10] . The methods described here will be used for modeling time series and computing the orders of their rational structures. Since the covariance structure of an underlying process exhibits features related to the orders of a model, it is possible to use Hankel determinants with these numerical algorithms to estimate the unknown orders from observations [11] and expectations [12] .
Regarding the paper's structure, we chose not to separate the preliminaries and the auxiliary part of the new results, since the reading and understanding of the text are more consistent if organized into three large blocks (the algorithms, ARMA models and TF models). This explains the criterion selected for the use of indices, that is to say, to maintain the general notation in each field: k, n for the algorithms in numerical analysis; i, j for the ARMA models and k, j for TF models. This notation is the same as in González-Concepción and Gil-Fariña [13] .
First, some known numerical methods for ARMA models (corner method, epsilon algorithm, r-s algorithm and q-d algorithm) and their statistical significance are briefly reviewed. These related techniques can be considered as "feasible alternatives" to each other and are useful for contrasting and/or confirming models which have been identified by other methods. We illustrate the role of the statistical significance of the r-s and q-d algorithms in a simulated ARMA model given in Berlinet and Francq [3] to compare results.
Next, we consider these techniques in the context of a causal TF model with one output and one or more inputs and we present our research on the statistical significance of the r-s and q-d algorithms. This original contribution is along similar lines as Tsay [14] for the corner method and González-Concepción et al. [8] for the epsilon-algorithm, from which we have selected some examples for the purpose of comparing results.
Finally, we present the empirical results and illustrate the main role of the statistical significance in both proposals for a TF model. A simulated TF model from Liu and Hanssens [7] is studied. Some economic applications are also considered using data from Box and Jenkins [16] (a single input TF model for the series M) and from González-Concepción and Gil-Fariña [13] and Gil-Fariña and Lorenzo-Alegría [18] (a multiple input TF model for financial data). Empirical findings shed light on the value of the statistical significance in a real data context.
Each of the numerical examples required costly simulation exercises and considerable computational resources.
This study concludes with a brief summary of our results, the most relevant conclusions and some open questions of interest.

MATERIALS AND METHODS
The Univariate case: Some methods of rational characterization in ARMA models: Theoretical characterization: Let us consider a minimal stationary and invertible ARMA model of order (p,q), defined as: is a sequence of independently and identically distributed random variables with mean zero and variance σ 2 a ; i.e., a noise term. It is assumed that Φ p (L) and Θ q (L) have no common factors.
Various methods have been proposed for identifying the polynomial orders p and q, generally starting from data autocorrelations and the application of criteria from PA theory. We can mention the C-table method from PA [9] and the corner method from the econometric literature [1] . Both methods use a similar rational characterization based on Hankel determinants, although the two approaches evolved independently.
Beguin et al. [1] have also studied the statistical significance of the C-table method. Their study was followed by Tsay [14] and Lii [2] , who proposed an estimator of the asymptotic variance written in terms of the partial derivatives of entries in the C-table.
The relationship between the Hankel determinant and the PA has stimulated the study of other algorithms for ARMA models. The epsilon algorithm, proposed by Berlinet and Francq [3,6] , is of particular interest. Its relationship to OP, the PA and the corner method is examined in Brezinski [10] and its use to characterize ARMA processes is described in Berlinet and Francq [3] . The latter pointed out that the table entries of the epsilon algorithm have also useful statistical properties, based on the methods and assumptions used by [1,2,14] for the corner method.
Finally, we also mention the r-s algorithm proposed by Gray et al. [15] and the q-d algorithm by Berlinet and Francq [3] , both of which have already been proposed for ARMA models.
Brezinski [10] proved that: n k 1,k n k 1 n k 2,k 1 1 n k 1,k n k n k 1,k ρ is the determinant of the Hankel matrix associated with the sequence Lρ. It follows that: The statistical significance of this algorithm can be shown by a method similar to that used by González-Concepción et al. [8] for the epsilon algorithm, namely by computing the Asymptotic Standard Deviation ρ ) represents the estimated asymptotic variance.
Following Tsay [14] and Berlinet and Francq [3] , the variance can be approximately represented by Note that the initial values are: i n i n 0 1 1 i n s ( ) 0 and r ( ) The q-d algorithm and its statistical significance: This algorithm is defined by the relations: It is not directly related to the PA or to OP Brezinski [10] , but it can be shown that: n k,k n k 2,k 1 n k n k 1,k n k 1,k 1 n k 1,k 1 n k,k 1 n k n k 1,k n k,k It follows that: In order to study the statistical significance of the q-d algorithm elements, we use the same statistical technique and a similar notation. Partial derivatives are computed according to the following iterative procedure: o t h e r w i s e The methods explained here have allowed us to obtain a tentative specification of the orders, or even several possible models which can be discriminated between by statistical methods and/or during the estimation stage. Other techniques can also be found, for example in Berlinet and Francq [3] .
Both of these proposals will now be illustrated. To this end, we follow Berlinet and Francq [3] and consider the simulated ARMA(1,1) model X t -0.7X t-1 = a t +0.5a t-1 , t ∀ ∈ where a t ∼N(0,1). The initial values of this series were set to zero and 200 values were generated, but only the last 100 values were considered in the analysis.
We applied the r-s algorithm to the ρ sequence using SCA software. The results concerning to the standard deviations of the i j 1 j r − + elements are shown in Table 1. Taking into account that a stationary process has a minimal ARMA (p,q) representation if [13] : and comparing these numerical entries to certain critical values, different patterns of probable orders can be selected ( Table 2). Note that the empirical results in Table 3 and 4 confirm the simulated model considering a given critical value.   The results obtained via the q-d algorithm relative to the standard deviation of the i j j q − elements and possible selected models are shown in Table 3 and 4 respectively. Again, these results take into account that a stationary process has a minimal ARMA(p,q) representation if [13] : These results suggest that both methods are efficient alternatives for reproducing the orders of the simulated model.
The multivariate case: Some methods of rational characterization in causal TF models: The PA can also be applied in the context of causal rational models to identify the proper VARMA model and is particularly useful in finding the more adequate TF model to available data [16] .
Theoretical characterization: Let us consider a VARMA (p,q) process defined as: Where: Φ p (L), Θ q (L) = Now matrix polynomials of dimension m and degrees p and q respectively L = The back-shift operator Z t = A multiple process u t = A vector of independent white noise components A structure of particular interest arises when: which can be expressed as: If φ is invertible, then Y t is given by: In this expression, which is usually referred to as the TF model, the output Y t is a function of both contemporary and delayed effects of the input variable X t . The existence of a dynamic one-way causal relation X t →Y t is assumed; that is, we assume that there is no "feedback" between output and input. We also assume the presence of a disturbance series described by N t = φ -1 (L)θ(L)a t , where a t is a Gaussian white noise process.
This type of causal model has a great number of practical applications, not only in economics but also in fields such as engineering, geography and business. It is especially useful where there is a correlative or causal structure between variables that are temporally or spatially related.
From the perspective of time series analysis, the Box-Jenkins approach deals with modeling input-output dynamic relations. Here we refer to TF models with only one output Y t ≡ y t and one or more inputs X t ≡ (x it ) I We wish to identify the values of b i , s i and r i and to obtain a satisfactory response of y t to each input. To this aim, several proposals based on algorithms related to the PA have already been considered in previous research. The use of PA offers consistent and reliable initial parameter estimations and no prior information about the noise structure is required.
We can write the relation in the following compact form: where v i (L), which transforms x it into y t , denotes the Impulse Response Function (IRF). The matrix covariance and the impulse response weights v ij for each input are computed first. This may be done either by using the ordinary least squares method or by maximizing the likelihood function in accordance with the following expression: This theoretical sequence of relative weights satisfies the following linear difference equation of order r i and rank b i +s i : which constitutes a characterization of the TF model.
Several methods have been proposed to provide a rational characterization of a TF model, among them the corner method [2,7,14] . This is a generalization of the corner method given in the univariate case, which is applied to each sequence η i . A study of its statistical significance in terms of the ASD can be found in Tsay [14] .
In this context we should also mention the epsilon algorithm of González-Concepción et al. [8] , which can be applied either to the sequence of relative weights or to a transformed sequence if necessary. A study of its statistical significance can be found in Berlinet and Francq [3] and González-Concepción et al. [8] .

The r-s algorithm and its statistical significance:
This iterative procedure has been proposed by González-Concepción and Gil-Fariña [13] for identifying a TF model, in accordance with the following result: Displaying these values in a tabular form, we obtain the block structures given in Table 5 and 6 for each input x it .
.C 1 In certain cases, for example in the epsilon algorithm, some transformation of the sequence of relative weights could be necessary to avoid computational instability.
The statistical significance of this algorithm can be analyzed in a manner similar to the ARMA models, by replacing the autocorrelation sequence (ρ) with the sequence of relative weights (η).

The q-d algorithm and its statistical significance:
This last algorithm was also proposed by González-Concepción and Gil-Fariña [13] and identifies a TF model according to the following characterization:

Theorem 2:
If v i (L) has a rational representation with orders (b i ,s i ,r i ), then one of the following statements holds: Displaying the elements q and d in tabular form, the structures in Table 7 and 8 can be obtained for each x it : The comments about statistical significance of the rs algorithm can be apply to this case.
≠ 0 0 If η i0 is significantly different from zero, then b i = 0

RESULTS AND DISCUSSION
Empirical results: To illustrate these identification methods, we consider the following simulated model with two inputs [7] : t 1,..,100 In this model a t is independent of c t and d t and c t and d t are contemporaneously correlated (with a correlation coefficient of 0.7).
Previous results for the corner method and the epsilon algorithm are given in Liu and Hanssens [7] and González-Concepción et al. [8] , respectively. They do not differ substantially from those given for the r-s algorithm and the q-d algorithm in Table 9-12 of this study.
The Impulse Response Function (IRF) is now computed using the Cochrane-Orcutt iterative method. Ordinary least squares and the other methods yield similar results.  The correct orders for the first input can be adequately identified. For the second input, the possible patterns are (b 2 ,s 2 ,r 2 ) = (2,0,1), (2,0,2), (2,1,1), (2,1,2). In this case, the tables of critical values are not included since the changes in the element values are negligible.

Applications:
Economic data: Now we consider some empirical results of these algorithms, using a data set from a leading sales indicator. These data are identified as "series M" in Box and Jenkins [16] and have also been studied by Tsay [14] . The data consist of 150 bivariate observations and were used by these authors to forecast the sales y t using the leading indicator x t .
The TF model found by Box and Jenkins [16] is: is the difference operator and a t and b t are white noise processes. The identification pattern is thus b = 3, s = 0 and r = 1, which is clearly confirmed by the corner method given in Box and Jenkins [16] . Tsay [14] carried the investigation further, however, studying the statistical significance of the elements in the corner method.     González-Concepción et al. [8] studied the statistical significance of null entries in the epsilon algorithm table, starting from the results of Berlinet and Francq [3] for ARMA models and adopting the approach of Tsay [14] to compute first-order approximations of the variances in the corner method.
We computed the sequence η and applied the epsilon algorithm to the transformed sequence {(-1) j η j } using SCA software and FORTRAN programming respectively, which yielded the standard deviations table and possible pattern orders in Table 13 and 14. These results confirm the model proposed by Box and Jenkins [16] , which was also obtained by Tsay [14] . Other possible patterns could be given, but they are less parsimonious. Applying the r-s algorithm instead yields the results in Table 15 and 16. Finally, Table 17 and 18 provide results obtained with the q-d algorithm.   Comparing the results of these four methods confirms that the pattern (3,0,1) is probably the best model.
Financial data: Volatility is one of the main variables in modern financial theory and particularly relevant for operators dealing with portfolio management, risk hedging and financial structure selection. Authors in the field of financial research have studied a variety of estimators in an effort to identify the one that provides the best information concerning real volatility. In particular, Canina and Figlewski [17] concluded that only a combination of historical volatility (historical data) and implied volatility (market perception) can explain the behavior of real volatility, which is defined as the standard deviation of the stock-pricing rate.
Next, empirical results provided by the algorithms under investigation are used to illustrate the dynamic relationship between real volatility (an endogenous variable) and the historical and implied volatilities (two exogenous variables). Our sample information (for the Spanish market) consists of 905 observations of the IBEX 35 daily futures price, taken from January 14, 1992, through January 26, 1996.
Due to the limitations of available software packages (SCA, Microfit) in dealing with long series, we decided to use Mathematical to estimate the IRF. We used the Ordinary Least Squares (OLS) approach to calculate the sequences η 1 and η 2 and standard covariance matrices as input data for algorithms programmed in FORTRAN.
We consider the model for y t (real volatility), x 1t (implied volatility) and x 2t (historical volatility) obtained by Gil-Fariña and Lorenzo-Alegría [18] and ∆ = 1-L is again the difference operator. The terms a t , b t and c t are all white noise processes.
Applying the r-s algorithm to identify the dynamic structure of both inputs, we obtain Table 19 and 20 for the transformed sequence {(-1) j η 1j }. These numerical results suggest a (0,1,2) model for the first input.
The results obtained for the sequence {η 2j } suggest a (0,0,0) pattern then, the second input has only a contemporaneous effect. Table 21 shows the ASD for the second input. Turning now to the q-d algorithm, we again identify a rational model for the first input. The suggested orders are now (0,1,2), (0,2,2) and (0,3,3), however, which are less parsimonious. The second input is again confirmed to have only a contemporaneous effect on the real volatility. The numerical results of this algorithm are shown in Table 22-24.   The results obtained confirm the models given in Gil-Fariña and Lorenzo-Alegría [18] and González-Concepción and Gil-Fariña [13] .

CONCLUSION
This study highlights the usefulness of several numerical methods (which are all closely related to the PA and OP) in identifying the rational structures associated with data series. We illustrate these methods in the context of causal time series models; that is, ARMA and TF Models.
Special emphasis is placed on the statistical significance of the r-s algorithm and q-d algorithm, in terms of their ASD, as a continuation of several previous contributions in the field of time series analysis.
Empirical findings emphasize the role of this statistical significance in determining the numerical values of the aforementioned algorithms. In general many different possible models may be obtained.
For future research topics, we point out that the correct generalization of these results to VARMA models is not evident. For the corner method, for example, consideration has to be given to the rank of the matrices and not the determinants. The generalization of the r-s and q-d algorithms to VARMA models has not yet been considered. It would also be of interest to consider the generalization of these numerical algorithms to non-causal contexts, including the estimated future values for the inputs of TF models.