Autoregressive approaches to import-export time series I: basic techniques

This work is the first part of a project dealing with an in-depth study of effective techniques used in econometrics in order to make accurate forecasts in the concrete framework of one of the major economies of the most productive Italian area, namely the province of Verona. In particular, we develop an approach mainly based on vector autoregressions, where lagged values of two or more variables are considered, Granger causality, and the stochastic trend approach useful to work with the cointegration phenomenon. Latter techniques constitute the core of the present paper, whereas in the second part of the project, we present how these approaches can be applied to economic data at our disposal in order to obtain concrete analysis of import--export behavior for the considered productive area of Verona.


Introduction
The analysis of time series data constitutes a key ingredient in econometric studies. Last years have been characterized by an increasing interest toward the study of econometric time series. Although various types of regression analysis and related forecast methods are rather old, the worldwide financial crisis experienced by markets starting from last months of 2007, and which is not yet finished, has put more attention on the subject. Moreover, analysis and forecast problems have become of great momentum even for medium and small enterprizes since their economic sustainability is strictly related to the propensity of a bank to give credits at reasonable conditions.
In particular, great efforts have been made to read economic data not as monads, but rather as constituting pieces of a whole. Namely, new techniques have been developed to study interconnections and dependencies between different factors characterizing the economic history of a certain market, a given firm, a specified industrial area, and so on. From this point of view, methods such as the vector autoregression, the cointegration approach, and the copula techniques have been benefitted by new research impulses.
A challenging problem is then to apply such instruments in concrete situations and the problem becomes even harder if we take into account the economies are hardly hit by the aforementioned crisis. A particularly important case study is constituted by a close analysis of import-export time series. In fact, such an information, spanning from countries to small firms, has the characteristic to provide highly interesting hints for people, for example, politicians or CEOs, to depict future economic scenarios and related investment plans for the markets in which they are involved.
Exploiting precious economic data that the Commerce Chamber of Verona Province has put at our disposal, we successfully applied some of the relevant approaches already cited to find dependencies between economic factor characterizing the Province economy and then to make effective forecasts, very close to the real behavior of studied markets.
For completeness, we have split our project into two parts, namely the present one, which aims at giving a self-contained introduction to the statistical techniques of interest, and the second one, where the Verona import-export case study have been treated in detail.
In what follows, we first recall univariate time series models, paying particular attention to the AR model, which relates a time series to its past values. We will explain how to make predictions, by using these models, how to choose the delays, for example, using the Akaike and Bayesian information crtiteria (AIC, resp. BIC), and how to behave in the presence of trends or structural breaks. Then we move to the vector autoregression (VAR) model, in which lagged values of two or more variables are used to forecast future values of these variables. Moreover, we present the Granger causality, and, in the last part, we return to the topic of stochastic trend introducing the phenomenon of cointegration.

Univariate time-series models
Univariate models have been widely used for short-run forecast (see, e.g., [6, Examples of Chapter 2]. In what follows, we recall some of these techniques, focusing ourselves particularly on the analysis of autoregressive (AR) processes, moving average (MA) processes, and a combination of both types, the so-called ARMA processes; for further details, see, for example, [3,2,8] and references therein.
The observation on the time-series variable Y made at date t is denoted by Y t , whereas T ∈ N + indicates the total number of observations. Moreover, we denote the jth lag of a time series {Y t } t=0,...,T by Y t−j (the value of the variable Y j periods ago); similarly, Y t+j denotes the value of Y j periods to the future, where, for any fixed t ∈ {0, . . . , T }, j is such that j ∈ N + , t − j ≥ 0, and t + j ≤ T . The jth autocovariance of a series Y t is the covariance between Y t and its jth lag, that is, autocovariance j = σ j := cov(Y t , Y t−j ), whereas the jth autocorrelation coefficient is the correlation between Y t and Y t−j , thats is, . When the average and variance of a variable are unknown, we can estimate them by taking a random sample of n observations. In a simple random sample, n objects are drawn at random from a population, and each object is equally likely to be drawn. The value of the random variable Y for the ith randomly drawn object is denoted Y i . Because each object is equally likely to be drawn and the distribution of Y i is the same for all i , the random variables Y 1 , . . . , Y n are independent and identically distributed (i.i.d.). Given a variable Y , we denote by Y its sample average with respect to the n observations Y 1 , . . . , Y n , thats is, The jth autocovariances, resp. autocorrelations, can be estimated by the jth sample autocovariances, resp. autocorrelations, as follows: where Y j+1,T denotes the sample average of Y t computed over the observations t = j + 1, . . . , T . Concerning forecast based on regression models that relates a time series variable to its past values, for completeness, we shall start with the first-order autoregressive process, namely the AR(1) model, which uses Y t−1 to forecast Y t . A systematic way to forecast is to estimate an ordinary least squares (OLS) regression. The OLS estimator chooses the regression coefficients so that the estimated regression line is as close as possible to the observed data, where the closeness is measured by the sum of the squared mistakes made in predicting Y t given Y t−1 . Hence, the AR(1) model for the series Y t is given by where β 0 and β 1 are the regression coefficients. In this case, the intercept β 0 is the value of the regression line when Y t−1 = 0, the slope β 1 represents the change in Y t associated with a unit change in Y t−1 , and u t denotes the error term whose nature will be later clarified. Let us assume that the value Y t0 of the time series Y t at initial time t 0 is given; then Y t0+1 = β 0 + β 1 Y t0 + u t0+1 , so that iterating relation (1) up to order τ > 0 , we get Hence, taking t = t 0 + τ with t 0 = 0, we obtain A time series Y t is called stationary if its probability distribution does not change over time, that is, if the joint distribution of (Y s+1 , Y s+2 , . . . , Y s+T ) does not depend on s; otherwise, Y t is said to be nonstationary. In (2), the process Y t consists of both time-dependent deterministic and stochastic parts, and, thus, it cannot be stationary. Formally, the process with stochastic initial conditions results from (2) if and only if |β 1 | < 1. It follows that if lim t0→−∞ Y t0 is bounded, then, as t 0 → −∞, we have where we have used that E[u t u s ] = 0 for t = s and |β 1 | < 1. Hence, both the mean and variance are constants, and thus the covariances are given by The previous AR(1) can be generalized by considering arbitrary but finite order p > 1.
In particular , an AR(p) process can be described by the equation where β 0 , . . . , β p are constants, whereas u t is the error term represented by a random variable with zero mean and variance σ 2 > 0. Using the lag operator, we can rewrite In such a framework, it is standard to assume that the following four properties hold (see, e.g., [7,Chap. 14.4]): • u t has conditional mean zero, given all the regressors, that is, E(u t |Y t−1 , Y t−2 , . . .) = 0, which implies that the best forecast of Y t is given by the AR(p) regression.
• Y i has a stationary distribution, and Y i , Y i−j are assumed to become independent as j gets large. If the time-series variables are nonstationary, then the forecast can be biased and inefficient, or conventional OLS-based statistical inferences can be misleading.
• All the variables have nonzero finite fourth moments.
• There is no perfect multicollinearity, namely it is not true that, given a certain regressor, it is a perfect linear function of the variables.

Forecasts
In this section, we show how the previously introduced class of models can be used to predict the future behavior of a certain quantity of interest. If Y t follows the AR(p) model and β 0 , β 1 , . . . , β p are unknown, then the forecast of Y T +1 is given by Forecasts must be based on estimates of the coefficients β i by using the OLS estimators based on historical data.
Then such a forecast refers to some data beyond the data set used to estimate the regression, so that the data on the actual value of the forecasted dependent variable are not in the sample used to estimate the regression. Forecasts and forecast error pertain to "out-of-sample" observations. The forecast error is the mistake made by the forecast; this is the difference between the value of Y T +1 that actually occurred and its forecasted value forecast error : The root mean squared forecast error RMSFE is a measure of the size of the , and it is characterized by two sources of error: the error arising because future values of u t are unknown and the error in estimating the coefficients β i . If the first source of error is much larger than the second, the RMSFE is approximately var(u t ), the standard deviation of the error u t , which is estimated by the standard error of regression (SER). One useful application used in time-series forecasting is to test whether the lags of one regressor have useful predictive content. The claim that a variable has no predictive content corresponds to the null hypothesis that the coefficients on all lags of that variable are zero. Such a hypothesis can be checked by the so-called Granger causality test (GCT), a type of F-statistic approach used to test joint hypothesis about regression coefficients. In particular, the GCT method tests the hypothesis that the coefficients of all the values of the variable in . . , Y t−p , are zero, and hence this null hypothesis implies that such regressors have no predictive content for Y t .

Lag length selection
Let us recall relevant statistical methods used to optimally choose the number of lags in an autoregression model; in particular, we focus our attention on the Bayes method (BIC) and on the Akaike method (AIC); for more details, see, for example, [7,Chap. 14.5]. The BIC method is specified by where SSR(p) is the sum of squared residuals of the estimated AR(p). The BIC estimator of p is the value that minimizes BIC(p) among all the possible choices. In the first term of Eq. (5), the sum of squared residuals necessarily decreases when adding a lag. In contrast, the second term is the number of estimated regression coefficients times the factor (ln T )/T , so this term increases when adding a lag. This implies that the BIC trades off these two aspects. The AIC approach is defined by and hence the main difference between the AIC and BIC is that the term ln(T ) in the BIC is replaced by 2 in the AIC , so the second term in the AIC is smaller. But the second term in the AIC is not large enough to assure choosing the correct length, so this estimator of p is not consistent. We recall that an estimator is consistent if, as the size of the sample increases, its probability distribution concentrates at the value of the parameter to be estimated. So, the BIC estimatorp of the lag length in an autoregression is correct in large samples, that is, Pr(p = p) → 1. This is not true for the AlC estimator, which can overestimate p even in large samples; for the proof, see, for example, [7, Appendix 14.5].

Trends
A further relevant topic in econometric analysis is constituted by nonstationarities that are due to trends and breaks. A trend is a persistent long-term movement of a variable over time. A time-series variable fluctuates around its trend. There are two types of trends, deterministic and stochastic. A deterministic trend is a nonrandom function of time. In contrast, a stochastic trend is characterized by a random behavior over time. Our treatment of trends in economic time series focuses on stochastic trend. One of the simplest models of time series with stochastic trend is the one-dimensional random walk defined by the where u t is the error term represented by a normally distributed random variable with zero mean and variance σ 2 > 0. In this case, the best forecast of tomorrow's value is its value today. A extension of the latter is the random walk with drift defined by where the best forecast is the value of the series today plus the drift β 0 . A random walk is nonstationary because the variance of a random walk increases over time, so the distribution of Y t changes over time. In fact, The random walk is a particular case of an AR(1) model with β 1 = 1. If |β 1 | < 1 and u t is stationary, then Y t is stationary. The condition for the stationarity of an AR(p) model is that the roots of 1 − β 1 z − β 2 z 2 − β 3 z 3 − · · · − β p z p = 0 are greater than one in absolute value.
If an AR(p) has a root equal to one, then we say that the series has a unit root and a stochastic trend. Stochastic trends usually bring many issues, for example, the autoregressive coefficients are biased toward zero. Because Y t is nonstationary, the assumptions for time-series regression do not hold, and we cannot rely on estimators and test statistics having their usual large-sample normal distributions; see, for example, [7,Chap. 3.2]. In fact, the OLS estimator of the autoregressive coefficientβ 1 is consistent, but it has a nonnormal distribution; then the asymptotic distribution ofβ 1 is shifted toward zero. Another problem caused by stochastic trend is the nonnormal distribution of the t-statistic, which means that conventional confidence intervals are not valid and hypothesis tests cannot be conducted as usual. The t-statistic is an important example of a test statistic, namely of a statistic used to perform a hypothesis test. A statistical hypothesis test can make two types of mistakes: a type I error, in which the null hypothesis is rejected when, in fact, it is true, and a type II error, in which the null hypothesis is not rejected when, in fact, it is false. The prespecified rejection probability of a statistical hypothesis test when the null hypothesis is true, that is, the prespecified probability of a type I error, is called the significance level of the test. The critical value of the test statistic is the value of the statistic for which the test just rejects the null hypothesis at the given significance level. The p-value is the probability of obtaining a test statistic, by random sampling variation, at least as adverse to the null hypothesis value as is the statistic actually observed, assuming that the null hypothesis is correct. Equivalently, the p-value is the smallest significance level at which you can reject the null hypothesis. The value of the t-statistic is t = estimator − hypothesized value standard error of the estimator and is well approximated by the standard normal distribution when n is large because of the central limit theorem (see, e.g., [ Then we assume that the following hypothesis test holds: For an AR(p) model, it is standard to use the augmented Dickey-Fuller test (ADF), which tests the null hypothesis H 0 : δ = 0 against the one-side alternative H 1 : δ < 0 in the regression under the null hypothesis. Let us note that since Y t has a stochastic trend, it follows that, under the alternative hypothesis, Y t is stationary. The ADF statistic is the OLS t-statistic testing δ = 0. If, instead, the alternative hypothesis is that Y t is stationary around a deterministic linear time trend, then this trend t must be added as an additional regressor. In this case, the Dickey-Fuller regression becomes and we test for δ = 0. The ADF statistic does not have a normal distribution, and hence different critical values have to be used.

Breaks
A second type of nonstationarity arises when the regression function changes over the course of the sample. In economics, this can occur for a variety of reasons, such as changes in economic policy, changes in the structure of the economy, or an invention that changes a specific industry. These breaks cannot be neglected by the regression model. A problem caused by breaks is that the OLS regression estimates over the full sample will estimate a relationship that holds "on average," in the sense that the estimate combines two different periods, and this leads to poor forecast. There are two types of testing for breaks: testing for a break at a known date and for a break at an unknown break date. We consider the first option for an AR(p) model. Let τ denote the hypothesized break date, and let D t (τ ) be the binary variable such that D t (τ ) = 0 if t > τ and D t (τ ) = 1 if t < τ . Then the regression including the binary break indicator and all interaction terms reads as follows: under the null hypothesis of no breaks, γ 0 = γ 1 = γ 2 = · · · = γ p = 0. Under the alternative hypothesis that there is a break, the regression function is different before and after the break date τ , and we can use the F-statistic performing the socalled the Chow test (see, e.g., [6, Chap. 5.3.3]). If we suspect a break between two dates τ 0 and τ 1 , the Chow test can be modified to test for breaks at all possible dates τ between τ 0 and τ 1 , then using the largest of the resulting F-statistics to test for a break at an unknown date. The latter technique is called the Quandt likelihood ratio statistic (QLR) (see, e.g., [7, Chap. 14.7]). Because the QLR statistic is the largest of many F-statistics, its distribution is not the same as that of an individual F-statistic; also, the critical values for the QLR statistic must be obtained from a special distribution.

MA and ARMA
In the following, we consider finite-order moving-average (MA) processes (see, e.g., [6, Chap. 2.2]). The moving-average process of order q, MA(q), is defined by Y t = α 0 + u t − α 1 u t−1 − α 2 u t−2 − · · · − α q u t−q ; equivalently, by using the lag operator we get Y t − α 0 = (1 − α 1 L − α 2 L 2 − · · · − α q L q )u t . Every finite MA(q) process is stationary, and we have Combining both an autoregressive (AR) term of order p and a moving-average (MA) term of order q, we can define the process denoted as ARMA(p, q) and represented by again, exploiting the lag operator, we can write

Vector autoregression
In what follows, we focus our study on the so-called vector autoregression (VAR) econometric model, also using some remarks on the relation between the univariate time series models described in the first part, and the set of simultaneous equations systems of traditional econometrics characterizing the VAR approach (see, e.g., [4,Chap. 2]).

Representation of the system
We have so far considered forecasting a single variable. However, it is often necessary to allow for a multidimensional statistical analysis if we want to forecast more than one-parameter dynamics. This section introduces a model for forecasting multiple variables, namely the vector autoregression (VAR) model, in which lagged values of two or more variables are used to forecast their future values. We start with the autoregressive representation in a VAR model of order p, denoted by VAR(p), where each component depends on its own lagged values up to p periods and on the lagged values of all other variables up to order p. It follows that the main idea behind the VAR model is to know how new information, appearing at a certain time point and concerning one of the observed variables, is processed in the system and which impact it has over time not only for this particular variable but also for the other system parameters. Hence, a VAR(p) model is a set of k time-series regressions (k ∈ N + ) in which the regressors are lagged values of all k series and the number of lags equals p for each equation. In the case of two time series variables, say, Y t and X t , the VAR(p) consists of two equations of the form where the βs and the γs are unknown coefficients, and u 1t and u 2t are error terms represented by normally distributed random variables with zero mean and variance σ 2 i > 0. The VAR assumptions are the same as those for the time-series regression defining AR models and applied to each equation; moreover, the coefficients of each VAR are estimated by means of the OLS approach. The reduced form of a vector autoregression of orderp is defined as where A i , i = 1, . . . , p, are k-dimensional quadratic matrices, U represents the kdimensional vector of residuals at time t, and δ is the vector of constant terms. System (6) can be rewritten compactly as A p (L)Z t = δ + U t , where A p (L) = t ] = σ uu , and E[U t U ′ s ] = 0 for t = s. Such a system is stable if and only if all included variables are stationary, that is, if all roots of the characteristic equation of the lag polynomial are outside the unit circle, namely det(I k − A 1 z − A 2 z − · · · − A p z) = 0 for |z| ≤ 1 (for details, see, e.g., [6,Chap. 4.1]). We use this condition because we saw in Section 2.3 that the condition for the stationarity of an AR(p) model is that the roots of 1 − β 1 z − β 2 z 2 − β 3 z 3 − · · · − β p z p = 0 are greater than one in absolute value. If an AR(p) has a root equal to one, we say that the series has a unit root and a stochastic trend. Moreover, the previous system can be rewritten by exploiting the MA representation as follows: The autocovariance matrices are defined as ; without loss of generality, we set δ = 0 and, therefore, µ = 0, whence we obtain and, for τ ≥ 0, Since the autocovariance matrix entries link a variable with both its delays and the remaining model variables, we have that if the autocovariance between X and Y is positive, then X tends to move accordingly with Y and vice versa, whereas if X and Y are independent, their autocovariance obviously equals zero.

Determining lag lengths in VARs
An appropriate method for the lag length selection of VAR is fundamental to determine properties of VAR and related estimates. There are two main approaches used for selecting or testing lag length in VAR models. The first consists of rules of thumb based on the periodicity of the data and past experience, and the second is based on formal information criteria. VAR models typically include enough lags to capture the full cycle of the data; for monthly data, this means that there is a minimum of 12 lags, but we will also expect that there is some seasonality that is carried over from year to year, so often lag lengths of 13-15 months are used (see, e.g., [4,Chap. 2.5]). For quarterly data, it is standard to use six lags. This captures the cyclical components in the year and any residual seasonal components in most cases. Usually, we decide to choose the number of lags not exceeding kp + 1 < T , where k is the number of endogenous variables, p is the lag length, and T is the total number of observations. We use this limitation because the estimate of all these coefficients increases the amount of forecast estimation errors, which can result in a deterioration of the accuracy of the forecast itself. The lag length in VAR can be formally determined using information criteria; letΣ uu be the estimate of the covariance matrix with the (i, j) element 1 T T t=1û itûjt , whereû it is the OLS residual from the jth equation. The BIC for the kth equation in a VAR model is whereas the AIC is computed using Eq. (7), modified by replacing the term ln T by 2.
Among a set of candidate values of p, the estimated lag lengthp is the value of p that minimizes BIC(p).

Multiperiod VAR forecast
Iterated multivariate forecasts are computed using a VAR in much the same way as univariate forecasts are computed using an autoregression. The main new feature of a multivariate forecast is that the forecast of one variable depends on the forecast of all variables in the VAR. To compute multiperiod VAR forecasts h periods ahead, it is necessary to compute forecast of all variables for all intervening periods between T and T +h. Then the following scheme applies: compute the one-period-ahead forecast of all the variables in the VAR, then use those forecasts to compute the two-periodahead forecasts, and repeat the previous stops until the desired forecast horizon. For example, the two-period-ahead forecast of Y T +2 based on the two-variable VAR(p) in Eq. (6) iŝ where the coefficients in (8) are the OLS estimates of the VAR coefficients.

Granger causality
An important question in multiple time series is to assign the value of individual variables to explain the remaining ones in the considered system of equations. An example is the value of a variable Y t for predicting another variable X t in a dynamic system of equations or understanding if the variable Y t is informative about future values of X t . The answer is based on the determination of the so-called Granger causality parameter for a time-series model (for details, see, e.g., [4,Chap. 2.5.4]).
To define the concept precisely, consider the bivariate VAR model for two variables (Y t , X t ) as in Eq. (6). Using this system of equations, Granger causality states that, for linear models, X t Granger causes Y t if the behavior of past Y t can better predict the behavior of X t than the past X t alone. For the model in system (6), if X t Granger causes Y t , then the coefficients for the past values of X t in the Y t equation are nonzero, that is, γ 1i = 0 for i = 1, 2, . . . , p. Similarly, if Y t Granger causes X t in the X t equation, then the coefficients for the past values of Y t are nonzero, that is, β 2i = 0 for i = 1, 2, . . . , p. The formal testing for Granger causality is then done by using an F test for the joint hypothesis that the possible causal variable does not cause the other variable. We can specify the null hypothesis for the Granger causality test as follows.
In the first model, we have γ 11 = 0, γ 12 = 0, . . . , γ 1p = 0, so the variable X t compares in the equation of Y t , namely the values of X t are useful to predict Y t . Instead, in the second model, γ 11 = γ 12 = · · · = γ 1p = 0, so X t does not Granger cause Y t . The test statistic has an F distribution with(p, T − 2p − 1) degrees of freedom: If this F statistic is greater than the critical value for a chosen level of significance, we reject the null hypothesis that X t has no effect on Y t and conclude that X t Granger causes Y t .

Cointegration
In Section 2.3, we introduced the model of random walk with drift as follows: If Y t follows Eq. (9), then it has an autoregressive root that equals 1. If we consider a random walk for the first difference of the trend, then we obtain Hence, if Y t follows Eq.(10), then ∆Y t follows a random walk, and accordingly ∆Y t − ∆Y t−1 is stationary; this is the second difference of Y t and is denoted ∆ 2 Y t . A series that has a random walk trend is said to be integrated of order one, or I(1); a series that has a trend of the form (10) is said to be integrated of order two, or I(2); and a series that has no stochastic trend and is stationary is said to be integrated of order zero, or I(0). The order of integration in the I(1) and I(2) terminology is the number of times that the series needs to be differenced for it to be stationary. If Y t is I(2), then ∆Y t is I(1), so ∆Y t has an autoregressive root that equals 1. If, however, Y t is I(1), then ∆Y t is stationary. Thus, the null hypothesis that Y t is I(2) can be tested against the alternative hypothesis that Y t is I(1) by testing whether ∆Y t has a unit autoregressive root. Sometimes, two or more series have the same stochastic trend in common. In this special case, referred to as cointegration, regression analysis can reveal long-run relationships among time series variables. One could think that a linear combination of two processes I(1) is a process I(1). However, this is not always true. Two or more series that have a common stochastic trend are said to be cointegrated. Suppose that X t and Y t are integrated of order one. If, for some coefficient θ, Y t − θX t is integrated of order zero, then X t and Y t are said to be cointegrated, and the coefficient θ is called the cointegrating coefficient. If X t and Y t are cointegrated, then they have a common stochastic trend that can be eliminated by computing the difference Y t − θX t , which eliminates this common stochastic trend. There are three ways to decide whether two variables can be plausibly modeled exploiting the cointegration approach, namely, by expert knowledge and economic theory, by a qualitative (graphical) analysis of the series checking for common stochastic trend, and by performing statistical tests for cointegration. In particular, there is a cointegration test when θ is unknown. Initially, the cointegrating coefficient θ is estimated by OLS estimation of the regression Y t = α + θX t + z t , and then we use the Dickey-Fuller test (see Section 2.3) to test for a unit root in z t ; this procedure is called the Engle-Granger augmented Dickey-Fuller test for cointegration (EG-ADF test); for details, see, for example, [6, Chap. 6.2] . The concepts covered so far can be extended to the case of more than two variables, for example, three variables, each of which is I(1), are said to be cointegrated if Y t −θ 1 X 1t −θ 2 X 2t is stationary. The Dickey-Fuller needs the use of different critical values (see Table 1), where the appropriate line depends on the number of regressors used in the first step of estimating the OLS cointegrating regression. A different estimator of the cointegrating coefficient is the dynamic OLS (DOLS) estimator, which is based on the equation In particular, from Eq. (12) we notice that DOLS includes past, present, and future values of the changes in X t . The DOLS estimator of θ is the OLS estimator of θ in Eq. (12). The DOLS estimator is efficient, and statistical inferences about θ and δs in Eq. (12) are valid. If we have cointegration in more than two variables, for example, three variable Y t , X 1t , X 2t , each of which is I(1), then they are cointegrated with cointegrating coefficients θ 1 and θ 2 if Y t − θ 1 X 1t − θ 2 X 2t is stationary. The EG-ADF procedure to test for a single cointegrating relationship among multiple variables is the same as for the case of two variables, except that the regression in Eq. (11) is modified so that both X 1t and X 2t are regressors. The DOLS estimator of a single cointegrating relationship among multiple Xs involves the level of each X along with lags of the first difference of each X.

Conclusion
In this first part of our ambitious project to use multivariate statistical techniques to study critic econometric data of one of the most influential economy in Italy, namely the Verona import-export time series, we have focused ourselves on a self-contained introduction to techniques of estimating OLS-type regressions, analysis of the correlations obtained between the different variables and various types of information criteria to check for the goodness of fit. A particular relevance has been devoted to the application of tests able to enlightening various types of nonstationarity for the considered time series, for example, the augmented Dickey-Fuller test (ADF) and the Quandt likelihood ratio statistic (QLR). Moreover, we have also exploited both the Granger causality test and the Engle-Granger augmented Dickey-Fuller test for cointegration (EG-ADF) in order to analyze if and how these variables are related to each other and to have a measure on how much a variable gives information on the other one. Such approaches constitute the core of the second part of our project, namely the aforementioned Verona case study.