{'fraction_non_alphanumeric': 0.043730016066684696, 'fraction_numerical': 0.014205494488013616, 'mean_word_length': 4.073769168684423, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 40, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Finite sample forecasting with estimated temporally aggregated linear processes

We propose a finite sample based predictor for estimated linear one dimensional time series models and compute the associated total forecasting error. The expression for the error that we present takes into account the estimation error. Unlike existing solutions in the literature, our formulas require neither assumptions on the second order stationarity of the sample nor Monte Carlo simulations for their evaluation. This result is used to prove the pertinence of a new hybrid scheme that we put forward for the forecast of linear temporal aggregates. This novel strategy consists of carrying out the parameter estimation based on disaggregated data and the prediction based on the corresponding aggregated model and data. We show that in some instances this scheme has a better performance than the"all-disaggregated"approach presented as optimal in the literature.


Introduction
The success of parametric time series models as a tool of choice in many research fields is due in part to their good performance when it comes to empirical forecasting based on historical samples. Once a data generating process (DGP) has been selected and estimated for the forecasting problem at hand, there is a variety of well studied forecasting procedures and algorithms available in the literature. The most widespread loss function used in the construction of predictors is the mean square forecasting error (MSFE); see the monographs [Bro06,BD02,Ham94,L05] and references therein for detailed presentations of the available MSFE minimization-based techniques. This is the approach to prediction that we follow in this work; the reader is referred to [Gra69] or Section 4.2 in [GN86] for forecasting techniques based on other optimality criteria.
The stochastic nature of the time series models that we consider implies that the forecasts produced with them, carry in their wake an error that cannot be minimized even if the parameters of the model are known with total precision; we refer to this as the characteristic error of the model. Additionally, all that is known in most applications is a historical sample of the variable that needs to be forecasted, out of which a model needs to be selected and estimated. There are well-known techniques to implement this, which are also stochastic in nature and that increase the total error committed when computing a forecast; we talk in that case of model selection error and estimation error. All these errors that one incurs in at the time of carrying out a forecasting task are of different nature and much effort has been dedicated in the literature in order to quantify them in the case of linear multivariate VARMA processes.
Most results obtained in this direction have to do with the combination of the estimation and the characteristic errors; this compound error is always studied assuming independence between the realizations of the model that are used for estimation and the ones used for prediction; we refer the reader to [Bai79,Rei80,Yam80,Yam81,Duf85,BS86,L87,SH88]. Explicit expressions for these errors in the VARMA context are available in the monograph [L05]. Indeed, if we assume that the sample out of which we want to forecast is a realization of the unique stationary solution of a VAR model, this error can be written down [L05, page 97] using the time-independent autocovariance of the process; the situation in the VARMA context is more complicated and the expression provided [L05, page 490] requires Monte Carlo simulations for its estimation.
The knowledge regarding the error associated to model selection is much more rudimentary and research in this subject seems to be in a more primitive state. A good description of the state of the art can be found in [L06, page 318] as well as in [L86b, page 89], and references therein. We do not consider this source of forecasting error in our work and hence in the sequel we will use the denomination total error to refer to the combination of the characteristic with the estimation errors.
In this paper we concentrate on one dimensional linear processes, a subclass of which is the ARMA family. The first contribution in this paper is the formulation of a MSFE based predictor that takes as ingredients a finite sample and the coefficients of a linear model estimated on it, as well as the computation of the corresponding total error. The main improvements that we provide with respect to preexisting work on this question are: • We make no hypothesis on the second order stationarity of the sample at hand; in other words, we do not assume that the sample is a realization of the stationary solution of the recursions that define the model. Such a hypothesis is extremely difficult to test in small and finite sample contexts and it is hence of much interest to be able to avoid it.
• The expression for the total forecasting error is completely explicit and does not require the use of Monte Carlo simulations.
The interplay between the characteristic error, the estimation errors, and the forecasting horizon is highly nonlinear and can produce surprising phenomena. For example, as it is well known, the characteristic error is an increasing function of the horizon, that is, the further into the future we forecast, the more error we are likely to commit. When we take into account the estimation error, the total error may decrease with the forecast horizon! We study this finite sample related phenomenon with the total error formula introduced in Theorem 3.3 and illustrate it with an example in Section 5.1. The characterization of the total forecasting error that we described serves as a basis for the second main theme of this paper, namely, the interplay between multistep forecasting, the prediction of temporal aggregates, and the use of temporal aggregation estimation based techniques to lower the total forecasting error. In this part of the paper we work strictly in the ARMA context. The temporal aggregation of ARMA processes is a venerable topic that is by now well understood [AW72, Tia72, Bre73, TW76, Wei79, SW86, Wei06, SV08] and has been extensively studied and exploited in the context of forecasting [Abr82, L86b, L86a, L87, L89a, L89b, RS95, L06, L09, L10] mainly by H. Lütkepohl. A recurrent question in this setup consists of determining the most efficient way to compute multistep forecasts or, more generally, predictions of linear temporal aggregates of a given time series. More specifically, given a sample and an underlying model, we can imagine at least two ways to construct a h time steps ahead forecast, or in general the one that is a linear combination of the h steps ahead values for the time series. First, we can simply compute the h time steps ahead forecasts of the time series out of the original disaggregated sample and to determine the needed aggregate prediction out of them; another possibility would be to temporally aggregate the sample and the time series model in such a way that the required forecast becomes a one time step ahead forecast for the new aggregated sample and model. If we do not take into consideration estimation errors and we only consider the characteristic error, there is a general result that states that the forecast based on high frequency disaggregate data has an associated error that is smaller or equal than the one associated to the aggregate sample and model (we will recall it in Proposition 4.2). In the VARMA context, H. Lütkepohl [L86b, L87, L09] has characterized the situations in which there is no loss of forecasting efficiency when working with temporally aggregated ingredients.
When estimation errors are taken into account, the inequality that we just described becomes strict [L86b, L87], that is, forecasts based on models estimated using the dissagregated high-frequency samples perform always better than those based on models estimated using aggregated data. This is so even in the situations described in [L86b, L87, L09] for which the characteristic errors associated to the use of the aggregated and the disaggregated models are identical; this is intuitively very reasonable due to the smaller sample size associated to the aggregated situation, which automatically causes an increase in the estimation error.
In Section 4.3 we propose a forecasting scheme that is a hybrid between the two strategies that we just described. We first use the high frequency data for estimating a model. Then, we temporally aggregate the data and the model and finally forecasting is carried out based on these two aggregated ingredients. We will show that this scheme presents two major advantages: • The model parameters are estimated using all the information available with the bigger sample size provided by the disaggregated data. Moreover, these parameters can be updated as new high frequency data becomes available.
• In some situations, the total error committed using this hybrid forecasting scheme is smaller than the one associated to the forecast based on the disaggregated data and model and hence our strategy becomes optimal. Examples in this direction for both stock and flow temporal aggregates are presented in Section 5. The increase in performance obtained with our method comes from minimizing the estimation error; given that the contribution of this error to the total one for univariate time series models is usually small for sizeable samples, the differences in forecasting performance that we will observe in practice are moderate. As we will show in a forthcoming work, this is likely to be different in the multivariate setup where in many cases, the estimation error is the main source of error.
To our knowledge, this forecasting scheme has not been previously investigated in the literature and the improvement stated in the last point seems to be the first substantial application of temporal aggregation techniques in the enhancing of forecasting efficiency.

Finite sample forecasting of linear processes
In this section we introduce notations and definitions used throughout the paper and describe the framework in which we work. Additionally, since we are interested in finite sample based forecasting, we spell out in detail the predictors as well as the information sets on which our constructions are based.

Linear processes
Let ε = {ε t } ∞ t=−∞ be a set of independent and identically distributed random variables with mean zero and variance σ 2 . We will write in short

Finite sample forecasting with estimated temporally aggregated linear processes
We say that X = {X t } ∞ t=−∞ is a linear causal process whenever it can be represented as is a set of real constants such that ∞ i=0 |ψ i | < ∞. Expression (2.1) can also be rewritten as where L is the backward shift operator and Ψ (z) is the power series and Π (L) can also be referred to as causal linear filter and invertible linear filter, respectively.

Finite sample forecasting of causal and invertible ARMA processes
Consider the causal and invertible ARMA(p, q) specification determined by the equivalent relations (2. 3) The innovations process ε = {ε t } can be either independent and identically distributed IID(0, σ 2 ) or white noise WN(0, σ 2 ). In this subsection we focus on how to forecast out of a finite sample ξ T = {x 1 , ..., x T } that satisfies the relations (2.3) and that has been generated out of a presample {x 1−p , ..., x 0 } and a preinnovations set {ε 1−q , ..., ε 0 }. A standard way to solve this problem [Bro06,BD02] consists of assuming that ξ T is a realization of the unique stationary process X that satisfies the ARMA relations (2.3) and to use its corresponding time independent autocovariance functions to formulate a linear system of equations whose solution provides the linear projection X T +h of the random variable X T +h onto {x 1 , ..., x T } using the L 2 norm; this projection X T +h minimizes the mean square error. We recall that writing the unique stationary solution of (2.3) usually requires knowledge about the infinite past history of the process. For example, for an AR(1) model of the form X t − φX t−1 = ε t , the unique stationary solution is given by Given that we are concentrating in the finite sample context, we prefer for this reason to avoid the stationarity hypothesis and the use of the corresponding autocovariance functions and to exploit in the forecast only the information that is strictly available, that is: (i) The model specification (2.3): we assume that the model parameters are known with certainty and we neglect estimation errors.
Let r = max {p, q} and define the enlarged preset I * as where: • if p > q: r = p and ε t := t+p−1 j=0 • if q = p: I = I * .
The enlarged preset I * is defined as a function of the elements in I. Consequently, the sigma-algebras σ(I) and σ(I * ) generated by I and I * , respectively, coincide, that is, σ(I) = σ(I * ). The following result is basically known (see for example [Ham94,L05]) but we include it in order to be explicit and self-contained about the forecasting scheme that we are using in the rest of the paper and also to spell out the peculiarities of the finite sample setup in which we are working. We include a brief proof in the appendix.
Proposition 2.1 In the conditions that we just described: (i) The information sets (sigma algebras) σ T := σ (I, ε 1 , ..., ε T ) and σ ξ T := σ (I, x 1 , ..., x T ) generated by the preset and the past histories of the innovations T := {ε 1 , . . . , ε T } and the sample If the innovations process is IID (respectively WN) then the optimal multistep forecast X T +h (respectively optimal linear forecast) based on σ (ξ T ) that minimizes the mean square forecasting error (MSFE) E X T +h − X T +h 2 is: The MSFE associated to this forecast is (2.6) (iv) For ARMA models, the forecasts constructed in (2.5) for different horizons with respect to the same information set F T := σ (ξ T ) = σ ( T ), satisfy the following recursive formula: Remark 2.2 Testing the stationarity of small or in general finite samples is a difficult task in practice. We emphasize that the prediction in Proposition 2.1 does not require any stationarity hypothesis. Moreover, we underline that the forecast (2.5) does not coincide in general neither with the standard linear forecast for second order stationary series that uses the corresponding time independent autocovariance function (see for example [Bro06], page 63), nor with the usual finite sample approximation to the optimal forecast (see [Ham94], page 85). The main difference with the latter consists of the fact that the innovations associated to the presample are not assumed to be equal to zero but they are reconstructed out of it so that there is no loss of information.
In the examples 2.3 and 2.4 below we show how our forecast allows us to construct a predictor that: (i) is different from the one obtained assuming stationarity; (ii) has a better performance in terms of characteristic forecasting error.
These statements do not generalize to arbitrary ARMA models; for example, for pure AR models, the predictor that we propose and those cited above coincide.
Consequently by (2.5), the forecast X 2 based on the information set F 1 = σ ({I, x 1 }), is given by and has the associated error MSFE( X 2 ) = σ 2 .
On the other hand, the forecast that assumes that x 1 is a realization of the unique stationary solution of (2.8) and that uses the corresponding autocovariance function [Bro06, page 63] is given by and has the associated error which shows that the forecast that we propose has a better performance than the one based on the stationarity hypothesis.
The better performance of the forecast that we propose in the preceding example can be in part due to the fact that we are using for X 2 additional information on the preinnovations that is not taken advantage of at the time of writing X S 2 . In the following example we consider an ARMA(1,1) model and we see that the statements of Remark 2.2 also hold, even though in this case, unlike in the MA(1) situation, the information sets on which the two forecasts considered are based are identical.
Example 2.4 Finite sample forecasting for the ARMA(1,1) process. Consider the model Then, • π 0 = 1, π j = (−1) j (φ + θ) θ j−1 , for any integer j ≥ 1, We consider the trivial sample x 1 generated by the preset I = {x 0 , ε 0 } = I * . Using Proposition 2.1, we have that the one-step ahead forecast X 2 based on the information set F 1 = σ ({I, x 1 }), is given by On the other hand, the forecast based on the stationarity hypothesis using the same information set, yields It is easy to check that the statement MSFE( X S 2 ) > MSFE( X 2 ) is equivalent to θ 4 (θ + φ) 2 > 0, which is always satisfied and shows that the forecast that we propose has a better performance than the one based on the stationarity hypothesis.

Forecasting with estimated linear processes
In Proposition 2.1 we studied forecasting when the parameters of the model are known with total precision. In this section we explore a more general situation in which the parameters are estimated out of a sample. Our goal is to quantify the joint mean square forecasting error that comes both from the stochastic nature of the model (characteristic error) and the estimation error; we will refer to this aggregation of errors as the total error. This problem has been extensively studied in the references cited in the introduction always using the following two main constituents: • Estimation and the forecasting are carried out using independent realizations of the time series model.
• The model parameter estimator is assumed to be asymptotically normal (for example, the maximum likelihood estimator); this hypothesis is combined with the use of the so called Delta Method [Ser80] in order to come up with precise expressions for the total error.
The most detailed formulas for the total error in the VARMA context can be found in [L05] where the Delta Method is applied to the forecast considered as a smooth function of the model parameters. If we assume that the sample out of which we want to forecast is a realization of the unique stationary solution of a VAR model, an explicit expression for this error can be written down by following this approach [L05, page 97] that involves the time-independent autocovariance of the process. In the VARMA setup, the situation is more complicated [L05, page 490] and the resulting formula requires the use of a Monte Carlo estimation.
In subsection 3.1 we start by obtaining a formula for the total error using a different approach at the time of invoking the Delta Method; our strategy uses this method at a more primitive level by considering the parameters of the linear representation of the process seen as a function of the ARMA coefficients. We show that discarding higher order terms on 1/ √ T , where T is the sample size used for estimation, the resulting formula for the total error can be approximated by a completely explicit expression that involves only the model parameters and the covariance matrix associated to the asymptotically normal estimator of the ARMA coefficients.
In subsection 3.2 we rederive the total error formula by H. Lütkepohl [L05, page 490] and show that it can be rewritten as explicitly as ours without using any stationarity hypothesis or Monte Carlo simulations. Moreover, we show that this formula coincides with the approximated one obtained in subsection 3.1 by discarding higher order terms on 1/ √ T .

The total error of finite sample based forecasting
Consider the causal and invertible ARMA(p, q) process {X t } determined by the equivalent relations and denote Ψ := {ψ 0 , ψ 1 , . . .}, Π := {π 0 , π 1 , . . .}. In Proposition 2.1 we studied forecasting for the process (3.1) when the parameters Ψ or Π of the model are known with total precision; in this section we suppose that these parameters are estimated by using a sample independent from the one that will be used for forecasting. A more preferable assumption would have been that the parameters Ψ are estimated based on the same sample that we intend to use for prediction, but exploiting only data up to the forecasting origin; Samaranayake [SH88] and Basu et al [BS86] have shown that many results obtained in the presence of the independence hypothesis remain valid under this more reasonable assumption. Under the independence hypothesis, the model coefficients Ψ or Π become random variables independent from the process X and the innovations ε. Moreover, we assume that these random variables are asymptotically normal, as for example in the case of maximum likelihood estimation of the ARMA coefficients.
For the sake of completeness, we start by recalling the Delta Method, that will be used profusely in the following pages. A proof and related asymptotic statements can be found in [Ser80].
Lemma 3.1 (Delta Method) Let β be an asymptotically normal estimator for the vector parameter β ∈ R n , that is, there exists a covariance matrix Σ such that where T is the sample size. Let f : R n → R m be a vector valued continuously differentiable function and let J f be its Jacobian matrix, that is, The next ingredient needed in the formulation of the main result of this section is the covariance matrix Σ Ξ P associated to the asymptotic normal character of the estimator Ξ P for the parameters Ξ P := (ψ 1 , . . . , ψ P , π 1 , . . . , π P ), for some integer P . This is spelled out in the following lemma whose proof is a straightforward combination of the Delta Method with the results in Section 8.8 of [Bro06].
Then, by the Delta Method we have that Details on how to algorithmically compute the Jacobian J Ξ P are provided in Appendix 7.2.
The next theorem is the main result in this section. Its proof can be found in the appendix.
. , x T } be a sample obtained as a realization of the causal and invertible ARMA(p,q) model in (3.1) using a preset I. In order to forecast out of this sample, we first estimate the parameters of the model Ψ = ψ 0 , ψ 1 , . . . , Π = { π 0 , π 1 , . . . } based on another sample ξ T that we assume to be independent of ξ T , using the maximum likelihood estimator β := Φ , Θ of the ARMA parameters. If we use the forecasting scheme introduced in Proposition 2.1, then: (i) The optimal multistep forecast X T +h for X T +h based on the information set F T generated by the sample ξ T and using the coefficients estimated on the independent sample ξ T is where r = max{p, q} andε t := t+r−1 j=0 π j x t−j .
(ii) The mean square forecasting error (MSFE) associated to this forecast is The first summand will be referred to as the characteristic forecasting error and the second one as the estimation based forecasting error. Notice that the characteristic error coincides with (2.6) and amounts to the forecasting error committed when the model parameters are known with the total precision.

On Lütkepohl's formula for the total forecasting error
As we already pointed out, H. Lütkepohl [L05, pages 97 and 490] has proposed formulas for VARMA models similar to the ones presented in Theorem 3.3 based on a different application of the Delta Method.
In this section, we rederive Lütkepohl's result in the ARMA context and show that it is identical to the approximated formula (3.5) presented in part (iii) of Theorem 3.3. In passing, this conveys that Lütkepohl's result can be explicitly formulated and computed using neither stationarity hypotheses nor Monte Carlo simulations.
The key idea behind Lütkepohl's formula is applying the Delta Method by thinking of the forecast X T +h in question as a differentiable function X T +h (β) of the model parameters β := (Φ , Θ ) . In order to develop further this idea, consider first the information sets F T := σ(ξ T ) and F T := σ(ξ T ) generated by two independent samples ξ T and ξ T of the same size. The sample ξ T is used for forecasting and hence F T determines the forecast X T +h (β) once the model parameters β have been specified. The sample ξ T is in turn used for model estimation and hence F T determines β. Consequently, the random variable X T +h ( β) is fully determined by F T and F T . In this setup, a straightforward application of the statement in Lemma 3.1 shows that which, as presented in the next result is enough to compute the total forecasting error.
Theorem 3.4 In the same setup as in Theorem 3.3, the total error associated to the forecast in (3.3) can be approximated by We refer to this expression as Lütkepohl's formula for the total forecasting error. Moreover: (i) Lutkepohl's formula coincides with the approximate expression for the total error stated in (3.5). In particular, the second summand in Lütkepohl's formula, which describes the contribution to the total error given by the estimation error, can be expressed as: (3.8) (ii) If we assume that the samples used for forecasting are second order stationary realizations of the model (3.1) and γ : Z → R is the corresponding time invariant autocovariance function, then the estimation error can be expressed as: (3.9)

Finite sample forecasting of temporally aggregated linear processes
The goal of this section is proposing a forecasting scheme for temporal aggregates based on using high frequency data for estimation purposes and the corresponding temporally aggregated model and data for the forecasting task. We show, using the formulas introduced in the previous section, that in some occasions this strategy can yield forecasts of superior quality than those based exclusively on high frequency data that are presented in the literature as the best performing option [L86b, L87, L09]. We start by recalling general statements about temporal aggregation that we need in the sequel. We then proceed by using various extensions of the results in Section 3 regarding the computation of total forecasting errors with estimated series in order to compare the performances of the schemes that we just indicated.

Temporal aggregation of time series
The linear temporal aggregation of time series requires the use of the elements provided in the following definition.
Definition 4.1 Given K ∈ N, X a time series, and w = (w 1 , ..., w K ) ∈ R K , define the K-period projection p K of X as .., X jK ∈ R K , and the corresponding temporally aggregated time series Y as The integer K is called the temporal aggregation period and the vector w the temporal aggregation vector. Notice that the aggregated time series Y is indexed using the time scale τ = mK, with m ∈ Z and its components are given by the X-aggregates X w t+K defined by Definition 4.1 can be reformulated in terms of the backward shift operator L as: and the indices of the components (Z j ) j∈Z of Z := There are four important particular cases covered by Definition 4.1, namely: (i) Stock aggregation (also called systematic sampling, skip-sampling, point-in-time sampling): it is obtained out of (4.1) or (4.2) by setting w = (0, 0, ..., 0, 1) .

Multistep approach to the forecasting of linear temporal aggregates
Let X be a time series and w = (w 1 , ..., w K ) an K-period aggregation vector. Given a finite time realization ξ T = {x 1 , ..., x T } of X such that T = M K with M ∈ N, we aim at forecasting the aggregate w 1 X T +1 + ... + w K X T +K . There are two obvious ways to carry this out; first, we can produce a multistep forecast X T +1 , ..., X T +K for X out of which we can obtain the forecast of the aggregate by setting X w T +K := w 1 X T +1 + ... + w K X T +K . Second, we can temporally aggregate X using (4.1) into the time series Y given by and produce a one-step forecast for Y . The following result recalls a well known comparison [AW72, L84, L86b, L89a] of the forecasting performances of the two schemes that we just described using the mean square characteristic error as an optimality criterion. In that setup, given an information set encoded as a σ-algebra, the optimal forecast is given by the conditional expectation with respect to it [Ham94, page 72]. Given a time series X, we will denote in what follows by σ X T the information set generated by a realization ξ T = {x 1 , . . . , x T } of length T of X and the preset I used to produce it; more specifically Proposition 4.2 Let X be a time series and w = (w 1 , ..., w K ) a K-period aggregation vector. Let Y = I w • p K (X) be the corresponding temporally aggregated time series. Let T = M K with M, T ∈ N and consider F T = σ X T , F M = σ Y M the information sets associated to two histories of X and Y of length T and M , respectively, related to each other by temporal aggregation. Then: Remark 4.3 The inequality (4.4) has been studied in detail in the VARMA context by H. Lütkepohl [L86b, L87, L09] who has fully characterized the situations in which the two predictors are identical and hence have exactly the same performance. This condition is stated and exploited in Section 5, where we illustrate with specific examples the performance of the forecasting scheme that we present in the following pages.
In the next two results we spell out the characteristic and the total errors associated to a multistep approach to the forecast of linear aggregates. The characteristic error is given in Proposition 4.4 and the total error is provided in Theorem 4.7 under the same independence hypothesis between the samples used for estimation and forecasting that were already invoked in Theorem 3.3.
Proposition 4.4 Let X be a time series model as in (3.1), r = max{p, q}, K a temporal aggregation period, w = (w 1 , ..., w K ) an aggregation vector, and F T := σ X T the information set generated by w i X T +i based on F T using the forecasting scheme in Proposition 2.1. Then: (i) The forecast X w T +K is given by: The corresponding mean square forecasting characteristic error is: (4.6) Example 4.5 Forecast of stock temporal aggregates. It is a particular case of the statement in Proposition 4.4 obtained by taking w = (0, ..., 0, 1) . In this case This shows that the forecast of the stock temporal aggregate coincides with the K-multistep forecast of the original time series. Consequently, it is easy to see by using (4.6) and (2.6) that MSFE X w T +K = MSFE X T +K . Example 4.6 Forecast of flow temporal aggregates. It is a particular case of the statement in Proposition 4.4 obtained by taking w = (1, ..., 1) . In this case Consequently, Theorem 4.7 (Multistep forecasting of linear temporal aggregates) Consider a sample ξ T = {x 1 , . . . , x T } obtained as the realization of a model of the type (3.1) using the preset I. In order to forecast out of this sample, we first estimate the parameters of the model Ψ = ψ 0 , ψ 1 , . . . and Π = { π 0 , π 1 , . . . } based on another sample ξ T of the same size that we assume to be independent of ξ T . Let w = (w 1 , . . . , w K ) be an aggregation vector and let X w T +K be the forecast of the aggregate X w T +K := K h=1 w h X T +h based on F T := σ (I ∪ ξ T ) using Proposition 4.4 and the estimated parameters Ψ, Π. Then: (i) The optimal forecast X w T +K for X w T +K given the samples ξ T and ξ T is where r = max{p, q} andε t = t+r−1 j=0 π j x t−j .
(ii) The mean square forecasting error associated to this forecast is where A, B, C are the matrices with components given by and σ 2 < w, A char w > is the characteristic forecasting error in part (ii) of Proposition 4.4.
(iii) Let Ξ P := ψ 1 , . . . , ψ P (K) , π 1 , . . . , π P (K) , with P (K) = T + K − 1 + r. Using the notation introduced in Lemmas 3.1 and 3.2 and discarding higher order terms in 1/ √ T , the MSFE in (4.8) can be approximated by Remark 4.8 In order to compute the total error in (4.12) it is necessary to determine the covariance matrix Σ Ξ P . By Lemma 3.2, it can be obtained out of the covariance Σ β matrix associated to the estimator of the ARMA parameters combined with the Jacobian J Ξ P . Details on how to algorithmically compute this Jacobian are provided in Appendix 7.2.
Remark 4.9 Notice that all the matrices involved in the statement of Theorem 4.7 are symmetric except for B and F .

A hybrid forecasting scheme using aggregated time series models
In the previous subsection we presented a forecasting method for linear temporal aggregates based exclusively on the use of high frequency data and models. The performance of this approach has been compared in the literature [L86b, L87] with the scheme that consists of using models estimated using the aggregated low frequency data; as it could be expected due to the resulting smaller sample size, this method yields a performance that is strictly inferior to the one based on working in the pure high frequency setup.
In this section, we introduce and compute the performance of a hybrid recipe that consists of estimating first the model using the high frequency data so that we can take advantage of larger sample sizes and of the possibility of updating the model as new high frequency data become available. This model and the data used to estimate it are subsequently aggregated and used for forecasting. We will refer to this approach as the hybrid forecasting scheme. The main goal in the following pages is writing down explicitly the total MSFE associated to this forecasting strategy so that we can compare it using Theorem 4.7 with the one obtained with the method based exclusively on the use of high frequency data and models. In the next section we use the resulting formulas in order to prove that there are situations in which the hybrid forecasting scheme provides optimal performance for various kinds of temporal aggregation.
The main tool at the time of computing the MSFE associated to the hybrid scheme is again the use of the Delta Method [Ser80] in order to establish the asymptotic normality of the estimation scheme resulting from the combination of high frequency data with the subsequent model temporal aggregation.
In order to make this more explicit, consider a time series model X determined by the parameters β X for which an asymptotically normal estimator β X is available, that is, there exists a covariance matrix Σ βX such that with T being the sample size on which the estimation is based. Now, let K ∈ N be an aggregation period, w ∈ R K an aggregation vector, and Y := I w • p K (X) the linear temporally aggregated process corresponding to X and w.
Proposition 4.10 In the setup that we just described, suppose that the temporally aggregated process Y is also a parametric time series model and that the parameters β Y that define it can be expressed as a C 1 function β Y (β X ) of the parameters β X that determine X. Using the estimator β X , we can construct an estimator β Y for β Y based on disaggregated X samples by setting where T is the disaggregated sample size and the Jacobian matrix corresponding to the function β Y (β X ).
Once the model temporal aggregation function and its Jacobian have been determined, this proposition can be used to formulate an analog of Theorem 4.7 for the hybrid forecasting scheme by mimicking the proof of Theorem 3.3; the only necessary modification consists of replacing the asymptotic covariance matrix Σ βX of the estimator for the disaggregated model by that of the aggregated model Σ βY obtained using Proposition 4.10.
We make this statement explicit in the following theorem and then describe how to compute the model aggregation function β Y (β X ) and its Jacobian J βY in order to make it fully functional. The construction of these objects is carried out in the ARMA context where the model aggregation question has already been fully studied. Even though all necessary details will be provided later on in the section, all we need to know at this stage in order to state the theorem is that the linear temporal aggregation of an ARMA(p,q) model is another ARMA(p, q * ) model where K is the temporal aggregation period, K * is the index of the first nonzero entry of the aggregation vector, and the symbol · denotes the integer part of its argument. We emphasize that if the innovations that drive the disaggregated model are independent with variance σ 2 , this is not necessarily the case for the resulting aggregated model, whose innovations may be only uncorrelated with a different variance σ 2 * , making it into a so called weak ARMA model. We forecast the value of the temporal aggregate Y M +1 = X w T +K out of the sample η M by first estimating the parameters β X of the model X using another disaggregated sample ξ T of the same size, that we assume to be independent of ξ T . Let β Y (β X ) be the function that relates the ARMA parameter values of the disaggregated and the aggregated model and let J βY be its Jacobian. Consider the ARMA(p, q * ) model, with q * as in (4.15), determined by the parameters β Y := β Y ( β X ). Then: (i) The optimal forecast Y M +1 of the temporal aggregate Y M +1 = X w T +K based on the information set F M := σ (I K ∪ η M ) using Proposition 4.4 and the estimated parameters β Y is given by
(ii) Let Ξ P := (ψ 1 , . . . , ψ P , π 1 , . . . , π P ) with P = T + r * . Discarding higher order terms in 1/ √ T , the MSFE corresponding to the forecast (4.16) can be approximated by where σ 2 * is the variance of the innovations of the aggregated ARMA(p, q * ) model Y and Σ Ξ P is the covariance matrix given by with J βY the Jacobian matrix corresponding to the function β Y (β X ) and J Ξ P the Jacobian of As we announced above, we conclude this section by describing in detail the parameters aggregation function β Y (β X ) and its Jacobian J βY , so that all the ingredients necessary to apply formula (4.17) are available to the reader. In order to provide explicit expressions regarding these two elements, we provide a brief review containing strictly the results on the temporal aggregation of ARMA processes that are necessary for our discussion; for more ample discussions about this topic we refer the reader to [AW72, Tia72, Bre73, TW76, Wei79, SW86, Wei06, SV08] and references therein.
In order to simplify the discussion and to avoid hidden periodicity phenomena, we place ourselves in a generic situation in which the autoregressive and moving average polynomials of the model that we want to temporally aggregate have no common roots and all roots are different (see [Wei06] for the general case). Consider now a K-period aggregation vector w = (w 1 , ..., w K ) . Our first goal is to find polynomials T (z) and Φ * (z) that satisfy with Π K : j∈Z R −→ j∈Z R as in (4.3). The intuition behind (4.18) is that for any time series X, its temporal aggregation In other words, the polynomial T (L), that we will call the temporal aggregation polynomial transforms the AR polynomial for X into an AR polynomial for Y in the aggregated time scale units. Let T (L) = t 0 + t 1 L + ... + t n L n and Φ * L K = 1 − φ * 1 L K − ... − φ * c L Kc be the unknown polynomials. Equation (4.18) can be written in matrix form as [Bre73]: where w = (w K , w K−1 , ..., w 1 ) is the reflection of w. We start by determining the unknown values n and c using the two following dimensional restrictions: • Since A is a matrix of size (n + p + 1) × (p + 1), Z and D are vectors of size p + 1 and cK + K, respectively, and we have AZ = D then, necessarily n + p + 1 = cK + K. (4.21) • The system AZ = D contains n + p + 1 equations that need to coincide with the number of unknowns, that is, the n + 1 + c coefficients (t 0 , t 1 , . . . , t n ) and (φ * 1 , . . . , φ * c ) of the polynomials T (z) and Φ * (z), respectively. Consequently, n + p + 1 = n + c + 1. (4.22) The conditions (4.21) and (4.22) yield c = p, n = pK + K − p − 1 = (p + 1) (K − 1) . (4.23) The first condition in (4.23) shows that the autoregressive order does not change under temporal aggregation.
Let now K * ≤ K be the index of the first nonzero component in the vector w. This implies that w is a vector of the form w = (w K , w K−1 , ..., w K * , 0, ..., 0) with w K , w K−1 , . . . , w K * = 0. Since (4.20) is a matrix representation of the polynomial equality in (4.18), we hence have that deg (T (L) Φ (L)) = deg (D (L)), where D (L) is the polynomial associated to the vector in the right hand side of (4.20). It is clear that deg (D (L)) = K (p + 1) − K * ; as deg (T (L) Φ (L)) = n + p, the degree n of T (L) is therefore n = K (p + 1) − p − K * . (4.24) Solving the polynomial equalities (4.20), we have found polynomials T (L) and Φ * L K such that the temporally aggregated time series Y satisfies Our goal now is showing the existence of a polynomial Θ * L K and a white noise {ε * } ∼ WN 0, σ 2 * such that T (L) Θ (L) ε lK = Θ * L K ε * lK for any l ∈ Z. (4.26) This equality, together with (4.25) shows that the temporally aggregated process Y out of the ARMA process X is a weak ARMA process, as it satisfies the relation is a polynomial of degree p, and Θ * (B) a polynomial of degree n + q K whose coefficients will be determined in the following paragraphs. We recall that the symbol · denotes the integer part of its argument. Indeed, by (4.24) we have that deg (T (L) Θ (L)) = K (p + 1)−p−K * +q = n + q. Additionally, by (4.25) The coefficients of the polynomial Θ * (B) are obtained by equating the autocovariance functions of the processes on both sides of (4.30) at lags 0, K, 2K, ..., q * K, which provide q * + 1 nonlinear equations that determine uniquely the q * + 1 unknowns corresponding to the coefficients θ * 1 , . . . , θ * q * of the polynomial Θ * (B) and the variance σ 2 * of the white noise of the aggregated model. In order to explicitly write down the equations that we just described, let us denote C (L) := T (L) Θ (L) as in (4.30) and set C (L) = n+q i=0 c i L i . Let now γ and Γ be the autocovariance functions of the MA (n + q) and MA (q * ) processes Consequently, the coefficients of the polynomial Θ * (B) and the variance σ 2 * are uniquely determined by the q * + 1 equations The equations (4.34) can be written down in matrix form, which is convenient later on at the time of spelling out the Jacobian of the aggregation function. Indeed, we can write: where the bars over the polynomials in the previous expressions denote the corresponding coefficient vectors, that is, given a polynomial q (x) = n i=1 a i x i , then q (x) = (a 0 , a 1 , . . . , a n ) . Additionally S i is the lower ith-shift matrix, that is, In conclusion, if we denote β X = (Φ, Θ) and β Y = (Φ * , Θ * ), the construction that we just examined shows that The function Φ * (Φ) is given by the solution of the polynomial equalities (4.20) and Θ * (Φ, Θ) by the coefficients (t 0 , t 1 , . . . , t n ) determined by (4.20) and the solutions of (4.37).
The Jacobian J βY of the temporal aggregation function β Y (β X ). The goal of the following paragraphs is the computation of the Jacobian J βY of the function β Y (β X ) = (Φ * (Φ) , Θ * (Φ, Θ)) in (4.38). We first compute the Jacobian of the function Φ * (Φ) by taking derivatives with respect to the components of the vector Φ on both sides of the equations (4.20) that determine Φ * (Φ). This results in the following p matrix equations that will be needed later on in the computation of the remaining blocks of the Jacobian. Given that there is no Θ dependence on the function Φ * (Φ), the (1, 2)-block of the Jacobian is a zero matrix of size p × q * . The remaining two blocks are computed by using the function Θ * (Φ, Θ) uniquely determined by the equations (4.37). Its derivatives are obtained out of a new set of equations resulting from the differentiation of both sides of this relation, namely, (4.41) j = 0, 1, . . . , q * , i = 1, . . . , p + q.
In the previous section we proposed a new hybrid scheme for the forecasting of temporal aggregates coming from ARMA processes. We recall that this strategy consists of first using high frequency disaggregated data for estimating a model; then we temporally aggregate both the data and the model, and finally we forecast using these two ingredients. As we announced in the introduction, there are situations in which our strategy is optimal with respect to the total error, that is, the predictor constructed following this procedure performs better than the one based exclusively on high frequency data and the underlying disaggregated model. In this section we give a few examples of specific models for which our scheme provides optimal efficiency of prediction. Before we proceed, we introduce abbreviations for the various predictors that we will be working with: (i) Temporally aggregated multistep predictor (TMS predictor): this is the denomination that we use for the forecast of the aggregate that is constructed out of the disaggregated data and the underlying disaggregated model estimated on them.
(ii) Temporally aggregated predictor (TA predictor): this is the forecast based on use of the temporally aggregated sample and a model estimated on it.
(iii) Hybrid predictor (H predictor): this is the predictor introduced in Section 4.3 whose performance is spelled out in Theorem 4.11. In this scheme, a first model is estimated on the disaggregated high frequency data sample, then the data and the model are temporally aggregated with an aggregation period that coincides with the forecasting horizon; finally, both the temporally aggregated model and the sample are used to produce a one-step ahead forecast that amounts to a prediction of the aggregate we are interested in.
(iv) Optimal hybrid predictor (OH predictor): this predictor is constructed by taking the multistep implementation of the H predictor that yields the smallest total error. More explicitly, suppose that the aggregate that we want to forecast involves K time steps; let {K 1 , . . . , K r } be the positive divisors of K and {C 1 , . . . , C r } the corresponding quotients, that is, K = K i C i for each i ∈ {1, . . . , r}. There are aggregation schemes (stock and flow for example) for which a K-temporal aggregate can be obtained as the aggregation of C i K i -temporal aggregates, for all i ∈ {1, . . . , r}. The total error associated to the forecasting of these aggregates using a multistep version of the H predictor obviously depends on the factor K i used. The OH predictor is the one associated to the factor K i that minimizes the total error.
As we already mentioned, the forecasting performance of the TMS predictor is always superior or equal than that of the TA predictor when we take into account only the characteristic error, and it is strictly superior when the total error is considered. In view of these results and given that the H and the OH predictors carry out the forecasting with temporally aggregated data, they are going to produce worse characteristic errors than their TMS counterpart; hence, the only way in which the H and OH predictors can be competitive in terms of total error performance is by sufficiently lowering the estimation error. In order to check that they indeed do so, we will place ourselves in situations that are particularly advantageous in this respect and will choose models for which the TMS and the TA predictors have identical characteristic errors and hence it is only the estimation error that makes a difference as to the total error. The linear models for which this coincidence of characteristic errors takes place have been identified in the works of H. Lütkepohl [L86b, L87, L09] via the following statement.
Theorem 5.1 (Lütkepohl) Let X t = ∞ i=0 ψ i ε t−i be a linear causal process and let w = (w 1 , . . . , w K ) ∈ R K be a K-period aggregation vector. Then the TMS and TA predictors for the K-temporal aggregate determined by w have identical associated characteristic errors if and only if the following identity holds: (5.1) The equality (5.1) is satisfied for both stock (w = (0, . . . , 0, 1) ) and flow aggregation (w = (1, . . . , 1) ) if {X t } is a purely seasonal process with period K, that is, Given a specific model we want to compare the performances of the H and the TMS predictors for a variety of forecasting horizons. Given that condition (5.1) is different for each aggregation period K and cannot be solved simultaneously for several of them, we will content ourselves either with approximate solutions that are likely to produce very close H and TMS characteristic errors for several periods K or with exact solutions that provide exactly equal errors for only a prescribed aggregation period. The following points describe how we have constructed examples following the lines that we just indicated: • We first choose the orders p and q of the disaggregated ARMA(p,q) model that we want to use as the basis for the example.
• We fix an aggregation period K and a number n of parameters ψ i for which the equation (5.1) will be solved. The choice of p and q imposes a minimal number n min = q − p + 1.
• We determine a vector Ψ * = (ψ 0 , ψ 1 , ..., ψ n−1 ) that consists of the n first components of the set Ψ = {ψ 0 , ψ 1 , ...} that satisfies condition (5.1). We emphasize that in general this condition does not determine uniquely the vector Ψ * and that arbitrary choices need to be made. The vector Ψ * is a truncation at order n − 1 of the MA representation of the ARMA process that we want to construct.
• We conclude the construction of the ARMA(p,q) model that we are after by designing either an AR(p) polynomial Φ consistent with causality or a MA(q) polynomial Θ consistent with invertibility. Then: -In the first case, the required model is given by -In the second case, the required model is given by In both cases, the MA and AR polynomials that are obtained in this way have to be checked regarding invertibility and causality, respectively. Additionally, the finite truncation of Ψ is likely to give rise to common roots between the AR and MA polynomials in (5.3) or in (5.4) which may make necessary a slight perturbation of the coefficients in order to be avoided.
• We emphasize that the resulting ARMA model satisfies (5.1) only approximately and hence the characteristic errors of the two predictors will be not identical but just close to each other for the specific aggregation period K used. For pure MA models no truncation is necessary and hence exact equality can be achieved.

Stock aggregation examples
In the particular case of stock temporal aggregation, condition (5.1) is written as: We now consider the truncated vector Ψ * with n components, that is, Ψ * = (ψ 0 , . . . , ψ n−1 ) . Then, the truncated version of (5.5) is: where the symbol · denotes the integer part of its argument. We now provide a few examples of models whose specification is obtained following the approach proposed in the previous subsection and the relation (5.6).
This system of nonlinear equations has many solutions. We choose one of them by setting ψ j = 0, for j = 1, . . . , 9, and ψ 10 = 0.3. This way we can trivially determine a MA(10) model which satisfies exactly the relation (5.5) by taking θ j = 0 for j = 1, . . . , 9 and θ 10 = 0.3. Figure 1 shows the values of the characteristic errors for different values of the forecasting horizon for the TMS predictor, the H predictor, and the OH predictor. For the horizon h = 2, the values of the characteristic errors of all the predictors coincide, which is a consequence of the fact that the model has been constructed using the relation (5.5) with K = 2. Moreover, it is easy to see by looking at (5.2), that the particular choice of MA coefficients that we have adopted ensures that the resulting model is seasonal for the periods 2, 5, and 10; this guarantees that (5.5) is also satisfied for the corresponding values of K and hence there is coincidence for the characteristic errors at those horizons too.
The total errors for a sample size of T = 50 are then computed using the formulas presented in sections 3 and 4. The corresponding results are also plotted in Figure 1 for the different forecasting schemes. This plot shows that for several forecasting horizons both the H and the OH predictors perform better than the TMS predictor.
A quick inspection of this plot reveals another interesting phenomenon consisting on the decrease of the total error associated to the three predictors as the forecasting horizon increases; this feature is due to the decrease of the estimation error using these forecasting schemes as the horizon becomes longer. The characteristic errors for the H and OH predictors do not increase monotonically with the forecasting horizon either; however, in this case, this is due to the fact that for each value of the forecasting horizon, these predictors are constructed using a different model since the aggregation period changes and hence so does the aggregated model used for forecasting.
In conclusion, in this particular example, both the H and the OH predictors exhibit a better forecasting performance than the TMS predictor and, additionally, the results regarding the OH predictor help in making a decision on what is the best possible aggregation period to work with in order to minimize the associated total forecasting error.  Figure 2 shows the errors with respect to the forecasting horizon for all the predictors as in the previous example. The H and the OH predictors perform better than the TMS predictor for h = 3, 6, 9, 10. Additionally the OH predictor performs better than the H predictor for h = 4. Taking into account the initial choice of K = 3 when constructing the example, it becomes clear why the characteristic errors associated to the H and the OH predictors are very close to those associated to the TMS predictor for horizons h that are multiples of 3.  Figure 3 shows that the H and the OH predictors have equal associated total errors and exhibit a better forecasting efficiency than the TMS predictor for all forecasting horizons except at h = 2. The initial choice of K = 4 at the model construction stage results in the fact that for h = 4 the values of the characteristic errors associated to the H and the OH predictors are very close to the one committed by the TMS predictor.

Flow aggregation examples
In the particular case of flow temporal aggregation, condition (5.1) can be written as: We now consider the truncation Ψ * of Ψ with n components, that is, Ψ * = (ψ 0 , . . . , ψ n−1 ). Then, the truncated version of (5.7) can be expressed as: where symbol · denotes the integer part of its argument. We now provide a few examples of models whose specification is obtained following the approach described in the beginning of the section and the relation (5.8).
Let p = 0, q = 10, n = n min = 11, and let K = 2. Then the expression (5.8) reads and consequently ψ 0 = 1, (1 + ψ 1 )(ψ i + ψ i+1 ) = ψ i+1 + ψ i+2 , for i = 1, . . . , 9, and ψ i = 0 for i ≥ n. (5.9) A solution for these equations is given by the choice ψ j = 0, for j = 1, . . . , 9, and ψ 10 = 0.3. Since the order of the AR polynomial is zero, the method that we proposed determines uniquely in this case the MA(10) polynomial that we are after with Θ = (0, 0, 0, 0, 0, 0, 0, 0, 0, 0.3) . The evolution of the forecasting errors versus the forecasting horizon in plotted in the Figure 4. Both the H and the OH predictors perform better than the TMS predictor. For h = 4 the OH predictor has the smallest total error among the three predictors.  Figure 5 shows the corresponding characteristic and total errors for the three predictors. The H and the OH predictors exhibit better performance than the TMS predictor for horizons h = 2, 4, 5, 6, 7. Remark 5.7 The model that we just presented can be used to illustrate the fact that the construction method that we presented in this sections is not the unique source of examples in which the H and the OH predictors perform better than the TMS scheme. Indeed, as it can be seen in Figure 6, the very same ARMA(3,10) model prescription used in the context of stock aggregation also shows this feature even though it has not been obtained by finding a solution of the equation (5.5).
Remark 5.8 Notice that when the forecasting horizon h equals one all predictors coincide because there is no temporal aggregation and hence they obviously have the same errors associated.

Conclusions
In this work we have carried out a detailed study of the total error committed when forecasting with one dimensional linear models by minimizing the mean square error. We have introduced a new hybrid scheme for the forecasting of linear temporal aggregates that in some situations shows optimal performance in comparison with other prediction strategies proposed in the literature. We work in a finite sample context. More specifically, the forecasting is based on the information set generated by a sample and a model whose parameters have been estimated on it and we avoid the use of second order stationarity hypotheses or the use of time independent autocovariance functions.
In this setup, we provide explicit expressions for the forecasting error that incorporate both the error incurred in due to the stochastic nature of the model (we call it characteristic error) as well as the one associated to the sample based estimation of the model parameters (estimation error). In order to derive these expressions we use certain independence and asymptotic normality hypotheses that are customary in the literature; our main contribution consists of providing expressions for the total error that do not require neither stationarity on the samples used nor Monte Carlo simulations to be evaluated. We subsequently use these formulas to evaluate the performance of a new forecasting strategy that we propose for the prediction of linear temporal aggregates; we call it hybrid scheme. This approach consists of using high frequency data for estimation purposes and the corresponding temporally aggregated data and model for forecasting. This scheme uses all the information available at the time of estimation by using the bigger sample size provided by the disaggregated data, and allows these parameters to be updated as new high frequency data become available. More importantly, as we illustrate with various examples, in some situations the total error committed using this scheme is smaller than the one associated to the forecast based on the disaggregated data and model; in those cases our strategy becomes optimal. As the increase in performance obtained with our method comes from minimizing the estimation error, we are persuaded that this approach to forecasting may prove very relevant in the multivariate setup where in many cases the estimation error is the main source of error. (ii) Suppose first that the innovations {ε t } are IID(0, σ 2 ). Then the forecast X T +h that minimizes the mean square forecasting error E X T +h − X T +h 2 is given by the conditional expectation (see for example [Ham94], page 72): as required.
When {ε t } is WN(0, σ 2 ) our goal is finding the linear combination Given that by (7.1), the elements X T −i can be written as a linear combination of the elements in T −i , this task is equivalent to finding the vector b = (b 0 , ..., b T −1+r ) that minimizes the function Hence, in order to minimize the function S (b 0 , ..., b T −1+r ) we compute the partial derivatives ∂S/∂b i and we set them to zero, which shows that the optimal values are attained when b i = ψ i+h . Consequently, the optimal linear forecast is given by X T +h = T −1+r i=0 ψ i+h ε T −i , as required.
(iii) We first compute X T +h − X T +h . By (2.5) and (7.1) we have (iv) Given the model Φ (L) X = Θ (L) ε, we have We first recall that by (2.4) we have that σ ξ T = σ T =: F T . We now project both sides of this equality onto the information set F T by thinking of this σ-algebra as σ ξ T for the left hand side projection and as σ T for the right hand side. We obtain: In the presence of white noise innovations, the conditional expectation in the previous equality should be replaced by a linear projection.
The causality and invertibility hypotheses on the ARMA process we are working with, guarantee that for any P ≥ max {p, q} there exist polynomials Ψ P (z), Π P (z) of order P uniquely determined by the relations: which are equivalent to These polynomial relations determine the functions Ψ P (Φ, Θ), Π P (Φ, Θ) needed in the computation of the Jacobian. We now rewrite (7.2) and (7.3) as Φ(z)Ψ P (Φ, Θ)(z) = Θ(z), (7.4) Φ(z) = Π P (Φ, Θ)(z)Θ(z). (7.5) If we take derivatives with respect to θ j and φ i , j ∈ {1, . . . , q}, i ∈ {1, . . . , p} on both sides of (7.4), we obtain a set of p + q polynomial equations: that determine uniquely the corresponding entries of the Jacobian due to the invertibility of Φ(z). At the same time, taking derivatives on both the right and left hand sides of (7.5) with respect to θ j and φ i , we obtain another set of p + q polynomial equations ∂Φ ∂θ j (z) = ∂Π P (φ, θ) ∂θ j (z) + Π P z j , j ∈ {1, . . . , q} , that determine uniquely the corresponding entries of the Jacobian due to the invertibility of Θ(z).

Proof of Theorem 3.3
(i) It is a straightforward consequence of the independence hypothesis between the samples ξ T and the ξ T , and part (ii) of Proposition 2.1.
(ii) By (3.1) and part (i) of Proposition 2.1 we have that In order to compute this error notice thatε T +h−i can be rewritten in terms of the original innovations asε Hence, (iii) By part (i) the forecast X T +h is given by According to the statement (3.2), both ψ i and π j can be asymptotically written as with r i and t j as Gaussian random variables of mean 0 and variances (Σ Ξ ) i,i and (Σ Ξ ) j+P,j+P , respectively. Consequently by (7.7) and (7.6) We now recall that and we eliminate in this expression the term that decays as 1/T; we hence approximate X T +h as Using this approximation we compute now the MSFE:

Proof of Proposition 4.4
First, we notice that by (4.2) we have F M ⊂ F T .
The same relation guarantees that X w T +k = Y M +1 . Hence the result is a consequence of the following general fact: Lemma 7.1 Let z be a random variable in the probability space (Ω, P, F). Let F * be a sub-sigma algebra of F , that is, F * ⊂ F. Then

Proof of Proposition 4.4
Part (i) is a straightforward consequence of (2.5). Regarding (ii), we first have that Consequently, 7.7 Proof of Theorem 4.7 (i) It is a straightforward consequence of part (i) in Theorem 3.3.
(iii) We recall that by (4.7), the forecast X w T +K is given by According to (3.2), both ψ i and π j can be asymptotically written as with r i and t j Gaussian random variables of mean 0 and variances (Σ Ξ P ) i,i and (Σ Ξ P ) j+P (K),j+P (K) , respectively. Consequently, (7.20) We now recall that and we eliminate in (7.20) the term that decays as 1/T. We hence approximate X w T +K as (7.21) Using this approximation we now compute the MSFE: Using Lemma 3.2, this expression can be asymptotically approximated by: where A char , D, F , G are matrices whose components are given by