Time Series: Cointegration Analysis

Unit root (ADF and KPSS) test

The ADF, available in the function adf.test() (in the package tseries) implements the t-test of \(H_0: \gamma = 0\) in the regression, below.

\[\begin{equation} \label{eqn:lagged-ts-regression} \Delta {{Y}_{t}}={{\beta }_{1}}+{{\beta }_{2}}t+\gamma {{Y}_{t-1}}+ \sum\limits_{i=1}^{m}{\delta_i \Delta {{Y}_{t-i}}+{{\varepsilon }_{t}}} \end{equation}\]

The null is therefore that series has a unit root. If only the series has a non-unit root, then it is stationary (rejection of null hypothesis).

The ADF test was parametrized with the alternative hypothesis of stationarity. This extends to following assumption in the model parameters;

\[ -2 \leq \gamma \leq 0\ \text{or } (-1 < 1+\phi < 1) \]

k in the function refers to the number of \(\delta\) lags, i.e., \(1, 2, 3, ...., m\) in the model equation.

The number of lags k defaults to trunc((length(x)-1)^(1/3)), where x is the series being tested. The default value of k corresponds to the suggested upper bound on the rate at which the number of lags, k, should be made to grow with the sample size for the general ARMA(p,q) setup citation(package = "tseries").

For a Dickey-Fueller test, only up to AR(1) time dependency in our stationary process is considered, we set k = 0. Hence we have no \(\delta\)s (lags) in our test.

The DF model can be written as:

\[ Y_t = \beta_1 + \beta_2 t + \phi Y_{t-1} + \varepsilon_t \]

It can be re-written so we can do a linear regression of \(\Delta Y_t\) against \(t\) and \(Y_{t-1}\) and test if \(\phi\) is different from 0. If only, \(\phi\) is not zero and assumption above (\(-1 < 1+\phi < 1\)) holds, the process is stationary. If \(\phi\) is straight up 0, then we have a random walk process – all white noise.

\[ \Delta {Y}_{t}=\beta_1+\beta_2 t+\gamma {Y}_{t-1} + \varepsilon_{t} \]

Alternative to above discussed tests, the Phillips-Perron test with its nonparametric correction for autocorrelation (essentially employing a HAC estimate of the long-run variance in a Dickey-Fuller-type test instead of parametric decorrelation) may be used. It is available in the function pp.test().

Unit root test based lag order differencing determination

An alternative to decomposition for removing trends is differencing [@woodward2017applied]. We define the order \(d\) difference operator as,

\[\begin{equation} \nabla^d x_t = (1-\mathbf{B})^d x_t, \label{eqn:order-d-difference-operator} \end{equation}\]

Where \(\mathbf{B}\) is the backshift operator (i.e., \(\mathbf{B}^k x_t = x_{t-k}\) for \(k \geq 1\)).

Applying the difference to a random walk, the most simple and widely used time series model, will yield a time series of Gaussian white noise errors \(\{w_t\}\):

\[\begin{equation} \begin{aligned} \nabla (x_t &= x_{t-1} + w_t) \\ x_t - x_{t-1} &= x_{t-1} - x_{t-1} + w_t \\ x_t - x_{t-1} &= w_t \end{aligned} \label{eqn:random-walk-series} \end{equation}\]

Vector autoregressive model (VAR)

VAR is a system regression model, i.e., there are more than one dependent variable. The regression is defined by a set of linear dynamic equations where each variable is specified as a function of an equal number of lags of itself and all other variables in the system. Any additional variable, adds to the modeling complexity by increasing an extra equation to be estimated.

The vector autoregression (VAR) model extends the idea of univariate autoregression to \(k\) time series regressions, where the lagged values of all \(k\) series appear as regressors. Put differently, in a VAR model we regress a vector of time series variables on lagged vectors of these variables. As for \(AR(p)\) models, the lag order is denoted by \(p\) so the \(VAR(p)\) model of two variables \(X_t\) and \(Y_t\) (\(k=2\)) is given by a vector of equations (Equation ??).

\[\begin{equation} \label{eqn:vector-regression-ts} \begin{split} \begin{aligned} Y_t =& \, \beta_{10} + \beta_{11} Y_{t-1} + \dots + \beta_{1p} Y_{t-p} + \gamma_{11} X_{t-1} + \dots + \gamma_{1p} X_{t-p} + u_{1t}, \\ X_t =& \, \beta_{20} + \beta_{21} Y_{t-1} + \dots + \beta_{2p} Y_{t-p} + \gamma_{21} X_{t-1} + \dots + \gamma_{2p} X_{t-p} + u_{2t}. \end{aligned} \end{split} \end{equation}\]

The \(\beta\)s and \(\gamma\)s can be estimated using OLS on each equation.

Simplifying this to a bivariate \(VAR(1)\), we can write the model in matrix form as:

\[\begin{equation} \label{eqn:matix-var1-model} Y_t = \beta_0 + \beta_1 Y_{t-1} + \mu_t \end{equation}\]

Where,

• \(Y_t, Y_{t-1}\) and \(\mu_t\) are (2 x 1) column vectors
• \(\beta_0\) is a (2 x 1) column vector
• \(\beta_1\) is a (2 x 2) matrix

also,

\[ Y_t = \begin{pmatrix} y_{1t} \\ y_{2t} \end{pmatrix},\ Y_{t-1} = \begin{pmatrix} y_{1t-1} \\ y_{2t-1} \end{pmatrix} \]

\[ \mu_t = \begin{pmatrix} \mu_{1t} \\ \mu_{2t} \end{pmatrix}, \beta_{0} = \begin{pmatrix} \beta_{10} \\ \beta_{20} \end{pmatrix}, \beta_{1} = \begin{pmatrix} \beta_{11} & \alpha_{11} \\ \alpha_{21} & \beta_{21} \end{pmatrix} \]

It is straightforward to estimate VAR models in R. A feasible approach is to simply use lm() for estimation of the individual equations. Furthermore, the vars package provides standard tools for estimation, diagnostic testing and prediction using this type of models.

Only when the assumptions presented below hold, the OLS estimators of the VAR coefficients are consistent and jointly normal in large samples so that the usual inferential methods such as confidence intervals and \(t\)-statistics can be used [@metcalfe2009introductory].

Two series \(w_{x,t}\) and \(w_{y,t}\) are bivariate white noise if they are stationary and their cross-covariances \(\gamma_{xy}(k) = Cov(w_{x,t}, w_{y, t+k})\) satisfies

\[ \gamma_{xx}(k) = \gamma_{yy}(k) = \gamma_{xy}(k) = 0\ \text{for all } k \neq 0 \]

The parameters of a var(p) model can be estimated using the ar function in R, which selects a best-fitting order \(p\) based on the smallest information criterion values.

##       [,1]  [,2]
## [1,] 1.009 0.782
## [2,] 0.782 0.947

##      [,1] [,2]
## [1,]  0.4  0.3
## [2,]  0.2  0.1

## [1]  1.86 26.86

##       x      y
## x 0.399 0.2813
## y 0.218 0.0419

The structure of VARs also allows to jointly test restrictions across multiple equations. For instance, it may be of interest to test whether the coefficients on all regressors of the lag \(p\) are zero. This corresponds to testing the null that the lag order \(p-1\) is correct. Large sample joint normality of the coefficient estimates is convenient because it implies that we may simply use an \(F\)-test for this testing problem. The explicit formula for such a test statistic is rather complicated but fortunately such computations are easily done using the ttcode("R") functions we work with in this chapter. Just as in the case of a single equation, for a multiple equation model we choose the specification which has the smallest \(BIC(p)\), where

\[ \begin{aligned} BIC(p) =& \, \log\left[\text{det}(\widehat{\Sigma}_u)\right] + k(kp+1) \frac{\log(T)}{T}. \end{aligned} \]

with \(\widehat{\Sigma}_u\) denoting the estimate of the \(k \times k\) covariance matrix of the VAR errors and \(\text{det}(\cdot)\) denotes the determinant.

As for univariate distributed lag models, one should think carefully about variables to include in a VAR, as adding unrelated variables reduces the forecast accuracy by increasing the estimation error. This is particularly important because the number of parameters to be estimated grows qudratically to the number of variables modeled by the VAR.

Residual based cointegration test

Since the food commodities are spatially linked, more of so because they occupy the same domestic market, it is obvious that factor affecting price of one inevitably affects other, especially that of same crop in a nearby market. Having evidence for nonstationarity, it is of interest to test for a common nonstationary component by means of a cointegration test (Non-stationarity is more valid for development regionwise price series).

A two step method proposed by @hylleberg1990seasonal can be used to test for cointegration. The procedure simply regressess one series on the other and performs a unit root test on the residuals. This test is often named after @phillips1990asymptotic. Specifically, po.test() performs a Phillips-Perron test using an auxiliary regression without a constant and linear trend and the Newey-West estimator for the required long-run variance.

The test computes the Phillips-Ouliaris test for the null hypothesis that series is not cointegrated [@R-tseries].

Note po.test does not handle missing values, so we fix them through imputation. It is implemented through tidyr::fill(..., .direction = "down").

The problem with this approach is that it treats both series in an asymmetric fashion, while the concept of cointegration demands that the treatment be symmetric.

The po.test() function is testing the cointegration with Phillip’s Z_alpha test, which is the second residual-based test described by @phillips1990asymptotic. Because the po.test() will use the series at the first position to derive the residual used in the test, results would be determined by the series on the most left-hand side¹.

The Phillips-Ouliaris test implemented in the ca.po() function from the urca package is different. In the ca.po() function, there are two cointegration tests implemented, namely “Pu” and “Pz” tests. Although both the ca.po() function and the po.test() function are supposed to do the Phillips-Ouliaris test，outcomes from both functions are completely different.

Similar to Phillip’s Z_alpha test, the Pu test also is not invariant to the position of each series and therefore would give different outcomes based upon the series on the most left-hand side. On the contrary, the multivariate trace statistic of Pz test has its appeal in that the outcome won’t change by the position of each series.

VAR based cointegration test

The standard tests proceeding in a symmetric manner stem from Johansen’s full-information maximum likelihood approach [@johansen1991estimation]. A general vector autoregressive model is similar to the AR(p) model except that each quantity is vector valued and matrices are used as the coefficients. The general form of the VAR(p) model, without drift, is given by:

\[\begin{equation} \label{eqn:var-general} {\bf y_t} = {\bf \mu} + A_1 {\bf y_{t-1}} + \ldots + A_j {\bf y_{t-j}} + {\bf \varepsilon_t} \end{equation}\]

Where \({\bf \mu}\) is the vector-valued mean of the series, \(A_i\) are the coefficient matrices for each lag and \({\bf \varepsilon_t}\) is a multivariate Gaussian noise term with mean zero.

At this stage we can form a Vector Error Correction Model (VECM) by differencing the series (Equation ??).

\[\begin{equation} \label{eqn:vecm-differenced} \Delta {\bf y_t} = {\bf \mu} + A {\bf y_{t-1}} + \Gamma_1 \Delta {\bf y_{t-1}} + \ldots + \Gamma_j \Delta {\bf y_{t-j}} + {\bf \varepsilon_t} \end{equation}\]

Where \(\Delta {\bf y_t} = {\bf y_t} - {\bf y_{t-1}}\) is the differencing operator, \(A\) is the coefficient matrix for the first lag and \(\Gamma_i\) are the matrices for each differenced lag.

For a \(p^{th}\) -order cointegrated vector autoregressive (VAR) model, the error correction form is (omitting deterministic components; both no intercept or trend in either cointegrating equation or test var), we may rewrite the VAR in the form of Equation ?? [@johansen1991estimation].

For more information refer to: - Eviews documentation - Quantstart post - Kevin kotze’s post

\[\begin{equation} \label{eqn:johansens} \Delta y_t = \Pi y_{t-1} + \sum_{j = 1}^{p-1} {\Gamma_j \Delta y_{t-j}} + \varepsilon_t \end{equation}\]

Where,

\[ \Pi = \sum^{p}_{i = 1}{A_{i}-I}; \Gamma = -\sum^{p}_{j = i + 1}{j} \]

(Although, for simplicity sake, we assume absence of deterministic trends, there are five popular scenarios of including such trends in a cointegration test. All of these are described in [@johansen1995identifying].)

Granger’s representation theorem asserts that if the coefficient matrix \(\Pi\) has reduced rank \(r < k\), then there exist \(kxr\) matrices \(\alpha\) and \(\beta\) each with rank \(k\) such that \(\Pi = \alpha \beta^{\prime}\) and \(\beta^{\prime}y_t\) is \(I(0)\).

To achieve this an eigenvalue decomposition of \(A\) is carried out. The rank of the matrix \(A\) is given by \(r\) and the Johansen test sequentially tests whether this rank \(r\) is equal to zero, equal to one, through to \(r=n-1\), where \(n\) is the number of time series under test.

The null hypothesis of \(r=0\) means that there is no cointegration at all. A rank \(r > 0\) implies a cointegrating relationship between two or possibly more time series.

The eigenvalue decomposition results in a set of eigenvectors. The components of the largest eigenvector admits the important property of forming the coefficients of a linear combination of time series to produce a stationary portfolio. Notice how this differs from the CADF test (often known as the Engle-Granger procedure) where it is necessary to ascertain the linear combination a priori via linear regression and ordinary least squares (OLS).

In summary, the test checks for the situation of no cointegration, which occurs when the matrix \(A=0\). So, starting with the base value of \(r\) (i.e., \(r=0\)), if the test statistic is greater than critical values of at the 10%, 5% and 1% levels, this would imply that we are able to reject the null of no cointegration. For the case r<=1, we if the calculated test statistic is below the critical values of, we are unable to reject the null, and the number of cointegrating vectors is between 0 and 1. The relevant tests are available in the function urca::ca.jo(). The basic version considers the eigenvalues of the matrix \(\Pi\) in the preceding equation.

Causality test

Causality test is VAR based approach to explain cause-effect relationship among endogenous variables. However, the Granger-causality [@granger1988causality] inference does not, of course, establish the real causation phenomena. If one of the variables is sufficiently correlated to the other so that forecast of former depends on the later to a considerable extent, then the first variable is granger-caused by the second one.

The Granger causality has been briefed to be useable test in certain cases of two series violating the stationarity assumption².

@papana2014identifying state that GC test can only to applied if both the series are stationary³. Same paper also cautioned that VAR(1) models of cointegrated endogenous series will fail to capture long-run relationships. Therefore, the authors suggest surplus lag Granger-causality test be used if the series are a nonstationary data.

Common statistical tests and econometric models are interpreted based on cut-off probability value 5%. Where model selection is desired, model terms minimizing Akaike or Bayesian Information Criterion, or where relevant both, were selected.