Nothing
# -------------------------------------------------------------------------------------------
#' Estimation of a \code{NAWRUmodel}
#'
#' @description Estimates a two-dimensional state-space model and performs filtering and
#' smoothing to obtain the NAWRU using either maximum likelihood estimation or bayesian
#' methods.
#'
#' @param model An object of class NAWRUmodel.
#' @param parRestr A list of matrices containing the parameter restrictions for the cycle,
#' trend, and the Phillip's curve. Each matrix contains the lower and upper bound of the
#' involved parameters. \code{NA} implies that no restriction is present. Autoregressive
#' parameters are automatically restricted to the stationary region unless box constraints
#' are specified. By default, \code{parRestr} is initialized by the function
#' \code{initializeRestr(model)}. Only used if \code{method = "MLE"}.
#' @param signalToNoise (Optional) signal to noise ratio. Only used if \code{method = "MLE"}.
#' @param method The estimation method. Options are maximum likelihood estimation \code{"MLE"}
#' and bayesian estimation \code{"bayesian"}. The default is \code{method = "MLE"}.
#' @param prior A list of matrices with parameters for the prior distribution and box
#' constraints. By default, \code{prior} is initialized by \code{initializePrior(model)}.
#' See details. Only used if \code{method = "bayesian"}.
#' @param R An integer specifying the number of MCMC draws. The default is \code{R = 10000}.
#' Only used if \code{method = "bayesian"}.
#' @param burnin An integer specifying the burn-in phase of the MCMC chain. The default is
#' \code{burnin = ceiling(R / 10)}. Only used if \code{method = "bayesian"}.
#' @param thin An integer specifying the thinning interval between consecutive draws. The
#' default is \code{thin = 1}, implying that no draws are dopped. For \code{thin = 2},
#' every second draw is dropped and so on. Only used if \code{method = "bayesian"}.
#' @param HPDIprob A numeric in the interval \code{(0,1)} specifying the target probability
#' of the highest posterior density intervals. The default is \code{HPDIprob = 0.9}. Only
#' used if \code{method = "bayesian"}.
#' @param pointEstimate Posterior distribution's statistic of central tendency. Possible
#' options are \code{"mean"} and \code{"median"}. The default is \code{pointEstimate = "mean"}.
#' Only used if \code{method = "bayesian"}.
#' @param MLEfit (Optional) An object of class \code{NAWRUfit} which is used for
#' initialization. Only used if \code{method = "bayesian"}.
#' @param control (Optional) A list of control arguments to be passed on to \code{optim}.
#' @param ... additional arguments to be passed to the methods functions.
#'
#' @details The list object \code{prior} contains three list elements \code{cycle},
#' \code{trend}, and \code{pcInd}. Each list element is a \code{4 x n} matrix where \code{n}
#' denotes the number of parameters involved in the respective equation. The upper two
#' elements specify the distribution, the lower two parameters specify box constraints.
#' \code{NA} denotes no constraints. Autoregressive parameters are automatically restricted
#' to the stationary region unless box constraints are specified. For instance,
#' \code{prior$cycle[, 1]} contains the mean, standard deviation, lower and upper bound for
#' the first variable, in that respective order.
#' @details The respective prior distributions are defined through their mean and standard
#' deviation.
#' @details The Gibbs sampling procedure is as follows. For each \eqn{r = 1, ..., R}
#' \itemize{\item The states are sampled by running the Kalman filter and smoother
#' conditional on the parameters of the previous step, \eqn{\theta_{r-1}}
#' \item Trend equation parameters \eqn{\theta_{trend}}: Conditional on the
#' states \eqn{\alpha_r}, a draw \eqn{\theta_{trend,k}} is obtained either
#' by a sequential Gibbs step, a Metropolis Hasting step, or by conjugacy,
#' depending on the trend model specification.
#' \item Cycle equation parameters \eqn{\theta_{cycle}}: Conditional on the
#' states \eqn{\alpha_r}, a draw \eqn{\theta_{cycle,k}} is obtained either
#' by a sequential Gibbs step, a Metropolis Hasting step, or by conjugacy,
#' depending on the cycle model specification.
#' \item Phillip's curve equation parameters \eqn{\theta_{pcInd}}: Conditional on the
#' states \eqn{\alpha_r}, a draw \eqn{\theta_{pcInd,k}} is obtained either
#' by a sequential Gibbs step, a Metropolis Hasting step, a combination
#' thereof, or by conjugacy, depending on the Phillip's curve equation specification.}
#'
#' @return For maximum likelihood estimation, an object of class \code{NAWRUfit} containing
#' the following components:
#' \item{model}{The input object of class \code{NAWRUmodel}.}
#' \item{SSMfit}{The estimation output from the function \code{fitSSM} from \code{KFAS}.}
#' \item{SSMout}{The filtering and smoothing output from the function \code{KFS} from
#' \code{KFAS}.}
#' \item{parameters}{A data frame containing the estimated parameters, including standard
#' errors, t-statistics, and p-values.}
#' \item{parRestr}{A list of matrices containing the enforced parameter constraints.}
#' \item{fit}{A list of model fit criteria (see below).}
#' \item{call}{Original call to the function. }
#' The list component \code{fit} contains the following model fit criteria:
#' \item{loglik}{Log-likelihood function values.}
#' \item{AIC}{Akaike information criterion.}
#' \item{BIC}{Bayesian information criterion.}
#' \item{AIC}{Hannan-Quinn information criterion.}
#' \item{RMSE}{root mean squared error of the Phillip's curve equation.}
#' \item{R2}{R squared of the Phillip's curve equation.}
#' \item{LjungBox}{Ljung-Box test output of the Phillip's curve equation.}
#' \item{signal-to-noise}{Signal-to-noise ratio.}
#' For bayesian estimation, an object of class \code{NAWRUfit} containing the following
#' components:
#' \item{model}{The input object of class \code{NAWRUmodel}.}
#' \item{tsl}{A list of time series containing the estimated states.}
#' \item{parameters}{A data frame containing the estimated parameters, including standard
#' errors, highest posterior density credible sets.}
#' \item{prior}{A list of matrices containing the used prior distributions.}
#' \item{fit}{A list of model fit criteria (see below).}
#' \item{call}{Original call to the function. }
#' The list component \code{fit} contains the following model fit criteria:
#' \item{R2}{R squared of the phillips curve equation,}
#' \item{signal-to-noise}{Signal-to-noise ratio.}
#'
#' @export
#' @family fitting methods
#' @examples
#' # define nawru model for France
#' data("gap")
#' country <- "France"
#' tsList <- amecoData2input(gap[[country]])
#' model <- NAWRUmodel(tsl = tsList)
#' # estimate nawru model via MLE
#' parRestr <- initializeRestr(model = model, type = "hp")
#' \donttest{
#' f <- fit(model = model, parRestr = parRestr)
#' }
#' # initialize priors and estimate model via Bayesian methods
#' prior <- initializePrior(model = model)
#' \donttest{
#' f <- fit(model = model, method = "bayesian", prior = prior, R = 5000, thin = 2)
#' }
fit.NAWRUmodel <- function(model, parRestr = initializeRestr(model = model), signalToNoise = NULL,
method = "MLE", control = NULL,
prior = initializePrior(model), R = 10000, burnin = ceiling(R / 10),
thin = 1, HPDIprob = 0.85, pointEstimate = "mean", MLEfit = NULL, ...) {
# save call
mc <- match.call(expand.dots = FALSE)
if (method == "MLE") {
f <- .MLEfitNAWRU(
model = model, parRestr = parRestr, signalToNoise = signalToNoise,
control = control
)
} else if (method == "bayesian") {
if (!is.null(signalToNoise)) {
warning("'signalToNoise' not applicable for method = 'bayesian'.")
}
f <- .BayesFitNAWRU(
model = model, prior = prior, R = R, burnin = burnin,
thin = thin, HPDIprob = HPDIprob, FUN = pointEstimate, MLEfit = MLEfit
)
} else {
stop("Invalid estimation method: please use either 'MLE' or 'bayesian'.")
}
f$call <- mc
print(f)
f
}
# -------------------------------------------------------------------------------------------
#' NAWRU model
#'
#' @description Creates a state space object object of class \code{NAWRUmodel} which can be
#' fitted using \code{fit}.
#'
#' @param tsl A list of time series objects, see details.
#' @param trend A character string specifying the trend model. \code{trend = "RW1"} denotes
#' a first order random walk, \code{trend = "RW2"} a second order random walk (local linear
#' trend) and \code{trend = "DT"} a damped trend model. The default is \code{trend = "RW2"}.
#' @param cycle A character string specifying the cycle model. \code{cycle = "AR1"} denotes
#' an AR(1) process, \code{cycle = "AR2"} an \code{AR(2)} process. The default is
#' \code{cycle = "AR2"}.
#' @param type A character string specifying the type of the Phillip's curve.
#' \code{type = "TKP"} denotes the traditional Keynesian Phillip's curve and
#' \code{type = "NKP"} the New Keynesian Phillip's curve, see details. The default is
#' \code{type = "TKP"}.
#' @param cycleLag A vector specifying the cycle lags that are included in the Phillip's
#' curve. The default is \code{cycleLag = 0}, see details.
#' @param pcErrorARMA A \code{2 x 1} vector with non-negative integers specifying the AR
#' and MA degree of the error term in the Phillip's curve equation. The default is
#' \code{pcErrorARMA = c(0, 0)}, see details.
#' @param exoType An optional \code{n x m x 2} array specifying the possible difference
#' and lag transformation for the variables. \code{exoType} can be initialized using the
#' function \code{inizializeExo}. The column names give the variable names.
#' \code{exoType[, , 1]} contains the difference transformations and \code{exoType[, , 2]}
#' the subsequent lag transformations, see details.
#' @param start (Optional) Start vector for the estimation, e.g. \code{c(1980, 1)}.
#' @param end (Optional) End vector for the estimation, e.g. \code{c(2020, 1)}.
#' @param anchor (Optional) Anchor value for the unemployment rate.
#' @param anchor.h (Optional) Anchor horizon in the frequency of the given time series.
#'
#' @details The list of time series \code{tsl} needs to have the following components:
#' \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"}.}
#' }
#' and optionally other variables included in \code{exoType}.
#' @details A \code{cycleLag} equal to \code{0} implies that only the contemporaneous cycle
#' is included in the Phillip's curve. A \code{cycleLag} equal to \code{0:1} implies that
#' the contemporaneous as well as the lagged cycle are included.
#' @details A \code{pcErrorARMA} equal to \code{c(0, 0)} implies that the error term in the
#' Phillip's curve is white noise. \code{pcErrorARMA = c(1, 0)} implies that the error is
#' an AR(1) process and for \code{pcErrorARMA = c(1, 2)} the error follows an ARMA(1, 2)
#' process.
#' @details For the New Keynesian Phillip's curve, the \code{cycleLag} cannot be chosen.
#' \code{cycleLag} will be set to \code{0} if \code{cycle = "AR1"} and to \code{1} if
#' \code{cycle = "AR2"}. In the latter case, the forward solution of the Phillip's curve
#' implies parameter restrictions for the lagged cycle on the Phillip's curve. Moreover,
#' exogenous variables will be ignored in the case of the New Keynesian Phillip's curve.
#' @details The array \code{exoType} consists of non-negative integers or \code{NA}s.
#' \code{exoType[, , 1] = c(NA,1)} and \code{exoType[, , 2] = c(NA,2)} implies that
#' the first variable is not included in the Phillip's curve whereas the second lag of
#' the first difference of the second variable is included.
#'
#' @return Object of class \code{NAWRUmodel}, which is a list with the following components:
#' \item{tsl}{A list of used time series.}
#' \item{SSModel}{An object of class SSModel specifying the state-space model.}
#' \item{loc}{A data frame containing information on each involved parameter, for instance
#' its corresponding system matrix, variable names, and parameter restrictions.}
#' \item{call}{Original call to the function.}
#' In addition, the object contains the following attributes:
#' \item{cycle}{Cycle specification.}
#' \item{trend}{Trend specification.}
#' \item{phillipsCurve}{A list containing the components \code{type, cycleLag, errorARMA, exoVariables}.}
#' \item{anchor}{A list containing the components \code{value, horizon}.}
#' \item{period}{A list containing the components \code{start, end, frequency}.}
#'
#' @export
#' @importFrom KFAS SSModel SSMregression SSMcustom
#' @importFrom stats start end window ts lag frequency time
#' @importFrom zoo na.trim
#' @examples
#' # load data for France
#' data("gap")
#' tsList <- amecoData2input(gap$France, alpha = 0.65)
#' # Traditional phillips curve
#' model <- NAWRUmodel(tsl = tsList, trend = "RW2", cycle = "AR2", type = "TKP", cycleLag = 0)
#'
#' # New-Keynesian Phillips curve
#' model <- NAWRUmodel(tsl = tsList, trend = "RW2", cycle = "AR2", type = "NKP", cycleLag = 0:1)
#'
#' # Traditional Phillips curve with 6 exogenous variables
#' # specify exogenous variable transformations
#' D <- matrix(c(2, 2, 2, 1, 1, 1), 2, 3, byrow = TRUE)
#' L <- matrix(c(0, 0, 0, 1, 1, 1), 2, 3, byrow = TRUE)
#' exoType <- initializeExo(varNames = c("tot", "prod","ws"), D = D, L = L)
#' model <- NAWRUmodel(tsl = tsList, cycleLag = 0:1, exoType = exoType)
NAWRUmodel <- function(tsl, trend = "RW2", cycle = "AR2", type = "TKP", cycleLag = 0, pcErrorARMA = c(0, 0),
exoType = NULL, start = NULL, end = NULL, anchor = NULL, anchor.h = NULL) {
# local variables
nPar <- nObs <- nState <- nStateV <- stateNames <- loc <- sys <- NULL
# save call
mc <- match.call(expand.dots = FALSE)
# ----- check input for consistency
exoNames <- colnames(exoType)
errorARMA <- pcErrorARMA
if (length(errorARMA) == 1) {
errorARMA <- c(errorARMA, 0)
}
list2env(.checkNawru(
tsl = tsl,
trend = trend,
cycle = cycle,
type = type,
cycleLag = cycleLag,
errorARMA = errorARMA,
exoNames = exoNames,
exoType = exoType,
start = start,
end = end,
anchor = anchor,
anchor.h = anchor.h
),
envir = environment())
# ----- preprocess data
# assign end and start
if (is.null(end)) end <- end(tsl[[1]])
if (is.null(start)) start <- start(tsl[[1]])
# collect used variables
varUsed <- c("ur")
nExo <- 0
exoNamesTmp <- NULL
# merge observation equation data
if (type == "TKP") {
tsl$pcInd <- diff(tsl$winfl) * 100
tsl$pcInd <- diff(diff(log(tsl$nulc))) * 100
varUsed <- c(varUsed, "pcInd")
obsNames <- c("ur", "ddlognulc")
} else {
tsl$pcInd <- diff(log(tsl$rulc)) * 100
tsl$pcIndl <- stats::lag(tsl$pcInd, k = -1)
varUsed <- c(varUsed, "pcInd", "pcIndl")
obsNames <- c("ur", "dlogrulc")
nExo <- 1
exoNamesTmp <- c(exoNamesTmp, "pcIndl")
}
# assign and compute exogenous variables
if (any(!is.na(exoType)) && type == "TKP") {
tsl[exoNames[!(substr(exoNames, 1, 5) == "dummy")]] <- lapply(tsl[exoNames[!(substr(exoNames, 1, 5) == "dummy")]], log)
for (jj in (1:dim(exoType)[2])) {
name <- exoNames[jj]
for (ii in (1:dim(exoType)[1])) {
if (!is.na(exoType[ii, jj, 1])) {
name_tmp <- paste0(
"pc", paste(rep("d", exoType[ii, jj, 1]), collapse = ""), name,
paste(rep("l", exoType[ii, jj, 2]), collapse = "")
)
tsl[[name_tmp]] <- tsl[[name]]
try(tsl[[name_tmp]] <- diff(tsl[[name_tmp]], differences = exoType[ii, jj, 1]), silent = TRUE)
tsl[[name_tmp]] <- stats::lag(tsl[[name_tmp]], k = -abs(exoType[ii, jj, 2]))
exoNamesTmp <- c(exoNamesTmp, name_tmp)
}
}
}
nExo <- length(exoNamesTmp)
varUsed <- c(varUsed, exoNamesTmp)
}
# list of used time series
tslUsed <- tsl[varUsed]
tslUsed <- lapply(tslUsed, function(x) suppressWarnings(window(x, start = start, end = end))) # warns if start or end are not changed.
tslUsed <- lapply(tslUsed, zoo::na.trim)
# update start and end date if necessary
start <- dateTsList(tslUsed, FUN1 = stats::start, FUN2 = base::max)
end <- dateTsList(tslUsed, FUN1 = stats::end, FUN2 = base::min)
tslUsed <- lapply(tslUsed, function(x) suppressWarnings(window(x, start = start, end = end))) # warns if start or end are not changed.
tslUsed <- lapply(tslUsed, zoo::na.trim)
# ----- system matrices pre processing
list2env(.SSSystem(
tsl = tslUsed,
cycle = cycle,
trend = trend,
cycleLag = cycleLag,
type = type,
errorARMA = errorARMA
),
envir = environment())
tslUsed <- lapply(tslUsed, function(x) suppressWarnings(window(x, start = start, end = end))) # warns if start or end are not changed.
# ----- state space model
if (is.null(exoNamesTmp)) {
modelSS <- SSModel(cbind(ur, pcInd) ~ -1
+ SSMcustom(
Z = sys$Zt, T = sys$Tt, R = sys$Rt, Q = sys$Qt,
a1 = sys$a1, P1 = sys$P1, P1inf = sys$P1inf, state_names = stateNames
),
H = sys$Ht,
data = tslUsed
)
} else {
tmp <- paste0("~ ", paste(exoNamesTmp, collapse = " + "))
modelSS <- SSModel(cbind(ur, pcInd) ~ -1
+ SSMregression(as.formula(tmp), index = 2, data = tslUsed, state_names = exoNamesTmp)
+ SSMcustom(
Z = sys$Zt, T = sys$Tt, R = sys$Rt, Q = sys$Qt,
a1 = sys$a1, P1 = sys$P1, P1inf = sys$P1inf, state_names = stateNames
),
H = sys$Ht,
data = tslUsed
)
}
# ----- nawru model
model <- list(
tsl = tslUsed,
SSModel = modelSS,
call = mc
)
lPc <- list(
type = type,
cycleLag = cycleLag,
errorARMA = errorARMA,
exoVariables = exoNamesTmp
)
lAnchor <- list(
value = anchor,
horizon = anchor.h
)
lPeriod <- list(
start = paste0(start(tslUsed$ur)[1], ifelse(frequency(tslUsed$ur) == 4, paste0(" Q", start(tslUsed$ur)[2]), "")),
end = paste0(end(tslUsed$ur)[1], ifelse(frequency(tslUsed$ur) == 4, paste0(" Q", end(tslUsed$ur)[2]), "")),
frequency = ifelse(frequency(tslUsed$ur) == 4, "quarterly", "annual")
)
class(model) <- c("NAWRUmodel", "model")
attr(model, "cycle") <- cycle
attr(model, "trend") <- trend
attr(model, "phillips curve") <- lPc
attr(model, "anchor") <- lAnchor
attr(model, "period") <- lPeriod
loc <- .initializeLoc(model)
if (type != "NKP") {
loc$variableRow[loc$sysMatrix == "exo"] <- loc$varName[loc$sysMatrix == "exo"]
}
model$loc <- loc
invisible(return(model))
}
# -------------------------------------------------------------------------------------------
#' \code{NAWRUodel} object check
#'
#' @description Tests whether the input object is a valid object of class \code{NAWRUmodel}.
#'
#' @param object An object to be tested.
#' @param return.logical If \code{return.logical = FALSE} (default), an error message is printed
#' if the object is not of class \code{NAWRUmodel}. If \code{return.logical = TRUE}, a logical
#' value is returned.
#'
#' @return A logical value or nothing, depending on the value of \code{return.logical}.
#'
#' @export
#' @importFrom KFAS is.SSModel
#' @examples
#' # load data for France
#' data("gap")
#' tsList <- amecoData2input(gap$France, alpha = 0.65)
#'
#' # Traditional phillips curve
#' model <- NAWRUmodel(tsl = tsList, trend = "RW2", cycle = "AR2", type = "NKP", cycleLag = 0:1)
#' is.NAWRUmodel(model, return.logical = TRUE)
#' attr(model, "phillips curve")$cycleLag <- 0
#' is.NAWRUmodel(model, return.logical = TRUE)
is.NAWRUmodel <- function(object, return.logical = FALSE) {
cycle <- attr(object, "cycle")
trend <- attr(object, "trend")
type <- attr(object, "phillips curve")$type
cycleLag <- attr(object, "phillips curve")$cycleLag
errorARMA <- attr(object, "phillips curve")$errorARMA
exoNames <- attr(object, "phillips curve")$exoVariables
anchor <- attr(object, "anchor")$value
anchor.h <- attr(object, "anchor")$horizon
components <- c("tsl", "SSModel", "loc", "call")
trendPossibilities <- c("RW1", "RW2", "DT")
cyclePossibilities <- c("AR1", "AR2")
typePossibilities <- c("TKP", "NKP")
cycleLagPossibilities <- c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
errorARPossibilities <- c(0, 1, 2)
errorMAPossibilities <- 0
obsNames <- c("ur", "pcInd")
exoNames <- exoNames[!grepl(obsNames[2], exoNames)]
if (return.logical) {
x <- inherits(object, "NAWRUmodel") &&
all(components %in% names(object)) &&
inherits(object$SSModel, "SSModel") &&
is.SSModel(object$SSModel, return.logical = TRUE) &&
all(trend %in% trendPossibilities) &&
all(cycle %in% cyclePossibilities) &&
all(type %in% typePossibilities) &&
all(cycleLag %in% cycleLagPossibilities) &&
all(errorARMA[1] %in% errorARPossibilities) &&
all(errorARMA[2] %in% errorMAPossibilities) &&
all(obsNames %in% names(object$tsl)) &&
!(all("NKP" %in% type) && all("AR1" %in% cycle) && max(cycleLag) != 0) &&
!(all("NKP" %in% type) && all("AR2" %in% cycle) && max(cycleLag) != 1) &&
!(!is.null(anchor) && is.null(anchor.h)) &&
!(all("NKP" %in% type) && !is.null(exoNames) && length(exoNames) != 0)
x
} else {
if (!inherits(object, "NAWRUmodel")) {
stop("Object is not of class 'NAWRUmodel'")
}
if (!all(components %in% names(object))) {
stop(paste0(
"Model is not a proper object of class 'NAWRUmodel'.",
"The following components are missing: ", paste0(components[!(components %in% names(object))], collapse = ", ")
))
}
if (!inherits(object$SSModel, "SSModel")) {
stop("Object component 'SSModel' is not of class 'SSModel'")
}
if (!is.SSModel(object$SSModel, return.logical = TRUE)) {
stop("Object component 'SSModel' is not a proper object of class 'SSModel'.")
}
if (!all(trend %in% trendPossibilities)) {
stop(paste0(
"Invalid trend specification. ",
"Valid trends: ", paste0(trendPossibilities, collapse = ", ")
))
}
if (!all(cycle %in% cyclePossibilities)) {
stop(paste0(
"Invalid cycle specification. ",
"Valid cycles: ", paste0(cyclePossibilities, collapse = ", ")
))
}
if (!all(type %in% typePossibilities)) {
stop(paste0(
"Invalid Phillips curve type specification. ",
"Valid types: ", paste0(typePossibilities, collapse = ", ")
))
}
if (!all(cycleLag %in% cycleLagPossibilities)) {
stop(paste0(
"Invalid cycleLag specification. ",
"Valid cycleLags: ", paste0(cycleLagPossibilities, collapse = ",")
))
}
if (!all(errorARMA[1] %in% errorARPossibilities)) {
stop(paste0(
"Invalid error AR specification. ",
"Valid error AR specification: ", paste0(errorARPossibilities, collapse = ",")
))
}
if (!all(errorARMA[2] %in% errorMAPossibilities)) {
stop(paste0(
"Invalid error MA specification. ",
"Valid error MA specification: ", paste0(errorMAPossibilities, collapse = ",")
))
}
if (!all(obsNames %in% names(object$tsl))) {
stop(paste0("Object component 'tsl' does not contain one or both of ", paste0(obsNames, collapse = ", ")))
}
if (all("NKP" %in% type) && all("AR1" %in% cycle) && max(cycleLag) != 0) {
stop(paste0(
"Invalid cycleLag for the New Keynesian Phillip's curve and an AR(1) cycle. ",
"Maximum cycleLag must be set to 0."
))
}
if (all("NKP" %in% type) && all("AR2" %in% cycle) && max(cycleLag) != 1) {
stop(paste0(
"Invalid cycleLag for the New Keynesian Phillip's curve and an AR(2) cycle. ",
"Maxmimum cycleLag must be set to 1."
))
}
if (!is.null(anchor) && is.null(anchor.h)) {
stop("Nawru anchor specified but no anchor horizon specified. ")
}
if (all("NKP" %in% type) && !is.null(exoNames) && length(exoNames) != 0) {
stop("Exogenous variables are not compatible with the NeW Keynesian Phillip's curve.")
}
}
}
# -------------------------------------------------------------------------------------------
#' \code{NAWRUfit} object check
#'
#' @description Tests whether the input object is a valid object of class \code{NAWRUfit}.
#'
#' @param object An object to be tested.
#' @param return.logical If \code{return.logical = FALSE} (default), an error message is printed
#' if the object is not of class \code{NAWRUfit}. If \code{return.logical = TRUE}, a logical
#' value is returned.
#'
#' @return A logical value or nothing, depending on the value of \code{return.logical}.
#'
#' @keywords internal
is.NAWRUfit <- function(object, return.logical = FALSE) {
components <- list()
components$MLE <- c("model", "tsl", "SSMfit", "SSMout", "parameters", "parRestr", "fit", "call")
components$bayesian <- c("model", "tsl", "parameters", "parametersSampled",
"statesSampled", "mcmc", "prior", "fit", "MLE", "call")
type <- attr(object, "method")
if (return.logical) {
x <- inherits(object, "NAWRUfit") &&
all(components[[type]] %in% names(object))
x
} else {
if (!inherits(object, "NAWRUfit")) {
stop("Object is not of class 'NAWRUfit'")
}
if (!all(components[[type]] %in% names(object))) {
stop(paste0(
"Model is not a proper object of class 'NAWRUfit'.",
"The following components are missing: ", paste0(components[!(components %in% names(object))], collapse = ", ")
))
}
}
}
# -------------------------------------------------------------------------------------------
#' Print \code{NAWRUmodel} object
#'
#' @description Prints the model specifications of an object of class \code{NAWRUmodel}.
#'
#' @param x An object of class \code{NAWRUmodel}.
#' @param call A logical. If \code{TRUE}, the call will be printed.
#' @param check A logical. If \code{TRUE}, the model class will be checked.
#' @param ... Ignored.
#' @return No return value, model information is printed.
#' @export
print.NAWRUmodel <- function(x, call = TRUE, check = TRUE, ...) {
.printSSModel(x = x, call = call, check = check)
}
# -------------------------------------------------------------------------------------------
#' Print \code{NAWRUfit} object
#'
#' Prints the model specifications and the estimation results of an object of class \code{NAWRUfit}.
#'
#' @param x An object of class \code{NAWRUfit}.
#' @param ... Ignored.
#' @return No return value, results are printed.
#' @export
print.NAWRUfit <- function(x, ...) {
.printSSModelFit(x = x, call = TRUE, check = FALSE, print.model = TRUE)
}
# -------------------------------------------------------------------------------------------
#' Plots for a \code{NAWRUfit} object
#'
#' @description Plots the NAWRU and the Phillip's curve and gives diagnostic plots based on
#' standardized residuals for objects of class \code{NAWRUfit}.
#'
#' @param x An object of class \code{NAWRUfit}.
#' @param alpha The significance level for the NAWRU (\code{alpha in [0,1]}). Only used if
#' \code{bounds = TRUE}.
#' @param bounds A logical indicating whether significance intervals should be plotted around
#' the nawru. The default is \code{bounds = TRUE}.
#' @param combine A logical indicating whether the diagnostic plots should be combined or not,
#' the default is \code{TRUE}.
#' @param posterior A logical indicating whether posterior diagnostics should be plotted. The
#' default is \code{FALSE}. Only applied in the case of bayesian estimation.
#' @inheritParams plot.gap
#'
#' @return No return value, plots are printed.
#'
#' @export
plot.NAWRUfit <- function(x, alpha = 0.05, bounds = TRUE, path = NULL, combine = TRUE, prefix = NULL,
posterior = FALSE, device = "png", width = 10, height = 3, ...) {
if (!is.NAWRUfit(x, return.logical = TRUE)) {
stop("x is no valid object of class 'NAWRUfit'.")
}
if (!is.null(path)) {
check <- dir.exists(paths = path)
if (!check) {
warning(paste0("The given file path '", path, "' does not exist. The plots will not be saved."))
path <- NULL
}
}
method <- attr(x, "method")
# check whether prediction is present
prediction <- "prediction" %in% names(attributes(x))
if (posterior == FALSE) {
if (!prediction) {
# ----- estimation
if (method == "MLE") {
# confidence bounds
tvalue <- -qnorm((alpha) / 2)
boundName <- paste0(100 * (1 - alpha), "% CI")
tslBounds <- list(
ub = (x$tsl$nawru + x$tsl$nawruSE * tvalue),
lb = (x$tsl$nawru - x$tsl$nawruSE * tvalue),
ub2 = (x$tsl$obsFitted[, 2] + x$tsl$obsFittedSE[, 2] * tvalue),
lb2 = (x$tsl$obsFitted[, 2] - x$tsl$obsFittedSE[, 2] * tvalue) )
# residuals
res <- x$tsl$obsResidualsRecursive[, "pcInd"]
} else { # Bayesian
# confidence bounds
HPDI <- attr(x, "HPDIprob")
boundName <- paste0(100 * HPDI, "% HPDCI")
tslBounds <- list(
ub = x$tsl$nawruSummary[, paste0(100 * HPDI, "% HPDI-UB")],
lb = x$tsl$nawruSummary[, paste0(100 * HPDI, "% HPDI-LB")],
ub2 = x$tsl$pcIndFittedSummary[, paste0(100 * HPDI, "% HPDI-UB")],
lb2 = x$tsl$pcIndFittedSummary[, paste0(100 * HPDI, "% HPDI-LB")]
)
# residuals
res <- NULL
}
# --- data
# nawru
tsl1 <- list(
trend = x$tsl$nawru,
# orig = x$model$tsl$ur,
orig = x$tsl$obs[, 1],
lb = tslBounds$lb,
ub = tslBounds$ub
)
if (!is.null(x$tsl$nawruAnchored)) {
tsl1 <- c(tsl1, list(anchor = x$tsl$nawruAnchored))
}
tsl1 <- do.call(cbind, tsl1)
# phillips curve
tsl2 <- do.call(cbind, list(
fitted = x$tsl$obsFitted[, 2],
# E2 = x$model$tsl$pcInd,
E2 = x$tsl$obs[, 2],
lb = tslBounds$lb2,
ub = tslBounds$ub2
))
# combine lists
tsl <- list(tsl1, tsl2)
# --- legends and titles and print names
legend <- list(
c("nawru", "unemployment rate", "anchored nawru"),
c("fitted", "Phillips curve indicator")
)
title <- list(
"Unemployment rate in %",
"Phillips curve",
"Phillips curve residuals"
)
namesPrint <- c("nawru", "phillips_curve")
if (!is.null(prefix)) namesPrint <- paste(prefix, namesPrint , sep = "_")
# plot
plotSSresults(
tsl = tsl, legend = legend, title = title,
boundName = boundName, res = res, namesPrint = namesPrint,
bounds = bounds, combine = combine, path = path, device = device,
width = width, height = height
)
} else {
# ----- prediction
n.ahead <- attr(x, "prediction")$n.ahead
if (method == "MLE") {
# confidence bounds
tvalue <- -qnorm((alpha) / 2)
boundName <- paste0(100 * (1 - alpha), "% CI")
tslBounds <- list(
lb = (x$tsl$nawru - x$tsl$nawruSE * tvalue),
ub = (x$tsl$nawru + x$tsl$nawruSE * tvalue),
lb1 = (x$tsl$obs[, 1] - x$tsl$obsSE[, 1] * tvalue),
ub1 = (x$tsl$obs[, 1] + x$tsl$obsSE[, 1] * tvalue),
lb2 = (x$tsl$obs[, 2] - x$tsl$obsSE[, 2] * tvalue),
ub2 = (x$tsl$obs[, 2] + x$tsl$obsSE[, 2] * tvalue),
lb3 = (x$tsl$stateSmoothed[, "cycle"] - x$tsl$stateSmoothedSE[, "cycle"] * tvalue),
ub3 = (x$tsl$stateSmoothed[, "cycle"] + x$tsl$stateSmoothedSE[, "cycle"] * tvalue) )
} else { # Bayesian
# confidence bounds
HPDI <- attr(x, "HPDIprob")
boundName <- paste0(100 * HPDI, "% HPDCI")
tslBounds <- list(
ub = x$tsl$nawruSummary[, paste0(100 * HPDI, "% HPDI-UB")],
lb = x$tsl$nawruSummary[, paste0(100 * HPDI, "% HPDI-LB")],
ub1 = x$tsl$obsSummary[, paste0("ur.", 100 * HPDI, "% HPDI-UB")],
lb1 = x$tsl$obsSummary[, paste0("ur.", 100 * HPDI, "% HPDI-LB")],
ub2 = x$tsl$obsSummary[, paste0("pcInd.", 100 * HPDI, "% HPDI-UB")],
lb2 = x$tsl$obsSummary[, paste0("pcInd.", 100 * HPDI, "% HPDI-LB")],
ub3 = x$tsl$stateSmoothedSummary[, paste0("cycle.", 100 * HPDI, "% HPDI-UB")],
lb3 = x$tsl$stateSmoothedSummary[, paste0("cycle.", 100 * HPDI, "% HPDI-LB")] )
}
# nawru
tsl1 <- list(
trend = x$tsl$nawru,
orig = x$tsl$obs[, 1],
lb = tslBounds$lb,
ub = tslBounds$ub,
lb_fc = tslBounds$lb1,
ub_fc = tslBounds$ub1
)
tsl1 <- do.call(cbind, tsl1)
# phillips curve
tsl2 <- do.call(cbind, list(
E2 = x$tsl$obs[, 2],
lb_fc = tslBounds$lb2,
ub_fc = tslBounds$ub2
))
# cycle
tsl3 <- do.call(cbind, list(
cycle = 100 * x$tsl$stateSmoothed[, "cycle"],
lb = tslBounds$lb3,
ub = tslBounds$ub3
))
# combine lists
tsl <- list(tsl1, tsl2, tsl3)
# --- legends and titles and print names
legend <- list(
c("nawru", "unemployment rate", rep(paste0(boundName, " (ur)"), 2)),
c("Phillips curve indicator", rep(paste0(boundName, ""), 2)),
c("unemployment gap")
)
title <- list(
"Unemployment rate in %",
"Phillips curve",
"Unemployment gap"
)
namesPrint <- c("nawru", "phillips_curve", "u_gap")
if (!is.null(prefix)) namesPrint <- paste(prefix, namesPrint , sep = "_")
# plot
plotSSprediction(
tsl = tsl, legend = legend, title = title, n.ahead = n.ahead,
boundName = boundName, res = NULL, namesPrint = namesPrint,
bounds = bounds, combine = combine, path = path, device = device,
width = width, height = height
)
}
# ----- bayesian estimation
} else if (posterior == TRUE & method != "MLE") {
R <- attr(x, "R")
burnin <- attr(x, "burnin")
thin <- attr(x, "thin")
HPDIprob <- attr(x, "HPDIprob")
FUN <- attr(x, "FUN")
loc <- x$model$loc
param <- x$parametersSampled
parNames <- colnames(param)
cycle <- attr(x$model, "cycle")
prior <- x$prior
# convert prior mean and standard deviation to actual parameters
# gamma and beta needs to be adjusted
distrPar <- .priorMSd2Parameter(
prior = cbind(prior$cycle, prior$trend, prior[[3]])[1:2, ],
restr = cbind(prior$cycle, prior$trend, prior[[3]])[3:4, ],
namesInvGammaDistr = loc$varName[loc$distribution == "invgamma"],
namesNormalDistr = loc$varName[loc$distribution == "normal"],
namesBetaDistr = loc$varName[loc$distribution == "beta"]
)
for (jj in parNames) {
distr <- loc$distribution[loc$varName == jj]
inputl <- list(
draws = param[, jj], burnin = burnin / thin, parName = gsub("E2", "pc", jj),
lb = distrPar[3, jj], ub = distrPar[4, jj]
)
switch(distr,
"invgamma" = do.call(plot_gibbs_output, args = c(path, inputl,
prefix = prefix, shape = distrPar[2, jj] / 2,
scale = distrPar[1, jj] / 2, device = device,
width = width, height = height
)),
"beta" = do.call(plot_gibbs_output, args = c(path, inputl,
prefix = prefix, shape1 = distrPar[1, jj],
shape2 = distrPar[2, jj], device = device,
width = width, height = height
)),
"normal" = do.call(plot_gibbs_output, args = c(path, inputl,
prefix = prefix, mu = distrPar[1, jj],
prec = 1 / distrPar[2, jj], device = device,
width = width, height = height
))
)
}
}
}
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.