Machado <- function(y, max.lag=5, tau=seq(.1, .9, by=.005), trim.prop=.025)
{
################################################################################
#                                                                              #
#         Machado Test Function by Breno de Andrade Pinheiro Néri              #
#                                                                              #
#            brenoneri@gmail.com        www.fgv.br/aluno/bneri                 #
#                                                                              #
#                                                    October, 2005 - Version 1 #
#                                                                              #
################################################################################
#                                                                              #
# This code computes the test statistics for the null hypothesis of the first  #
# q lags are jointly non statiscally significant for every quantil, following  #
# the procedure described in Koenker, Roger and José A.F. Machado, 1999,       #
# "Goodness of Fit and Related Inference Processes for Quantile Regression",   #
# Journal of American Statistical Association 94 (448), 1296-1310.             #
#                                                                              #
# You enter a time series (y), the maximum number of lags to be analyzed       #
# (max.lag), a vector of quantiles (tau) and a trim value (trim.prop).         #
#                                                                              #
# Critical Value at  5% Significance Level: 9.31                               #
# Critical Value at 10% Significance Level: 7.36                               #
#                                                                              #
# These critical values are tabulated in Andrews, D.W.K., 1993, "Tests for     #
# Parameter Instability and Structural Change With Unknown Change Point",      #
# Econometrica 61, 821-856.                                                    #
#                                                                              #
# Defaults: max.lag=5, tau=seq(.1, .9, by=.005) and trim.prop=.025             #
#                                                                              #
################################################################################
#                                                                              #
# You may use this code only if you accept the conditions:                     #
# (1) I am not liable for any problem caused by this code (i.e. you use it at  #
# your own risk);                                                              #
# (2) You must give me credit in your papers where this code has been used.    #
#                                                                              #
################################################################################

y <- ts(y)
library(quantreg)
qrq1 <-function(s1, a)
{
    r <- s1$sol[1,]
    q <- s1$sol[2,]
    q <- c(q[1], q)
    J <- length(r)
    r <- c(0, (r[-J]+r[2:J])/2, 1)
    u <- 0
    for(k in 1:length(a))
    {
        w <- (a[k]-r[i<-sum(r<a[k])])/(r[i+1]-r[i])
        u[k] <- w*q[i+1]+(1-w)*q[i]
    }
    return(u)
}
for(q in max.lag:1)
{
    Y <- y
    for(i in 1:q) Y <- ts.intersect(Y,lag(y,-i))
    fit.QAR.1 <- rq(Y[,1]~Y[,2:(q+1)], tau=-1)
    h <- .810576*((dnorm(x0<-qnorm(tau))^2/(2*x0^2+1))^2/nrow(Y))^0.2
    trim <- tau-h<trim.prop|tau+h>1-trim.prop
    f.F.hat <- 2*h[!trim]/(qrq1(fit.QAR.1, (tau+h)[!trim])-qrq1(fit.QAR.1, (tau-h)[!trim]))
    L <- 0
    for(quantil in 1:length(tau))
    {
        u <- residuals(rq(Y[,1]~Y[,2:(q+1)], tau=tau[quantil]))
        V.hat <- sum(u*(tau[quantil]-(u<0)))
        u <- residuals(rq(Y[,1]~Y[,2:q], tau=tau[quantil]))
        V.tild <- sum(u*(tau[quantil]-(u<0)))
        L[quantil] <- 2*(V.tild-V.hat)*f.F.hat[quantil]/(tau[quantil]*(1-tau[quantil]))
    }
    message("For q = ", q, ", Sup(L) = ", max(abs(L)))
}
}