R/nonparametric.R

In midasr: Mixed Data Sampling Regression

##' Estimates non-parametric MIDAS regression
##'
##' Estimates non-parametric MIDAS regression accodring Breitung et al.
##' 
##' @title Estimate non-parametric MIDAS regression
##' @param formula formula specifying MIDAS regression 
##' @param data a named list containing data with mixed frequencies
##' @param lambda smoothing parameter, defaults to \code{NULL}, which means that it is chosen by minimising AIC.
##' @return a \code{midas_r_np} object
##' @author Vaidotas Zemlys
##' @references Breitung J, Roling C, Elengikal S (2013). \emph{Forecasting inflation rates using daily data: A nonparametric MIDAS approach} Working paper, URL http://www.ect.uni-bonn.de/mitarbeiter/joerg-breitung/npmidas.
##' @export
##' @import Matrix
##' @importFrom stats optimize
##' @examples
##' data("USunempr")
##' data("USrealgdp")
##' y <- diff(log(USrealgdp))
##' x <- window(diff(USunempr),start=1949)
##' trend <- 1:length(y)
##' midas_r_np(y~trend+fmls(x,12,12))
midas_r_np <- function(formula, data, lambda=NULL) {
    Zenv <- new.env(parent=environment(formula))
   
    if(missing(data)) {
        ee <- NULL
    }
    else {
        ee <- data_to_env(data)
        parent.env(ee) <- parent.frame()
    }
    
    assign("ee",ee,Zenv)
    formula <- as.formula(formula)
    cl <- match.call()    
    mf <- match.call(expand.dots = FALSE)
    mf$formula <- formula
    m <- match(c("formula", "data"), names(mf), 0L)
    mf <- mf[c(1L, m)]
    mf[[1L]] <- as.name("model.frame")
    mf[[3L]] <- as.name("ee")   
    mf[[4L]] <- as.name("na.omit")
    names(mf)[c(2,3,4)] <- c("formula","data","na.action")
    
    mf <- eval(mf,Zenv)
    mt <- attr(mf, "terms")

    terms.lhs <- as.list(attr(mt,"variables"))[-2:-1]
    term.labels <- attr(mt,"term.labels") 

    rfd <- vector("list",length(terms.lhs))
    yname <- all.vars(mt[[2]])

    for(i in 1:length(rfd)) {
        fr <- terms.lhs[[i]]
        #This a behaviour of R we rely on. It might be non-standard one.
        fun <- as.character(fr)[1]
        rfd[[i]] <- if(fun %in% c("fmls","mls","dmls")) {
            lags <- eval(fr[[3]],Zenv)
            nm <-as.character(fr[[2]])            
            nol <- switch(fun,
                              fmls = lags+1,
                              dmls = lags+1,
                              mls = length(lags)
                          )
            lagstruct <- switch(fun,
                                fmls = 0:lags,
                                dmls = 0:lags,
                                mls = lags
                          )
                                
            if(nol<3 & nm!=yname)stop("For nonparametric MIDAS you need at least 3 high frequency lags")          
           
            wlab <- ifelse(nm==yname,"",nm)
            list(weight=function(p)p,
                 term_name = nm,
                 gradient = NULL,
                 start = NULL,
                 weight_name = "non-parametric weight",
                 frequency = eval(fr[[4]],Zenv),
                 lag_structure = lagstruct
                )
        } else {
            list(weight=function(p)p,
                 term_name = term.labels[i],
                 gradient = NULL,
                 start = NULL,
                 weight_name = "",
                 frequency = 1,
                 lag_structure = 0
                 )
        }
    }
    
    if(attr(mt,"intercept")==1) {
        rfd <- c(list(list(weight = function(p)p,
                           term_name ="(Intercept)",
                           gradient = NULL,
                           start = NULL,
                           weight_name = "",
                           frequency = 1,
                           lag_structure = 0)),rfd)
    }


    names(rfd) <- sapply(rfd,"[[","term_name")

    ##Note this is a bit of misnomer. Variable weight_names is actualy a vector of term names which have MIDAS weights.
    ##It *is not* the same as actual name of weight function. This is a left-over from the old code. Grabbed from prepmidas_r 

    weight_names <- sapply(rfd,"[[","weight_name")
    weight_inds <- which(weight_names!="")
    weight_names <- names(rfd)[weight_names!=""]

    lengths <- sapply(lapply(rfd,"[[","lag_structure"),length)
    
    build_indices <- function(ci,nm) {
        inds <- cbind(c(1,ci[-length(ci)]+1),ci)
        inds <- apply(inds,1,function(x)list(x[1]:x[2]))
        inds <- lapply(inds,function(x)x[[1]])
        names(inds) <- nm
        inds
    }
    
    pinds <- build_indices(cumsum(lengths),names(rfd))
        
    if(length(weight_names)>1)stop("Only one non-autoregressive mixed frequency term is currently supported")
        
    resplace <- pinds[[weight_names]][1]     
    rno <- length(rfd[[weight_names]]$lag_structure)
       
    y <- model.response(mf, "numeric")
    X <- model.matrix(mt, mf)    

    k <- ncol(X)
    if(k < nrow(X)) {
        unrestricted <- lm(y~.-1,data=data.frame(cbind(y,X),check.names=FALSE))
    }else unrestricted <-  NULL
    
    D <- bandSparse(rno-2,k,resplace-1+c(0,1,2),diagonals=list(rep(1,rno-2),rep(-2,rno-2),rep(1,rno-2)))
  
    DD <- crossprod(D)
    ol <- opt_lambda(y,X,DD,lambda)

    fit <- X%*%ol$beta
    res <- y-fit
    
    cf <- as.numeric(ol$beta)
    names(cf) <- names(unlist(pinds))

    term_info <- rfd
    names(term_info) <- sapply(term_info,"[[","term_name")
    term_info <- mapply(function(term,pind,xind){
        term$start <- NULL
        term$coef_index <- pind
        term$midas_coef_index <- pind
        term
    },term_info,pinds[names(term_info)],SIMPLIFY=FALSE)

    
    out <- list(coefficients=cf,
                midas_coefficients=cf,
                model=cbind(y,X),
                unrestricted = unrestricted,
                call=cl,        
                terms=mt,
                fitted.values=as.numeric(fit),
                residuals=as.numeric(res),
                term_info = term_info,
                lambda=ol$lambda,
                klambda=ol$klambda,
                AIC=ol$AIC,
                opt=ol$opt,
                Zenv=Zenv
                )    
    class(out) <- "midas_r_np"
    out
}

##' @importFrom stats AIC
##' @export
##' @method AIC midas_r_np
AIC.midas_r_np <- function(object, ..., k) {
    object$AIC(object$lambda)
}

##' @importFrom stats BIC
##' @export
##' @method BIC midas_r_np
BIC.midas_r_np <- function(object, ..., k) {
    #for compatibility
    object$AIC(object$lambda)
}

##' @export
##' @method forecast midas_r_np
forecast.midas_r_np <- forecast.midas_r

##' @export
##' @method predict midas_r_np
predict.midas_r_np <- predict.midas_r

##' @export
##' @method print midas_r_np
print.midas_r_np <- function(x,...) {
    cat("Nonparametric MIDAS regression model",paste0("(",nrow(x$model)), "low frequency observations)")
    cat("\nFormula: ", deparse(formula(terms(x))))    
    cat("\nThe smoothing parameter: ", x$lambda)
    cat("\nThe effective number of parameters:", x$klambda)
    cat("\nAIC of the model: ",AIC(x))
    cat("\nRoot mean squared error: ", sqrt(mean(residuals(x)^2)),"\n")
}

##' @export
##' @method summary midas_r_np
summary.midas_r_np <- function(object,...) {
  print(object,...)
}

opt_lambda <- function(y,X,DD,lambda) {
    n <- length(y)
    XX <- crossprod(X)
    Xy <- crossprod(X,y)
    tX <- t(X)
    AIC <- function(lambda) {
            Qlambda <- XX+lambda*n*DD        
            klambda <- sum(X*t(solve(Qlambda,tX)))
            beta <- solve(Qlambda,Xy)
            res <- y-X%*%beta
            log(sum(res^2))+2*(klambda+1)/(n-klambda-2)
        }
    if(is.null(lambda)) {       
        opt <- optimize(AIC,c(0,100))     
        lambda <- opt$minimum
    }
    else opt <- NULL
    Qlambda <- XX+lambda*n*DD
    klambda <- sum(X*t(solve(Qlambda,tX)))
    beta <- solve(Qlambda,Xy)    
    AICl <- AIC(lambda)
    list(beta=beta,klambda=klambda,lambda=lambda,AIC=AIC,opt=opt)
    
}