Nothing
# -------------------------------------------------------------------------------------------
#' Draws from the posterior the autoregressive parameter of a stationary AR(1) process
#' without starting values.
#'
#' @param Y a \code{Tn x 1} vector.
#' @param phi a scalar containing the last draw of the autoregressive parameter \eqn{\phi}.
#' @param sigma a scalar containing the innovation variance.
#' @param phi0 a scalar containing the prior mean for \code{phi}.
#' @param Q0 a scalar containing the prior precision for \code{phi}.
#' @param lb (optional) lower bound for \code{phi}.
#' @param ub (optional) upper bound for \code{phi}.
#'
#' @details The corresponding model is given by \eqn{Y_t = \phi Y_{t-1} + e_t}, where
#' \eqn{e_t ~ N(0, \sigma)} with prior distribution
#' \eqn{p(\phi) = N(\phi_0, 1/Q_0 )}.
#' @details Conditional on the variance \eqn{\sigma}, the posterior is normal and known.
#' @details Stationarity and box constraints are enforced. If the stationarity constraint
#' is not fulfilled, the last draw is returned.
#'
#' @importFrom stats runif
#' @keywords internal
.postAR1 <- function(Y, phi, phi0, Q0, sigma, lb = -Inf, ub = Inf) {
# gibbs step for phi of AR(1) process: Yt = phi Yt-1 + et, et ~ N(0, sigma)
# draws for phi | sigma, Y
# dimension
Tn <- length(Y)
# posterior parameters
Qstar <- (sum(Y[2:Tn]^2) - Y[1]^2) * solve(sigma) + Q0
phistar <- solve(Qstar) * (Y[2:Tn] %*% Y[1:(Tn - 1)] * solve(sigma) + phi0 * Q0)
# draw from posterior with box constraints
tmp <- pnorm(lb, mean = phistar, sd = sqrt(solve(Qstar))) +
runif(1, 0, 1) * (pnorm(ub, mean = phistar, sd = sqrt(solve(Qstar))) - pnorm(lb, mean = phistar, sd = sqrt(solve(Qstar))))
phiProp <- qnorm(tmp, mean = phistar, sd = sqrt(solve(Qstar)))
# stationarity
if (!all(abs(polyroot(z = c(1, -phiProp))) > 1)) {
phiProp <- phi
}
names(phiProp) <- "phi"
return(phiProp)
}
# -------------------------------------------------------------------------------------------
#' Draws from the posterior of the autoregressive paramteres of a stationary AR(p),
#' \eqn{p > 1} process without starting values.
#'
#' @param Y A \code{Tn x 1} vector with the time series.
#' @param phi a \code{1 x p} vector containing the last draw of the autoregressive parameters
#' \eqn{\phi}. \code{p} has to be larger than one.
#' @param sigma a scalar containing the innovation variance.
#' @param phi0 a \code{1 x p} vector containing the prior mean for \code{phi}.
#' @param Q0 a \code{p x p} matrix containing the prior precision for \code{phi}.
#' @param lb (optional) \code{1 x p} vector with lower bounds for \code{phi}.
#' @param ub (optional) \code{1 x p} vector with upper bounds for \code{phi}.
#'
#' @details The corresponding model is given by
#' \eqn{Y_t = \phi_1 Y_{t-1} + ... + \phi_p Y_{t-p} + e_t}, where
#' \eqn{e_t ~ N(0, \sigma)} with prior distribution \eqn{p(\phi) = N(\phi_0, 1/Q_0 )}.
#' @details The posterior draw is obtained via a Metropolis Hastings step with proposal density
#' \eqn{ q = \prod_{t=p+1}^Tn p(Y_t, \phi, \sigma, Y_{t-1}, ..., Y_{t-p} )} which is
#' known due to conjugacy. The acceptance probability is given by
#' \eqn{ \alpha = \min{1, p(Y_1, ... Y_p | \phi_r, \sigma) / p(Y_1, ... Y_p | \phi_{r-1}, \sigma)}}
#' where the subscript \eqn{r} denotes the r-th draw. \eqn{p(Y_1, ... Y_p | \phi_r, \sigma)}
#' is itself normal.
#' @details Stationarity and box constraints are enforced. If the constraints
#' are not fulfilled, the last draw is returned.
#'
#' @importFrom stats runif
#' @keywords internal
.postARp <- function(Y, phi, phi0, Q0, sigma, lb = -Inf, ub = Inf) {
phi <- drop(phi)
phi0 <- drop(phi0)
# dimension
Tn <- length(Y)
p <- length(phi)
X <- sapply(1:p, function(x) {
Y[seq((p + 1 - x), Tn - x)]
})
y <- Y[(p + 1):Tn]
Qstar <- t(X) %*% X / sigma + Q0
phistar <- solve(Qstar) %*% (t(X) %*% y / sigma + Q0 %*% phi0)
# Draw theta from proposal density (procedure such that ub and lb hold)
phiProp <- .mvrnorm(mu = phistar, sigma = solve(Qstar))
# variance covariance of c_1 and c_2
covCLast <- .covAR(k = p, phi = phi, sigma = sigma)
covCProp <- .covAR(k = p, phi = phiProp, sigma = sigma)
# Compute acceptance probability
alpha_ic <- min(1, det(covCLast)^(1 / 2) / det(covCProp)^(1 / 2) *
exp(1 / 2 * Y[1:p] %*% (solve(covCLast) - solve(covCProp)) %*% Y[1:p]))
# stationarity and box constraints for phi
if (!all(abs(polyroot(z = c(1, -phiProp))) > 1) | !(all(phiProp > lb) & all(phiProp < ub))) {
alpha_ic <- 0
}
# Accept-Reject the current draw
if (runif(1, 0, 1) < alpha_ic) {
phir <- phiProp
covCstar <- covCProp
} else {
phir <- phi
covCstar <- covCLast
}
res <- list(
phi = phir,
sigmaY12 = covCstar,
aProb = alpha_ic
)
return(res)
}
# -------------------------------------------------------------------------------------------
#' Draws from the posterior of the variance parameter of a random walk or a random walk with
#' constant or stochastic drift.
#'
#' @param Y A \code{Tn x 1} vector with the dependent variable.
#' @param sigmaDistr A \code{1 x k} matrix with prior distribution and box constraints for
#' the innovation variance. The first two entries contain the prior hyperparameters and
#' the last two entries the upper and lower bound.
#' @param sigmaLast A scalar containing the last draw of the innovation variance.
#' @param muDistr A \code{k x 1} matrix with prior distribution and box constraints for
#' the constant trend. The first two entries contain the prior hyperparameters and
#' the last two entries the upper and lower bound.
#'
#' @details If the process follows a random walk with constant drift, the two parameters are
#' drawn sequentially (conditional on the other parameter). The constant is drawn from a
#' normal posterior given by conjugacy.
#' @details The innovation variance is drawn from a Inverse-Gamma posterior given by
#' conjugacy.
#'
#' @importFrom stats runif pgamma
#' @keywords internal
.postRW <- function(Y, sigmaDistr, sigmaLast = NULL, muDistr = NULL) {
# RW1 or RW2
parr <- NULL
diff <- 1
if (is.null(muDistr) | length(muDistr) == 0) { # RW2
diff <- 2
muProp <- 0
}
# dimension
Tn <- length(Y)
# difference series
y <- diff(Y, lag = 1, differences = diff)
if (diff == 1) {
# ----- constant
# prior parameters
mu0 <- muDistr[1, ]
V0 <- muDistr[2, ]
lb <- muDistr[3, ]
ub <- muDistr[4, ]
# posterior parameters
Vstar <- 1 / ((Tn - 1) / sigmaLast + 1 / V0)
mustar <- Vstar * (sum(y) / sigmaLast + mu0 / V0)
# draw from posterior with box constraints
tmp <- pnorm(lb, mean = mustar, sd = sqrt(Vstar)) +
runif(1, 0, 1) * (pnorm(ub, mean = mustar, sd = sqrt(Vstar)) - pnorm(lb, mean = mustar, sd = sqrt(Vstar)))
muProp <- qnorm(tmp, mean = mustar, sd = sqrt(Vstar))
parr <- muProp
names(parr) <- colnames(muDistr)
}
# ----- sigma
# prior parameters
s0 <- sigmaDistr[1, ]
nu0 <- sigmaDistr[2, ]
lb <- max(0, sigmaDistr[3, ])
ub <- sigmaDistr[4, ]
# posterior parameters
nustar <- nu0 + Tn - diff
sstar <- s0 + sum((y - muProp)^2) # seems low
# draw from posterior with box constraints
# x in [lb, ub] <=> 1/x in [1/ub, 1/lb]
tmp <- pgamma(1 / ub, shape = nustar / 2, rate = sstar / 2) +
runif(1, 0, 1) * (pgamma(1 / lb, shape = nustar / 2, rate = sstar / 2) - pgamma(1 / ub, shape = nustar / 2, rate = sstar / 2))
sigmar <- 1 / qgamma(tmp, shape = nustar / 2, rate = sstar / 2)
names(sigmar) <- colnames(sigmaDistr)
parr <- c(sigmar, parr)
res <- list("par" = parr)
return(res)
}
# -------------------------------------------------------------------------------------------
#' Draws from the posterior of parameters of a normal model with normal inverse gamma prior.
#'
#' @param Y A \code{Tn x 1} vector with the dependent variable.
#' @param X A \code{Tn x k} matrix with the \code{k} explanatory variables.
#' @param betaLast A \code{k x 1} vector containing the last draw.
#' @param betaDistr A \code{4 x k} matrix with prior distribution and box constraints for
#' the parameters of each variable. In each column, the first two entries contain the
#' prior hyperparameters and the last two entries the upper and lower bound.
#' @param sigmaDistr A \code{1 x k} matrix with prior distribution and box constraints for
#' the innovation variance. The first two entries contain the prior hyperparameters and
#' the last two entries the upper and lower bound.
#'
#' @details Draws from the posterior are obtained by conjugacy of the Normal-Inverse-Gamma
#' distribution.
#'
#' @return A list with a draws for beta and the innovation variance.
#'
#' @importFrom stats runif pgamma
#' @keywords internal
.postNIG <- function(Y, X, betaLast, betaDistr, sigmaDistr) {
betaLast <- drop(betaLast)
# dimension
Tn <- length(Y)
# ----- beta and gamma
# prior parameters
betaDistr <- as.matrix(betaDistr)
beta0 <- betaDistr[1, ]
Q0 <- solve(diag(betaDistr[2, ]))
s0 <- sigmaDistr[1]
nu0 <- sigmaDistr[2]
# posterior parameters
Qstar <- Q0 + t(X) %*% X
betastar <- solve(Qstar) %*% (Q0 %*% beta0 + t(X) %*% Y)
nustar <- nu0 + Tn
sstar <- s0 + t(Y) %*% Y + t(beta0) %*% Q0 %*% beta0 - t(betastar) %*% Qstar %*% betastar
# sstar <- s0 + t(Y - X %*% beta0) %*% solve( diag(Tn) + X %*% solve(Q0) %*% t(X) ) %*% (Y - X %*% beta0)
# (is the same)
# dram from inverse gamma enacting constraints if present
lb <- max(0, sigmaDistr[3])
ub <- sigmaDistr[4]
# x in [lb, ub] <=> 1/x in [1/ub, 1/lb]
tmp <- pgamma(1 / ub, shape = nustar / 2, rate = sstar / 2) +
runif(1, 0, 1) * (pgamma(1 / lb, shape = nustar / 2, rate = sstar / 2) - pgamma(1 / ub, shape = nustar / 2, rate = sstar / 2))
sigmar <- 1 / qgamma(tmp, shape = nustar / 2, rate = sstar / 2)
# Get draw for mean from normal distribution conditional on draw for gamma
lb <- betaDistr[3, ]
ub <- betaDistr[4, ]
betar <- .mvrnorm(mu = betastar, sigma = solve(Qstar) * sigmar)
if (!(all(betar > lb) & all(betar < ub))) {
betar <- betaLast
}
# result
res <- list(
"beta" = betar,
"sigma" = sigmar
)
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.