2. MODEL CLASS

In this paper we consider ARMAX processes that are defined as follows: Given a sequence of observed exogenous inputs

and a martingale difference sequence

(for detailed assumptions see Section 4) the process

is an ARMAX process if there exist matrices

and (possibly random) initial conditions y_t,t = 1 − p,…,0,z_t,t = 1 − r,…,0,ε_t,t = 1 − q,…,0 such that

solves the vector difference equation (VDE)

for

. The integers p,q,r define the complexity of the model.

Equivalently every ARMAX process

can be written as the solution to the state space equations

for

. Here

are real matrices. The recursions are started using the initial state

, which is assumed to be a random variable with finite second moments uncorrelated with

. The n-dimensional process

is the so-called state process. The state x_t is not observed. For explicit expressions of the matrices (A,B,C,D,K) for given ARMAX system see Appendix A, where some relations between ARMAX and state space representations also are provided. A more detailed exposition can be found in Chapters 1 and 2 of Hannan and Deistler (1988). In the following discussion we will deal exclusively with the state space framework. To make the paper accessible for readers unfamiliar with these concepts, the main ideas needed in this paper are discussed in more detail in Appendix A, where the corresponding notions for ARMAX systems also are given for illustration purposes.

The solution to the state space equations is easily seen to be (cf. Hannan and Deistler, 1988, p. 14)

where the last equation defines

. Therefore

denotes the solution to the initial conditions x₁ = 0. Denoting K₀ := I_s,L₀ := D,K_j := CA^j−1K,L_j := CA^j−1B,j > 0 the transfer functions

are both rational functions in

. These equations define the mapping π : (A,B,C,D,K) [map ] (k(z),l(z)) = (I_s + zC(I_n − zA)⁻¹K,D + zC(I_n − zA)⁻¹B). The transfer functions converge at least for

in a neighborhood of the origin. The state space system (A,B,C,D,K) is called state space representation of the pair (k(z),l(z)) of transfer functions. From (3) it follows that y_t can be additively decomposed into

, the solution corresponding to zero initial state, and the effects of a nonzero initial state, i.e., CA^t−1x₁. From the representation

it is clear that

depends only on the pair (k(z),l(z)) and not on the particular state space representation (A,B,C,D,K) thereof. Thus if the effects of the initial conditions are modeled separately it may be assumed without loss of generality that the system (A,B,C,D,K) is minimal (as defined in Hannan and Deistler, 1988, p. 47; see also Appendix A). For a given pair of transfer functions (k(z),l(z)) a minimal state space representation uses the smallest possible state dimension to represent (k(z),l(z)). The dimension of the state in a minimal state space representation is called the order of the system. The order measures the complexity of the model.

With the exception of Section 6 we will also impose the stability assumption λ_max(A) < 1 (where λ_max(A) denotes the maximum modulus of the eigenvalues of A), implying that k(z) and l(z) converge absolutely for |z| ≤ 1,

. This implies that for stationary input

the initial state x₁ can be chosen such that

is stationary. Furthermore the so-called strict minimum-phase assumption λ_max(A − KC) < 1 is imposed throughout, implying that k(z) is stably invertible on |z| ≤ 1. Let M_n denote the set of pairs of rational transfer functions (k(z),l(z)) having a state space representation of order n where (k(0),l(0)) = (I_s,0), which are analytic on the closed unit disc and where k(z)⁻¹ is analytic on the closed unit disc. The set of pairs M_n is equipped with the so-called pointwise topology (cf. Hannan and Deistler, 1988, p. 54).

Minimal state space systems are not identified. Changing coordinates in the state space (i.e., using

for nonsingular

as the state) leads to a different state space representation. It is possible to restrict sufficiently many entries of the system matrices (A,B,C,D,K) to zero or one such that the remaining unrestricted entries are identifiable parameters for pairs of transfer functions (k(z),l(z)) in a small enough neighborhood of any given pair (k₀(z),l₀(z)). In this paper we will use the parameterization based on the overlapping echelon forms described in Section 2.6 of Hannan and Deistler (1988) using such local identification conditions. In the following discussion we will only establish the relevant notation. For details on the construction of these forms see Appendix B. It should also be noted that the parameterization is only introduced as a means to conveniently compare the estimation accuracy of the estimators. Most other results and in particular the description of the algorithms will not use a particular coordinate system for M_n.

The parameterization consists of a finite collection of mappings φ_α indexed by a multi-index

. For each given α the mapping φ_α is defined on a subset M_n,α ⊂ M_n such that

attaches the parameter vector θ collecting all free entries to each pair of transfer functions in M_n,α. Theorems 2.6.4 and 2.6.5 of Hannan and Deistler (1988) imply the following properties: The parameter space Θ_α is an open subset of

for each α. For each multi-index α it holds that M_n,α is open and dense in M_n and φ_α is continuous on M_n,α. It follows that the specification of α once the order n has been specified is less critical. Finally denoting with

the order of the parameterized pairs of transfer functions, we have ∪_α,|α|=n M_n,α = M_n.

From this discussion and Appendix B it is seen that the question of finding a parameterization of ARMAX systems (or for that matter state space systems) for multivariate systems is not trivial. The underlying technical reason is that there does not exist a continuous parameterization for M_n (cf. the discussion in Hannan and Deistler, 1988, Sect. 2.6.). It might seem to be the case that using the maximal lag lengths (p,q,r) as the measure of complexity rather than the system order n simplifies the situation. However, as discussed on p. 77 of Hannan and Deistler (1988), such a parameterization has theoretical disadvantages. Therefore we will use the relatively complicated parameterization given previously. This carries the advantage that for echelon forms the relations between ARMAX and state space parameters are very simple (see Hannan and Deistler, 1988, Theorem 2.6.5, p. 72).

3. DESCRIPTION OF THE ALGORITHM

The basic idea underlying the state approach to subspace algorithms is the fact that if the state were known the state space system equations (2) could be used to obtain estimates of the system matrices using simple regressions. The main problem then is to obtain a good estimate of the state.

In this respect the properties of the state for prediction purposes are central. Consider the problem of predicting y_t+j,j ≥ 0 based on the knowledge of

. Let y_t+j|t denote the best mean square prediction. Analogous to (3) for j ≥ 0 we have the equation

because

is assumed to be a martingale difference sequence. Note that y_t+j|t depends on z_l,y_l,l < t only via x_t and that furthermore y_t+j|t is a linear function of the state x_t and the future inputs z_t+i,i ≥ 0 for all j ≥ 0. On the other hand it follows also from the system equations that

where for A := A − KC,B := B − KD we have

and W_t,p⁻ := [w_t−1′,…,w_t−p′]′ denoting with w_t the joint process, i.e., w_t := [y_t′,z_t′]′. Therefore it follows that x_t is a linear function of the data obtained up to time t and the initial state x₁. Under the strict minimum-phase assumption, (A − KC)^p will be small for large p, and hence a reasonable approximation of x_t is given by

. Therefore an estimate of the state can be obtained by replacing

with an estimate

in this approximation.

These two facts build the basis for the estimation of the state. Consider the prediction equation (4) written jointly for j = 0,…,f − 1 for some integer f (where t > p is assumed):

where Y_t,f⁺ := [y_t′,…,y_t+f−1′]′ collects the predicted variables, Z_t,f⁺ := [z_t′,…,z_t+f−1′]′ is the future of the exogenous inputs, and N_t,f⁺ denotes the prediction errors. Finally

is sometimes called the truncated observability matrix. Expressions for

are easy to derive but not essential for our purposes and hence are omitted.

Equation (5) lies at the core of the subspace idea. Neglecting

under the assumption A^p ≈ 0 the matrices

can be estimated using regression because W_t,p⁻ and Z_t,f⁺ are observed and the residuals N_t,f⁺ are orthogonal to the regressors. Moreover for f sufficiently large the fact that

, where n denotes the order of the system, can be exploited. This motivates the following procedure.

2. Estimate the equation Y_t,f⁺ = β₁W_t,p⁻ + β₂ Z_t,f⁺ + N_t,⊥,t = p + 1,…,T − f + 1 under the restriction rank(β₁) = n. This leads to the estimates

of β₁ (rank n). In this step also the order n of the estimated system has to be specified.

3. Estimate the state

.

4. Regress y_t onto

to obtain estimates

and residuals

.

5. Regress

for t = p + 1,…,T to obtain estimates

to compute the estimates

.

6. Check the minimum-phase assumption. If the estimated system is not minimum-phase the estimated spectrum

where

(and .* denotes conjugate transposition) has to be refactored to ensure that the minimum-phase assumption holds (cf. Hannan and Deistler, 1988, Theorem 1.3.3, p. 26).

The central step in the algorithm is the reduced rank regression. Let

denote the least squares estimator ignoring the rank restriction on β₁. Then the rank n approximation

is obtained using the singular value decomposition (SVD)

where

denotes the matrix whose columns are the left singular vectors of

corresponding to the n dominating singular values

denotes the matrix of the corresponding right singular vectors. The approximation error arising from using a system order n is denoted as

. Finally

denote two weighting matrices that for the moment only are assumed to be nonsingular a.s. More specific conditions will be provided subsequently. This results in estimates

. If a different factorization into

for some nonsingular matrix

is used the resulting estimates will be observationally equivalent; i.e., the estimated pair of transfer functions will be unchanged. The asymptotic properties of the state space representation

on the contrary depend on the choice of the factorization.

The term canonical correlation analysis (CCA) has its origin in the approximation step. Let

for two vector sequences

somewhat sloppily using the same symbol for the processes and their realizations at time t. Introduce the notation

such that Y_t,f^+,Π denotes the residuals of a regression of Y_t,f⁺ onto Z_t,f⁺ and correspondingly W_t,p^−,Π are the residuals of the regression of W_t,p⁻ onto Z_t,f⁺ . Then the original version of CCA uses

where X^1/2 = (X^1/2)′,(X^1/2)² = X denotes the symmetric square root of the positive definite matrix X. Further the unrestricted estimate

. Therefore in the case of no inputs the weightings used in CCA imply that the SVD amounts to the estimation of the canonical correlations and vectors of the two subspaces spanned by the components of Y_t,f⁺ and W_t,p⁻, respectively (because there is no input here Y_t,f^+,Π = Y_t,f⁺,W_t,p^−,Π = W_t,p⁻).

The choice of weightings according to (6) corresponds to the maximum-likelihood solution to the reduced rank regression Y_t,f⁺ = β₁W_t,p⁻ + β₂ Z_t,f⁺ + N_t,⊥ such that rank(β₁) = n based on Gaussanity (and uncorrelatedness over time) assumptions on N_t,⊥. The least squares solution is obtained by using

. With respect to

only the choice given previously will be considered throughout. In the following discussion the algorithm using the weightings (6) will be referred to as CCA. Algorithms using a different weighting scheme but the same outline in the literature have been called Larimore approach or state approach. We will use the second term in the following discussion.

In the reduced rank regression both the outputs and the regressors are lags of y_t,z_t and w_t, respectively, and hence the problem has a similar structure as autoregressions. As in autoregressions there are two possible implementations that in the stationary case possess identical asymptotics. Either only time instants are used in the regression where all observations are available, i.e., the range of t is limited to t = p + 1,…,T − f + 1, or alternatively, the regression can be based on the estimated second moments as in the Yule–Walker equations, essentially imposing zero initial and end conditions (cf. Hannan and Deistler, 1988, p. 212). The first approach might seem to be favorable for small sample sizes and relatively large values of f and p as the imposition of zero initial conditions introduces avoidable bias. On the other hand the usage of covariance estimates makes the procedure easily adaptable to data sets including missing values. Using robust estimators of the covariances leads to (statistically) robust versions of the state approach. A recursive (in sample size) formulation of the regression problem leads to recursive estimators.

Finally the term subspace algorithm needs explanation. It originates from the reduced rank regression in the second step as the most popular methods in the class of subspace algorithms, that is, MOESP (Verhaegen, 1994) and N4SID (Van Overschee and DeMoor, 1994), use the matrix

in subsequent steps rather than

. Because the factorization into

is of no importance for the estimated pair of transfer functions, the reduced rank regression problem essentially estimates the subspace spanned by the columns of the matrix

, hence the name. MOESP continues to obtain estimates for (A,B,C,D) but does not estimate K, which implies that only l(z) is estimated. Frequently in econometric analysis one is primarily interested in the noise dynamics, and therefore MOESP does not seem to be a suitable algorithm for these purposes. Therefore it is not considered further in this survey. The same remark is true for one variant of N4SID. An alternative version of N4SID and a recently proposed algorithm by Chiuso and Picci (2003) use a more complicated method to estimate the shifted state

. Because the analysis of these methods is still in its infancy we will not consider these algorithms further. We refer to Chiuso and Picci (2003) for details.

The description of the state approach makes it obvious that there are a couple of design parameters involved. First f and p have to be chosen; second the weighting matrices

and the order n have to be specified. A large part of the literature on subspace methods deals with the effects of these choices on the asymptotic properties of the corresponding estimators. Section 4 will survey the main results in this respect for the stationary case, and Section 6 will deal with the case of integrated processes, whereas Section 5 discusses procedures to estimate the order n.

5. ORDER ESTIMATION

Up to now it was always assumed that the true order of the system is known. This section will discuss methods to specify the order.

The first step in the subspace procedure is a reduced rank regression in the equation

under the constraint rank(β₁) = n. The first idea for order estimation hence would be to use the reduced rank regression framework. Under the assumption that N_t,⊥ is i.i.d. Gaussian the likelihood ratio between using rank(β₁) = n and rank(β₁) = n + 1 is equal to

where

denotes the ith largest singular value of 〈Y_t,f^+,Π,Y_t,f^+,Π〉^−1/2〈Y_t,f^+,Π,W_t,p^−,Π〉〈W_t,p^−,Π,W_t,p^−,Π〉^−1/2 (Anderson, 1951). But

where

Therefore (N_t,⊥)_{t=p+1,…,T−f+1} is not uncorrelated over time and moreover depends on p. Consequently the distribution of the pseudo-likelihood ratio does not follow from standard theory. To the best of my knowledge the asymptotic distribution of the likelihood ratio is not known.

An alternative idea then is to formally use the information criteria defined as

where

denotes the number of nonzero singular values. Here d(n) = 2sn + sm + nm equals the number of parameters included in a model of order n and C_T > 0,C_T /T → 0 denotes the penalty term punishing complex models. Here SVC stands for singular value criterion. It should be stressed that this criterion is not an information criterion being defined in a misspecified case for the same reasons that the likelihood ratio does not follow standard asymptotics. This criterion has been discussed by Camba-Mendez and Kapetanios (2001). Peternell (1995) introduced a similar criterion using

, the Frobenius norm of the approximation error. Bauer (2000) proposes

, i.e., the two norm of the approximation error, based on simulations. The advantage of the last two criteria is that they can be used for general weighting matrices

, whereas the first criterion uses the CCA choices to guarantee

. In the following discussion we will only deal with these three choices.

The order is estimated as the minimizing integer of SVC(n;C_T) resulting in an estimator

. The rationale behind these estimates is that

where the matrix U_n0(f,p)Σ_n0(f,p) × V_n0(f,p)′ contains at most n₀ nonzero singular values, the remaining ones being zero. Recall that it is the rank of β₁ that we want to infer from the noisy estimate

. Letting the penalty factor C_T > 0 tend to infinity such that C_T /T converges more slowly to zero than the maximal estimation error in the singular values then ensures that asymptotically the penalty dominates and asymptotically the order is not overestimated, implying that

. In fact we have the following condition.

Convergence of Order Estimation Using SVC(n;C_T). Let the processes

fulfill Assumption 1 with homoskedastic innovations

where the true order is denoted as n₀. Additionally it is assumed that for some constants 0 < c < ∞ and 0 < ρ < 1 it holds that

. Assume that f and p are either fixed and finite or tend to infinity, in which case f ≤ p,p = O((log T)^a),a < ∞ is assumed. Let U_n0(∞)Σ_n0(∞)V_n0(∞)′ denote the limit of

, which for f → ∞ or p → ∞ (a.s.) has to be interpreted in an appropriate sense (for details see Bauer, 2000). Then C_T > 0,C_T /T → 0,C_T /(fp log log T) → ∞ is sufficient for

to converge to rank(U_n0(∞)Σ_n0(∞)V_n0(∞)′).

Note that U_n0(∞)Σ_n0(∞)V_n0(∞)′ depends on the choice of f and p. Jansson and Wahlberg (1998) show by example that f ≥ n₀ + 1,p ≥ n₀ is not sufficient for σ_n0(f,p) > 0. However, violation of this condition is rare because Bauer and Jansson (2000) show that for fixed pseudostationary input sequence

that is persistently exciting of order f + p − 1 (i.e., fulfilling conditions (4.1.1.), (4.1.14), and (4.2.21) of Hannan and Deistler, 1988, the latter with f + p − 1 inserted for sp; these assumptions hold for a particular realization of a process

fulfilling Assumption 1) for fixed f and p the subset of pairs of transfer functions in M_n0 with σ_n0(f,p) = 0 is the complement of an open and dense subset. This implies that for any given f,p it is unlikely that the true system happens to have the property σ_n0(f,p) = 0. Nevertheless the problem is relevant because for systems close (in an appropriate sense) to these exceptional systems the smallest singular value will be very small and hence a large number of observations will be needed to differentiate these singular values from noise. Chui (1997) gives sharp conditions on f and p for n₀ = rank(U_n0(f,p) × Σ_n0(f,p)V_n0(f,p)′) that imply the sufficient conditions min(f,p) ≥ 3n₀. Another set of sufficient conditions is f ≥ n₀ and p → ∞.

Bauer (2001, Theorem 3) states consistency of order estimation using SVC(n;C_T) for p → ∞ under the assumptions on the input given previously. Only the Frobenius norm and the two norm of

in the definition of SVC(n;C_T) are considered in that paper. There convergence of

for fixed and finite f and p also is discussed (allowing a harmonic process to be added to the process fulfilling the assumptions just stated). Because again mainly the uniform convergence of the sample covariances is used in the proof of Theorem 3 of Bauer (2001) the consistency proof for SVC(n;C_T)-based order estimation given there can be straightforwardly generalized to yield the result given earlier.

. In this case the condition C_T /(fp log log T) → ∞ has to hold a.s. Again the sufficient conditions are by no means necessary.

Simulation evidence shows that for

the sufficient condition on the penalty term seems to be too pessimistic, leading to a high risk of underestimating the order in finite samples (see the simulations in Bauer, 2001). The criterion based on the two norm seems to be less sensitive with respect to the choice of f and p than the method based on the Frobenius norm. No thorough motivation for any penalty term has been found.

A second approach to the problem of order specification is to note that the problem of selecting the order is equivalent to inferring the rank of the noisy matrix

. Different procedures using sequential testing are developed in Sorelius (1999) in the case of no exogenous inputs included (m = 0) and Camba-Mendez and Kapetanios (2001). All these procedures center on the assessment of the rank of 〈Y_t,f⁺,W_t,f⁻〉 and proceed by sequentially testing for singularity of the matrix and increasing the size by one in case of a rejection until the test is not rejected for the first time.

The Hankel structure of the matrix 〈Y_t,f⁺,W_t,f⁻〉 complicates the derivation of the test statistic considerably. Additionally one runs into a sequential testing problem. Having said that, it should also be mentioned that this approach is statistically much sounder than the ad hoc definition of criteria such as SVC(n;C_T). At present the development of this strand of literature is deemed to be premature, however. We will not go into detail in this respect, and we refer the reader to Sorelius (1999) and Camba-Mendez and Kapetanios (2001) for details.

A third approach to order estimation uses analogies to traditional information criteria. In the context of pseudo-maximum-likelihood estimation the order is estimated using information criteria such as the Akaike information criterion (AIC) or BIC (for details see Hannan and Deistler, 1988, Ch. 5). Assume that the Gaussian likelihood is maximized over M_n for n = 0,1,2,…,H_T for some upper bound H_T such as (log T)^a, 0 < a < ∞, resulting in an estimate

for each n. Then an estimate of the one step ahead prediction errors

is obtained from Kalman filtering (for details see Hannan and Deistler, 1988, Ch. 3), leading to an estimated one step ahead prediction variance

. Then the information criterion estimate is taken as a minimizing integer of the criterion

where again d(n) = 2sn + sm + nm denotes the number of parameters included in a model of order n and C_T > 0, C_T /T → 0 denotes the penalty term punishing complex models. BIC uses C_T = log T and AIC uses C_T = 2.

For subspace methods there are two obvious generalizations of this approach. The pseudo-maximum-likelihood estimator

could simply be replaced by the subspace estimator

in the definition of

. Alternatively the estimation of the equation

immediately leads to an estimate of the residual variance

that can be inserted into

. Here the dependence of

on the usage of the n-dimensional state

is emphasized. These two alternatives are fundamentally different. Whereas

takes the full system dynamics into account,

is calculated on the basis of the n-dimensional estimated state, which in the limit need not correspond to an n-dimensional system.

To illustrate the point assume for the moment that we are selecting the model order based on the knowledge of the true state

corresponding to a state space representation such that C = [I_s,0] where n₀ > s is assumed. Then let x_t(n),1 ≤ n ≤ n₀ denote the vector of the first n components of x_t. Now order estimation based on

, i.e., the residual variance of the equation

, will result in an order estimator converging to s (for suitable choice of the penalty term C_T) because the remaining coordinates of x_t do not contribute to the prediction of y_t. The order estimate based on

, however, will lead to consistent order estimates because in this case for each order n a new estimate of the state also is obtained from the Kalman filter, which for n < n₀ will not coincide with x_t(n).

For the subspace methods the estimated state

may not correspond to a true system where the last columns of C are zero for any pair (k₀(z),l₀(z)) because the limiting system is determined implicitly in the SVD. But problems for the order estimation procedure are expected also in situations where the norm of the last column is small. An example in this respect is provided in Bauer (2001), where some asymptotic properties of the order estimate obtained by minimizing

also are derived. Here IVC stands for innovation variance criterion. Because of these problems the use of the criterion is strongly discouraged.

Finally in this section the asymptotic properties of the order selection based on

where

denotes the CCA estimator using f = p,p ≥ −d log T/(2 log ρ₀), d > 1,p = O((log T)^a),p chosen deterministically, are investigated. Additionally assume that

so that no inputs are included. In this case it can be shown that

and therefore C_T > 0,C_T /T → 0,C_T /log log T → ∞ is sufficient for weak consistency of the order estimate obtained by minimizing

over 0 ≤ n ≤ ps. The expression for n ≠ n₀ is immediate because

minimizes

. For n = n₀ the expression is a consequence of

. The details of the proof are omitted.

In the situation without inputs z_t the criterion

has the tremendous advantage compared to

that once the state has been estimated for a maximal order the subspace estimators for a range of different orders can be obtained extremely numerically efficiently (in fact, only one regression has to be estimated; for details, see Bauer, 2001, p. 1564). Therefore, Kalman filtering to obtain

constitutes the only main computational cost. This idea also has another interesting feature. As stated in Section 4, CCA produces asymptotically optimal estimators in the correctly specified case for Gaussian innovations. For misspecified cases no such result is known, and therefore the subspace estimator might result in worse values of the criterion function than the pseudo-maximum-likelihood estimator that minimizes

. Hence order estimation based on

could be expected to improve on BIC order estimation. It should be stressed, however, that the performance in cases where there exists a true model might not be the only criterion to assess the quality of a model selection procedure (cf. the discussion in Shibata, 1981).

7. SIMULATIONS

It has been argued that the main advantage of subspace identification lies in the numerical simplicity and reliability as an initial guess. These claims have been confirmed in a number of simulation studies. The statistical properties of CCA have been highlighted by Peternell (1995), Deistler et al. (1995), Peternell et al. (1996), and Bauer (1998). Bauer (2004c) also contains a comparison of the numerical load. In the integrated case Bauer and Wagner (2002a) demonstrate the usage of CCA for estimating the cointegrating rank, and Wagner (2004) examines the finite sample properties of the estimation of the cointegrating space and the predictions based on the CCA estimates. In this section a small simulation study will be used to demonstrate the claims of the paper. The simulation is not intended to mirror any particular application, and the simulated systems are totally arbitrary. The results reflect the findings in the literature and therefore can be seen to be typical. The main point to be made in the simulations, however, is that there are situations in which subspace methods appear to be a useful alternative to the classical estimation methods.

In the simulation experiment CCA is compared to the algorithm PEM of the MATLAB system identification toolbox. This routine uses an initial estimate obtained by a version of N4SID to afterward maximize the Gaussian likelihood in a Gauss–Newton type method. The method has been chosen because of its widespread use in the engineering community. The basic setup is the following: For each experiment three integers (s,n,T) are given where s denotes the output dimension, n the system order, and T the sample size. Based on these integers an ARMA system is drawn randomly and 100 time series of dimension s and sample size T + 1,000 are generated using this one system and Gaussian white noise of zero mean and variance equal to the identity matrix. The initial conditions are set to zero, which might be seen as unrealistic. Then CCA setting

and estimating the order via SVC(n;log T) based on the two norm

is used to estimate the system on the first T data points. After the completion of this estimation another system is estimated on the same data set using the function PEM and the order estimated in CCA. For both estimated systems the one step ahead prediction error on the remaining 1,000 time instants is calculated. As a measure of accuracy

is chosen where

denotes the estimated innovation variance on the validation sample for the two procedures. For the true system

asymptotically is normally distributed with mean zero and variance 2s. For each experiment three numbers are reported for each of the two algorithms.

2. The average of the 100 accuracy measures

where i = 1 for CCA and i = 2 for PEM.

3. The average over the subsample where PEM is not obviously trapped in a local minimum: In some cases the estimation error for PEM is many orders of magnitude too high compared to the error attained by CCA. Formally these cases are defined somewhat arbitrarily by

.

These three numbers are reported for each setting of integers (s,n,T) in Table 1. The main findings are shown.

• CCA is always numerically cheaper than PEM. The difference becomes remarkable for large sample sizes and high-dimensional data sets.
• The specification of the initial guess is an issue. If one does not account for local optima the performance of PEM is worse than the performance of CCA in most cases.
• In medium sample sizes (here T = 1,000) PEM achieves slightly better performance.
• For large sample size PEM does not seem to be trapped in local optima. At the same time the performance is not differentiable from the performance of CCA.
• For small sample sizes SVC picks overly small orders. This is especially the case for high-dimensional data sets. The consequence is a sizable bias.

These findings demonstrate the claims based on asymptotic theory raised in the literature.

Numerical and statistical properties of CCA and PEM