R/midas_sp_methods.R

In midasr: Mixed Data Sampling Regression

##' Fitted values for semi-parametric MIDAS regression model
##'
##' Returns the fitted values of a fitted semi-parametric MIDAS regression object
##' @param object a \code{\link{midas_r}} object
##' @param ... currently nothing
##' @return the vector of fitted values
##' @author Virmantas Kvedaras, Vaidotas Zemlys
##' @export
##' @method fitted midas_sp
fitted.midas_sp <- function(object, ...) {
    fit <- midas_sp_fit(object, Xeval = NULL, Zeval = NULL)
    as.numeric(fit$fitted.values)
}

##' @export
##' @method print midas_sp
print.midas_sp <- function(x, digits=max(3,getOption("digits")-3),...) {
    #Code adapted from dynln:::print.dynlm code
    model_string <-  "\nSemi-parametric MIDAS regression model with \""
    cat(paste(model_string, class(x$lhs)[1], 
              "\" data:\n", sep = ""))
    cat(paste("Start = ", x$lhs_start, 
              ", End = ", x$lhs_end, 
              "\n", sep = ""))
    cat(" model:", deparse(x$call$formula),"\n")
    print(coef(x),digits = digits, ...)
    cat("\n")
    cat("Function", x$argmap_opt$Ofunction, "was used for fitting\n")
    invisible(x)
}


##' @export
##' @importFrom stats deviance pt pnorm residuals printCoefmat cor
##' @method summary midas_sp
summary.midas_sp <- function(object, df=NULL, ...) {
    r <- as.vector(residuals(object))
    param <- coef(object)
    pnames <- names(param)
    n <- length(r)
    p <- length(param)
    rdf <- n - p
    
    resvar <- deviance(object)/rdf
    
    dg <- jacobian(object$rhs, param)
    V <- resvar * ginv(crossprod(dg)/nrow(dg))/n
    
    colnames(V) <- pnames
    rownames(V) <- pnames
   
    se <- sqrt(diag(V)) 
       
    f <- as.vector(object$fitted.values)
    mss <- sum((f - mean(f))^2)
    rss <- sum(r^2)
    
    n <- length(r)
    p <- length(coef(object))
    rdf <- n-p
    df.int <- 0L
    
    r_squared <- R2.np(object$model[,1], f) 
    tval <- param/se
    
    #Code stolen from coeftest function
    if(is.null(df)) {
        df <- rdf
    }
    if (is.finite(df) && df > 0) {
        pval <- 2 * pt(abs(tval), df = df, lower.tail = FALSE)
        cnames <- c("Estimate", "Std. Error", "t value", "Pr(>|t|)")
        mthd <- "t"
    }
    else {
        pval <- 2 * pnorm(abs(tval), lower.tail = FALSE)
        cnames <- c("Estimate", "Std. Error", "z value", "Pr(>|z|)")
        mthd <- "z"
    }
    
    param <- cbind(param,se,tval,pval)
    dimnames(param) <- list(pnames, c("Estimate", "Std. Error", 
                                      "t value", "Pr(>|t|)"))
    ans <- list(formula = formula(object$call$formula), residuals=r, sigma=sqrt(resvar),
                df=c(p,rdf), cov.unscaled = V/resvar, call=object$call,
                coefficients=param,
                r_squared = r_squared, lhs_start = object$lhs_start, lhs_end = object$lhs_end, class_lhs = class(object$lhs))
    class(ans) <- "summary.midas_sp"
    ans
}

##' Semi-parametric MIDAS regression model deviance
##'
##' Returns the deviance of a fitted MIDAS regression object
##' @param object a \code{\link{midas_r}} object
##' @param ... currently nothing
##' @return The sum of squared residuals
##' @author Virmantas Kvedaras, Vaidotas Zemlys
##' @rdname deviance.midas_sp
##' @method deviance midas_sp
##' @export
deviance.midas_sp <- function(object,...) {
    sum(residuals(object)^2,na.rm=TRUE)
}

##' @export
##' @method print summary.midas_sp
print.summary.midas_sp <- function(x, digits=max(3, getOption("digits") - 3 ), signif.stars = getOption("show.signif.stars"), ...) {
    cat(paste("\nSemi-parametric MIDAS regression model with \"", x$class_lhs[1], 
              "\" data:\n", sep = ""))
    cat(paste("Start = ", x$lhs_start, 
              ", End = ", x$lhs_end, 
              "\n", sep = ""))
    cat("\n Formula:", deparse(x$formula),"\n")
    df <- x$df
    rdf <- df[2L]
    cat("\n Parameters:\n")
    printCoefmat(x$coefficients,digits=digits,signif.stars=signif.stars,...)
    cat("\n Residual standard error:", format(signif(x$sigma,digits)), "on", rdf , "degrees of freedom\n")
    #   cat(" Multiple R-squared: ", formatC(x$r_squared, digits = digits))
    #   cat(",\tAdjusted R-squared: ", formatC(x$adj_r_squared, digits = digits),"\n")
    invisible(x)
}


##' @include midas_r_methods.R
setMethod("extract", signature = className("midas_sp", "midasr"), definition = extract.midas_r)

##' @method update midas_sp
##' @importFrom stats getCall setNames
##' @export
update.midas_sp <- function(object, formula.,..., evaluate = TRUE) {
    if (is.null(call <- getCall(object))) 
        stop("need an object with call component")
    extras <- match.call(expand.dots = FALSE)$...
    if (!missing(formula.)) 
        call$formula <- update(Formula(object), formula.)
    
    if (length(extras)) {
        existing <- !is.na(match(names(extras), names(call))) 
        for (a in names(extras)[existing]) call[[a]] <- extras[[a]]
        if (any(!existing)) {
            call <- c(as.list(call), extras[!existing])
            call <- as.call(call)
        }        
    }
    
    ##1. If no start is supplied update the start from the call
    ##2. If start is supplied intersect it with already fitted values.
    
    cf <- coef(object)
    ustart <- lapply(object$term_info,function(x)cf[x[["coef_index"]]])
    
    redo <- FALSE
    if(!("start" %in% names(extras))) {        
        if(!("start" %in% names(call) && is.null(call$start))) {
            call$start <- ustart
            object$start_opt <- cf
        }
    } else {
        cstart <- eval(call$start,object$Zenv)
        ustart[names(cstart)] <- cstart
        call$start <- ustart
        object$start_opt <- unlist(ustart)
    } 
    
    if (evaluate) {
        if(!missing(formula.) || "data" %in% names(extras)  || "weight_gradients" %in% names(extras) || redo) {
            eval(call, parent.frame())
        } else {
            ##If we got here, we assume that we do not need to reevaluate terms.
            if(!is.null(extras$Ofunction)) {
                Ofunction <- eval(extras$Ofunction)
                extras$Ofunction <- NULL
            } else Ofunction <- object$argmap_opt$Ofunction            
            dotargnm <- names(extras)
            if (length(dotargnm) > 0) {
                offending <- dotargnm[!dotargnm %in% names(formals(Ofunction))]
                if (length(offending) > 0) {
                    stop(paste("The function ", Ofunction, " does not have the following arguments: ", 
                               paste(offending, collapse = ", "), sep = ""))
                }
            }
            else {
                extras <- NULL
            }
            if (Ofunction != object$argmap_opt$Ofunction) {
                argmap <- c(list(Ofunction = Ofunction), extras)
            }
            else {
                argmap <- object$argmap_opt
                argmap$Ofunction <- NULL
                argnm <- union(names(argmap), names(extras))
                marg <- vector("list", length(argnm))
                names(marg) <- argnm
                marg[names(extras)] <- extras
                oldarg <- setdiff(names(argmap), names(extras))
                marg[oldarg] <- argmap[oldarg]
                argmap <- c(list(Ofunction = Ofunction), marg)
            }
            object$call <- call
            object$argmap_opt <- argmap
            midas_nlpr.fit(object)
        }
    }
    else call
}

##' Predicted values based on \code{midas_sp} object.
##'
##' \code{predict.midas_sp} produces predicted values, obtained by evaluating regression function in the frame \code{newdata}. This means that the appropriate model matrix is constructed using only the data in \code{newdata}. This makes this function not very convenient for forecasting purposes. If you want to supply the new data for forecasting horizon only use the function \link{forecast.midas_r}. Also this function produces only static predictions, if you want dynamic forecasts use the \link{forecast.midas_r}.
##' 
##' @title Predict method for semi-parametric MIDAS regression fit
##' @param object \code{\link{midas_nlpr}} object
##' @param newdata a named list containing data for mixed frequencies. If omitted, the in-sample values are used.
##' @param na.action function determining what should be done with missing values in \code{newdata}. The most likely cause of missing values is the insufficient data for the lagged variables. The default is to omit such missing values.
##' @param ... additional arguments, not used
##' @return a vector of predicted values
##' @author Virmantas Kvedaras, Vaidotas Zemlys-Balevičius
##' @method predict midas_sp
##' @rdname predict.midas_sp
##' @export
predict.midas_sp <- function(object, newdata, na.action = na.omit, ... ) {
    Zenv <- new.env(parent=parent.frame())
    
    if(missing(newdata))
        return(as.vector(fitted(object)))
    else {
        ee <- data_to_env(newdata)    
        ZZ <- object$Zenv    
        parent.env(ee) <- ZZ
    }
    
    assign("ee",ee,Zenv)
    cll <- match.call()
    mf <- match.call(expand.dots = FALSE)
    m <- match(c("object", "newdata"), names(mf), 0L)
    mf <- mf[c(1L, m)]
    mf[[1L]] <- as.name("model.frame")
    
    formula <- Formula(formula(object))
    Terms <- delete.response(formula)
    mf[[2L]] <- Terms    
    mf[[3L]] <- as.name("ee")   
    mf[[4L]] <- na.action
    names(mf)[c(2,3,4)] <- c("formula","data","na.action")
    
    mf <- eval(mf,Zenv)
    mt <- attr(mf, "terms")
    
  

    Xeval <- model.matrix(formula, data = mf, rhs = 1)
    if (length(attr(formula, "rhs")) > 1) {
        Zeval <- model.matrix(formula, data = mf, rhs = 2) 
        if (attr(object$terms2,"intercept") == 1)  {
            Zeval <- Zeval[, -1, drop = FALSE]
        }
        if (attr(object$terms,"intercept") == 1)  {
            Xeval <- Xeval[, -1, drop = FALSE]
        }
    }
    else { 
        if (attr(object$terms2,"intercept") == 1)  {
            Zeval <- Xeval[, -1, drop = FALSE]
        }
        Xeval <- NULL
    }
    
    as.vector(midas_sp_fit(object, Xeval, Zeval)$fitted.values)
}

##' Extracts various coefficients of MIDAS regression
##'
##' MIDAS regression has two sets of cofficients. The first set is the coefficients associated with the parameters
##' of weight functions associated with MIDAS regression terms. These are the coefficients of the NLS problem associated with MIDAS regression.
##' The second is the coefficients of the linear model, i.e  the values of weight
##' functions of terms, or so called MIDAS coefficients. By default the function returns the first set of the coefficients.
##' 
##' @title Extract coefficients of MIDAS regression
##' @param object \code{midas_nlpr} object
##' @param type one of plain, midas, or nlpr. Returns appropriate coefficients.
##' @param term_names a character vector with term names. Default is \code{NULL}, which means that coefficients of all the terms are returned
##' @param ... not used currently
##' @return a vector with coefficients
##' @author Vaidotas Zemlys
##' @method coef midas_sp
##' @rdname coef.midas_sp
##' @export
coef.midas_sp <- function(object, type = c("plain", "midas", "bw"), term_names = NULL, ...) {
    type <- match.arg(type)
    if (is.null(term_names) & type == "midas") stop("Please provide a term name to get midas type coefficients")
    if (is.null(term_names)) {
        if(type == "bw") {
            object$coefficients[1:length(object$bws)]
        }
        else return(object$coefficients)
    } else {
        if (type == "bw") {
            warning("Bandwith setting is the same for all terms") 
            coef(object, type = "bw", term_names = NULL)
        }
        if (length(setdiff(term_names,names(object$term_info))) > 0) stop("Some of the term names are not present in estimated MIDAS regression")
        if (type == "plain") {
            res <- lapply(object$term_info[term_names], function(x) {
                if (!is.null(x$coef_index)) object$coefficients[x$coef_index]
            })
        }  
        if (type == "midas") {
            res <- lapply(object$term_info[term_names], function(x) {
                if (!is.null(x$coef_index))  x$weight(object$coefficients[x$coef_index])
            })
        }
        names(res) <- NULL
        if (length(res) == 1) return(res[[1]])
        else return(unlist(res))
        
    }
}

##' @include midas_r_methods.R
##' @include midas_nlpr_methods.R
##' @method plot_midas_coef midas_sp
##' @export
plot_midas_coef.midas_sp <- plot_midas_coef.midas_nlpr

##' Plot non-parametric part of the single index MIDAS regression
##'
##' Plot non-parametric part of the single index MIDAS regression
##' of unrestricted MIDAS regression
##' @param x \code{midas_r} object
##' @param term_name the term name for which the coefficients are plotted. Default is \code{NULL}, which selects the first MIDAS term
##' @param title the title string of the graph. The default is \code{NULL} for the default title.
##' @param compare the parameters for weight function to compare with the model, default is NULL
##' @param ... not used
##' @return a data frame with restricted MIDAS coefficients, unrestricted MIDAS coefficients and lower and upper confidence interval limits. The data
##' frame is returned invisibly.
##' @author Virmantas Kvedaras, Vaidotas Zemlys
##' @importFrom graphics plot points
##' @importFrom numDeriv jacobian
##' @importFrom stats na.omit approx density
##' @export
plot_sp <- function(x, term_name, title = NULL,  compare = NULL, ... ) {
    
    if (length(x$bws) > 1) stop("Plotting is currently supported for univariate non-parametric functions only")
    ti <- x$term_info[[term_name]]
    
    if (ti$term_type != "Z") stop("The term name supplied does not correspond to non-parametric term") 
    
    gfun <- function(p) {
        if(is.null(x$model_matrix_map$X)) { 
            xi <- 0 
        } else {
            xi <- x$model[, x$model_matrix_map$X] %*% x$coef_X(p)
        }
        Z <- x$model[, x$model_matrix_map$Z, drop = FALSE]
        if (is.null(ti$coef_index)) zi <- Z[, ti$midas_coef_index, drop = FALSE]
        else zi <- Z[, ti$midas_coef_index] %*% ti$weight(p[ti$coef_index])
        if (ncol(zi) == 1) zi <- as.numeric(zi)
        u <- x$model[, 1] - xi
        list(xi = xi, zi = zi, u = u, g = cv_np(u, zi, p[1:length(x$bws)], x$degree))
    }
    
    
    gg <- gfun(coef(x))
    
    kappa <- 0.5/pi # for standard Gaussian kernel
    dns <- density(gg$zi, na.rm = T) # choose a better bandwith here !!!
    f_z <- approx(dns$x, dns$y, xout = gg$zi)$y
    se_np <- sqrt(kappa*deviance(x)/f_z/(exp(coef(x)[1])*x$nobs))
    
    
    ozi <- order(gg$zi)
    
    se_npo <- se_np[ozi] 
    xio <- gg$zi[ozi]
    trmo <-  gg$g[ozi]
    pd <- data.frame(xi = xio,  term = trmo, compare = NA, lower = trmo - 1.96*se_npo, upper = trmo + 1.96*se_npo)
    
    if (!is.null(compare)) {
       pd$compare <- compare(pd$xi)
    }
    
    ylim <- range(na.omit(c(pd$term,pd$compare, pd$lower, pd$upper)))
    
    plot(pd$xi, pd$term, col = "blue", ylab = "Estimated non-parametric function", xlab = "MIDAS aggregate", ylim = ylim, type = "l")
    
    points(pd$xi, pd$compare, type = "l", col="black")
    points(pd$xi, pd$lower, type = "l", col = "grey", lty = 2)
    points(pd$xi, pd$upper, type = "l", col = "grey", lty = 2)
    
    if (is.null(title)) {
        title(main = paste0("Non-parametric estimate for term ",term_name," with weight: ",ti$weight_name))
    } else title(main = title) 
    
    invisible(pd)
    
}