Time Series – Choosing the Correct Dickey-Fuller Test for Models with Intercept and Linear Trend

augmented-dickey-fullerstationaritytime seriesunit root

Short version:

I have a time series of climate data that I'm testing for stationarity. Based on previous research, I expect the model underlying (or "generating", so to speak) the data to have an intercept term and a positive linear time trend. To test these data for stationarity, should I use the Dickey-Fuller test that includes an intercept and time trend, i.e. equation #3?

$\nabla y_t = \alpha_0+\alpha_1t+\delta y_{t-1}+u_t$

Or, should I use the DF test that only includes an intercept because the first difference of the equation I believe underlies the model has only an intercept?

Long version:

As stated above, I have a time series of climate data that I'm testing for stationarity. Based on previous research, I expect the model underlying the data to have an intercept term, a positive linear time trend, and some normally distributed error term. In other words, I expect the underlying model to look something like this:

$y_t = a_0 + a_1t + \beta y_{t-1} + u_t $

where $u_t$ is normally distributed. Since I'm assuming the underlying model has both an intercept and a linear time trend, I tested for a unit root with equation #3 of the simple Dickey-Fuller test, as shown:

$\nabla y_t = \alpha_0+\alpha_1t+\delta y_{t-1}+u_t$

This test returns a critical value that would lead me to reject the null hypothesis and conclude that the underlying model is non-stationary. However, I question if I'm applying this correctly, since even though the underlying model is assumed to have an intercept and a time trend, this does not imply that the first difference $\nabla y_t$ will as well. Quite the opposite, in fact, if my math is correct.

Calculating the first difference based on the equation of the assumed underlying model gives:
$\nabla y_t = y_t – y_{t-1} = [a_0 + a_1t + \beta y_{t-1} + u_t] – [a_0 + a_1(t-1) + \beta y_{t-2} + u_{t-1}]$

$\nabla y_t = [a_0 – a_0] + [a_1t – a_t(t-1)] + \beta[y_{t-1} – y_{t-2}] + [u_t – u_{t-1}]$

$\nabla y_t = a_1 + \beta \cdot \nabla y_{t-1} + u_t – u_{t-1}$

Therefore, the first difference $\nabla y_t$ appears to only have an intercept, not a time trend.

I think my question is similar to this one, except I'm not sure how to apply that answer to my question.

Sample data:

Here is some of the sample temperature data that I'm working with.

Best Answer

You need to consider the drift and (parametric/linear) trend in the levels of the time series in order to specify the deterministic terms in the augmented Dickey-Fuller regression which is in terms of the first differences of the time series. The confusion arises exactly from deriving the first-differences equation in the way that you have done.

(Augmented) Dickey-Fuller regression model

Suppose that the levels of the series include a drift and trend term $$ Y_t = \beta_{0,l} + \beta_{1,l} t + \beta_{2, l}Y_{t-1} + \varepsilon_{t} $$ The null hypothesis of nonstationarity in this case would be $\mathfrak{H}_0{}:{}\beta_{2, l} = 1$.

One equation for the first differences implied by this data-generating process [DGP], is the one that you have derived $$ \Delta Y_t = \beta_{1,l} + \beta_{2, l}\Delta Y_{t-1} + \Delta \varepsilon_{t} $$ However, this is not the (augmented) Dickey Fuller regression as used in the test.

Instead, the correct version can be had by subtracting $Y_{t-1}$ from both sides of the first equation resulting in $$ \begin{align} \Delta Y_t &= \beta_{0,l} + \beta_{1,l} t + (\beta_{2, l}-1)Y_{t-1} + \varepsilon_{t} \\ &\equiv \beta_{0,d} + \beta_{1,d}t + \beta_{2,d}Y_{t-1} + \varepsilon_{t} \end{align} $$ This is the (augmented) Dickey-Fuller regression, and the equivalent version of the null hypothesis of nonstationarity is the test $\mathfrak{H}_0{}:{}\beta_{2, d}=0$ which is just a t-test using the OLS estimate of $\beta_{2, d}$ in the regression above. Note that the drift and trend come through to this specification unchanged.

An additional point to note is that if you are not certain about the presence of the linear trend in the levels of the time series, then you can jointly test for the linear trend and unit root, that is, $\mathfrak{H}_0{}:{}[\beta_{2, d}, \beta_{1,l}]' = [0, 0]'$, which can be tested using an F-test with appropriate critical values. These tests and critical values are produced by the R function ur.df in the urca package.

Let us consider some examples in detail.

Examples

1. Using the US investment series

The first example uses the US investment series which is discussed in Lutkepohl and Kratzig (2005, pg. 9). The plot of the series and its first difference are given below.

enter image description here

From the levels of the series, it appears that it has a non-zero mean, but does not appear to have a linear trend. So, we proceed with an augmented Dickey Fuller regression with an intercept, and also three lags of the dependent variable to account for serial correlation, that is: $$ \Delta Y_t = \beta_{0,d} + \beta_{2,d}Y_{t-1} + \sum_{j=1}^3 \gamma_j \Delta Y_{t-j} + \varepsilon_{t} $$ Note the key point that I have looked at the levels to specify the regression equation in differences.

The R code to do this is given below:

    library(urca)
    library(foreign)
    library(zoo)

    tsInv <- as.zoo(ts(as.data.frame(read.table(
      "http://www.jmulti.de/download/datasets/US_investment.dat", skip=8, header=TRUE)), 
                       frequency=4, start=1947+2/4))
    png("USinvPlot.png", width=6,
        height=7, units="in", res=100)
    par(mfrow=c(2, 1))
    plot(tsInv$USinvestment)
    plot(diff(tsInv$USinvestment))
    dev.off()

    # ADF with intercept
    adfIntercept <- ur.df(tsInv$USinvestment, lags = 3, type= 'drift')
    summary(adfIntercept)

The results indicate that the the null hypothesis of nonstationarity can be rejected for this series using the t-test based on the estimated coefficient. The joint F-test of the intercept and the slope coefficient ($\mathfrak{H}{}:{}[\beta_{2, d}, \beta_{0,l}]' = [0, 0]'$) also rejects the null hypothesis that there is a unit root in the series.

2. Using German (log) consumption series

The second example is using the German quarterly seasonally adjusted time series of (log) consumption. The plot of the series and its differences are given below.

enter image description here

From the levels of the series, it is clear that the series has a trend, so we include the trend in the augmented Dickey-Fuller regression together with four lags of the first differences to account for the serial correlation, that is $$ \Delta Y_t = \beta_{0,d} + \beta_{1,d}t + \beta_{2,d}Y_{t-1} + \sum_{j=1}^4 \gamma_j \Delta Y_{t-j} + \varepsilon_{t} $$

The R code to do this is

# using the (log) consumption series
tsConsump <- zoo(read.dta("http://www.stata-press.com/data/r12/lutkepohl2.dta"), frequency=1)
png("logConsPlot.png", width=6,
    height=7, units="in", res=100)
par(mfrow=c(2, 1))
plot(tsConsump$ln_consump)
plot(diff(tsConsump$ln_consump))
dev.off()

# ADF with trend
adfTrend <- ur.df(tsConsump$ln_consump, lags = 4, type = 'trend')
summary(adfTrend)

The results indicate that the null of nonstationarity cannot be rejected using the t-test based on the estimated coefficient. The joint F-test of the linear trend coefficient and the slope coefficient ($\mathfrak{H}{}:{}[\beta_{2, d}, \beta_{1,l}]' = [0, 0]'$) also indicates that the null of nonstationarity cannot be rejected.

3. Using given temperature data

Now we can assess the properties of your data. The usual plots in levels and first differences are given below.

enter image description here

These indicate that your data has an intercept and a trend, so we perform the ADF test (with no lagged first difference terms), using the following R code

# using the given data
tsTemp <- read.table(textConnection("temp 
64.19749  
65.19011  
64.03281  
64.99111  
65.43837  
65.51817  
65.22061  
65.43191  
65.0221  
65.44038  
64.41756  
64.65764  
64.7486  
65.11544  
64.12437  
64.49148  
64.89215  
64.72688  
64.97553  
64.6361  
64.29038  
65.31076  
64.2114  
65.37864  
65.49637  
65.3289  
65.38394  
65.39384  
65.0984  
65.32695  
65.28  
64.31041  
65.20193  
65.78063  
65.17604  
66.16412  
65.85091  
65.46718  
65.75551  
65.39994  
66.36175  
65.37125  
65.77763  
65.48623  
64.62135  
65.77237  
65.84289  
65.80289  
66.78865  
65.56931  
65.29913  
64.85516  
65.56866  
64.75768  
65.95956  
65.64745  
64.77283  
65.64165  
66.64309  
65.84163  
66.2946  
66.10482  
65.72736  
65.56701  
65.11096  
66.0006  
66.71783  
65.35595  
66.44798  
65.74924  
65.4501  
65.97633  
65.32825  
65.7741  
65.76783  
65.88689  
65.88939  
65.16927  
64.95984  
66.02226  
66.79225  
66.75573  
65.74074  
66.14969  
66.15687  
65.81199  
66.13094  
66.13194  
65.82172  
66.14661  
65.32756  
66.3979  
65.84383  
65.55329  
65.68398  
66.42857  
65.82402  
66.01003  
66.25157  
65.82142  
66.08791  
65.78863  
66.2764  
66.00948  
66.26236  
65.40246  
65.40166  
65.37064  
65.73147  
65.32708  
65.84894  
65.82043  
64.91447  
65.81062  
66.42228  
66.0316  
65.35361  
66.46407  
66.41045  
65.81548  
65.06059  
66.25414  
65.69747  
65.15275  
65.50985  
66.66216  
66.88095  
65.81281  
66.15546  
66.40939  
65.94115  
65.98144  
66.13243  
66.89761  
66.95423  
65.63435  
66.05837  
66.71114"), header=T)
tsTemp <- as.zoo(ts(tsTemp, frequency=1))

png("tempPlot.png", width=6,
    height=7, units="in", res=100)
par(mfrow=c(2, 1))
plot(tsTemp$temp)
plot(diff(tsTemp$temp))
dev.off()

# ADF with trend
adfTrend <- ur.df(tsTemp$temp, type = 'trend')
summary(adfTrend)

The results for both the t-test and the F-test indicate that the null of nonstationarity can be rejected for the temperature series. I hope that clarifies matter somewhat.

Related Solutions

Augmented Dickey-Fuller Test – Difference Among None, Drift, and Trend

The Wikipedia page states the following:

The testing procedure for the ADF test is the same as for the Dickey–Fuller test but it is applied to the model $$ \Delta y_t = \alpha + \beta t + \gamma y_{t-1} + \delta_1 \Delta y_{t-1} + \cdots + \delta_{p-1} \Delta y_{t-p+1} + \varepsilon_t $$

As you very well note, there are variations of the test, which involve restricting $\alpha$ and/or $\beta$ equal to 0. Imposing the restriction on $\alpha$ corresponds to omitting a constant while restricting $\beta$ corresponds to omitting a time trend.

To understand what you're doing when using the adf.test() function from the tseries package in R, we should first consult the documentation provided by the package authors. To do this, we execute ?adf.test in the R console. Doing this will provide us details about the function; what it does, how we can use it, etc. For present purposes, we just need to be aware that the documentation states:

The general regression equation which incorporates a constant and a linear trend is used and the t-statistic for a first order autoregressive coefficient equals one is computed.

(Do we need more information than that?)

Coupled with that fact, if we look at the usage of the function; namely,

adf.test(x, alternative = c("stationary", "explosive"),
         k = trunc((length(x)-1)^(1/3)))

one begins to think that the function has limited capabilities with regard to the restricted variations of the ADF test. Reading all of the documentation seems to make it clear that the function only runs one variation of the test; the unrestricted version, which includes both a constant and a trend.

(Do we need more information than that?)

Since you're using R, we don't have to be left wondering if the function somehow imposes the restrictions internally without us knowing! To really be sure what's going on behind the scenes, we can look at the source code of the adf.test() function. Below, I step through the code, which I have shortened, and I hope it's instructive to you.

# Import some toy data
data(sunspots)

# Set arguments that are normally function inputs
x           <- sunspots
alternative <- "stationary"
k           <- trunc((length(x) - 1)^(1/3))

# Let the function go to work! (short version)
k <- k + 1          # Number of lagged differenced terms
y <- diff(x)        # First differences
n <- length(y)      # Length of first differenced series
z <- embed(y, k)    # Used for creating lagged series

# Things get interesting here as variables are prepared for the regression
yt  <- z[, 1]       # First differences
xt1 <- x[k:n]       # Series in levels - the first k-1 observations are dropped
tt  <- k:n          # Time-trend
yt1 <- z[, 2:k]     # Lagged differenced series - there are k-1 of them

# Next, the key pieces of code.

# Regression 1: if k > 0
# The augmented Dickey-Fuller test (with constant and time-trend)
res <- lm(yt ~ xt1 + 1 + tt + yt1) 

# Regression 2: if k = 0
# The standard Dickey-Fuller test (with constant and time-trend)
res <- lm(yt ~ xt1 + 1 + tt)

By my count, the adf.test() function is, in fact, made up of 57 lines of code, which I encourage you to inspect. The rest of the function code is not important in the context of this question. All that needs to be known is that the function does do what it says on the tin. Importantly, there does not seem to be a high level way of using the function to run a restricted variation of the ADF test and retrieve the associated critical values.

What to do? Your first instinct should be to check out the CRAN Task View: Time Series Analysis page. In doing so, you'll learn that the urca package provides an alternative implementation of the ADF test. Indeed, as I mentioned in the comments, the ur.df() function should be able to meet your needs. Inspecting the function usage is quite informative!

ur.df(y, type = c("none", "drift", "trend"), lags = 1, 
      selectlags = c("Fixed", "AIC", "BIC"))

The urca package can be found here and I recommend consulting the package documentation and the source code if you need to. I suspect that you should be able to use the function and not worry about issues regarding critical values; the authors of the package will have taken care of that so you can concentrate on using it as a high-level function and doing your research.

In terms of applying the ADF test (knowing which tests to run and in which order), I would suggest the Dolado et al. procedure. The reference is:

Dolado, J. J., Jenkinson, T., and Sosvilla-Rivera, S. (1990). Cointegration and unit roots, Journal of Economic Surveys, 4, 249-273.

Final note on matching the R code to the mathematical equation. You can basically think of it as follows (strictly speaking, the parameters should be omitted, but...):

yt = $\Delta y_{t}$

xt = $\gamma y_{t-1}$

+ 1 = $\alpha$

tt = $\beta t$

yt1 = $\delta_1 \Delta y_{t-1} + \cdots + \delta_{p-1} \Delta y_{t-p+1}$

Solved – Augmented Dickey-Fuller Trend and Intercept

library(tseries)
library(urca)
data2<-na.remove(data1$SANDP) #assumed one dataset
    data2
    c<-adf.test(data2,k=2)
    c
    summary(c$statistic)

useful command for you

summary(ur.df(data2,type=c("drift"),lags=2,selectlags =c("Fixed")))

Short version:

Long version:

Sample data:

Best Answer

(Augmented) Dickey-Fuller regression model

Examples

1. Using the US investment series

2. Using German (log) consumption series

3. Using given temperature data

Related Solutions

Augmented Dickey-Fuller Test – Difference Among None, Drift, and Trend

Solved – Augmented Dickey-Fuller Trend and Intercept

Related Question