# Can I trust the R^2 from an OLS regression of two cointegrated series?

Suppose $x_t$ is a non-stationary time series and $v_t$ is white noise. $$x_t = \Sigma{v_t}$$

set.seed(123)
x <- cumsum(rnorm(100,10,10))


Also suppose $y_t$ is co-integrated with $x_t$. Where $x_t$ is integrated of order 1 and $\epsilon_t$ fulfills all the classic OLS assumptions. $$y_t = \beta_0 + \beta_1x_t + \epsilon_t$$

y <- 50 + 5*x + rnorm(100, 0, 100)


Augmented Dickey-Fuller reveals that both $x_t$ and $y_t$ are non-stationary.

tseries::adf.test(x)


However, I know that $x_t$ and $y_t$ are co-integrated, so I fit an OLS model.

lm.mod.coint <- lm(y ~ x)
summary(lm.mod.coint)

Call:
lm(formula = y ~ x)

Residuals:
Min      1Q  Median      3Q     Max
-189.98  -66.85   -8.78   56.51  331.14

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 24.65687   19.47269   1.266    0.208
x           5.02754    0.03189 157.636   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 96.82 on 98 degrees of freedom
Multiple R-squared:  0.9961,    Adjusted R-squared:  0.996
F-statistic: 2.485e+04 on 1 and 98 DF,  p-value: < 2.2e-16


A Phillips–Ouliaris Cointegration Test confirms that the model is indeed co-integrated.

tseries::po.test(lm.mod.coint$model)  I have several questions: 1. Did I specify the functional form correctly for cointegration? $$y_t = \beta_0 + \beta_1x_t + \epsilon_t$$ 2. Can I trust the$R^2$value from the OLS output? The$R^2$is very high. 3. Can I trust the p-values from the OLS output? When I fit a model of first differences, the model fit looks horrible. Is this a result of overdifferencing? df.coint <- lm.mod.coint$model

df.coint.diff <- as.data.frame(diff(as.matrix(df.coint)))  ## calc first-diffs

names(df.coint.diff) <- paste0('diff_', names(df.coint))

lm.diff.coint <- lm(diff_y ~ diff_x - 1, data=df.coint.diff)

summary(lm.diff.coint)

• Google for super-consistency and cointegration. And as $T \rightarrow \infty$, what do you think the estimated $R^2$ will converge on? – Matthew Gunn May 20 '17 at 9:04
• thank you. I found these slides very helpful. econ.ku.dk/metrics/econometrics2_05_ii/slides/… – William Chiu May 20 '17 at 16:52
• Those slides look quite good. Speaking loosely, the variation from the non-stationary, I(1) component is infinite. The variation from the stationary $\epsilon_t$ is finite. So "unexplained" variation as a fraction of total variation goes to $\frac{\sigma^2_\epsilon}{\infty} = 0$ and the "explained" fraction of variation (i.e. estimated $R^2$) goes to 1. – Matthew Gunn May 20 '17 at 18:12

1: Yes. Your data are generated exactly according to this equation. You cannot find a so perfect model in the real life.

2: Yes. Plot a scatterplot of Y vs X, you will find the points are on the straight line. So it is not strange to have so high R square.

3: Yes. Why not?

BTW, OLS is an algorithm or estimator, not a model.

• A scatter plot of $x$ and $y$ would also show a straight line if $x$ and $y$ were not cointegrated. If there was no cointegration, the correlation is spurious. set.seed(123); x <- cumsum(rnorm(100,10,10)); set.seed(321); y <- cumsum(rnorm(100,5,5)); cor(x,y)^2; plot(x,y) – William Chiu May 20 '17 at 8:07
• In your model, there are 3, and only 3 parameters (50, 5, 100). OLS estimator successfully guessed them based on the provided data. What else you want to do by using statistics? – user158565 May 20 '17 at 13:29
• Based on these slides (econ.ku.dk/metrics/econometrics2_05_ii/slides/…) as T approaches infinity, the OLS estimate of the intercept converges to 0. It also appears that although the estimate of the slope is super-consistent, the p-values are not reliable. – William Chiu May 20 '17 at 16:55
• Did you read the slide #2? If yes, I am speechless. – user158565 May 21 '17 at 17:31
• slide two says the results (including $R^2$) are nonsense when the residuals are non-stationary. However, in my case, the residuals are stationary (as demonstrated by the po.test). – William Chiu May 22 '17 at 19:07

The answer by a_statistician answers your particular questions, but I find that when you're working with theoretically tractable models like this, it is worth having a look at the underlying relationships to see why you are getting the simulation results that occur.

Analysis of initial series: In order to see why you are getting the OLS results of your regression simulation, it is worth looking at the theoretical model underneath, and the true relationships between the observable series. Letting $v_t \sim \text{N}(\mu, \sigma^2)$ be your white noise, your Brownian motion process $\{ x_t \}$ has mean $\mathbb{E}(x_t) = t \mu$ and variance $\mathbb{V}(x_t) = t \sigma^2$ which gives the covariance result:

\begin{aligned} \mathbb{C}(y_t, x_t) = \mathbb{E}(y_t x_t) -\mathbb{E}(y_t) \mathbb{E}(x_t) &= \mathbb{E}((\beta_0 + \beta_1 x_t + \epsilon_t) x_t) -\mathbb{E}(\beta_0 + \beta_1 x_t + \epsilon_t) \mathbb{E}(x_t) \\[6pt] &= \mathbb{E}(\beta_0 x_t + \beta_1 x_t^2 + \epsilon_t x_t) -\mathbb{E}(\beta_0 + \beta_1 x_t + \epsilon_t) \mathbb{E}(x_t) \\[6pt] &= \beta_0 \mathbb{E}(x_t) + \beta_1 \mathbb{E}(x_t^2) + \mathbb{E}(\epsilon_t x_t) - (\beta_0 + \beta_1 \mathbb{E}(x_t) + \mathbb{E}(\epsilon_t)) \mathbb{E}(x_t) \\[6pt] &= \beta_0 t \mu + \beta_1 (t \sigma^2 + t^2 \mu^2) + 0 - (\beta_0 + \beta_1 t \mu + 0) t \mu \\[6pt] &= \beta_0 t \mu + \beta_1 t \sigma^2 + \beta_1 t^2 \mu^2 - \beta_0 t \mu - \beta_1 t^2 \mu^2 \\[6pt] &= \beta_1 t \sigma^2. \\[6pt] \end{aligned}

Since $\mathbb{V}(y_t) = \mathbb{V}(\beta_0 + \beta_1 x_t + \epsilon_t) = \beta_1^2 t \sigma^2 + \sigma_\epsilon^2$ we then have squared-correlation:

\begin{aligned} \mathbb{Corr}(y_t, x_t)^2 = \frac{\mathbb{C}(y_t, x_t)^2}{\mathbb{V}(y_t)\mathbb{V}(x_t)} &= \frac{\beta_1^2 t^2 \sigma^4}{(\beta_1^2 t \sigma^2 + \sigma_\epsilon^2) t \sigma^2} \\[6pt] &= \frac{\beta_1^2 t \sigma^2}{\beta_1^2 t \sigma^2 + \sigma_\epsilon^2} = \frac{\phi_t}{1 + \phi_t}, \\[6pt] \end{aligned}

where $\phi_t \equiv \beta_1^2 t \cdot \sigma^2 / \sigma_\epsilon^2$. Now, in your analysis you have $\beta_1 = 5$, $\sigma^2 = 10$ and $\sigma_\epsilon^2 = 100$, which gives you $\phi_t = 2.5t$. Hence, for even moderately large $t$, the correlation between the co-integrated series is high. It is unsurprising that simulation of this model gives a high coefficient-of-determination. You can lower this correlation by lowering $\beta_1$ or $\sigma$, or by raising $\sigma_\epsilon$, and when you repeat your simulations this ought to lower the coefficient-of-determination.

Analysis of first-differences: When you take first-differences you obtain the model:

$$\Delta y_t = \beta_1 \Delta x_t + \Delta \epsilon_t = \beta_1 v_t + \nu_t,$$

where $v_t \sim \text{N}(0, \sigma_v^2)$ and $\nu_t \sim \text{N}(0, 2 \sigma_\epsilon^2)$ are independent white noise processes. In this case the true covariance is:

$$\mathbb{C}(\Delta y_t, \Delta x_t) = \mathbb{C}(\beta_1 v_t + \nu_t, v_t) = \beta_1 \sigma^2.$$

Since $\mathbb{V}(\Delta y_t) = \mathbb{V}(\beta_1 v_t + \nu_t) = \beta_1^2 \sigma^2 + 2 \sigma_\epsilon^2$ we then have squared-correlation:

\begin{aligned} \mathbb{Corr}(\Delta y_t, \Delta x_t)^2 = \frac{\mathbb{C}(\Delta y_t, \Delta x_t)^2}{\mathbb{V}(\Delta y_t)\mathbb{V}(\Delta x_t)} &= \frac{\beta_1^2 \sigma^4}{(\beta_1^2 \sigma^2 + 2 \sigma_\epsilon^2) \sigma^2} \\[6pt] &= \frac{\beta_1^2 \sigma^2}{\beta_1^2 \sigma^2 + 2 \sigma_\epsilon^2} = \frac{\phi_1}{2 + \phi_1}. \\[6pt] \end{aligned}

You can see that this gives a lower squared-correlation than for the initial series. With your values this gives the true squared-correlation of 56%. You should still be getting a moderate coefficient-of-determination under this model.