Nothing
tvmvar <- function(data, # n x p data matrix
type, # p vector indicating the type of each variable
level, # p vector indivating the levels of each variable
timepoints,
estpoints, # vector of estimation points in [0,1]
bandwidth, # choice of bandwidth
... # arguments passed to mvar
)
{
# -------------------- Input Checks -------------------
args <- list(...)
lags <- args$lags
# ----- Fill in Defaults -----
if(is.null(args$lambdaSeq)) args$lambdaSeq <- NULL
if(is.null(args$lambdaSel)) args$lambdaSel <- 'EBIC'
if(is.null(args$lambdaFolds)) args$lambdaFolds <- 10
if(is.null(args$lambdaGam)) args$lambdaGam <- .25
if(is.null(args$alphaSeq)) args$alphaSeq <- 1
if(is.null(args$alphaSel)) args$alphaSel <- 'CV'
if(is.null(args$alphaFolds)) args$alphaFolds <- 10
if(is.null(args$alphaGam)) args$alphaGam <- .25
if(is.null(args$threshold)) args$threshold <- 'HW'
if(is.null(args$method)) args$method <- 'glm'
if(is.null(args$binarySign)) args$binarySign <- FALSE
if(is.null(args$scale)) args$scale <- TRUE
if(is.null(args$consec)) args$consec <- NULL
if(is.null(args$verbatim)) args$verbatim <- FALSE
if(is.null(args$warnings)) args$warnings <- FALSE
if(is.null(args$saveModels)) args$saveModels <- TRUE
if(is.null(args$saveData)) args$saveData <- FALSE
if(is.null(args$pbar)) args$pbar <- pbar <- TRUE
if(is.null(args$overparameterize)) args$overparameterize <- FALSE
if(is.null(args$signInfo)) args$signInfo <- TRUE
if(is.null(args$weights)) args$weights <- rep(1, nrow(data))
if(missing(level)) level <- NULL
if(args$verbatim) args$pbar <- FALSE
if(args$verbatim) args$warnings <- FALSE
# Switch all warnings off
if(!args$warnings) {
oldw <- getOption("warn")
options(warn = -1)
}
# ----- Create Output Object -----
tvmvar_object <- list('call' = NULL,
'wadj' = NULL,
'signs' = NULL,
'edgecolor' = NULL,
'intercepts' = NULL,
'tvmodels' = list())
# ----- Copy the Call -----
tvmvar_object$call <- list('data' = NULL,
'type' = type,
'level' = level,
'timepoints' = NULL,
'estpoints' = estpoints,
'estpointsNorm' = NULL,
'bandwidth' = bandwidth,
'lags' = args$lags,
'consec' = args$consec,
'dayvar' = args$dayvar,
'beepvar' = args$beepvar,
'weights' = args$weights,
'lambdaSeq' = args$lambdaSeq,
'lambdaSel' = args$lambdaSel,
'lambdaFolds' = args$lambdaFolds,
'lambdaGam' = args$lambdaGam,
'alphaSeq' = args$alphaSeq,
'alphaSel' = args$alphaSel,
'alphaFolds' = args$alphaFolds,
'alphaGam' = args$alphaGam,
'threshold' = args$threshold,
'method' = args$method,
'binarySign' = args$binarySign,
'scale' = args$scale,
'verbatim' = args$verbatim,
'warnings' = args$warnings,
'saveModels' = args$saveModels,
'saveData' = args$saveData,
'pbar' = args$pbar,
'overparameterize' = args$overparameterize,
'signInfo' = args$signInfo)
if(args$saveData) tvmvar_object$call$data <- data
# ----- Input Checks & compute Aux Variables -----
if(is.null(lags)) stop('No lags specified. See ?mvar.')
n <- nrow(data)
p <- ncol(data)
n_var <- n - max(lags) # this is how many rows there are after transforming the data
if(any(estpoints < 0 | estpoints > 1)) stop('Estimation points have to be specified on the unit interval [0,1].')
if(bandwidth <= 0) stop('The bandwidth parameter has to be strictly positive')
# -------------------- Compute Weightings -------------------
# Define time vector: if not provided, assume equally spaced time points
if(missing(timepoints)) {
timevec <- seq(0, 1, length = n_var)
tvmvar_object$call$timepoints <- seq(0, 1, length = n) # included for the case when reusing call in resample()
} else {
tvmvar_object$call$timepoints <- timepoints
# normalize to [0,1]
timepoints <- timepoints[-(1:max(lags))] # delete first x rows that have to be exluded by definition of VAR model
timepoints <- timepoints - min(timepoints)
timevec <- timepoints / max(timepoints)
}
tvmvar_object$call$timepoints_cut <- timevec
# Normalize time estimation points to interval [0,1]
estpoints_norm <- estpoints
tvmvar_object$call$estpointsNorm <- estpoints_norm
no_estpoints <- length(estpoints_norm)
# Compute weights for each estimation point
l_weights <- list()
for(i in 1:no_estpoints) {
l_weights[[i]] <- dnorm(timevec, mean = estpoints_norm[i], sd = bandwidth)
# Normalize to [x,1]
l_weights[[i]] <- l_weights[[i]] / max(l_weights[[i]])
# If tvmvar is used within bwSelect: set weights to zero for test samples
if(!is.null(args$zero_weights)) {
# Sanity
if(length(l_weights[[i]]) != length(args$zero_weights)) stop('Zero weights vector does not have same length as weights vector!')
# Set to zero
l_weights[[i]] <- l_weights[[i]] * args$zero_weights
}
# Add zero weights (necessary since I changed $included in lagData() s.t. it is defined over the entire input data matrix)
l_weights[[i]] <- c(rep(0, max(lags)),
l_weights[[i]])
} # end for:i (estpoints)
# -------------------- Loop over Estpoints -------------------
# Progress bar
if(args$pbar==TRUE) pb <- txtProgressBar(min = 0, max = no_estpoints, initial = 0, char="-", style = 3)
# Storage
l_mvar_models <- list()
for(i in 1:no_estpoints) {
l_mvar_models[[i]] <- mvar(data = data,
type = type,
level = level,
lambdaSeq = args$lambdaSeq,
lambdaSel = args$lambdaSel,
lambdaFolds = args$lambdaFolds,
lambdaGam = args$lambdaGam,
alphaSeq = args$alphaSeq,
alphaSel = args$alphaSel,
alphaFolds = args$alphaFolds,
alphaGam = args$alphaGam,
lags = args$lags,
consec = args$consec,
beepvar = args$beepvar,
dayvar = args$dayvar,
weights = l_weights[[i]],
threshold = args$threshold,
method = args$method,
binarySign = args$binarySign,
scale = args$scale,
verbatim = args$verbatim,
pbar = FALSE,
warnings = args$warnings,
saveModels = args$saveModels,
saveData = args$saveData,
overparameterize = args$overparameterize,
signInfo = FALSE, # to avoid msg for each model
bootstrap = args$bootstrap,
boot_ind = args$boot_ind)
# Update Progress Bar
if(args$pbar==TRUE) setTxtProgressBar(pb, i)
} # End for: timepoints
# Save into output list
tvmvar_object$tvmodels <- l_mvar_models
# -------------------- Process Output -------------------
# Restructure output to make time-varying model easier accessible
n_lags <- length(lags)
a_wadj <- a_signs <- a_edgecolor <- array(dim=c(p, p, n_lags, no_estpoints))
for(lag in 1:n_lags) {
for(i in 1:no_estpoints) {
a_wadj[ , , lag, i] <- l_mvar_models[[i]]$wadj[,,lag]
a_signs[ , , lag, i] <- l_mvar_models[[i]]$signs[,,lag]
a_edgecolor[ , , lag, i] <- l_mvar_models[[i]]$edgecolor[,,lag]
}
}
# Save intercepts
l_intercepts <- list()
for(i in 1:no_estpoints) l_intercepts[[i]] <- l_mvar_models[[i]]$intercepts
tvmvar_object$intercepts <- l_intercepts
# Save in output list
tvmvar_object$wadj <- a_wadj
tvmvar_object$signs <- a_signs
tvmvar_object$edgecolor <- a_edgecolor
# Compute effectively used Sample size (relative to n)
Ne <- lapply(l_weights, sum)
tvmvar_object$Ne <- unlist(Ne)
# Copy lagged Data if saveData == TRUE (take from first time point, because always the same)
tvmvar_object$call$data_lagged <- l_mvar_models[[1]]$call$data_lagged
# -------------------- Output -------------------
class(tvmvar_object) <- c('mgm', 'tvmvar')
if(args$pbar) {
if(args$signInfo) cat('\nNote that the sign of parameter estimates is stored separately; see ?tvmvar')
} else {
if(args$signInfo) cat('Note that the sign of parameter estimates is stored separately; see ?tvmvar')
}
return(tvmvar_object)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.