Nothing
#' @title \code{l1}-regularised Yule-Walker estimation for VAR processes
#' @description Estimates the VAR parameter matrices via \code{l1}-regularised Yule-Walker estimation
#' and innovation covariance matrix via constrained \code{l1}-minimisation.
#' @details Further information can be found in Barigozzi, Cho and Owens (2024+).
#' @param x input time series each column representing a time series variable; it is coerced into a \link[stats]{ts} object
#' @param center whether to de-mean the input \code{x}
#' @param method a string specifying the method to be adopted for VAR process estimation; possible values are:
#' \describe{
#' \item{\code{"lasso"}}{ Lasso-type \code{l1}-regularised \code{M}-estimation}
#' \item{\code{"ds"}}{ Dantzig Selector-type constrained \code{l1}-minimisation}
#' }
#' @param var.order order of the VAR process; if a vector of integers is supplied, the order is chosen via \code{tuning}
#' @param lambda \code{l1}-regularisation parameter; if \code{lambda = NULL}, \code{tuning} is employed to select the parameter
#' @param tuning.args a list specifying arguments for \code{tuning}
#' for selecting the regularisation parameter (and VAR order). It contains:
#' \describe{
#' \item{\code{tuning}}{a string specifying the selection procedure for \code{var.order} and \code{lambda}; possible values are:
#' \code{"cv"} for cross validation, and \code{"bic"} for information criterion}
#' \item{\code{n.folds}}{ if \code{tuning = "cv"}, positive integer number of folds}
#' \item{\code{penalty}}{ if \code{tuning = "bic"}, penalty multiplier between 0 and 1; if \code{penalty = NULL}, it is set to \code{1/(1+exp(dim(x)[1])/dim(x)[2]))}} by default
#' \item{\code{path.length}}{ positive integer number of regularisation parameter values to consider; a sequence is generated automatically based in this value}
#' }
#' @param do.threshold whether to perform adaptive thresholding of VAR parameter estimator with \link[fnets]{threshold}
#' @param n.iter maximum number of descent steps, by default depends on \code{var.order}; applicable when \code{method = "lasso"}
#' @param tol numerical tolerance for increases in the loss function; applicable when \code{method = "lasso"}
#' @param n.cores number of cores to use for parallel computing, see \link[parallel]{makePSOCKcluster}; applicable when \code{method = "ds"}
#' @return a list which contains the following fields:
#' \item{beta}{ estimate of VAR parameter matrix; each column contains parameter estimates for the regression model for a given variable}
#' \item{Gamma}{ estimate of the innovation covariance matrix}
#' \item{lambda}{ \code{l1}-regularisation parameter}
#' \item{var.order}{ VAR order}
#' \item{mean.x}{ if \code{center = TRUE}, returns a vector containing row-wise sample means of \code{x}; if \code{center = FALSE}, returns a vector of zeros}
#' @example R/examples/var_ex.R
#' @importFrom parallel detectCores
#' @importFrom stats as.ts
#' @references Barigozzi, M., Cho, H. & Owens, D. (2024+) FNETS: Factor-adjusted network estimation and forecasting for high-dimensional time series. Journal of Business & Economic Statistics (to appear).
#' @references Owens, D., Cho, H. & Barigozzi, M. (2024+) fnets: An R Package for Network Estimation and Forecasting via Factor-Adjusted VAR Modelling. The R Journal (to appear).
#' @export
fnets.var <- function(x,
center = TRUE,
method = c("lasso", "ds"),
lambda = NULL,
var.order = 1,
tuning.args = list(
tuning = c("cv", "bic"),
n.folds = 1,
penalty = NULL,
path.length = 10
),
do.threshold = FALSE,
n.iter = NULL,
tol = 0,
n.cores = 1) {
x <- t(as.ts(x))
p <- dim(x)[1]
n <- dim(x)[2]
if(!is.null(lambda)) lambda <- max(0, tol)
if(!is.null(n.iter)) n.iter <- posint(n.iter)
tol <- max(0, tol)
n.cores <- posint(n.cores)
tuning.args <- check.list.arg(tuning.args)
method <- match.arg(method, c("lasso", "ds"))
tuning <- match.arg(tuning.args$tuning, c("cv", "bic"))
args <- as.list(environment())
args$x <- t(args$x)
ifelse(center,mean.x <- apply(x, 1, mean), mean.x <- rep(0, p))
xx <- x - mean.x
dpca <- dyn.pca(xx, q = 0)
acv <- dpca$acv
ive <- fnets.var.internal(xx, acv, method = method,
lambda = lambda,
var.order = var.order,
tuning.args = tuning.args,
do.threshold = do.threshold,
n.iter = n.iter,
tol = tol,
n.cores = n.cores)
ive$mean.x <- mean.x
args$do.lrpc <- ive$do.lrpc <- FALSE
ive$q <- 0
attr(ive, "class") <- "fnets"
attr(ive, "factor") <- "none"
attr(ive, "args") <- args
return(ive)
}
#' @title internal function for \code{fnets.var}
#' @keywords internal
fnets.var.internal <- function(xx,
acv,
method = c("lasso", "ds"),
lambda = NULL,
var.order = 1,
tuning.args = list(
tuning = c("cv", "bic"),
n.folds = 1,
penalty = NULL,
path.length = 10
),
do.threshold = FALSE,
n.iter = NULL,
tol = 0,
n.cores = 1){
for (ii in 1:length(var.order)) var.order[ii] <- posint(var.order[ii])
method <- match.arg(method, c("lasso", "ds"))
tuning <- match.arg(tuning.args$tuning, c("cv", "bic"))
if(tuning == "cv") {
icv <- yw.cv(
xx,
method = method,
lambda.max = NULL,
var.order = var.order,
n.folds = tuning.args$n.folds,
path.length = tuning.args$path.length,
q = 0,
kern.bw = NULL,
n.cores = n.cores
)
}
if(tuning == "bic") {
icv <- yw.ic(
xx,
method = method,
lambda.max = NULL,
var.order = var.order,
penalty = tuning.args$penalty,
path.length = tuning.args$path.length,
q = 0,
kern.bw = NULL
)
}
mg <- make.gg(acv$Gamma_i, icv$order.min)
gg <- mg$gg
GG <- mg$GG
if(method == "lasso"){
if(is.null(n.iter)) n.iter <- icv$order.min*100
ive <-
var.lasso(
GG,
gg,
lambda = icv$lambda,
symmetric = "min",
n.iter = n.iter,
tol = tol
)
}
if(method == "ds")
ive <-
var.dantzig(
GG,
gg,
lambda = icv$lambda,
symmetric = "min",
n.cores = n.cores
)
ive$var.order <- icv$order.min
if(do.threshold)
ive$beta <- threshold(ive$beta)$thr.mat
attr(ive, "data") <- icv
return(ive)
}
#' @title Lasso-type estimator of VAR processes via \code{l1}-regularised \code{M}-estimation
#' @keywords internal
var.lasso <-
function(GG,
gg,
lambda,
symmetric = "min",
n.iter = 100,
tol = 0) {
backtracking <- TRUE
p <- ncol(gg)
d <- nrow(gg) / ncol(gg)
ii <- 0
t.new <- t <- 1
x <- gg * 0
x.new <- y <- x
diff.val <- tol - 1
if(backtracking) {
L <- norm(GG, "F") / 5
gamma <- 2
} else {
L <- norm(GG, "F")
}
obj.val <- rel.err <- c()
while (ii < n.iter & diff.val < tol) {
ii <- ii + 1
if(backtracking) {
L.bar <- L
found <- FALSE
while (!found) {
prox <- prox.func(y, lambda, L = 2 * L.bar, GG, gg)
ifelse(f.func(GG, gg, prox) <= Q.func(prox, y, L.bar, GG, gg),
found <- TRUE,
L.bar <- L.bar * gamma
)
}
L <- L.bar
} else {
prox <- prox.func(y, lambda, L = 2 * L, GG, gg)
}
x <- x.new
x.new <- prox
t <- t.new
t.new <- (1 + sqrt(1 + 4 * t^2)) / 2
y <- x.new + (t - 1) / t.new * (x.new - y)
obj.val <-
c(obj.val, f.func(GG, gg, x.new) + lambda * sum(abs(x.new)))
if(ii > 1)
diff.val <- obj.val[ii] - obj.val[ii - 1]
}
A <- t(x.new)
Gamma <- GG[1:p, 1:p]
for (ll in 1:d)
Gamma <-
Gamma - A[, (ll - 1) * p + 1:p] %*% gg[(ll - 1) * p + 1:p, ]
Gamma <- make.symmetric(Gamma, symmetric)
out <-
list(
beta = x.new,
Gamma = Gamma,
lambda = lambda
)
return(out)
}
#' @title Dantzig selector-type estimator of VAR processes via constrained \code{l1}-minimisation
#' @importFrom parallel makePSOCKcluster stopCluster detectCores
#' @importFrom doParallel registerDoParallel
#' @importFrom foreach foreach %dopar%
#' @importFrom lpSolve lp
#' @keywords internal
var.dantzig <-
function(GG,
gg,
lambda,
symmetric = "min",
n.cores = 1) {
p <- dim(gg)[2]
d <- dim(gg)[1] / dim(gg)[2]
beta <- gg * 0
f.obj <- rep(1, 2 * p * d)
f.con <- rbind(-GG, GG)
f.con <- cbind(f.con, -f.con)
f.dir <- rep("<=", 2 * p * d)
cl <- parallel::makePSOCKcluster(n.cores)
doParallel::registerDoParallel(cl)
ii <- 1
beta <-
foreach::foreach(
ii = 1:p,
.combine = "cbind",
.multicombine = TRUE,
.export = c("lp")
) %dopar% {
b1 <- rep(lambda, p * d) - gg[, ii]
b2 <- rep(lambda, p * d) + gg[, ii]
f.rhs <- c(b1, b2)
lpout <- lpSolve::lp("min", f.obj, f.con, f.dir, f.rhs)
lpout$solution[1:(p * d)] - lpout$solution[-(1:(p * d))]
}
parallel::stopCluster(cl)
A <- t(beta)
Gamma <- GG[1:p, 1:p]
for (ll in 1:d)
Gamma <-
Gamma - A[, (ll - 1) * p + 1:p] %*% gg[(ll - 1) * p + 1:p, ]
Gamma <- make.symmetric(Gamma, symmetric)
out <- list(beta = beta,
Gamma = Gamma,
lambda = lambda)
return(out)
}
#' @title Cross validation for factor-adjusted VAR estimation
#' @keywords internal
yw.cv <- function(xx,
method = c("lasso", "ds"),
lambda.max = NULL,
var.order = 1,
n.folds = 1,
path.length = 10,
q = 0,
kern.bw = NULL,
n.cores = 1) {
n <- ncol(xx)
p <- nrow(xx)
if(is.null(kern.bw))
kern.bw <- 4 * floor((n / log(n))^(1 / 3))
if(is.null(lambda.max))
lambda.max <- max(abs(xx %*% t(xx) / n)) * 1
lambda.path <-
round(exp(seq(
log(lambda.max), log(lambda.max * .0001), length.out = path.length
)), digits = 10)
cv.err.mat <-
matrix(0, nrow = path.length, ncol = length(var.order))
dimnames(cv.err.mat)[[1]] <- lambda.path
dimnames(cv.err.mat)[[2]] <- var.order
ind.list <- split(1:n, ceiling(n.folds * (1:n) / n))
for (fold in 1:n.folds) {
train.ind <- 1:ceiling(length(ind.list[[fold]]) * .5)
train.x <- xx[, ind.list[[fold]][train.ind]]
test.x <- xx[, ind.list[[fold]][-train.ind]]
train.acv <-
dyn.pca(train.x,
q = q,
kern.bw = kern.bw,
mm = max(var.order))$acv$Gamma_i
test.acv <-
dyn.pca(test.x,
q = q,
kern.bw = kern.bw,
mm = max(var.order))$acv$Gamma_i
for (jj in 1:length(var.order)) {
mg <- make.gg(train.acv, var.order[jj])
gg <- mg$gg
GG <- mg$GG
mg <- make.gg(test.acv, var.order[jj])
test.gg <- mg$gg
test.GG <- mg$GG
for (ii in 1:path.length) {
if(method == "ds")
train.beta <-
var.dantzig(GG, gg, lambda = lambda.path[ii], n.cores = n.cores)$beta
if(method == "lasso")
train.beta <-
var.lasso(GG, gg, lambda = lambda.path[ii])$beta
beta.gg <- t(train.beta) %*% test.gg
cv.err.mat[ii, jj] <- cv.err.mat[ii, jj] +
sum(diag(
test.acv[, , 1] - beta.gg - t(beta.gg) + t(train.beta) %*% test.GG %*% (train.beta)
))
}
}
}
cv.err.mat[cv.err.mat < 0] <- Inf
lambda.min <-
min(lambda.path[apply(cv.err.mat, 1, min) == min(apply(cv.err.mat, 1, min))])
order.min <-
min(var.order[apply(cv.err.mat, 2, min) == min(apply(cv.err.mat, 2, min))])
out <-
list(
lambda = lambda.min,
order.min = order.min,
error = cv.err.mat,
lambda.path = lambda.path,
var.order = var.order
)
return(out)
}
# ebic
#' @title logarithmic factorial of `n`
#' @keywords internal
logfactorial <- function(n)
sum(log(1:max(n, 1)))
#' @title full likelihood
#' @keywords internal
f.func.full <- function(GG, gg, A) {
return(0.5 * sum(diag(
GG[1:ncol(A), 1:ncol(A)] + t(A) %*% GG %*% A - t(A) %*% gg - t(gg) %*% A
)))
}
f.func.mat <- function(GG, gg, A) {
return(0.5 * (diag(
GG[1:ncol(A), 1:ncol(A)] + t(A) %*% GG %*% A - t(A) %*% gg - t(gg) %*% A
)))
}
#' @title extended Bayesian Information Criterion
#' @keywords internal
ebic <- function(object, n, penalty = 0) {
penalty <- max( min(1,penalty) , 0)
beta <- object$idio.var$beta
p <- ncol(beta)
d <- nrow(beta) / ncol(beta)
sparsity <- sum(beta != 0)
mg <- make.gg(object$acv$Gamma_i, d)
gg <- mg$gg
GG <- mg$GG
n / 2 * log(2 * f.func.full(GG, gg, beta)) + sparsity * log(n) +
2 * penalty * (logfactorial(p^2 * d) - logfactorial(sparsity) - logfactorial(p ^
2 * d - sparsity))
}
#' @title Information criterion for factor-adjusted VAR estimation
#' @keywords internal
yw.ic <- function(xx,
method = c("lasso", "ds"),
lambda.max = NULL,
var.order = 1,
penalty = NULL,
path.length = 10,
q = 0,
kern.bw = NULL) {
n <- ncol(xx)
p <- nrow(xx)
if(is.null(kern.bw))
kern.bw <- 4 * floor((n / log(n))^(1 / 3))
if(is.null(penalty))
penalty <- 1 / (1 + exp(5 - p / n))
if(is.null(lambda.max))
lambda.max <- max(abs(xx %*% t(xx) / n)) * 1
lambda.path <-
round(exp(seq(
log(lambda.max), log(lambda.max * .0001), length.out = path.length
)), digits = 10)
ic.err.mat <-
matrix(0, nrow = path.length, ncol = length(var.order))
dimnames(ic.err.mat)[[1]] <- lambda.path
dimnames(ic.err.mat)[[2]] <- var.order
acv <-
dyn.pca(xx,
q = q,
kern.bw = kern.bw,
mm = max(var.order))$acv$Gamma_i
for (jj in 1:length(var.order)) {
mg <- make.gg(acv, var.order[jj])
gg <- mg$gg
GG <- mg$GG
test.gg <- mg$gg
test.GG <- mg$GG
for (ii in 1:path.length) {
if(method == "ds")
beta <- var.dantzig(GG, gg, lambda = lambda.path[ii])$beta
if(method == "lasso")
beta <- var.lasso(GG, gg, lambda = lambda.path[ii])$beta
beta <- threshold(beta)$thr.mat
sparsity <- sum(beta[, 1] != 0)
ifelse(sparsity == 0,pen <- 0,
pen <-2 * penalty * (logfactorial(var.order[jj] * p) - log(sparsity) - log(var.order[jj] * p - sparsity)))
ic.err.mat[ii, jj] <-
ic.err.mat[ii, jj] + n / 2 * log(2 * f.func.mat(GG, gg, beta)[1]) + sparsity * log(n) + pen
}
}
#ic.err.mat[ic.err.mat < 0] <- Inf
ic.err.mat[is.nan(ic.err.mat)] <- Inf
lambda.min <-
min(lambda.path[apply(ic.err.mat, 1, min) == min(apply(ic.err.mat, 1, min))])
order.min <-
min(var.order[apply(ic.err.mat, 2, min) == min(apply(ic.err.mat, 2, min))])
out <-
list(
lambda = lambda.min,
order.min = order.min,
error = ic.err.mat,
lambda.path = lambda.path,
var.order = var.order
)
return(out)
}
#' @title Plotting output for tuning parameter selection in fnets
#' @description Tuning plots for S3 objects of class \code{fnets}.
#' Produces up to two plots visualising CV and IC procedures for selecting tuning parameters and the VAR order.
#' @details See Owens, Cho and Barigozzi (2024+) for further details.
#' @param x \code{fnets} object
#' @param ... additional arguments
#' @return CV/IC plot for the VAR component, and CV plot for the lrpc component (when \code{x$do.lrpc = TRUE}).
#' @seealso \link[fnets]{fnets}
#' @references Owens, D., Cho, H. & Barigozzi, M. (2024+) fnets: An R Package for Network Estimation and Forecasting via Factor-Adjusted VAR Modelling. The R Journal (to appear).
#' @importFrom graphics par abline box axis legend matplot
#' @keywords internal
tuning_plot <- function(x, ...){
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
data <- attr(x, "data")
lambda.min <- data$lambda
order.min <- data$order.min
error <- data$error
lambda.path <- data$lambda.path
var.order <- data$var.order
args <- attr(x, "args")
par(mfrow = c(1,1+x$do.lrpc)) ## check
if(args$tuning == "bic") {
ylab = "IC"
main = "IC for VAR parameter estimation"
} else if(args$tuning == "cv") {
ylab = "CV"
main = "CV for VAR parameter estimation"
}
par(xpd=FALSE)
error.plot <- error
if(any(is.infinite(error.plot))) {
error.plot[which(is.infinite(error.plot))] <- NA
keep <- colSums(!is.na(error.plot)) > 0
error.plot <- error.plot[,keep]
} else keep <- rep(TRUE, length(var.order))
matplot(
lambda.path,
error.plot,
type = "b",
col = (2:(max(var.order) + 1))[keep],
pch = (2:(max(var.order) + 1))[keep],
log = "x",
xlab = "lambda (log scale)",
ylab = ylab,
main = main
)
abline(v = lambda.min)
legend(
"topleft",
legend = var.order[keep],
col = (2:(max(var.order) + 1))[keep],
pch = (2:(max(var.order) + 1))[keep],
lty = 1
)
if (x$do.lrpc) {
data <- attr(x$lrpc, "data")
plot(
data$eta.path,
data$cv.err,
type = "b",
col = 2,
pch = 2,
log = "x",
xlab = "eta (log scale)",
ylab = "CV error",
main = "CV for (LR)PC matrix estimation"
)
abline(v = data$eta)
}
}
#' @title Forecasting idiosyncratic VAR process
#' @description Produces forecasts of the idiosyncratic VAR process
#' for a given forecasting horizon by estimating the best linear predictors
#' @param object \code{fnets} object
#' @param x input time series, with each row representing a variable
#' @param cpre output of \link[fnets]{common.predict}
#' @param n.ahead forecast horizon
#' @return a list containing
#' \item{is}{ in-sample estimator of the idiosyncratic component (with each column representing a variable)}
#' \item{fc}{ forecasts of the idiosyncratic component for a given forecasting horizon \code{h} (with each column representing a variable)}
#' \item{n.ahead}{ forecast horizon}
#' @examples
#' \dontrun{
#' out <- fnets(data.unrestricted,
#' do.lrpc = FALSE, var.args = list(n.cores = 2))
#' cpre <- common.predict(out)
#' ipre <- idio.predict(out, cpre)
#' }
#' @importFrom stats as.ts
#' @keywords internal
idio.predict <- function(object, x, cpre, n.ahead = 1) {
p <- dim(x)[1]
n <- dim(x)[2]
# if(attr(object, 'factor') == 'dynamic'){
xx <- x - object$mean.x
beta <- object$idio.var$beta
if(is.null(beta)) beta <- object$beta
d <- dim(beta)[1] / p
A <- t(beta)
is <- xx - t(cpre$is)
if(n.ahead >= 1) {
fc <- matrix(0, nrow = p, ncol = n.ahead)
for (ii in 1:n.ahead) {
for (ll in 1:d)
fc[, ii] <-
fc[, ii] + A[, p * (ll - 1) + 1:p] %*% is[, n + ii - ll]
is <- cbind(is, fc[, ii])
}
} else {
fc <- NA
}
out <- list(is = as.ts(t(is[, 1:n])), fc = as.ts(t(fc)), n.ahead = n.ahead)
return(out)
}
#' @keywords internal
make.gg <- function(acv, d) {
p <- dim(acv)[1]
gg <- matrix(0, nrow = p * d, ncol = p)
GG <- matrix(0, p * d, p * d)
for (ll in 1:d) {
gg[(ll - 1) * p + 1:p, ] <- acv[, , ll + 1]
for (lll in ll:d) {
GG[(ll - 1) * p + 1:p, (lll - 1) * p + 1:p] <-
t(acv[, , 1 + lll - ll])
GG[(lll - 1) * p + 1:p, (ll - 1) * p + 1:p] <-
acv[, , 1 + lll - ll]
}
}
out <- list(gg = gg, GG = GG)
return(out)
}
#' @keywords internal
f.func <- function(GG, gg, A) {
return(.5 * sum(diag(t(A) %*% GG %*% A - 2 * t(A) %*% gg)))
}
#' @keywords internal
gradf.func <- function(GG, gg, A) {
return(GG %*% (A) - gg)
}
#' @keywords internal
Q.func <- function(A, A.up, L, GG, gg) {
Adiff <- (A - A.up)
return(f.func(GG, gg, A.up) + sum(Adiff * gradf.func(GG, gg, A.up)) + 0.5 * L * norm(Adiff, "F") ^
2)
}
#' @keywords internal
prox.func <- function(B, lambda, L, GG, gg) {
b <- B - (1 / L) * gradf.func(GG, gg, B)
sgn <- sign(b)
ab <- abs(b)
sub <- ab - 2 * lambda / L
sub[sub < 0] <- 0
out <- sub * sgn
return(as.matrix(out))
}
#' @title Threshold the entries of the input matrix at a data-driven level
#' @description Threshold the entries of the input matrix at a data-driven level.
#' This can be used to perform edge selection for VAR parameter, inverse innovation covariance, and long-run partial correlation networks.
#' @details See Owens, Cho & Barigozzi (2024+) for more information on the threshold selection process
#' @param mat input parameter matrix
#' @param path.length number of candidate thresholds
#' @return an S3 object of class \code{threshold}, which contains the following fields:
#' \item{threshold}{ data-driven threshold}
#' \item{thr.mat}{ thresholded input matrix}
#' @seealso \link[fnets]{plot.threshold}, \link[fnets]{print.threshold}
#' @example R/examples/thresh_ex.R
#' @importFrom graphics par
#' @references Owens, D., Cho, H. & Barigozzi, M. (2024+) fnets: An R Package for Network Estimation and Forecasting via Factor-Adjusted VAR Modelling. The R Journal (to appear).
#' @export
threshold <- function(mat,
path.length = 500) {
path.length <- posint(path.length)
p <- dim(mat)[1]
M <- max(abs(mat), 1e-3)
m <- max(min(abs(mat)), M * .01, 1e-4)
rseq <-
round(exp(seq(log(M), log(m), length.out = path.length)), digits = 10)
cusum <- ratio <- rseq * 0
for (ii in 1:path.length) {
A <- mat
A[abs(A) < rseq[ii]] <- 0
edges <- sum(A != 0)
ratio[ii] <- edges / (prod(dim(mat)) - edges)
}
dif <- diff(ratio) / diff(rseq)
for (ii in 2:(path.length - 1))
cusum[ii] <-
(mean(dif[2:ii - 1]) - mean(dif[ii:(path.length - 1)])) * (ii / sqrt(path.length)) * (1 - ii / path.length)
thr <- rseq[which.max(abs(cusum))]
A <- mat
A[abs(A) < thr] <- 0
out <- list(threshold = thr, thr.mat = A)
attr(out, "seqs") <- list(rseq = rseq, ratio = ratio, dif = dif, cusum = cusum)
attr(out, "class") <- "threshold"
return(out)
}
#' @title Plotting the thresholding procedure
#' @method plot threshold
#' @description Plotting method for S3 objects of class \code{threshold}.
#' Produces a plot visualising three diagnostics for the thresholding procedure, with threshold values t_k (x axis) against
#' (i) Ratio_k, the ratio of the number of non-zero to zero entries in the matrix, as the threshold varies
#' (ii) Diff_k, the first difference of \code{Ratio_k}
#' (iii) |CUSUM_k|, the absolute scaled cumulative sums of \code{Diff_k}
#' @details See Owens, Cho and Barigozzi (2024+) for further details.
#' @param x \code{threshold} object
#' @param plots logical vector, which plots to use (Ratio, Diff, CUSUM respectively)
#' @param ... additional arguments
#' @return A network plot produced as per the input arguments
#' @references Owens, D., Cho, H. & Barigozzi, M. (2024+) fnets: An R Package for Network Estimation and Forecasting via Factor-Adjusted VAR Modelling. The R Journal (to appear).
#' @seealso \link[fnets]{threshold}
#' @example R/examples/thresh_ex.R
#' @export
plot.threshold <- function(x, plots = c(TRUE, FALSE, TRUE), ...){
seqs <- attr(x, "seqs")
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
if(sum(plots) == 0 | length(plots) != 3){
warning("setting all plots to TRUE")
plots <- c(TRUE, TRUE, TRUE)
}
par(mfrow = c(1, sum(plots)))
if(plots[1]){
plot(seqs$rseq, seqs$ratio, type = "l", xlab = "Threshold", ylab = "Ratio")
abline(v = x$threshold)
}
if(plots[2]){
plot(seqs$rseq[-1], seqs$dif, type = "l", xlab = "Threshold", ylab = "Diff")
abline(v = x$threshold)
}
if(plots[3]){
plot(seqs$rseq, abs(seqs$cusum), type = "l", xlab = "Threshold", ylab = "|Cusum|")
abline(v = x$threshold)
}
}
#' @title Print threshold
#' @method print threshold
#' @description Prints a summary of a \code{threshold} object
#' @param x \code{threshold} object
#' @param ... not used
#' @return NULL, printed to console
#' @seealso \link[fnets]{threshold}
#' @example R/examples/thresh_ex.R
#' @export
print.threshold <- function(x,
...){
cat(paste("Thresholded matrix \n"))
cat(paste("Threshold: ", x$threshold, "\n", sep = ""))
cat(paste("Non-zero entries: ", sum(x$thr.mat != 0), "/", prod(dim(x$thr.mat)), "\n", sep = ""))
}
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.