Nothing
# -------------------------------------------------------------------------------------------
#' Fit best production function model
#'
#' @description Finds the most suitable model for the NAWRU and the TFP trend according to
#' the BIC or the RMSE. The function computes the output gap based on the chosen models.
#'
#' @param tsl A list of time series objects, see details.
#' @param type The variance restriction type. Possible options are \code{"basic"},
#' \code{"hp"}, see \code{initializeRestr}. The default is \code{type = "hp"}.
#' @param q Quantile for the Inverse Gamma distribution (only used if \code{type = "hp"}),
#' see \code{initializeRestr}. The default is \code{q = 0.01}.
#' @param method The estimation method. Options are maximum likelihood estimation \code{"MLE"}
#' and bayesian estimation \code{"bayesian"}. If \code{method = c("MLE", "bayesian")} the
#' NAWRU is fitted by MLE and the TFP trend by Bayesian methods. The default is
#' \code{method = "MLE"}.
#' @param criterion Model selection criterion. Options are the Bayesian information criterion
#' \code{"BIC"} and the root mean squared error \code{"RMSE"}, both computed for the second
#' observation equation. The default is \code{criterion = "BIC"}. For Bayesian estimation
#' and \code{criterion = "RMSE"}, the mean RMSE is used.
#' @param fast Boolean, indicating whether a "fast" procedure should be used, see details.
#' @param nModels Integer, the maximum number of models for each unobserved component model.
#' @param nawruPoss List with possible model specifications for the NAWRU, see details.
#' @param tfpPoss List with possible model specifications for the NAWRU, see details.
#' @param auto If \code{auto = "NAWRU"} or \code{auto = "TFP"}, the function only
#' finds the most suitable NAWRU or TFP model, respectively. The default is
#' \code{auto = "gap"}.
#'
#' @details For \code{fast = TRUE}, the function pre-selects suitable models by applying the
#' following procedure: A HP-filtered trend is computed based on which the best trend and
#' cycle models are chosen according to the BIC. Also based on the HP trend, a variety of
#' different specifications for the second observation equation are estimated in a
#' univariate regression and the best models are selected via the BIC. The \code{nModels}
#' best models are subsequently estimated in the usual bivariate unobserved component
#' model. For \code{fast = FALSE}, a variety of models is estimated in the usual bivariate
#' unobserved component framework.
#' @details The input component \code{nawruPoss} is a list containing a (sub-) set of the
#' following components:
#' \describe{
#' \item{maxCycleLag}{Maximum cycle lag included in the second observation equation.}
#' \item{trend}{Trend model specification.}
#' \item{cycle}{Cycle model specification.}
#' \item{errorARmax}{Maximum autoregressive order of the error term in the second
#' observation equation.}
#' \item{errorMAmax}{Maximum moving average order of the error term in the second
#' observation equation.}
#' \item{type}{Type of Phillip's curve.}
#' \item{exoNames}{Names of the exogenous variables potentially included in the Phillip's
#' curve (need to be included in the list of time series \code{tsl}).}
#' \item{signal-to-noise}{Signal-to-noise ratio.}
#' }
#' @details The input component \code{tfpPoss} is a list containing a (sub-) set of the
#' following components:
#' \describe{
#' \item{maxCycleLag}{Maximum cycle lag included in the second observation equation.}
#' \item{trend}{Trend model specification.}
#' \item{cycle}{Cycle model specification.}
#' \item{cubsARmax}{Maximum CUBS autoregressive order.}
#' \item{errorARmax}{Maximum autoregressive order of the error term in the second
#' observation equation.}
#' \item{errorMAmax}{Maximum moving average order of the error term in the second
#' observation equation.}
#' \item{signal-to-noise}{Signal-to-noise ratio.}
#' }
#' @details The list of time series \code{tsl} needs to have the following components
#' (plus those series included in the list component \code{exoNames} in \code{nawruPoss}):
#' \describe{
#' \item{ur}{Unemployment rate.}
#' \item{nulc}{Nominal Unit labor costs, if \code{type = "TKP"}.}
#' \item{rulc}{Real unit labor costs, if \code{type = "NKP"}.}
#' \item{tfp}{Total factor productivity.}
#' \item{cubs}{Capacity utilization economic sentiment indicator.}
#' \item{lfnd}{Labor force non-domestic (unit: 1000 persons).}
#' \item{parts}{Participation rate.}
#' \item{ahours}{Average hours worked (unit: hours).}
#' \item{gdp}{Gross domestic product at constant prices (unit: bn National currency, code: OVGD).}
#' \item{k}{Net capital stock at constant prices: total economy (unit: bn National currency, code: OKND).}
#' \item{popw}{Population: 15 to 64 years (unit: 1000 persons, code: NPAN).}
#' }
#' @details The set of tested models is extensive but not exhaustive. The best model is
#' solely based on convergence and the chosen criterion (RMSE or BIC). A manual check
#' of the results is highly recommended.
#' @details In some cases, more than \code{nModels} are checked. For instance, if a
#' re-parametrized and regular AR(2) process are options for the cycle.
#'
#' @return A list containing three components: \code{gap} (the best model of class \code{"gap"}),
#' \code{tfp} (a nested list of TFP models, fitted objects and model fit criteria), \code{nawru}
#' (a nested list of NAWRU models, fitted objects and model fit criteria). The lists \code{nawru}
#' and \code{tfp} contain a list of models, a list of fitted objects and a dataframe \code{info},
#' which contains
#' \item{loglik}{log-likelihood function at optimum}
#' \item{AIC}{Akaike information criterion}
#' \item{BIC}{Bayesian information criterion}
#' \item{HQC}{Hannan-Quinn information criterion}
#' \item{RMSE}{Root mean squared error}
#' \item{R2}{Coefficient of determination (R squared)}
#' \item{signal-to-noise}{Signal-to-noise ratio}
#' \item{LjungBox}{p-value of Ljung-Box test for autocorrelation (H0 = no autocorrelation)}
#' \item{convergence}{0 indicates convergence of the optimization}
#' \item{rrange}{relative range of trend series w.r.t original series}
#' \item{neg}{1 indicates that negative values are present in the trend series}
#' \item{rev}{relative excess volatility w.r.t. original series (stationary series)}
#' \item{rsd}{relative standard deviation w.r.t. original series (stationary series)}
#' \item{cor}{correlation between trend and original series (stationary series)}
#' \item{msdtg}{mean standardized deviation (stationary trend)}
#' \item{magtg}{mean absolute growth of trend (stationaty trend)}
#' \item{drop}{1 indicates the model should be dropped}
#'
#' @export
autoGapProd <- function(tsl,
type = "hp",
q = 0.01,
method = "MLE",
criterion = "BIC",
fast = TRUE,
nModels = 5,
nawruPoss = list(maxCycleLag = 2,
trend = c("RW2", "DT"),
cycle = c("AR1", "AR2"),
errorARmax = 1,
errorMAmax = 0,
type = c("TKP", "NKP"),
exoNames = c("ws", "prod", "tot"),
signalToNoise = NULL),
tfpPoss = list(maxCycleLag = 2,
trend = c("RW2", "DT"),
cycle = c("AR1", "AR2", "RAR2"),
cubsARmax = 0,
errorARmax = 1,
errorMAmax = 0,
signalToNoise = NULL),
auto = "gap") {
nawruPossBase <- list(maxCycleLag = 5,
trend = c("RW2", "DT"),
cycle = c("AR1", "AR2"),
errorARmax = 2,
errorMAmax = 0,
type = c("TKP", "NKP"),
exoNames = c("ws", "prod", "tot"),
signalToNoise = NULL)
tfpPossBase <- list(maxCycleLag = 5,
trend = c("RW2", "DT"),
cycle = c("AR1", "AR2", "RAR2"),
cubsARmax = 2,
errorARmax = 2,
errorMAmax = 0,
signalToNoise = NULL)
nawruPoss <- c(nawruPoss, nawruPossBase[!(names(nawruPossBase) %in% names(nawruPoss))])
tfpPoss <- c(tfpPoss, tfpPossBase[!(names(tfpPossBase) %in% names(tfpPoss))])
if (length(method) == 1) {
method <- rep(method, 2)
}
result <- list()
# NAWRU ------------------------------------
if (auto == "gap" || auto == "NAWRU") {
if (fast == TRUE) {
comb <- autoNAWRUmodel(tsl = tsl, poss = nawruPoss, nModels = nModels)
comb <- lapply(comb, function(x) {
y <- names(x)
index <- which(names(x) == "BIC")
x <- x[-index]
x
})
} else {
# nawru possibilities
possN <- list()
possN$cycleLag <- lapply(0:nawruPoss$maxCycleLag, function(x) (0:10)[0:x + 1])
possN$trend <- nawruPoss$trend
possN$cycle <- nawruPoss$cycle
possN$errorAR <- c(0:nawruPoss$errorARmax)
possN$errorMA <- c(0:nawruPoss$errorMAmax)
possN$type <- nawruPoss$type
possN$exoType[[1]] <- NULL
if (!is.null(nawruPoss$exoNames)) {
possN$exoType[[2]] <- initializeExo(varNames = nawruPoss$exoNames)
possN$exoType[[2]][1, , "difference"] <- 2
possN$exoType[[2]][2, , "difference"] <- 1
possN$exoType[[2]][1, , "lag"] <- 0
possN$exoType[[2]][2, , "lag"] <- 1
# if (length(nawruPoss$exoNames) > 1) {
# for (k in 1:length(nawruPoss$exoNames)) {
# possN$exoType[[2 + k]] <- initializeExo(varNames = nawruPoss$exoNames[k])
# possN$exoType[[2 + k]][1, , "difference"] <- 2
# possN$exoType[[2 + k]][2, , "difference"] <- 1
# possN$exoType[[2 + k]][1, , "lag"] <- 0
# possN$exoType[[2 + k]][2, , "lag"] <- 1
# }
# }
}
# number of models
comb <- expand.grid(lapply(possN, function(x) 1:length(x)))
comb <- lapply(1:nrow(comb), function(x) {
tmp <- lapply(1:ncol(comb), function(y) {
(possN[[y]][[comb[x,y]]] )
})
names(tmp) <- names(possN)
tmp
})
}
# define and fit models
res <- helper_model_fit(FUNmodel = NAWRUmodel,
FUNfit = fit.NAWRUmodel,
tsl = tsl,
comb = comb,
poss = nawruPoss,
method = method[1],
type = type,
q = q,
modelName = "NAWRU")
# compute criteria for model selection
crit <- helper_fit_comparison(fit = res$fit,
E1name = "ur",
E1Trendname = "nawru",
fitBayes = res$fitBayes)
# eliminate models
info <- cbind(helper_model_comparison(models = res$model), crit$info)
info$drop <- 1 * (info$convergence != 0
| info$neg == 1
| info$R2 < 0
| info$LjungBox < 0.1)
info <- info[order(info$drop, info[[criterion]]),]
# order and save results
index <- as.numeric(rownames(info))
nawru <- list()
nawru$model <- res$model[index]
nawru$fit <- res$fit[index]
nawru$info <- info
if (method[1] == "bayesian") {
# eliminate models
info <- cbind(helper_model_comparison(models = res$model), crit$infoBayes)
info$drop <- 1 * (info$R2 < 0
| info$neg == 2)
criterion_bayes <- grep(criterion, c("MRMSE", "R2"), value = TRUE)
if (length(criterion_bayes) == 0) { criterion_bayes <- "MRMSE" }
info <- info[order(info$drop, info[[criterion_bayes]]),]
# order and save results
index <- as.numeric(rownames(info))
nawru$modelBayes <- res$model[index]
nawru$fitBayes <- res$fitBayes[index]
nawru$infoBayes <- info
}
result$nawru <- nawru
}
# TFP ------------------------------------
if (auto == "gap" || auto == "TFP") {
if (fast == TRUE) {
comb <- autoTFPmodel(tsl = tsl, poss = tfpPoss, nModels = nModels)
comb <- lapply(comb, function(x) {
y <- names(x)
index <- which(names(x) == "BIC")
x <- x[-index]
x
})
} else {
possT <- list()
possT$cycleLag <- lapply(0:tfpPoss$maxCycleLag, function(x) (0:10)[0:x + 1])
possT$trend <- tfpPoss$trend
possT$cycle <- tfpPoss$cycle
possT$cubsAR <- tfpPoss$cubsARmax
possT$errorAR <- c(0:tfpPoss$errorARmax)
possT$errorMA <- c(0:tfpPoss$errorMAmax)
comb <- expand.grid(lapply(possN, function(x) 1:length(x)))
comb <- lapply(1:nrow(comb), function(x) {
tmp <- lapply(1:ncol(comb), function(y) {
(possN[[y]][[comb[x,y]]] )
})
names(tmp) <- names(possN)
tmp })
}
# define and fit models
res <- helper_model_fit(FUNmodel = TFPmodel,
FUNfit = fit.TFPmodel,
tsl = tsl,
comb = comb,
poss = tfpPoss,
method = method[2],
type = type,
q = q,
modelName = "TFP")
# compute criteria for model selection
crit <- helper_fit_comparison(fit = res$fit,
E1name = "logtfp",
E1Trendname = "tfpTrend",
E1trans = exp,
fitBayes = res$fitBayes)
# eliminate models
info <- cbind(helper_model_comparison(models = res$model), crit$info)
info$drop <- 1 * (info$convergence != 0
| info$neg == 1
| info$R2 < 0
| info$LjungBox < 0.1)
info <- info[order(info$drop, info[[criterion]]),]
# order and save results
index <- as.numeric(rownames(info))
tfp <- list()
tfp$model <- res$model[index]
tfp$fit <- res$fit[index]
tfp$info <- info
if (method[2] == "bayesian") {
# eliminate models
info <- cbind(helper_model_comparison(models = res$model), crit$infoBayes)
info$drop <- 1 * (info$R2 < 0
| info$neg == 1)
criterion_bayes <- grep(criterion, c("MRMSE", "R2"), value = TRUE)
if (length(criterion_bayes) == 0) { criterion_bayes <- "MRMSE" }
info <- info[order(info$drop, info[[criterion_bayes]]),]
# order and save results
index <- as.numeric(rownames(info))
tfp$modelBayes <- res$model[index]
tfp$fitBayes <- res$fitBayes[index]
tfp$infoBayes <- info
}
result$tfp <- tfp
}
# gap --------------------------------------
if (auto == "gap") {
if (nrow(tfp$info) >= 1 & nrow(nawru$info) >= 1) {
if (method[1] == "MLE") {
nawrufit <- nawru$fit[[1]]
} else {
nawrufit <- nawru$fitBayes[[1]]
}
if (method[2] == "MLE") {
tfpfit <- tfp$fit[[1]]
} else {
tfpfit <- tfp$fitBayes[[1]]
}
bestgap <- gapProd(tsl = tsl, NAWRUfit = nawrufit, TFPfit = tfpfit)
} else {
warning("No valid model left.")
bestgap <- NULL
}
result$gap <- bestgap
}
return(result)
}
# -------------------------------------------------------------------------------------------
#' NAWRU model suggestion
#'
#' @description Finds the most suitable NAWRU models.
#'
#' @param poss A list with the characteristics of possible models (see \code{autoGapProd)}
#' @inheritParams autoGapProd
#'
#' @return A nested list with one model specification per list entry.
#' @keywords internal
autoNAWRUmodel <- function(tsl, poss, nModels = 10) {
# trend and cycle
model <- NAWRUmodel(tsl = tsl)
trend <- trendOptim(x = model$tsl$ur, opt = poss$trend)
cycle <- cycleOptim(x = model$tsl$ur, opt = poss$cycle)
# second observation equation
typel <-list()
for (k in poss$type) {
model <- suppressWarnings(NAWRUmodel(tsl = tsl, type = k))
xexo <- switch(k,
"TKP" = tsl[poss$exoNames],
"NKP" = NULL)
typel[[k]] <- obs2Optim(x1 = model$tsl$ur,
x2 = model$tsl$pcInd,
xexo = xexo,
errorARmax = poss$errorARmax,
errorMAmax = poss$errorMAmax,
maxCycleLag = poss$maxCycleLag,
maxAR = 0,
nModels = nModels)
}
modell <- c(typel[[1]],typel[[2]])
# sort and discard models
bic <- unlist(lapply(modell, "[[", "BIC"))
index <- order(bic)[1:nModels]
modell <- modell[index]
for (k in 1:nModels) {
# create exoType input
if (!is.null(modell[[k]]$exo) && !is.na(modell[[k]]$exo)) {
exoType <- initializeExo(varNames = unique(modell[[k]]$exo$var))
for (k1 in unique(modell[[k]]$exo$var)) {
if (k1 %in% modell[[k]]$exo$var) {
exo_tmp <- modell[[k]]$exo[modell[[k]]$exo$var == k1, ]
for (k2 in 1:nrow(exo_tmp)) {
exoType[k2, k1, "difference"] <- exo_tmp$diff[k2]
exoType[k2, k1, "lag"] <- exo_tmp$lag[k2]
}
}
}
modell[[k]]$exoType <- exoType
modell[[k]]$exo <- NULL
}
# assign type
modell[[k]]$type <- rep(poss$type, each = nModels)[index][k]
}
# assign trend, cycle
nModels <- length(modell)
modell <- rep(modell, length(cycle))
for (k in 1:length(cycle)) {
for (j in 1:nModels) {
modell[[j]]$cycle <- cycle[k]
}
}
nModels <- length(modell)
modell <- rep(modell, length(trend))
for (k in 1:length(trend)) {
for (j in 1:nModels) {
modell[[j]]$trend <- trend[k]
}
}
modell <- lapply(modell, function(x) {
index <- which(names(x) == "errorARMA")
names(x)[index] <- "pcErrorARMA"
index <- which(names(x) == "ar")
x <- x[-index]
x
})
return(modell)
}
# -------------------------------------------------------------------------------------------
#' TFP model suggestion
#'
#' @description Finds the most suitable TFP models.
#'
#' @param poss A list with the characteristics of possible models (see \code{autoGapProd)}
#' @inheritParams autoGapProd
#'
#' @return A nested list with one model specification per list entry.
#' @keywords internal
autoTFPmodel <- function(tsl, poss, nModels = 10) {
# trend and cycle
model <- TFPmodel(tsl = tsl)
trend <- trendOptim(x = model$tsl$logtfp, opt = poss$trend)
cycle <- cycleOptim(x = model$tsl$logtfp, opt = poss$cycle)
# second observation equation
modell <- obs2Optim(x1 = model$tsl$logtfp,
x2 = model$tsl$cubs,
xexo = NULL,
errorARmax = poss$errorARmax,
errorMAmax = poss$errorMAmax,
maxCycleLag = poss$maxCycleLag,
maxAR = poss$cubsARmax,
nModels = nModels)
# sort by bic
bic <- unlist(lapply(modell, "[[", "BIC"))
# assign trend and cycle
nModels <- length(modell)
modell <- rep(modell, length(cycle))
for (k in 1:length(cycle)) {
for (j in 1:nModels) {
modell[[j]]$cycle <- cycle[k]
}
}
nModels <- length(modell)
modell <- rep(modell, length(trend))
for (k in 1:length(trend)) {
for (j in 1:nModels) {
modell[[j]]$trend <- trend[k]
}
}
modell <- lapply(modell, function(x) {
index <- which(names(x) == "errorARMA")
names(x)[index] <- "cubsErrorARMA"
index <- which(names(x) == "ar")
names(x)[index] <- "cubsAR"
index <- which(names(x) == "exo")
x <- x[-index]
x
})
return(modell)
}
# -------------------------------------------------------------------------------------------
#' model selection helper function
#'
#' @description Defines and fits models given the input parameters
#'
#' @param FUNmodel The model definition function.
#' @param FUNmodel The model fitting function.
#' @param poss A list with the characteristics of possible models (see \code{autoGapProd)}
#' @param comb A nested list with one model specification per list entry.
#' @param modelName Name of the model, i.e., NAWRU or TFP.
#' @inheritParams autoGapProd
#'
#' @return A nested list with the models and fitted objects.
#' @keywords internal
helper_model_fit <- function(FUNmodel, FUNfit, tsl, comb, poss, method, type, q, modelName) {
n_poss <- length(comb)
# define models
model <- f <- fBayes <- list()
valid <- rep(NA, n_poss)
message(paste0("\nDefining ", n_poss," ", modelName, " models ..."))
pb <- utils::txtProgressBar(min = 0, max = n_poss, style = 3)
for (k in 1:n_poss) {
if (k %% (n_poss/1) == 0) {
utils::setTxtProgressBar(pb, k)
}
model_tmp <- tryCatch(
{
tmp <- suppressWarnings(do.call(FUNmodel, c(list(tsl = tsl), comb[[k]])))
tmp
},
error=function(cond) {
return(NULL)
}
)
model[[k]] <- model_tmp
valid[k] <- !is.null(model_tmp)
}
model <- model[valid]
comb <- comb[valid]
n_poss <- length(comb)
# find and delete duplicates
valid <- rep(NA, n_poss)
valid[n_poss] <- TRUE
for (k in 1:(n_poss-1)) {
valid_tmp <- NULL
for (k2 in (k+1):n_poss) {
valid_tmp[k2-k] <- tryCatch(
{
!all(unlist(attributes(model[[k]])) == unlist(attributes(model[[k2]])))
},
error=function(cond) {
return(TRUE)
},
warning=function(cond) {
return(TRUE)
}
)
}
valid_tmp
valid[k] <- all(valid_tmp)
}
model <- model[valid]
comb <- comb[valid]
n_poss <- length(comb)
message(paste0("\n", n_poss," valid ", modelName, " models remain."))
# fit models
message(paste0("\nFitting ", n_poss," ", modelName, " models ..."))
pb <- utils::txtProgressBar(min = 0, max = n_poss, style = 3)
for (k in 1:n_poss) {
if (k %% 1 == 0) {
utils::setTxtProgressBar(pb, k)
cat("\n")
}
parRestr <- initializeRestr(model = model[[k]], type = type, q = q)
f[[k]] <- FUNfit(model = model[[k]],
parRestr = parRestr,
signalToNoise = poss$signalToNoise,
control = list(maxit = 1000))
if (method == "bayesian") {
prior <- initializePrior(model = model[[k]])
fBayes[[k]] <- tryCatch(
{ FUNfit(model = model[[k]],
parRestr = parRestr,
prior = prior,
signalToNoise = poss$signalToNoise,
method = method,
R = 1000,
MLEfit = f[[k]])
},
error=function(cond) {
return(NA)
}
)
}
}
res <- list(model = model,
fit = f,
fitBayes = fBayes)
return(res)
}
# -------------------------------------------------------------------------------------------
#' model selection fit comparison helper function
#'
#' @description Defines and fits models given the input parameters
#'
#' @param fit Fitted object.
#' @param E1name Name of first observation equation.
#' @param E1Trendname Name of trend of first observation equation.
#' @param E1trans Transformation function for the first observation equation..
#' @param fitBayes Fitted Bayesian object.
#'
#' @return A data frame with information criteria, other goodness-of-fit measures,
#' convergence status, trend volatility measures.
#' @keywords internal
helper_fit_comparison <- function(fit, E1name, E1Trendname, E1trans = identity, fitBayes) {
n_poss <- length(fit)
# fit criteria
crit <- names(fit[[1]]$fit)
crit <- c(crit["LjungBox" != crit])
info <- data.frame(matrix(NA, length(fit), length(crit)))
colnames(info) <- crit
for (k in 1:(length(crit))) {
info[, k] <- unlist(lapply(lapply(fit, "[[", "fit"), "[[", crit[k]))
}
info$LjungBox <- unlist(lapply(lapply(lapply(fit, "[[", "fit"), "[[", "LjungBox"), "[[","p.value"))
info$convergence <- unlist(lapply(lapply(lapply(fit, "[[", "SSMfit"), "[[", "optim.out"), "[[", "convergence"))
# trend volatility measures
df <- data.frame()
for (k in 1:n_poss) {
nDiff <- switch(attr(fit[[k]]$model, "trend"),
"RW1" = 1,
"RW2" = 2,
"DT" = 1)
tmp <- trendVolaMeasures(tsOriginal = E1trans(fit[[k]]$model$tsl[[E1name]]),
tsTrend = fit[[k]]$tsl[[E1Trendname]],
nDiff = nDiff)
tmp <- as.data.frame(tmp)
df <- rbind(df, tmp)
}
info <- cbind(info, df)
info_mle <- info
info <- NULL
if (length(fitBayes)!=0) {
nFit <- length(fitBayes)
# index_na <- c(1:length(fitBayes))[unlist(lapply(fitBayes, function(x) all(is.na(x))))]
index_na_invers <- c(1:length(fitBayes))[unlist(lapply(fitBayes, function(x) !all(is.na(x))))]
fitBayes <- fitBayes[unlist(lapply(fitBayes, function(x) !all(is.na(x))))]
# fit criteria
crit <- names(fitBayes[[1]]$fit)
info <- data.frame(matrix(NA, length(fitBayes), length(crit)))
colnames(info) <- crit
for (k in 1:(length(crit))) {
info[, k] <- unlist(lapply(lapply(fitBayes, "[[", "fit"), "[[", crit[k]))
}
df <- data.frame()
for (k in 1:length(fitBayes)) {
nDiff <- switch(attr(fitBayes[[k]]$model, "trend"),
"RW1" = 1,
"RW2" = 2,
"DT" = 1)
tmp <- trendVolaMeasures(tsOriginal = fitBayes[[k]]$model$tsl[[E1name]],
tsTrend = fitBayes[[k]]$tsl[[E1Trendname]],
nDiff = nDiff)
tmp <- as.data.frame(tmp)
df <- rbind(df, tmp)
}
info <- cbind(info, df)
# add NA rows for error runs.
info_na <- data.frame(matrix(NA, nFit, ncol(info)))
info_na[index_na_invers, ] <- info
colnames(info_na) <- colnames(info)
info <- info_na
}
res <- list(info = info_mle,
infoBayes = info)
return(res)
}
# -------------------------------------------------------------------------------------------
#' model comparison helper function
#'
#' @description Gets model attributes
#'
#' @param models A list of models.
#'
#' @return A data frame with model attributes
#' @keywords internal
helper_model_comparison <- function(models) {
model_info <- lapply(models, function(x) {
att <- attributes(x)
xclass <- att$class[1]
E2name <- ifelse(xclass == "NAWRUmodel", "phillips curve", "cubs")
df <- data.frame(cycle = att$cycle,
trend = att$trend,
E2cycleLag = paste0(att[[E2name]]$cycleLag, collapse = ","),
E2exoVariables = paste0(att[[E2name]]$exoVariables, collapse = ", "),
E2errorAR = att[[E2name]]$errorARMA[1],
E2errorMA = att[[E2name]]$errorARMA[2])
if (xclass == "NAWRUmodel") {
df$E2type <- att[[E2name]]$type
}
df
})
do.call(rbind, model_info)
}
# -------------------------------------------------------------------------------------------
#' Trend volatility measures
#'
#' @description Computes trend volatility measures.
#'
#' @param tsOriginal The original time series.
#' @param tsTrend The trend time series.
#' @param nDiff Integer indicating the order of differencing applied to the input series.
#'
#' @return A list containing the different measures.
#' @importFrom stats cor
#' @keywords internal
trendVolaMeasures <- function(tsOriginal, tsTrend, nDiff) {
result <- list()
# negative values
result$neg <- 1 * any(tsTrend < 0)
# relative range
result$rrange <- (max(tsTrend) - min(tsTrend)) / (max(tsOriginal) - min(tsOriginal))
# transform to stationary series
tsOriginal <- diff(tsOriginal, differences = nDiff)
tsTrend <- diff(tsTrend, differences = nDiff)
# rev: relative excess volatility (result close of above 1 -> excess vola)
result$rev <- (max(tsTrend) - min(tsTrend)) / (max(tsOriginal) - min(tsOriginal))
# rstd: relative standard deviations
result$rsd <- sqrt(var(tsTrend, na.rm = TRUE)) / sqrt(var(tsOriginal, na.rm = TRUE))
# corr: correlation
result$corr <- cor(tsTrend, tsOriginal)
# normalized mean absolute deviation of trend growth (diff) (the smaller the smoother the trend)
result$msdtg <- mean(abs(tsTrend - mean(tsTrend))) / mean(abs(tsTrend))
# magtg: mean absolute growth of trend growth (diff)
result$magtg <- mean(abs(growth(tsTrend)))
result
}
# -------------------------------------------------------------------------------------------
#' Find suitable trend specification
#'
#' @description Finds the most suitable trend model according to the BIC.
#'
#' @param x A time series.
#' @param opt A character vector with the trend models to be tested. The default is
#' \code{opt = c("RW1", "RW2", "DT")}.
#'
#' @return A character string with the chosen trend model.
#' @importFrom stats BIC arima frequency
#' @importFrom zoo na.trim
#' @keywords internal
trendOptim <- function(x, opt = c("RW1", "RW2", "DT")) {
x <- na.trim(x)
freq <- frequency(x)
lambda <- 1600 / ((4 / freq)^4)
trend <- hpfilter(x, lambda = lambda)
bic <- NULL
if ("RW1" %in% opt) {
bic["RW1"] <- BIC(arima(x = diff(trend), order = c(0, 0, 0)))
}
if ("RW2" %in% opt) {
bic["RW2"] <- BIC(arima(x = diff(trend, differences = 2), order = c(0, 0, 0)))
}
if ("DT" %in% opt) {
bic["DT"] <- BIC(arima(x = diff(trend), order = c(1, 0, 0)))
}
res <- which(min(bic) == bic)
return(names(res))
}
# -------------------------------------------------------------------------------------------
#' Find suitable cycle specification
#'
#' @description Finds the most suitable cycle model according to the BIC.
#'
#' @param x A time series.
#' @param opt A character vector with the cycle models to be tested. The default is
#' \code{opt = c("AR1", "AR2", "RAR2")}.
#'
#' @return A character string with the chosen cycle model.
#' @importFrom stats BIC arima frequency
#' @importFrom zoo na.trim
#' @keywords internal
cycleOptim <- function(x, opt = c("AR1", "AR2", "RAR2")) {
x <- na.trim(x)
freq <- frequency(x)
lambda <- 1600 / ((4 / freq)^4)
trend <- hpfilter(x, lambda = lambda)
cycle <- x - trend
bic <- NULL
if ("AR1" %in% opt) {
bic["AR1"] <- BIC(arima(cycle, order = c(1, 0, 0), method = "ML", include.mean = FALSE))
}
if ("AR2" %in% opt | "RAR2" %in% opt) {
bic["AR2"] <- BIC(arima(cycle, order = c(2, 0, 0), method = "ML", include.mean = FALSE))
}
res <- which(min(bic) == bic)
if (names(res) == "AR2" & "RAR2" %in% opt) {
res <- c("AR2", "RAR2")
} else {
res <- names(res)
}
return(res)
}
# -------------------------------------------------------------------------------------------
#' Find suitable 2nd observation specification
#'
#' @description Finds the most suitable model for the second observation equation according
#' to the BIC.
#'
#' @param x1 A time series, the first observation equation.
#' @param x1 A time series, the second observation equation.
#' @param xexo (Optional) A (multiple) time series with exogenous variables.
#' @param errorARmax Integer, maximal AR order of the error process of the 2nd observation
#' equation.
#' @param errorMAmax Integer, maximal MA order of the error process of the 2nd observation
#' equation.
#' @param maxCycleLag Integer, maximal cycle lag included in the 2nd observation
#' equation.
#' @param maxAR Integer, maximal AR order of the time series \code{x2} in the 2nd observation
#' equation. \code{0} means that no lag is included.
#' @param nModels Integer, maximum number of models chosen to be fitted.
#'
#' @return A list containing the chosen parameters: \code{errorAR, errorMA, cycleLag, ar, exo}.
#' @keywords internal
#' @importFrom stats BIC predict arima frequency window lag
obs2Optim <- function(x1, x2, xexo = NULL, errorARmax = 2, errorMAmax = 2, maxCycleLag = 2, maxAR = 2, nModels = 1) {
maxExo <- 10
freq <- frequency(x1)
lambda <- 1600 / ((4 / freq)^4)
trend <- hpfilter(x1, lambda = lambda)
cycle <- list()
cycle[1:(maxCycleLag + 1)] <- lapply(0:maxCycleLag, function(x) {
window(stats::lag(x1 - trend, -x), end = end(x1 - trend))
})
xregAR <- regExo <- list()
if (maxAR > 0) {
xregAR <- list()
xregAR[1:(maxAR)] <- lapply(1:maxAR, function(x) {
window(stats::lag(x2, -x), start = start(x2), end = end(x1 - trend))
})
}
xreg_base <- cycle[[1]]
bic_base <- BIC(arima(x2, order = c(0, 0, 0), xreg = xreg_base, method="ML"))
# xexo pre selection
tmp <- NULL
bic <- NULL
par_max <- length(x2) - 2 - 1
exo_lag_max <- min(ifelse(freq == 1, 2, 8), par_max)
if (!is.null(xexo)) {
count <- 0
for (k in 1:length(xexo)) {
for (p1 in 0:exo_lag_max) { # lag
for (p2 in 1:2) { # diff
count <- count + 1
xregExo <- cbind(xreg_base, window(diff(stats::lag(xexo[[k]], -p1), differences = p2), start = start(x2), end = end(x2)))
mod <- arima(x2, order = c(0, 0, 0), xreg = xregExo, method="ML")
bic[count] <- BIC(mod)
}
}
}
tmp <- cbind(bic, expand.grid(lapply(list( 1:2, 0:exo_lag_max, 1:length(xexo)), function(x) 1:length(x))))
tmp <- tmp[tmp[,1] < bic_base, ]
tmp <- tmp[order(tmp[, 1]), ][1:min(maxExo, nrow(tmp)), 2:4]
names(tmp) <- c("diff", "lag", "var")
tmp$lag <- tmp$lag-1
tmp$var <- names(xexo)[tmp$var]
tmp
# only 2 lags/diff combs, otherwise there can be linearly dependent combinations
tmpl <- list()
for (k in unique(tmp$var)) {
tmpl[[k]] <- tmp[tmp$var==k, ]
tmpl[[k]] <- tmpl[[k]][1:min(nrow(tmpl[[k]]), 2), ]
}
tmp <- do.call(rbind, tmpl)
tmp
regExo <- list()
for (k in 1:nrow(tmp)) {
regExo[[k]] <- window(diff(lag(xexo[[tmp$var[k]]], -tmp$lag[k]), differences = tmp$diff[k]), start = start(x2), end = end(x2))
}
}
# ARMA errors
par_cycle_max <- floor(length(x2) * 3 / 4) - nrow(tmp) - errorARmax - errorMAmax - maxAR
poss <- list(cycle = 0:(length(cycle)-1),
xregAR = 0,
AR = 0:errorARmax,
MA = 0:errorMAmax)
if (maxAR > 0) {
poss$xregAR = 0:maxAR
}
if (!is.null(xexo)) {
poss <- c(poss, lapply(1:length(regExo), function(x) 1:2))
names(poss)[length(poss) - (length(regExo)-1):0] <- paste0("exo", 1:length(regExo))
}
combinations <- expand.grid(lapply(poss, function(x) 1:length(x)))
comb <- as.data.frame(sapply(1:length(poss), function(x) poss[[x]][combinations[,x]]))
colnames(comb) <- colnames(combinations)
bic <- ljungp <- mae_ratio <- rep(NA, nrow(comb))
h <- min(10, length(x2) / 5) # see https://robjhyndman.com/hyndsight/ljung-box-test/
for (p in 1:nrow(comb)) {
pAR <- comb$AR[p]
pMA <- comb$MA[p]
xreg <- cbind(do.call(cbind, cycle[1:(comb$cycle[p] + 1)]),
do.call(cbind, regExo[comb[p, grepl("exo", names(poss))] == 1]))
if (comb$xregAR[p] > 0) {
xreg <- cbind(xreg, do.call(cbind, xregAR[1:comb$xregAR[p]]))
}
mod <- arima(x2, order = c(pAR, 0, pMA), xreg = xreg, method="ML")
bic[p] <- BIC(mod)
ljungp[p] <- Box.test(mod$residuals, type = "Ljung-Box", lag = h)$p.value
# out of sample performance
n_par <- ncol(xreg) + comb$AR[p] + comb$MA[p] + 1
index <- 1:max(floor(length(x2) * 3 / 4), n_par + 1)
mae_ratio[p] <- tryCatch(
{
mod <- arima(x2[index], order = c(comb$AR[p], 0, comb$MA[p]), xreg = as.matrix(xreg)[index,], method="ML")
pred <- stats::predict(mod, newxreg = as.matrix(xreg)[-index,])
mae_in <- 100 / length(index) * sum( abs( mod$residuals ), na.rm = TRUE )
mae_out <- 100 / length(pred$pred) * sum( abs( x2[-index] - pred$pred ), na.rm = TRUE )
mae_out / mae_in
},
error=function(cond) {
return(Inf)
},
warning=function(cond) {
return(Inf)
}
)
}
if (sum(ljungp < 0.1) < nrow(comb)) { bic[ljungp < 0.1] <- Inf }
if (sum(mae_ratio > 1.5 | ljungp < 0.1) < nrow(comb)) { bic[mae_ratio > 1.5] <- Inf }
nModels <- min(nModels, length(bic))
index <- order(bic)[1:nModels]
bic_min <- bic[index]
errorAR <- comb$AR[index]
errorMA <- comb$MA[index]
cycleLag <- comb$cycle[index]
ar <- comb$xregAR[index]
exo <- as.list(rep(NA, nModels))
if (!is.null(xexo)) {
exo <- list()
for (k in 1:nModels)
exo[[k]] <- tmp[comb[index[k], grepl("exo", names(poss))] == 1, ]
}
res <- list(
errorARMA = as.list(as.data.frame(rbind(errorAR, errorMA))),
cycleLag = lapply(cycleLag, function(x) 0:x),
ar = as.list(ar),
exo = exo,
BIC = as.list(bic_min)
)
res <- lapply(1:nModels, function(x) lapply(res, "[[", x))
return(res)
}
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.