Several quantities may be recorded at each time: prices of related goods, temperatures at nearby stations, or measurements from several sensors. We model such observations by vectors \(X_t\in\R^d\). Their mean is a vector; their covariance records both variation within coordinates and variation between them. For random vectors \(U,V\) with finite second moments, define \[
\Cov(U,V):=\E[(U-\E[U])(V-\E[V])^\top].
\] Its \((i,j)\) entry is \(\Cov(U_i,V_j)\).
The definition of weak stationarity for vectors is directly analogous to that for scalars.
Definition 10.1 (Vector weak stationarity). A process \((X_t)\) in \(\R^d\) is weakly stationary if it has finite second moments, a constant mean vector \(\mu\), and \[
\Cov(X_{t+h},X_t)=\Gamma(h)
\] for every admissible \(t,h\), where the matrix \(\Gamma(h)\) depends only on the lag. The function \(\Gamma\) is its matrix autocovariance function.
The entry \(\Gamma_{ij}(h)\) compares coordinate \(i\) at time \(t+h\) with coordinate \(j\) at time \(t\). An important difference between the scalar and vector definitions:Transposing the covariance gives \(\Gamma(-h)=\Gamma(h)^\top\). Thus \(\Gamma(0)\) is symmetric and positive semidefinite; \(\Gamma(h)\) need not be symmetric at a nonzero lag (note that this means that it generalizes the scalar definition regardless). Stationarity of each coordinate checks only the diagonal entries. Joint stationarity also requires the cross-covariances to depend only on lag.
Filters Across Time and Coordinates
The backshift operator \(B\) moves an observation one step backward in time for vectors all the same. For a vector series, it shifts every coordinate together: \[
BX_t:=X_{t-1}=(X_{1,t-1},\ldots,X_{d,t-1})^\top,
\qquad B^jX_t=X_{t-j}.
\] A scalar filter combines observations at different times. With vectors, matrix coefficients can also combine coordinates.
Definition 10.2 (Finite matrix filter). For fixed real \(r\times d\) matrices \(C_0,\ldots,C_m\), the filtered process is \[
W_t:=\sum_{j=0}^{m}C_jB^jX_t
=\sum_{j=0}^{m}C_jX_{t-j}\in\R^r.
\] The output dimension \(r\) can differ from the input dimension \(d\).
The difference filter uses \(C_0=I_d\), \(C_1=-I_d\) and produces \(\nabla X_t=(1-B)X_t=X_t-X_{t-1}\), where \(I_d\) is the identity matrix. It differences each coordinate separately.
Proposition 10.3 (Weak stationarity under finite matrix filtering). Let \((X_t)\) be a weakly stationary process in \(\R^d\) with mean \(\mu\) and matrix ACVF \(\Gamma\). For any fixed real \(r\times d\) matrices \(C_0,\ldots,C_m\), the process \((W_t)\) in Definition 10.2 is weakly stationary, with mean \[
\mu_W=\sum_{j=0}^{m}C_j\mu
\] and matrix ACVF \[
\Gamma_W(h)=\sum_{j=0}^{m}\sum_{k=0}^{m}
C_j\Gamma(h-j+k)C_k^\top.
\] These formulas hold wherever all required observations are defined.
NoteProof
Each coordinate of \(W_t\) is a finite linear combination of random variables with finite second moments, so it also has a finite second moment. By linearity, \[
\E[W_t]=\sum_{j=0}^{m}C_j\E[X_{t-j}]
=\sum_{j=0}^{m}C_j\mu=\mu_W,
\] which does not depend on \(t\).
For the covariance, the lag between \(t+h-j\) and \(t-k\) is \(h-j+k\). Thus \[
\begin{aligned}
\Cov(W_{t+h},W_t)
&=\sum_{j=0}^{m}\sum_{k=0}^{m}
C_j\Cov(X_{t+h-j},X_{t-k})C_k^\top\\
&=\sum_{j=0}^{m}\sum_{k=0}^{m}
C_j\Gamma(h-j+k)C_k^\top.
\end{aligned}
\] This depends only on \(h\), so \((W_t)\) is weakly stationary. \(\square\)
In a similar way as for scalars, we can define \(I(1)\) and \(I(0)\) for vectors.
Definition 10.4 (Vector integration orders). A process \((X_t)\) in \(\R^d\) is \(I(0)\) if it is jointly weakly stationary, its matrix ACVF satisfies \(\sum_{h\in\Z}|\Gamma_{ij}(h)|<\infty\) for every \(i,j\), and its long-run covariance matrix \[
\Omega:=\sum_{h\in\Z}\Gamma(h)
\] has positive diagonal entries. We allow \(\Omega\) to be singular. It is \(I(1)\) if \[
X_t=X_0+\sum_{j=1}^{t}V_j,\qquad t\ge1,
\] where \((V_t)\) is a mean-zero \(I(0)\) vector process and \(X_0\) has finite second moments.
Just like with scalars, first differencing makes \(I(1)\) vector processes into \(I(0)\) vector processes. However, in the multivariate setting, there is a further possibility. Even when every coordinate of \(X_t\) is \(I(1)\), a fixed linear combination \(\beta^\top X_t\) can be \(I(0)\). This filter combines coordinates at the same time: \(m=0\) and \(C_0=\beta^\top\). Before we proceed with differencing and developing the theory of vector AR and MA models, we will explore this additional possibility to obtain stationary processes.
Definition 10.5 (Cointegration). Suppose every coordinate of \(X_t\in\R^d\) is \(I(1)\). The coordinates are cointegrated if there is a fixed nonzero vector \(\beta\in\R^d\) such that \(\beta^\top X_t\) is \(I(0)\). We call \(\beta\) a cointegrating vector.
The scalar combination is called a spread. Multiplying \(\beta\) by a nonzero constant describes the same relation. Subtracting a constant from the spread only changes its mean.
Example 10.6 (Cancelling a stochastic trend). Let \((Z_t)\), \((U_t)\), and \((V_t)\) be mutually independent \(I(0)\) processes, indexed by \(t\in\Z\), with \(\E[Z_t]=0\). Write their positive long-run variances as \(\tau_Z^2\), \(\tau_U^2\), and \(\tau_V^2\), respectively. Define \[
S_0:=0,\qquad S_t:=S_{t-1}+Z_t,\quad t\geq1,
\] and consider \[
X_t:=S_t+U_t,\qquad Y_t:=2S_t+V_t.
\] Their differences are \[
\nabla X_t=Z_t+\nabla U_t,\qquad
\nabla Y_t=2Z_t+\nabla V_t.
\] These are mean-zero, weakly stationary processes with absolutely summable autocovariances. Their long-run variances are \(\tau_Z^2\) and \(4\tau_Z^2\), respectively, because \(\nabla U_t\) and \(\nabla V_t\) have zero long-run variance. Both differences are therefore \(I(0)\), so both levels are \(I(1)\). Yet \[
Y_t-2X_t=V_t-2U_t
\] is also \(I(0)\), with long-run variance \(\tau_V^2+4\tau_U^2>0\). The common accumulated component cancels, so the coordinates are cointegrated with cointegrating vector \((-2,1)^\top\).
The gap can be serially correlated, so its history may help predict it… \(\diamond\)
The requirement for \((Y_t,X_t)\) to be \(I(1)\) in each coordinate is crucial for cointegration. There are many processes that posses a common trend that are not cointegrated, and their behavior is fundamentally different.
Example 10.7 (Cancelling a deterministic trend). Now replace \(S_t\) by \(t\): \[
X_t:=t+U_t,\qquad Y_t:=2t+V_t.
\] Neither process is stationary, because its mean changes with time. Their combination is again \[
Y_t-2X_t=V_t-2U_t,
\] which is \(I(0)\). Nevertheless, the coordinates are not cointegrated. Neither is \(I(1)\): their variances remain constant, whereas an \(I(1)\) process has variance growing asymptotically linearly. A stationary combination of nonstationary processes is therefore not sufficient for cointegration. \(\diamond\)
The mathematical difference between the two examples is subtle: both processes result in a stationary. However, in the first case, we can draw a specific conclusion about a statistical relationship between the two processes. In Example 10.6, the relation between the levels reflects a shared stochastic component that is (potentially) unobserved: both processes contain \(S_t\). They are \(I(1)\), and subtracting \(2X_t\) from \(Y_t\) cancels that component, leaving an \(I(0)\) gap.
In Example 10.7, a regression can fit the levels closely because both means increase with time, even though the processes are independent. Interpreting this fit as evidence of some kind of dependence (a statistical relationship) would be wrong! More generally, these so called spurious regressions refers to misleading evidence from regression when series happen to have a common trend. It is not defined by whether the gap is stationary: here the gap is \(I(0)\), but the coordinates are not \(I(1)\) and hence are not cointegrated: there was no common stochastic component to cancel, the processes just happened to travel in “at the same rate” over time.
Given observations, how might we find a candidate cointegrating coefficient? A levels regression fits the undifferenced observations \(y_t\) on an intercept and \(x_t\): \[
y_t\approx a+cx_t.
\] Least squares chooses the coefficients to make the observed squared gaps small. This criterion establishes neither that the original processes are \(I(1)\) nor that \(Y_t-a-cX_t\) is \(I(0)\). Trends and persistence can produce high \(R^2\) values and apparently significant regression coefficients even for independent processes. This misleading regression evidence is called spurious regression.
Spurious Regression
Australian air passengers and rice production in Guinea make an unusual regression pair. Hyndman and Athanasopoulos (2021, Section 7.3) use these 42 annual World Bank observations, 1970–2011, to illustrate the effect of trends on a levels regression, which uses the undifferenced values of the series. Here we regress each year’s passenger numbers on an intercept and that year’s rice production. Regressing their year-to-year changes on each other would instead be a regression in differences. Passengers are measured in millions and rice production in million metric tons.
Load the 42 observations and fit the line
air_rice <-read.csv(text ='"Year","Passengers","Production"1970,7.3187,0.3114051971,7.3266,0.325181972,7.7956,0.3395651973,9.3846,0.354581974,10.6647,0.3702651975,11.0551,0.386641976,10.8643,0.403741977,11.3065,0.42161978,12.1223,0.440251979,13.0225,0.459721980,13.6488,0.4800551981,13.2195,0.501291982,13.1879,0.523461983,12.6015,0.546611984,13.2368,0.570791985,14.4121,0.5960351986,15.4973,0.6223951987,16.8802,0.6499251988,18.8163,0.678671989,15.1143,0.708691990,17.5534,0.7400351991,21.8601,0.7727651992,23.8866,0.8069451993,26.9293,0.8426351994,26.8885,0.8799051995,28.8314,0.9188251996,30.0751,0.959461997,30.9535,1.00191998,30.1857,1.046221999,31.5797,1.092492000,32.577569,1.140812001,33.477398,1.033522002,39.021581,1.088672003,41.386432,1.146762004,41.596552,1.207962005,44.657324,1.272422006,46.951775,1.340312007,48.728837,1.401592008,51.488427,1.534092009,50.026967,1.455932010,60.640913,1.498962011,63.3603103378,1.5')fit_air_rice <-lm(Passengers ~ Production, data = air_rice)
Figure 1: Real data: Australian air passengers and Guinea rice production, 1970–2011, World Bank via Hyndman–Athanasopoulos. Both series rise over time, producing a high-fit levels regression between the two measurements.
The fitted line has \(R^2=\) 0.958: it accounts for about 95.8% of the variation in passenger numbers in this sample. Both series rise over time, so later years tend to have both more rice production and more air passengers. A line can therefore fit the observations closely even when the gap around that line is nonstationary. A high \(R^2\) also does not establish that rice production causes air travel.
Spurious regression occurs when a regression gives misleading evidence of a relationship because it ignores trends or persistence. Even independent persistent series can yield a high \(R^2\) and apparently significant coefficients: conventional regression tests can reject far more often than their stated significance levels suggest. The following example uses independent random walks, for which we can check directly whether any fixed nonzero linear combination is stationary.
Example 10.8 (Independent trends). Consider two independent random walks starting at \(A_0=B_0=0\): \[
A_t=A_{t-1}+\varepsilon_t,\qquad
B_t=B_{t-1}+\eta_t,\qquad t\ge1,
\] where \((\varepsilon_t)\) and \((\eta_t)\) are independent sequences of independent \(N(0,1)\) variables. First differencing recovers the shocks: \(\nabla A_t=\varepsilon_t\) and \(\nabla B_t=\eta_t\), so both differences are stationary white noise.
Each level is the sum of \(t\) independent shocks of variance one, so \(\Var(A_t)=\Var(B_t)=t\). Independence also gives \(\Cov(A_t,B_t)=0\). Thus, for any fixed coefficients \((c_1,c_2)\ne(0,0)\), \[
\Var(c_1A_t+c_2B_t)
=c_1^2\Var(A_t)+c_2^2\Var(B_t)
=t(c_1^2+c_2^2).
\] This variance grows with \(t\), whereas weak stationarity requires a constant variance. No nonzero fixed linear combination of the levels is therefore stationary. In particular, there are no fixed constants \(a,b\) for which the gap \(B_t-a-bA_t\) is stationary; subtracting an intercept does not change the variance.
Over a finite record, the independent paths can nevertheless happen to move together, giving a high \(R^2\) when \(b_t\) is regressed on an intercept and \(a_t\) (a levels regression). The usual regression \(t\)-test for a zero slope can then reject far more often than its stated significance level suggests: the nonstationary regression error violates the test’s assumptions. Exercise 10.2 investigates how often independent random walks produce a high \(R^2\) over repeated records. \(\diamond\)
Cointegration and Common Trends
Recall that an \(I(0)\) process is weakly stationary, has absolutely summable autocovariances, and has positive long-run variance. An \(I(1)\) process accumulates mean-zero \(I(0)\) increments: \[
X_t=X_0+\sum_{j=1}^{t}V_j,\qquad t\ge1,
\] where \((V_t)\) is mean-zero and \(I(0)\), and \(X_0\) has finite variance. First differencing recovers these increments: \(\nabla X_t=V_t\).
Here we allow the initial value \(X_0\) to have a nonzero mean. Adding a constant \(c\) to \(X_0\) shifts every \(X_t\) by the same amount; it changes neither the differences nor the variances.
Let \(v_V:=\sum_{h\in\Z}\gamma_V(h)>0\) be the long-run variance of the increments, where \(\gamma_V\) is their ACVF. The partial-sum variance result gives \[
\frac{\Var(X_t)}{t}\longrightarrow v_V.
\] Thus \(\Var(X_t)\) is approximately \(t v_V\) for large \(t\). Its growth rules out weak stationarity, which requires a constant variance. Under this convention, an \(I(1)\) process is therefore not itself \(I(0)\).
Exercise 10.3 (Can stationary series be cointegrated?). Let \((U_t)\) and \((V_t)\) be independent sequences of independent \(N(0,1)\) variables, indexed by \(t\in\mathbb Z\). Define \[
X_t=U_t,\qquad Y_t=2U_t+V_t.
\] A student claims that \(X_t\) and \(Y_t\) are cointegrated because \(Y_t-2X_t=V_t\) is stationary.
Determine whether \((X_t)\), \((Y_t)\), and \((Y_t-2X_t)\) are \(I(0)\). Justify each answer using the definition.
Assess the student’s claim under Definition 10.5.
Show that \((Y_t-bX_t)\) is \(I(0)\) for every fixed \(b\in\mathbb R\). Can stationarity identify the coefficient \(2\)? What distinguishes \(2\) from the other coefficients?
Example 10.9 (Modeling the cointegrating spread). We specialize the construction of Example 10.6. Let \((\eta_t)\), \((\varepsilon_t)\), and \((\xi_t)\) be mutually independent sequences, each i.i.d. with mean zero and positive finite variance \(\sigma_\eta^2\), \(\sigma_\varepsilon^2\), and \(\sigma_\xi^2\), respectively. Set \(S_0=0\) and \(S_t=S_{t-1}+\eta_t\) for \(t\ge1\). Independently of \(S\), let \(U_t\) and \(V_t\) be the stationary causal solutions of \[
U_t=\phi_1U_{t-1}+\varepsilon_t,\qquad
V_t=\phi_2V_{t-1}+\xi_t,\qquad |\phi_1|,|\phi_2|<1.
\] For constants \(a\) and \(b\ne0\), define \[
X_t=S_t+U_t,\qquad Y_t=a+bS_t+V_t.
\] We observe \(X_t\) and \(Y_t\); the common trend and the stationary components are not separately recorded. Neither coordinate equals \(S_t\).
Write \[
q_U:=\Var(U_t)=\frac{\sigma_\varepsilon^2}{1-\phi_1^2},\qquad
q_V:=\Var(V_t)=\frac{\sigma_\xi^2}{1-\phi_2^2}.
\] The level variances are \(\Var(X_t)=t\sigma_\eta^2+q_U\) and \(\Var(Y_t)=b^2t\sigma_\eta^2+q_V\). Their differences are \(I(0)\) with long-run variances \(\sigma_\eta^2\) and \(b^2\sigma_\eta^2\), respectively: the differenced stationary terms \(\nabla U_t\) and \(\nabla V_t\) contribute zero long-run variance. Both levels are therefore \(I(1)\).
For known \(a,b\), the observations give the mean-zero spread \[
E_t:=Y_t-a-bX_t=V_t-bU_t.
\] Independence of \(U\) and \(V\) gives \[
\gamma_E(h)=q_V\phi_2^{|h|}+b^2q_U\phi_1^{|h|},\qquad
\tau_E^2=\frac{\sigma_\xi^2}{(1-\phi_2)^2}
+\frac{b^2\sigma_\varepsilon^2}{(1-\phi_1)^2}>0.
\] Thus \(E_t\) is \(I(0)\), and \((-b,1)^\top\) is a cointegrating vector for \((X_t,Y_t)^\top\).
Compare a fixed regression line in this model with one between the independent random walks of Example 10.8. For a trial slope \(c\), the residual here is \(Y_t-a-cX_t=(b-c)S_t+V_t-cU_t\). The residual variances are \[
\begin{aligned}
\Var(B_t-a-cA_t)&=t(1+c^2),\\
\Var(Y_t-a-cX_t)&=t(b-c)^2\sigma_\eta^2+q_V+c^2q_U.
\end{aligned}
\] For the independent walks, every fixed slope leaves accumulating error. In the shared-trend model, choosing \(c=b\) cancels that accumulation and leaves the stationary error \(E_t\). A high sample \(R^2\) can occur in either example; cointegration requires a fixed population relation with an \(I(0)\) error.
The simulation below uses normal innovations and independent stationary initial values \(U_0\sim N(0,q_U)\) and \(V_0\sim N(0,q_V)\), independent of future innovations.
R code to simulate the common trend and stationary components
R code for the observed levels and their cancelling spread
wb_par(columns =2, simulated =TRUE)matplot(cbind(y_shared, a_shared + b_shared * x_shared), type ="l",col =c(wb_blue, wb_copper), lty =1, lwd =c(1.1, 1.4),xlab ="Time index", ylab ="Level (simulated units)",main ="The common trend in the levels")legend("topleft", c("Y", "a + bX"), col =c(wb_blue, wb_copper),lty =1, lwd =1.5, bty ="n", cex = .85)plot(e_shared, type ="l", col = wb_blue, lwd =1.3,xlab ="Time index", ylab ="Spread (simulated units)",main ="The known spread: Y − a − bX")abline(h =0, col ="gray60", lty =3)par(mfrow =c(1, 1))
Figure 2: Simulated data: 200 observations, a = 2, b = 1.3, phi_1 = 0.35, phi_2 = 0.8, normal innovation SDs 0.7, 0.3, and 0.3, seed 5421005. Both levels contain a random walk and a stationary component. Their known spread is V minus bU.
This model has several filters leading to stationary outputs. For example, both \[
\begin{pmatrix}\nabla X_t\\ \nabla Y_t\end{pmatrix}
=\begin{pmatrix}\eta_t+\nabla U_t\\b\eta_t+\nabla V_t\end{pmatrix}
\quad\text{and}\quad
\begin{pmatrix}\nabla X_t\\Y_t-a-bX_t\end{pmatrix}
=\begin{pmatrix}\eta_t+\nabla U_t\\E_t\end{pmatrix}
\] are jointly weakly stationary. The second records the current gap as one coordinate. Further differencing that gap gives \(\nabla E_t\), which is stationary but has zero long-run variance and hence is outside our \(I(0)\) class. The choice of a stationarity-producing filter affects which quantities are directly available for modeling. Exercise 10.1 gives a three-coordinate example with two independent stationary relations. \(\diamond\)
Prediction from the Stationary Gap
Retain the known model and parameters of Example 10.9. We can compute \(E_t=Y_t-a-bX_t\) from the observed levels. To forecast \(E_{t+1}\) using the current gap alone, consider linear predictors \(rE_t\). Since the spread has mean zero, an intercept cannot reduce the mean squared error.
Expanding the squared error gives \[
\E[(E_{t+1}-rE_t)^2]
=\gamma_E(0)-2r\gamma_E(1)+r^2\gamma_E(0).
\] The minimizing coefficient and forecast are \[
\kappa:=\frac{\gamma_E(1)}{\gamma_E(0)}
=\frac{q_V\phi_2+b^2q_U\phi_1}{q_V+b^2q_U},\qquad
\widehat E^{\mathrm{lin}}_{t+1\mid t}=\kappa E_t,
\] with mean squared error \[
\E[(E_{t+1}-\widehat E^{\mathrm{lin}}_{t+1\mid t})^2]
=\gamma_E(0)-\frac{\gamma_E(1)^2}{\gamma_E(0)}.
\] This is the best linear forecast based on \(E_t\) alone. It is the one-observation case of Proposition 3.8 in the workbook on prediction.
The coefficient \(\kappa\) is a weighted average of \(\phi_1\) and \(\phi_2\), so \(|\kappa|<1\). The predicted change is \((\kappa-1)E_t\): a positive gap has a negative predicted change, and a negative gap has a positive one. Compared with repeating the current gap, the mean squared error is smaller by \((1-\kappa)^2\gamma_E(0)\).
For the parameters used in the simulation, \(\kappa=\) 0.616. The model mean squared error is 0.263, compared with 0.325 for repeating the current gap.
When \(\phi_1\ne\phi_2\), older gaps improve the linear forecast. Indeed, \[
\begin{aligned}
\Cov(E_{t+1}-\kappa E_t,E_{t-1})
&=\gamma_E(2)-\frac{\gamma_E(1)^2}{\gamma_E(0)}\\
&=\frac{b^2q_Uq_V(\phi_1-\phi_2)^2}{q_V+b^2q_U}>0.
\end{aligned}
\] The forecast error is still correlated with an earlier observation. In particular, the spread is not an AR(1) process: its ACVF does not satisfy \(\gamma_E(2)\gamma_E(0)=\gamma_E(1)^2\). If \(\phi_1=\phi_2=\phi\), however, \(E_t=\phi E_{t-1}+\xi_t-b\varepsilon_t\), and the current gap suffices.
For a forecast using several past gaps, apply Proposition 3.8 with the ACVF in Example 10.9. Joint models for \((\nabla X_t,E_t)^\top\) can also use changes in \(X_t\); the workbooks on vector autoregression and matrix filters develop such models.
To forecast \(X_{t+1}\) and \(Y_{t+1}\) separately, we also need a forecast of one coordinate’s change. We predict \(\nabla X_{t+1}\) linearly from the same current gap \(E_t\). The model gives \[
\nabla X_{t+1}=(\phi_1-1)U_t+\eta_{t+1}+\varepsilon_{t+1}.
\] The new innovations are independent of \(E_t\), and \(\Cov(U_t,E_t)=-bq_U\). Hence the best linear forecast of this change from \(E_t\) is \[
\widehat{\nabla X}^{\mathrm{lin}}_{t+1\mid t}
=\frac{\Cov(\nabla X_{t+1},E_t)}{\Var(E_t)}E_t
=\frac{b(1-\phi_1)q_U}{q_V+b^2q_U}E_t.
\]
Adding the predicted change to \(X_t\) gives a forecast of \(X_{t+1}\). We then use \(Y_{t+1}=a+bX_{t+1}+E_{t+1}\) to obtain \[
\begin{aligned}
\widehat X^{\mathrm{lin}}_{t+1\mid t}
&=X_t+\frac{b(1-\phi_1)q_U}{q_V+b^2q_U}E_t,\\
\widehat Y^{\mathrm{lin}}_{t+1\mid t}
&=a+b\widehat X^{\mathrm{lin}}_{t+1\mid t}+\kappa E_t\\
&=Y_t-\frac{(1-\phi_2)q_V}{q_V+b^2q_U}E_t.
\end{aligned}
\] These minimize mean squared error among forecasts of the forms \(X_t+rE_t\) and \(Y_t+sE_t\), respectively: each uses only the current gap to predict the next change.
For \(b>0\), a positive gap predicts an increase in \(X\) and a decrease in \(Y\); the signs reverse for a negative gap. The forecasted levels satisfy \(\widehat Y^{\mathrm{lin}}_{t+1\mid t}-a-b\widehat X^{\mathrm{lin}}_{t+1\mid t}
=\kappa E_t\), agreeing with the gap forecast above.
Estimating and Assessing a Relation
The population spread uses unknown coefficients in an application. Given observations \((x_t,y_t)\), choose a response and estimate its relation by ordinary least squares: \[
(\widehat a,\widehat b)
:=\operatorname*{arg\,min}_{a,b}\sum_{t=1}^T(y_t-a-bx_t)^2,
\qquad \widehat u_t:=y_t-\widehat a-\widehat b x_t.
\] This fixes the coefficient of \(Y\) at one; rescaling all coefficients would describe the same population direction. In R, the calculation is fit <- lm(y ~ x), followed by resid(fit).
Fit the relation in the simulated record
fit_shared <-lm(y_shared ~ x_shared)
In the model of Example 10.9, OLS consistently estimates \(a\) and \(b\). A fixed incorrect slope leaves a multiple of \(S_t\) in the residual; the correct slope cancels it. The estimated slope in the simulated record is 1.291, compared with the true value 1.3. Consistency does not give conventional regression standard errors their usual interpretation.
R code to compare the fitted spread with the known spread
Figure 3: The same simulated record as the common-trend figure. The thick line is the known stationary spread; the thin line uses the estimated intercept and slope. Fitting the relation changes the spread.
For a coefficient \(b_0\) specified before looking at the observations, the spread \(Y_t-b_0X_t\) is an observed scalar series. The level-KPSS procedure tests whether it is \(I(0)\) around an unknown constant mean, subject to the same dependence and bandwidth conditions as before. Rejection rules against this particular coefficient choice; it does not rule out every possible cointegrating relation.
The OLS spread has instead been chosen from the observations. Even under Example 10.9, its fitted residuals satisfy \[
\widehat u_t
=e_t-(\widehat a-a)-(\widehat b-b)x_t,\qquad e_t:=y_t-a-bx_t.
\] Its sample mean is zero because an intercept was fitted. Neither that constraint nor a small plotted spread establishes a stationary population relation. Ordinary KPSS critical values do not account for estimating the slope from the same record. We use OLS to estimate a candidate relation and examine its later behavior; formal testing of the estimated relation requires a separate calibration.
Oil and Gasoline
The astsa::oil series records weekly West Texas Intermediate crude-oil spot prices in dollars per barrel; astsa::gas records New York Harbor regular gasoline spot prices in cents per gallon. These 545 matched EIA observations run from 2000 to mid-2010. Crude oil is an input to gasoline production, which motivates examining their relative prices. Write \(O_t\) and \(G_t\) for the recorded oil and gasoline prices. We specify the unit-slope log relation first: \[
R_t:=\log G_t-\log O_t.
\] Changing the recorded units adds a constant to \(R_t\), leaving its level-stationarity property unchanged.
Use weeks 1–350 for assessment and fitting, leaving weeks 351–545 for a later check. On the first period, apply level KPSS to \(R_t\), using the Bartlett estimate with bandwidth \(H=8\) and the supplied 5% cutoff \(0.463\). The null is that this fixed log ratio is \(I(0)\) around a constant mean. If both log prices are \(I(1)\), this null specifies a particular cointegrating relation.
Load the prices, test the fixed ratio, and fit the general relation
# Level KPSS: estimate a constant mean, use divisor n at every lag, and# include lags 1,...,H with Bartlett weights 1-h/(H+1). Here n corresponds# to the mathematical sample size T.kpss_level <-function(x, H =8L) {stopifnot(is.numeric(x), is.null(dim(x)), length(x) >=4L,all(is.finite(x)), length(H) ==1L, is.finite(H), H >=0L, H ==floor(H), H <length(x)) n <-length(x) e <-as.numeric(x) -mean(x) gamma <-vapply(0:H, function(h) {sum(e[seq_len(n - h)] * e[h +seq_len(n - h)]) / n }, numeric(1)) weights <-1- (0:H) / (H +1) lrv <-2*sum(weights * gamma) - gamma[1L]if (!is.finite(lrv) || lrv <=0) {stop("KPSS requires a positive estimated long-run variance.") } partial_sums <-cumsum(e) statistic <-sum(partial_sums^2) / (n^2* lrv)list(statistic = statistic, lrv = lrv, bandwidth = H, n = n,critical_5pct =0.463, reject_5pct = statistic >0.463,residuals = e, partial_sums = partial_sums,autocovariances = gamma)}data("oil", package ="astsa")data("gas", package ="astsa")stopifnot(length(oil) ==545L, length(gas) ==length(oil),identical(tsp(oil), tsp(gas)), all(oil >0), all(gas >0))log_oil <-log(as.numeric(oil))log_gas <-log(as.numeric(gas))fuel_time <-as.numeric(time(oil))formation_end <-350Lformation <-seq_len(formation_end)later <- (formation_end +1L):length(log_oil)# The coefficient one is prespecified; KPSS is applied only to this ratio.fixed_spread <- log_gas - log_oilfixed_kpss <-kpss_level(fixed_spread[formation], H =8L)# Fit on weeks 1-350 and retain the same coefficients for every later week.fuel_fit <-lm(log_gas ~ log_oil, subset = formation)fuel_a <-unname(coef(fuel_fit)[1])fuel_b <-unname(coef(fuel_fit)[2])fuel_spread <- log_gas - fuel_a - fuel_b * log_oilspread_summary <-function(z) {c(Mean =mean(z), SD =sd(z),`Lag-one ACF`=as.numeric(acf(z, lag.max =1, plot =FALSE)$acf[2]))}fuel_windows <-rbind(`Weeks 1-350: fitting`=spread_summary(fuel_spread[formation]),`Weeks 351-545: later`=spread_summary(fuel_spread[later]))
The KPSS statistic is 0.085, below \(0.463\). We do not reject stationarity of the prespecified log ratio at 5%. This is evidence compatible with the proposed relation over these weeks, subject to the KPSS assumptions; nonrejection does ofcourse not establish the null.
Now fit \(\log G_t=a+b\log O_t+U_t\) by OLS on the same 350 weeks. The estimates are \(\widehat a=\) 1.051 and \(\widehat b=\) 0.995. The estimated slope is close to the prespecified value one. For every later observation, keep these estimates fixed and compute \(\log G_t-\widehat a-\widehat b\log O_t\).
R code for the frozen relation and its later spread
wb_par(rows =2)matplot(fuel_time, cbind(log_gas, fuel_a + fuel_b * log_oil),type ="l", col =c(wb_blue, wb_copper), lty =1, lwd =1.3,xlab ="", ylab ="Log price (cents/gallon)",main ="Gasoline and the relation fitted to oil")abline(v = fuel_time[formation_end], lty =3, col = wb_gray)legend("topleft", c("Observed log gasoline", "Frozen oil-based relation"),col =c(wb_blue, wb_copper), lty =1, bty ="n", cex = .8)fixed_centered <- fixed_spread -mean(fixed_spread[formation])matplot(fuel_time, cbind(fuel_spread, fixed_centered), type ="l",col =c(wb_blue, wb_copper), lty =c(1, 2), lwd =c(1.5, 1.0),xlab ="Year (dataset index)", ylab ="Log spread",main ="The same coefficients in the later period")abline(v = fuel_time[formation_end], lty =3, col = wb_gray)abline(h =0, lty =3, col ="gray65")legend("topleft", c("Frozen OLS spread", "Fixed ratio, centered"),col =c(wb_blue, wb_copper), lty =c(1, 2), bty ="n", cex = .8,horiz =TRUE)par(mfrow =c(1, 1))
Figure 4: Real data: weekly WTI crude oil and New York Harbor gasoline spot prices, EIA via astsa. The vertical line ends the 350-week fitting period. The lower panel compares the frozen OLS spread with the prespecified log ratio centered using its first-period mean; the two curves are close because the fitted slope is close to one.
Mean
SD
Lag-one ACF
Weeks 1-350: fitting
0.000
0.074
0.832
Weeks 351-545: later
-0.052
0.097
0.898
The first-period fitted spread has mean zero by construction. Its later mean is -0.052, and its later SD is 1.3 times its fitting-period SD. These changes can reflect sampling variation, coefficient error, or a changing relation. The comparison describes how the frozen relation behaves beyond its fitting period; it does not separate those explanations.
The fixed-ratio test, the fitted slope, and the later spread answer different questions. The first checks one proposed stationary relation; the second estimates a more general relation; the third examines its behavior on later observations. Prediction additionally requires a model for how the gap evolves. A trading strategy would also require positions, timing, and transaction costs.
Exercises
Exercise 10.1 (Independent stationary relations). Let \(Q_0=0\) and \(Q_t=Q_{t-1}+\eta_t\) for \(t\ge1\). Suppose \((\eta_t)\), \((U_t)\), and \((V_t)\) are mutually independent sequences of independent \(N(0,1)\) variables. Define \[
X_t=(Q_t+U_t,\ Q_t+V_t,\ Q_t)^\top.
\]
Find every nonzero \(c\in\mathbb R^3\) for which \(c^\top X_t\) is weakly stationary. Give its variance and ACVF.
Find two linearly independent coefficient vectors with this property and express every stationary combination using their two spreads. Does \(X_{1t}-X_{2t}\) supply a third independent stationary relation?
Exercise 10.2 (What does a high \(R^2\) establish?). Compare two models for a pair \((X_t,Y_t)\). In both, \(X_0=0\) and \(X_t=X_{t-1}+\eta_t\). In Model A, \(Y_0=0\) and \(Y_t=Y_{t-1}+\zeta_t\). In Model B, \(Y_t=1+1.5X_t+U_t\). The sequences \((\eta_t)\), \((\zeta_t)\), and \((U_t)\) are mutually independent, with independent \(N(0,1)\) variables.
Using seed 542102, simulate 1,000 independent records of length \(T=200\). Within each record, use the same \(X_t\) for both models and fit each \(y_t\) on an intercept and \(x_t\). For the first record, report both fitted slopes and \(R^2\) values, and plot both residual series; overlay \(U_t\) on the Model B plot.
For each model, report the median fitted slope, median \(R^2\), and fraction \(\hat p\) of records with \(R^2>1/2\). Give the Monte Carlo standard error \(\sqrt{\hat p(1-\hat p)/1000}\) of each fraction. Explain why a high \(R^2\) does not establish cointegration. In Model B, distinguish the known stationary spread from the fitted residuals.
Exercise 10.4 (Can a small extra trend destroy cointegration?). Let \(S_0=R_0=0\) and \[
S_t=S_{t-1}+\eta_t,\qquad R_t=R_{t-1}+\xi_t,\qquad t\ge1.
\] Suppose \((\eta_t)\), \((\xi_t)\), and \((U_t)\) are mutually independent sequences of independent \(N(0,1)\) variables, with \(U_t\) defined also at \(t=0\). For a fixed \(\delta\ge0\), define \[
X_t=S_t,\qquad Y_t=2S_t+U_t+\delta R_t.
\]
Set \(\delta=0\). Find a cointegrating vector for \((X_t,Y_t)^\top\). Explain which accumulating component cancels in its spread.
Now suppose \(\delta>0\). For an arbitrary fixed slope \(b\), calculate \(\Var(Y_t-bX_t)\). Can any slope make this gap stationary? Are the coordinates cointegrated?
Using seed 542104, simulate one record of length \(T=5000\) for \(\delta=0\) and \(\delta=0.02\), using the same realizations of \(S_t\), \(R_t\), and \(U_t\) in both cases. Plot the gaps \(Y_t-2X_t\) over the first 200 observations and over the full record, using the same vertical scale. Describe what changes when the record is extended. Does a stationary-looking gap establish cointegration? Use (b) to explain why a small nonzero \(\delta\) can be difficult to detect in a short record.
Kwiatkowski, Phillips, Schmidt, and Shin (1992) give the level-KPSS reference used for the fixed spread. The folded R code implements the course’s divisor-\(T\) sample covariances and Bartlett weights \(1-h/(H+1)\). The oil and gasoline observations are distributed with astsa, with original spot-price data from the U.S. Energy Information Administration. Shumway and Stoffer use the same pair in their time-series regression examples.
