R/midas_sp.R

In midasr: Mixed Data Sampling Regression

##' Semi-parametric MIDAS regression
##'
##' Estimate semi-parametric MIDAS regression using non-linear least squares.
##'
##' @param formula formula for restricted MIDAS regression or \code{midas_r} object. Formula must include \code{\link{fmls}} function
##' @param data a named list containing data with mixed frequencies
##' @param bws a bandwith specification. Note you need to supply logarithm value of the bandwith. 
##' @param start the starting values for optimisation. Must be a list with named elements.
##' @param degree the degree of local polynomial. 0 corresponds to local-constant, 1 local-linear. For univariate models higher values can be provided.
##' @param Ofunction the list with information which R function to use for optimisation. The list must have element named \code{Ofunction} which contains character string of chosen 
##' R function. Other elements of the list are the arguments passed to this function.  The default optimisation function is \code{\link{optim}} with arguments
##'  \code{method="Nelder-Mead"} and \code{control=list(maxit=5000)}. Other supported functions are \code{\link{nls}}, \code{\link{optimx}}.
##' 
##' @param ... additional arguments supplied to optimisation function
##' @return a \code{midas_sp} object which is the list with the following elements:
##' 
##' \item{coefficients}{the estimates of parameters of restrictions}
##' \item{midas_coefficients}{the estimates of MIDAS coefficients of MIDAS regression}
##' \item{model}{model data}
##' \item{unrestricted}{unrestricted regression estimated using \code{\link{midas_u}}}
##' \item{term_info}{the named list. Each element is a list with the information about the term, such as its frequency, function for weights, gradient function of weights, etc.}
##' \item{fn0}{optimisation function for non-linear least squares problem solved in restricted MIDAS regression}
##' \item{rhs}{the function which evaluates the right-hand side of the MIDAS regression}
##' \item{gen_midas_coef}{the function which generates the MIDAS coefficients of MIDAS regression}
##' \item{opt}{the output of optimisation procedure}
##' \item{argmap_opt}{the list containing the name of optimisation function together with arguments for optimisation function}
##' \item{start_opt}{the starting values used in optimisation}
##' \item{start_list}{the starting values as a list}
##' \item{call}{the call to the function}
##' \item{terms}{terms object}
##' \item{gradient}{gradient of NLS objective function}
##' \item{hessian}{hessian of NLS objective function}
##' \item{gradD}{gradient function of MIDAS weight functions} 
##' \item{Zenv}{the environment in which data is placed}
##' \item{nobs}{the number of effective observations}
##' \item{convergence}{the convergence message}
##' \item{fitted.values}{the fitted values of MIDAS regression}
##' \item{residuals}{the residuals of MIDAS regression}
##' 
##' @author Virmantas Kvedaras, Vaidotas Zemlys-Balevičius
##' @rdname midas_sp
##' @details Given MIDAS regression:
##'
##' \deqn{y_t = \sum_{j=1}^p\alpha_jy_{t-j} +\sum_{i=0}^{k}\sum_{j=0}^{l_i}\beta_{j}^{(i)}x_{tm_i-j}^{(i)} + u_t,}
##' 
##' estimate the parameters of the restriction
##'
##' \deqn{\beta_j^{(i)}=g^{(i)}(j,\lambda).}
##'
##' Such model is a generalisation of so called ADL-MIDAS regression. It is not required that all the coefficients should be restricted, i.e the function \eqn{g^{(i)}}
##' might be an identity function. The regressors \eqn{x_\tau^{(i)}} must be of higher
##' (or of the same) frequency as the dependent variable \eqn{y_t}. 
##'
##' @importFrom stats as.formula formula model.matrix model.response terms lsfit time
##' @importFrom zoo index index2char
##' @importFrom Formula Formula
##' @export
midas_sp <- function(formula, data, bws, start, degree = 1,  Ofunction="optim", ...) {
    
    Zenv <- new.env(parent = environment(formula))
    
    if (missing(data)) {
        ee <- NULL
    }
    else {
        ee <- data_to_env(data)
        parent.env(ee) <- parent.frame()
    }
    
    if (missing(start)) {
        stop("Please supply starting values.")
    } 
    
    assign("ee",ee,Zenv)
    formula <- Formula(formula)
    environment(formula) <- Zenv
    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")
    args <- list(...)
    y <- as.numeric(model.response(mf, "numeric"))
    X <- model.matrix(formula, data = mf, rhs = 1)
    if(length(attr(formula,"rhs")) > 1) {
        Z <- model.matrix(formula, data = mf, rhs = 2)
    } else Z <- NULL
    
    #Save ts/zoo information
    if(is.null(ee)) { 
        yy <- eval(formula[[2]], Zenv)
    }else {
        yy <- eval(formula[[2]], ee)
    }
    
    y_index <- 1:length(yy) 
    if(!is.null(attr(mf, "na.action"))) {
        y_index <- y_index[-attr(mf, "na.action")]
    }
    if(length(y_index)>1) {
        if(sum(abs(diff(y_index) - 1))>0) warning("There are NAs in the middle of the time series")                
    }
    
    ysave <- yy[y_index]
    
    if(inherits(yy, "ts")) {
        class(ysave) <- class(yy)
        attr(ysave, "tsp") <- c(time(yy)[range(y_index)], frequency(yy))
    }
    
    if(inherits(yy,c("zoo","ts"))) {
        y_start <- index2char(index(ysave)[1], frequency(ysave))
        y_end <- index2char(index(ysave)[length(ysave)], frequency(ysave))
    } else {
        y_start <- y_index[1]
        y_end <- y_index[length(y_index)]
    }
    
    prepmd <- prep_midas_sp(y, X, Z, bws, degree, formula, Zenv,cl,args,start,Ofunction)
    
    prepmd <- c(prepmd, list(lhs = ysave, lhs_start = y_start, lhs_end = y_end))
    
    class(prepmd) <- "midas_sp"
    
    midas_nlpr.fit(prepmd)    
}

## Prepare necessary objects for fitting of the non-linear parametric MIDAS regression
##
## y the response
## X the model matrix
## mt the terms of the formula
## Zenv the environment to evaluate the formula
## cl call of the function
## args additional argument
## start starting values
## Ofunction the optimisation function
## weight_gradients a list of gradient functions for weights
## lagsTable the lagstable from checkARstar
## unrestricted the unrestricted model
## guess_start if TRUE, get the initial values for non-MIDAS terms via OLS, if FALSE, initialize them with zero.
prep_midas_sp <- function(y, X, Z, bws, degree, f, Zenv, cl, args, start, Ofunction,  guess_start = TRUE) {
    
    
    start <- start[!sapply(start,is.null)]
    
    if (!is.null(args$guess_start)) {
        guess_start <- args$guess_start
        args$guess_start <- NULL
    }  
    if(is.null(Z)) {
        ##We are having pure SI model.
        Z <- X
        X <- 0
        mt1 <- NULL
        mt2 <- terms(f, rhs = 1)
    } else {
        mt1 <- terms(f, rhs = 1)
        mt2 <- terms(f, rhs = 2)
        
        if (attr(mt1,"intercept") == 1)  {
            X <- X[, -1, drop = FALSE]
        }
        
    }

    if (attr(mt2,"intercept") == 1)  {
        Z <- Z[, -1, drop = FALSE]
    }
    
    bws_length <- length(bws)
    
    if(!is.null(mt1)) {
        terms.rhs1 <- as.list(attr(mt1,"variables"))[-2:-1]
        rfd1 <- lapply(terms.rhs1, dterm_nlpr,  Zenv = Zenv)
        term_names1 <- sapply(rfd1,"[[","term_name")
        names(rfd1) <- term_names1
        rf1 <- lapply(rfd1,"[[","weight")
   
        start_default1 <- lapply(rfd1,"[[","start")
        start_default1[intersect(names(start), term_names1)] <- start[intersect(names(start), term_names1)]
        np1 <- cumsum(sapply(start_default1,length))
        pinds1 <- build_indices(np1,names(start_default1))
        fake_x_coef <- all_coef_full(unlist(start_default1), pinds1, rf1) 
        npx <- cumsum(sapply(fake_x_coef,length))
        xinds <- build_indices(npx,names(start_default1))
        pinds1 <- lapply(pinds1, function(x) x + bws_length)
        cf1 <- function(p) unlist(all_coef_full(p, pinds1, rf1))
        p2_offset <- max(pinds1[[length(pinds1)]])
    } else {
        cf1 <- function(p) return(0)
        p2_offset <- bws_length
        start_default1 <- NULL
    }
    
    terms.rhs2 <- as.list(attr(mt2,"variables"))[-2:-1]
    rfd2 <- lapply(terms.rhs2, dterm_nlpr,  Zenv = Zenv)
    term_names2 <- sapply(rfd2,"[[","term_name")
    names(rfd2) <- term_names2
    rf2 <- lapply(rfd2,"[[","weight")
    weight_names2 <- sapply(rfd2,"[[","weight_name")
    weight_inds2 <- which(weight_names2 != "")
    weight_names2 <- names(rf2)[weight_names2 != ""]
    
    start_fake2 <- lapply(rfd2,"[[","start")
    start_fake2[intersect(names(start), weight_names2)] <- start[intersect(names(start), weight_names2)]
    np_fake2 <- cumsum(sapply(start_fake2,length))
    pinds_fake2 <- build_indices(np_fake2,names(start_fake2)) 
    fake_z_coef <- all_coef_full(unlist(start_fake2), pinds_fake2, rf2) 
    npz <- cumsum(sapply(fake_z_coef,length))
    zinds <- build_indices(npz,names(start_fake2))
   
    start_default2 <- start_fake2[weight_names2]
    if (length(start_default2) > 0) {
        np2 <- cumsum(sapply(start_default2,length))
        pinds2 <- build_indices(np2,names(start_default2))
        pinds2 <- lapply(pinds2, function(x) x + p2_offset)
    } else start_default2 <- NULL
    
    cf0 <- function(p) p[1:bws_length]
   
    rhs_cv <- function(p) {
        h <- cf0(p)
        xi <- as.numeric(X %*% cf1(p))
        zi <- do.call("rbind", lapply(term_names2, function(t2) {
            if (t2 %in% weight_names2) Z[, zinds[[t2]]] %*% rf2[[t2]](p[pinds2[[t2]]])
            else Z[, zinds[[t2]], drop = FALSE]
        }))
        if (ncol(zi) == 1) zi <- as.numeric(zi)
        u <- y - xi
        xi + cv_np(u, zi, h, degree)
    }
    
    fn0 <- function(p,...) {
        if(any(cf0(p)<= 0)) return(NA) 
        else {
            r <- y - rhs_cv(p)
            return(sum(r^2))
        }
    }
    
    hess <- function(x)numDeriv::hessian(fn0,x)
    
    control <- c(list(Ofunction = Ofunction),args)
    ##The default method is "Nelder-Mead" and number of maximum iterations is 5000
    if (!("method" %in% names(control)) & Ofunction == "optim") {        
        control$method <- "Nelder-Mead"
        if (is.null(control$control$maxit)) control$control <- list(maxit=5000)
    }
    if (bws_length > 1) names(bws) <- paste0("bw", 1:bws_length)
    else names(bws) <- "bw"
    
    if (is.null(mt1)) {
        starto <- c(bws, unlist(start_default2)) 
        term_info1 <- NULL   
        model_matrix_map = list(y = 1,
                                X = NULL,
                                Z = 1 + 1:ncol(Z))
        model <- cbind(y, Z)
    } else {
        starto <- c(bws, unlist(start_default1), unlist(start_default2)) 
        term_info1 <- mapply(function(l, pind, xind) {
            l$coef_index <- pind
            l$midas_coef_index <- xind
            l$term_type <- "X"
            l
        },rfd1, pinds1, xinds, SIMPLIFY = FALSE)
        model_matrix_map = list(y = 1,
                                X = 1 + 1:ncol(X),
                                Z = ncol(X) + 1 + 1:ncol(Z))
        model <- cbind(y, X, Z)
    }
    
    term_info2 <- lapply(term_names2, function(t2) {
        res <- rfd2[[t2]]
        if(t2 %in% weight_names2) {
            res$coef_index <- pinds2[[t2]]
            res$midas_coef_index <- zinds[[t2]]
        } else {
            res$midas_coef_index <- zinds[[t2]]
        }
        res$term_type <- "Z"
        res
    })
    
    names(term_info2) <- term_names2
    term_info <- c(term_info1, term_info2)
    
    
    list(coefficients = starto,
         model = model, 
         fn0 = fn0,
         rhs = rhs_cv,
         coef_X  = cf1,
         model_matrix_map = model_matrix_map,
         opt = NULL,
         argmap_opt = control,
         start_opt = starto,
         start_list = c(list(bws = bws), start1 = start_default1, start2 = start_default2),
         call = cl,
         terms = mt1,
         terms2 = mt2,
         formula = f,
         bws = bws,
         term_info = term_info,
         hessian = hess,
         Zenv = Zenv,
         nobs = nrow(Z),
         degree = degree)   
}


##' MIDAS Single index regression
##'
##' Function for fitting SI MIDAS regression without the formula interface
##' @param y model response
##' @param X prepared matrix of high frequency variable lags for MMM term
##' @param p.ar length of AR part
##' @param weight the weight function
##' @param degree the degree of local polynomial
##' @param start_bws the starting values bandwith
##' @param start_x the starting values for weight function
##' @param start_ar the starting values for AR part. Should be the same length as \code{p}
##' @param method a method passed to \link{optim}, defaults to Nelder-Mead
##' @param ... additional parameters to \link{optim}
##' @return an object similar to \code{midas_r} object
##' @author Virmantas Kvedaras, Vaidotas Zemlys
##' 
##' @importFrom stats na.omit optim
##'
##' @export
##' 
midas_si_plain <- function(y, X, p.ar = NULL, weight, degree = 1, start_bws, start_x, start_ar = NULL, method = "Nelder-Mead", ...) {
    d <- ncol(X)
    
    yy <- NULL 
    if(!is.null(p.ar)) {
        p.ar <- as.integer(p.ar)
        if (length(start_ar)!=p.ar) stop("The length of starting values vector for AR terms must the same as the number of AR terms")
        yy <- as.matrix(mls(y, 1:p.ar, 1))
    } 
    
    model <- na.omit(cbind(y,X,yy))
    
    y <- model[,1]
    X <- model[, 2:(ncol(X) + 1)]
    
    if (is.null(yy)) { 
        yy <- 0
    } else {
        yy <- model[, (ncol(X) + 2):ncol(model)]
    }
    n <- nrow(model)
    
    sx <- length(start_x)
    
    rhs_cv <- function(p) {
        h <- p[1]
        pr <- p[2:(1 + sx)]
        if (is.null(yy)) pyy <- 0 else {
            pyy <- p[(2 + sx):length(p)]
        }
        yar <- yy %*% pyy
        u <- y - yar
        yar + cv_np(u, as.vector(X %*% weight(pr, ncol(X))), h, degree)
    }
    
    fn0 <- function(p) {
        sum((y - rhs_cv(p))^2)
    }
    
    start <- c(start_bws, start_x, start_ar)
    opt <- optim(start, fn0, method = method,...)
    par <- opt$par
    call <- match.call()
    
    rhs <- function(p) {
        h <- p[1]
        pr <- p[2:(1 + sx)]
        if (is.null(yy)) pyy <- 0 else {
            pyy <- p[(2 + sx):length(p)]
        }
        yar <- yy %*% pyy
        u <- y - yar
        xi <-  as.vector(X %*% weight(pr, ncol(X)))
        np <- g_np(u, xi,  xeval = xi, h, degree)
        list(fitted = yar + np, xi = yar, z = xi, y = y, g = np)
    }
    fitted.values <- as.vector(rhs(par)$fitted)
    names(par) <- c("h", paste0("x", 1:length(start_x)), paste0("y", 1:length(start_ar)))

    list(coefficients = par,
         midas_coefficients = weight(par[2:(1 + sx)],ncol(X)),
         bws = par[1],
         model = model,
         weights = weight,
         fn0 = fn0,
         rhs_cv = rhs_cv,
         rhs = rhs,
         opt = opt,
         call = call,
         hessian = function(x)numDeriv::hessian(fn0,x),
         fitted.values = fitted.values,
         residuals = as.vector(y - fitted.values),
         start = start)
    
}

##' MIDAS Partialy linear non-parametric regression
##'
##' Function for fitting PL MIDAS regression without the formula interface
##' @param y model response
##' @param X prepared matrix of high frequency variable lags for MMM term
##' @param z a vector, data for the non-parametric part
##' @param p.ar length of AR part
##' @param weight the weight function
##' @param degree the degree of local polynomial
##' @param start_bws the starting values bandwith
##' @param start_x the starting values for weight function
##' @param start_ar the starting values for AR part. Should be the same length as \code{p}
##' @param method a method passed to \link{optim}
##' @param ... additional parameters to \link{optim}
##' @return an object similar to \code{midas_r} object
##' @author Virmantas Kvedaras, Vaidotas Zemlys
##' 
##' @importFrom stats na.omit 
##'
##' @export
##' 
midas_pl_plain <- function(y, X, z, p.ar = NULL, weight, degree = 1, start_bws, start_x, start_ar = NULL, method = c("Nelder-Mead"), ...) {
    d <- ncol(X)
    
    yy <- NULL 
    
    if(!is.null(p.ar)) {
        p.ar <- as.integer(p.ar)
        if (length(start_ar)!=p.ar) stop("The length of starting values vector for AR terms must the same as the number of AR terms")
        yy <- as.matrix(mls(y, 1:p.ar, 1))
    } 
    z <- as.numeric(z)
    y <- as.numeric(y)
    if ((length(z) != length(y)) &  (length(y) != nrow(X))) stop("The dimensions of the data do not match, z and y must be vectors of equal length, which must coincide with number of rows of X")
    
    model <- na.omit(cbind(y, X, z, yy ))
    
    y <- model[,1]
    X <- model[, 2:(ncol(X) + 1)]
    z <- model[, ncol(X) + 2]
   
    if (is.null(yy)) { 
        yy <- 0
    } else {
        yy <- model[, (ncol(X) + 3):ncol(model)]
    }
    n <- nrow(model)
    
    sx <- length(start_x)
    
    rhs_cv <- function(p) {
        h <- p[1]
        pr <- p[2:(1 + sx)]
        if (is.null(yy)) pyy <- 0 else {
            pyy <- p[(2 + sx):length(p)]
        }
        yar <- yy %*% pyy
        xi <- as.vector(X %*% weight(pr, ncol(X)))
        u <- y - yar - xi
        yar + xi + cv_np(u, z, h, degree)
    }
    
    fn0 <- function(p) {
        sum((y - rhs_cv(p))^2)
    }
    start <- c(start_bws, start_x, start_ar)
    opt <- optim(start, fn0, method = method, ...)
    par <- opt$par
    call <- match.call()
    
    rhs <- function(p) {
        h <- p[1]
        pr <- p[2:(1 + sx)]
        if (is.null(yy)) pyy <- 0 else {
            pyy <- p[(2 + sx):length(p)]
        }
        yar <- yy %*% pyy
        xi <-  as.vector(X %*% weight(pr, ncol(X)))
        u <- y - yar  - xi
        
        np <- g_np(u, z, xeval = z, h, degree)
        list(fitted = yar + xi + np, xi = xi + yar, z = z, g = np, y = y)
    }
    fitted.values <- as.vector(rhs(par)$fitted)
    names(par) <- c("h", paste0("x", 1:length(start_x)), paste0("y", 1:length(start_ar)))
    
    list(coefficients = par,
         midas_coefficients = weight(par[2:(1 + sx)],ncol(X)),
         bws = par[1],
         model = model,
         weights = weight,
         fn0 = fn0,
         rhs_cv = rhs_cv,
         rhs = rhs,
         opt = opt,
         call = call,
         hessian = function(x)numDeriv::hessian(fn0,x),
         fitted.values = fitted.values,
         residuals = as.vector(y - fitted.values),
         start = start)
    
}

cv_np <- function(y, x, h, degree = 1) {
    cvg <- rep(NA, length(y))
    for (i in 1:length(y)) {
        if (is.null(dim(x))) w <- kfun(x[i], x, h)
        else  w <- kfun(x[i, ], x, h)
        if (degree > 0) {
            if (is.null(dim(x))) {
                xc <-  matrix(x - x[i], ncol = 1)
            } else xc <- sweep(x, 2, x[i, ], "-")
            if(degree > 1) xc <- poly(xc, degree = degree, raw = TRUE, simple = TRUE)
            cc <- suppressWarnings(lsfit(xc[-i,], y[-i], wt = w[-i], intercept = TRUE))
            cvg[i] <- coef(cc)[1]
        } else {
            cvg[i] <- sum(w[-i]*y[-i])/sum(w[-i])
        }
    }
    cvg
}

g_np <- function(y, x, xeval, h, degree = 1) {
    res <- rep(NA, length(xeval))
    for (i in 1:length(xeval)) {
        w <- kfun(xeval[i], x, h)
        if(degree > 0) {
            xc <- poly(x - xeval[i], degree = degree, raw = TRUE, simple = TRUE)
            cc <- lsfit(xc, y, wt = w, intercept = TRUE)
            res[i] <- coef(cc)[1]
        } else {
            res[i] <- sum(w*y)/sum(w)
        }
    }
    res
}

g_np_mv <- function(u, z, zeval, h, degree = 1) {
    res <- rep(NA, nrow(zeval))
    for (i in 1:nrow(zeval)) {
        w <- kfun(zeval[i, ], z, h)
        if(degree > 0) {
            zc <- sweep(z, 2, zeval[i, ], "-")
            if(degree > 1) zc <- poly(zc, degree = degree, raw = TRUE, simple = TRUE)
            cc <- lsfit(zc, u, wt = w, intercept = TRUE)
            res[i] <- coef(cc)[1]
        } else {
            res[i] <- sum(w*u)/sum(w)
        }
    }
    res
}

midas_sp_fit <- function(object, Xeval = NULL, Zeval = NULL) {
    
    p <- coef(object)
    
    y <- object$model[, 1]
    h <- p[1:length(object$bws)]
    Z <- object$model[, object$model_matrix_map$Z, drop = FALSE]

    zterms <- sapply(object$term_info, "[[", "term_type") == "Z"
    
    ti <- object$term_info[zterms]
    
    calc_zi <- function(term, ZZ) {
        if (is.null(term$coef_index)) zi <- ZZ[, term$midas_coef_index, drop = FALSE]
        else zi <- ZZ[, term$midas_coef_index] %*% term$weight(p[term$coef_index])
    }
    
    zi <- do.call("rbind", lapply(ti, calc_zi, ZZ = Z))
    if (is.null(Zeval)) zi_eval <- zi
    else zi_eval <- do.call("rbind", lapply(ti, calc_zi, ZZ = Zeval))
   
    if (is.null(object$model_matrix_map$X)) {
        xi <- 0  
    } else {
        xi <- object$model[, object$model_matrix_map$X] %*% object$coef_X(p)
    } 
    
    u <- y - xi
    ghat <- g_np_mv(u, zi, zi_eval, h, degree = object$degree)
    
    if(is.null(Xeval)) xi_eval <- xi
    else xi_eval <- Xeval %*% object$coef_X(p)
    
    fit <- xi_eval + ghat
    list(fitted.values = fit, xi = xi_eval, zi = zi_eval, g = ghat, y = y)
}


#' @importFrom stats dnorm
kfun <- function(z0, z, h) {
    if (is.null(ncol(z))) {
        dnorm((z - z0)/h)
    } else {
        z0 <- as.numeric(z0)
        h <- as.numeric(h)
        apply(dnorm(scale(z, center = as.numeric(z0), scale = h)), 1, prod)
    }
}

all_coef_full <- function(p, pinds, rf) {              
    pp <- lapply(pinds,function(x)p[x])     
    res <- mapply(function(fun,param)fun(param),rf ,pp,SIMPLIFY=FALSE)
    return(res)
}