Re: [R] Dickey-Fuller Test

From: Clark Allan <Allan_at_stats.uct.ac.za>
Date: Mon 23 May 2005 - 20:49:13 EST

hello

i ammended the code. it was in some package. dont know which? try tseries. this might be of help.

adf.test.1<-function (x,kind=3,k = trunc((length(x)- 1)^(1/3))) {

#kind = the kind of test undertaken
#kind = 1 ==> No constant no trend
#kind = 2 ==> Constant
#kind = 3 ==> Constant and trend

#the null is ALWAYS non stationarity

    if (NCOL(x) > 1)

        stop("x is not a vector or univariate time series")

    if (any(is.na(x)))

        stop("NAs in x")

    if (k < 0)

        stop("k negative")

    DNAME <- deparse(substitute(x))

    k <- k + 1
    y <- diff(x)
    n <- length(y)
    z <- embed(y, k)

    yt <- z[,1]
    xt1 <- x[k:n]
    tt <- k:n

    if (kind==1)
    {
    table <- cbind(c(2.66, 2.62, 2.6, 2.58, 2.58, 2.58), c(2.26,

        2.25, 2.24, 2.23, 2.23, 2.23), c(1.95, 1.95, 1.95, 1.95, 
        1.95, 1.95), c(1.60, 1.61, 1.61, 1.62, 1.62, 1.62), c(0.92, 
        0.91, 0.90, 0.89, 0.89, 0.89), c(1.33, 1.31, 1.29, 1.29, 
        1.28, 1.28), c(1.70, 1.66, 1.64, 1.63, 1.62, 1.62), c(2.16, 
        2.08, 2.03, 2.01, 2.00, 2.00))

    if (k > 1)
    {

        yt1 <- z[,2:k]
        res <- lm(yt ~ xt1 - 1 + yt1)

    }
    else res <- lm(yt ~ xt1-1)
    res.sum <- summary(res)
    STAT <- res.sum$coefficients[1,1]/res.sum$coefficients[1,2]     

    }

    if (kind==2)
    {
    table <- cbind(c(3.75, 3.58, 3.51, 3.46, 3.44, 3.43), c(3.33,

        3.22, 3.17, 3.14, 3.13, 3.12), c(3.00, 2.93, 2.89, 2.88, 
        2.87, 2.86), c(2.62, 2.60, 2.58, 2.57, 2.57, 2.57), c(0.37, 
        0.40, 0.42, 0.42, 0.43, 0.44), c(0.00, 0.03, 0.05, 0.06, 
        0.07, 0.07), c(0.34, 0.29, 0.26, 0.24, 0.24, 0.23), c(0.72, 
        0.66, 0.63, 0.62, 0.61, 0.60))

    if (k > 1)
    {

        yt1 <- z[,2:k]
        res <- lm(yt ~ xt1 + 1 + yt1)

    }
    else res <- lm(yt ~ xt1 + 1)
    res.sum <- summary(res)
    STAT <- res.sum$coefficients[2,1]/res.sum$coefficients[2,2]

    }

    if (kind==3)
    {
    table <- cbind(c(4.38, 4.15, 4.04, 3.99, 3.98, 3.96), c(3.95,

        3.8, 3.73, 3.69, 3.68, 3.66), c(3.6, 3.5, 3.45, 3.43, 
        3.42, 3.41), c(3.24, 3.18, 3.15, 3.13, 3.13, 3.12), c(1.14, 
        1.19, 1.22, 1.23, 1.24, 1.25), c(0.8, 0.87, 0.9, 0.92, 
        0.93, 0.94), c(0.5, 0.58, 0.62, 0.64, 0.65, 0.66), c(0.15, 
        0.24, 0.28, 0.31, 0.32, 0.33))

    if (k > 1)
    {

        yt1 <- z[,2:k]
        res <- lm(yt ~ xt1 + 1 + tt + yt1)
    }
    else res <- lm(yt ~ xt1 + 1 + tt)
    res.sum <- summary(res)
    STAT <- res.sum$coefficients[2,1]/res.sum$coefficients[2,2]

    }     

    table <- -table

    tablen <- dim(table)[2]
    tableT <- c(25, 50, 100, 250, 500, 1e+05)
    tablep <- c(0.01, 0.025, 0.05, 0.1, 0.9, 0.95, 0.975, 0.99)
    tableipl <- numeric(tablen)

    for (i in (1:tablen)) tableipl[i] <- approx(tableT, table[,i], n, rule = 2)$y

    interpol <- approx(tableipl, tablep, STAT, rule = 2)$y

    if (is.na(approx(tableipl, tablep, STAT, rule = 1)$y))

        if (interpol == min(tablep)) 
            warning("p-value smaller than printed p-value")
        else warning("p-value greater than printed p-value")

    PVAL <- interpol

    PARAMETER <- k - 1
    METHOD <- "Augmented Dickey-Fuller Test"     names(STAT) <- "Dickey-Fuller"
    names(PARAMETER) <- "Lag order"

    structure(list(statistic = STAT, parameter = PARAMETER,alternative = "The series is stationary",

    p.value = PVAL, method = METHOD, data.name = DNAME),class = "htest")

}

Amir Safari wrote:
>
>
>
> Hi All ,
> Could you please tell using which library ,Dickey-Fuller Test can be run?
> Thanks a lot
>
> __________________________________________________
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Mon May 23 20:53:43 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:32:00 EST