## functions.R
## David Ardia
## bayesGARCH
## __input__
## y : time series
## mu.alpha = c(0,0) : prior mean for alpha
## v.alpha = 1000 * diag(1,2) : prior covariance matrix for alpha
## mu.beta = 0 : prior mean for beta
## v.beta = 1000 : prior variance for beta
## c.nu = 0.01
## d.nu = 2 : prior mean of 102
## control = list() : a list with
## n.chain = 1: number of chain
## l.chain = 1000 : length of each chain
## start.val = c(0.01, 0.1, 0.7, 20) : starting valuesn
## addPriorConditions : a function
## digits = 4 : digits for the output
## refresh = 10 : refreshing frequency
## __output__
## r : a list of the class MCMC.list of length n.chain containing the MCMC outputs
"bayesGARCH" <- function(y, mu.alpha = c(0,0), Sigma.alpha = 1000 * diag(1,2),
mu.beta = 0, Sigma.beta = 1000,
lambda = 0.01, delta = 2, control = list()){
if (missing(y))
stop ("'y' is missing")
if (!is.vector(y))
stop ("'y' must be a vector")
if (length(y) < 2)
stop ("'y' must be a longer time series")
if (any(is.na(y)))
stop ("'y' contains 'NA' values")
if (!is.vector(mu.alpha))
stop ("'mu.alpha' must be a vector")
if (length(mu.alpha) != 2)
stop ("'mu.alpha' is not of appropriate size")
if (!is.matrix(Sigma.alpha))
stop ("'Sigma.alpha' must be a matrix")
if (!all(dim(Sigma.alpha) == c(2, 2)))
stop ("'Sigma.alpha' is not of appropriate size")
if (Sigma.alpha[1,2] != Sigma.alpha[2,1])
stop ("'Sigma.alpha' is not symmetric")
if (any(eigen(Sigma.alpha)$values <= 0))
stop ("'Sigma.alpha' is not positive definite")
if (!is.vector(mu.beta) || length(mu.beta) != 1)
stop ("'mu.beta' must be a scalar")
if (!is.vector(Sigma.beta) || length(Sigma.beta) != 1)
stop ("'Sigma.beta' must be a scalar")
if (Sigma.beta <= 0)
stop ("'Sigma.beta' must be positive")
if (!is.vector(lambda) || length(lambda) != 1)
stop ("'lambda' must be a scalar")
if (lambda <= 0)
stop ("'lambda' must be positive")
if (!is.vector(delta) || length(delta) != 1)
stop ("'delta' must be a scalar")
if (delta < 2)
stop ("'delta' cannot be lower than '2'")
if (!is.list(control))
stop ("'control' must be a list")
if (!is.null(control$addPriorConditions))
if (!is.function(control$addPriorConditions))
stop ("'control$addPriorConditions' must be NULL or a function")
fn.bayesGARCH(y, mu.alpha, solve(Sigma.alpha), mu.beta, 1/Sigma.beta, lambda, delta, control)
}
## sub function fn.bayesGARCH
"fn.bayesGARCH" <- function(y, mu.alpha, iv.alpha, mu.beta, iv.beta, c.nu, d.nu, control){
con <- list(n.chain = NULL, l.chain = 10000, addPriorConditions = NULL, digits = 4, refresh = 10,
start.val = matrix(c(0.01,0.1, 0.7, 100), 1, 4,
dimnames = list("chain1", c("alpha0", "alpha1", "beta", "nu"))),
hypers = list(mu.alpha = mu.alpha, iv.alpha = iv.alpha,
mu.beta = mu.beta, iv.beta = iv.beta, c.nu = c.nu, d.nu = d.nu))
con[names(control)] <- control
if (is.vector(con$start.val))
con$start.val <- matrix(con$start.val, nrow = 1, byrow = TRUE)
if (is.null(con$n.chain))
con$n.chain <- nrow(con$start.val)
else con$start.val <- matrix(rep(con$start.val, con$n.chain),
con$n.chain, byrow = TRUE)
if (is.null(colnames(con$start.val)))
colnames(con$start.val) <- c("alpha0", "alpha1", "beta", "nu")
if (is.null(rownames(con$start.val)))
rownames(con$start.val) <- paste(rep("chain", con$n.chain), 1:con$n.chain, sep = "")
if (is.null(con$addPriorConditions)) {con$addPriorConditions <- function(psi){TRUE}}
r <- list()
for (i in 1:con$n.chain) {
k <- NULL
chain <- fn.initializeChain(con$l.chain, con$start.val[i,])
for (j in 2:con$l.chain) {
k <- j - 1
chain[j, ] <- fn.block(y, chain[k, 1:2], chain[k,3], chain[k, 4], con$hypers$mu.alpha, con$hypers$iv.alpha,
con$hypers$mu.beta, con$hypers$iv.beta, con$hypers$c.nu,
con$hypers$d.nu, con$addPriorConditions)
if (con$refresh > 0 & j%%con$refresh == 0)
cat("chain: ", i, " iteration: ", j, " parameters: ",
round(chain[j, ], con$digits), "\n")
}
r[[i]] <- chain
}
names(r) <- paste(rep("chain", con$n.chain), 1:con$n.chain, sep = "")
attr(r, "class") <- "mcmc.list"
r
}
## cvGARCH
## __input__
## y : Tx1 time series
## theta : vector of scedastic's function parameter, theta := (alpha, beta)
## CSC : covariance statinarity condition
## __ouput__
## h : (T+1)x1 vector of conditional variances
"cvGARCH" <- function(y, theta, CSC = FALSE){
if (missing(y))
stop ("'y' is missing")
if (!is.vector(y))
stop ("'y' should be a vector")
if (length(y) < 2)
stop ("'y' sould be a longer time series")
if (any(is.na(y)))
stop ("'y' contains 'NA' values")
if (missing(theta))
stop ("'theta' is missing")
if (!is.vector(theta))
stop ("'theta' should be a vector")
if (length(theta) != 3)
stop ("not appropriate size for 'theta'")
if (any(theta <= 0))
stop ("some components of 'theta' are not positive")
if (CSC && theta[2] + theta[3] >=1)
warning ("the covariance stationarity condition 'alpha1 + beta < 1' is violated")
fn.cvGARCH(y, theta)
}
"fn.cvGARCH" <- function(y, theta){
n <- length(y)
.C('fnGarchC',
n = as.integer(n),
order = as.integer(c(1,2)),
alpha = as.double(theta[1:2]),
delta = as.double(0),
beta = as.double(theta[3]),
h = vector('double', n+1),
u = as.double(y),
sim = vector('integer',n+1),
PACKAGE = 'bayesGARCH')$h
}
## fn.initializeChain
## __input__
## l.chain : length of the chain
## start.val : starting values
## __output__
## mcmc matrix
"fn.initializeChain" <- function (l.chain, start.val){
r <- matrix(NA, l.chain, length(start.val), dimnames = list(1:l.chain, names(start.val)))
r[1,] <- start.val
coda::mcmc(r, start = 1)
}
## fn.block
## __input__
## y : data
## alpha, beta, nu : parameters
## alpha0, iv.alpha0, beta0, iv.beta0, c.nu, d.nu : hyperparameters
## __output__
## block of parameters
"fn.block" <- function(y, alpha, beta, nu, alpha0, iv.alpha0, beta0, iv.beta0, c.nu, d.nu, addPriorConditions){
w.new <- fn.w.full(y, alpha, beta, nu,
alpha0, iv.alpha0, beta0, iv.beta0, c.nu, d.nu)
alpha.new <- fn.alpha.full(y, alpha, beta, nu, w.new,
alpha0, iv.alpha0, beta0, iv.beta0, c.nu, d.nu)
beta.new <- fn.beta.full(y, alpha.new, beta, nu, w.new,
alpha0, iv.alpha0, beta0, iv.beta0, c.nu, d.nu)
nu.new <- fn.nu.full(y, alpha.new, beta.new, nu, w.new,
alpha0, iv.alpha0, beta0, iv.beta0, c.nu, d.nu)
psi.new <- c(alpha.new, beta.new, nu.new)
if (addPriorConditions(psi.new)) psi.new else c(alpha,beta,nu)
}
## fn.alpha.full
## __input__
## y : data
## alpha.old, beta, nu, w : parameters
## alpha0, iv.alpha0, beta0, iv.beta0, c.nu, d.nu : hyperparameters
## __output__
## MH draw for alpha
"fn.alpha.full" <- function(y, alpha.old, beta, nu, w,
alpha0, iv.alpha0, beta0, iv.beta0, c.nu, d.nu){
r <- alpha.old
n <- length(y)
h <- fn.cvGARCH(y, c(alpha.old, beta))[1:n]
post.old <- fn.post.garch(y, w*h, alpha.old, beta, nu,
alpha0, iv.alpha0, beta0, iv.beta0, c.nu, d.nu)
tau <- 2 * h^2
filter.alpha <- fn.filterAlpha(y, beta)
Dd.D <- fn.Dd.D(y^2/w, filter.alpha, tau,
alpha0, iv.alpha0)
alpha.new <- t( mvtnorm::rmvnorm(n = 1, mean = Dd.D$Dd, sigma = Dd.D$D) )
if (all(alpha.new>0)){
h.new <- fn.cvGARCH(y, c(alpha.new, beta))[1:n]
tau.new <- 2 * h.new^2
post.new <- fn.post.garch(y, w*h.new, alpha.new, beta, nu,
alpha0, iv.alpha0, beta0, iv.beta0, c.nu, d.nu)
prop.new <- mvtnorm::dmvnorm(t(alpha.new), mean = Dd.D$Dd, sigma = Dd.D$D, log = TRUE)
Dd.D <- fn.Dd.D(y^2/w, filter.alpha, tau.new,
alpha0, iv.alpha0)
prop.old <- dmvnorm(alpha.old, mean = Dd.D$Dd, sigma = Dd.D$D, log = TRUE)
ratio <- post.new - post.old + prop.old - prop.new
r.u <- runif(1)
if (r.u < min(1, exp(ratio)))
r <- alpha.new
}
r
}
## fn.beta.full
## __input__
## y : data set
## alpha, beta.old, nu, w : parameters
## alpha0, iv.alpha0, beta0, iv.beta0, c.nu, d.nu : hyperparameters
## __output__
## MH draw for beta
"fn.beta.full" <- function(y, alpha, beta.old, nu, w,
alpha0, iv.alpha0, beta0, iv.beta0, c.nu, d.nu){
n <- length(y)
r <- beta.old
h <- fn.cvGARCH(y, c(alpha, beta.old))[1:n]
post.old <- fn.post.garch(y, w*h, alpha, beta.old, nu,
alpha0, iv.alpha0, beta0, iv.beta0, c.nu, d.nu)
tau <- 2 * h^2
beta.nls <- beta.old
w1 <- fn.filterW(y, alpha, beta.nls)
W <- fn.W(y^2-w1, -beta.nls)
z <- w1 + W * beta.nls
z <- z + y^2 * (1/w - 1)
Dd.D <- fn.Dd.D(z, W, tau, beta0, iv.beta0)
beta.new <- rnorm(n = 1, mean = Dd.D$Dd, sd = sqrt(Dd.D$D))
if (beta.new>0){
h.new <- fn.cvGARCH(y, c(alpha, beta.new))[1:n]
tau.new <- 2 * h.new^2
post.new <- fn.post.garch(y, w*h.new, alpha, beta.new, nu,
alpha0, iv.alpha0, beta0, iv.beta0, c.nu, d.nu)
prop.new <- dnorm(beta.new, mean = Dd.D$Dd, sd = sqrt(Dd.D$D), log = TRUE)
w1 <- fn.filterW(y, alpha, beta.new)
W <- fn.W(y^2-w1, -beta.new)
z <- w1 + W * beta.new
z <- z + y^2 * (1/w - 1)
Dd.D <- fn.Dd.D(z, W, tau.new, beta0, iv.beta0)
prop.old <- dnorm(beta.old, mean = Dd.D$Dd, sd = sqrt(Dd.D$D), log = TRUE)
ratio <- post.new - post.old + prop.old - prop.new
r.u <- runif(1)
if (r.u < min(1, exp(ratio)))
r <- beta.new
}
r
}
## fn.w.full
## __input__
## u : data set
## alpha, beta, nu : parameters
## alpha0, iv.alpha0, beta0, iv.beta0, c.nu, d.nu : hyperparameters
## __output__
## draw for w
"fn.w.full" <- function(u, alpha, beta, nu,
alpha0, iv.alpha0, beta0, iv.beta0, c.nu, d.nu){
n <- length(u)
h <- fn.cvGARCH(u, c(alpha, beta))[1:n]
r <- rgamma(n = n, shape = 0.5 * rep((nu+1),n),
rate = 0.5 * ((u^2/h) + (nu-2)) )
1/r
}
## fn.nu.full
## __input__
## u : data set
## alpha, beta, nu.old, w : parameters
## alpha0, iv.alpha0, beta0, iv.beta0, c.nu, d.nu : hyperparameters
## alpha.min, alpha.max, k.max : simulation parameters
## __output__
## simulated nu
"fn.nu.full" <- function(u, alpha, beta, nu.old, w,
alpha0, iv.alpha0, beta0, iv.beta0, c.nu, d.nu,
alpha.min = 0.00001, alpha.max = 500, k.max = 50000){
n <- length(u)
r <- nu.old
phi <- 0.5 * sum(log(w) + 1/w) + c.nu
optim.alpha <- uniroot(fn.neg.alpha,
interval = c(alpha.min, alpha.max),
tol = .Machine$double.eps,
maxiter = 1000,
n = n,
d = d.nu,
phi = phi)
alpha <- optim.alpha$root
value <- optim.alpha$f.root
iter <- optim.alpha$iter
if ( (alpha>alpha.min) && (alpha<alpha.max) ){
k <- 0
while (k<k.max){
k <- k+1
nu.new <- rexp(n = 1, rate = alpha)
nu.new <- nu.new + d.nu
p <- fn.accept.proba(n, alpha, d.nu, nu.new, phi)
r.u <- runif(1)
if(r.u < p){
r = nu.new
k = k.max
}
}
}
r
}
## fn.post.garch
## __input__
## u : data set
## h : conditional variance process
## alpha, beta, nu : parameters
## alpha0, iv.alpha0, beta0, iv.beta0, c.nu, d.nu : parameters
## __output__
## posterior density
"fn.post.garch" <- function(u, h, alpha, beta, nu,
alpha0, iv.alpha0, beta0, iv.beta0, c.nu, d.nu){
r <- 0
r <- -0.5 * sum(log(h))
r <- r - 0.5 * sum(u^2/h)
r <- r - 0.5 * t(alpha-alpha0) %*%
iv.alpha0 %*% (alpha-alpha0)
r <- r - 0.5 * (beta-beta0)^2 * iv.beta0
r <- r - c.nu * nu
r
}
## fn.neg.alpha
"fn.neg.alpha" <- function(alpha, n, d, phi){
r <- 0
r <- log( (1+alpha*(d-2))/(2*alpha) )
r <- r + (1+alpha*d)/(1+alpha*(d-2))
r <- r - digamma( (1+alpha*d)/(2*alpha) )
r <- 0.5 * n * r + alpha - phi
r
}
## fn.accept.proba
"fn.accept.proba" <- function(n, alpha, d, nu, phi){
u <- (1+alpha*d) / (2*alpha)
v <- nu/2
w <- (1+alpha*(d-2)) / (2*alpha)
r.ln <- 0
r.ln <- n * (lgamma(u) - lgamma(v))
r.ln <- r.ln + n * v * log((nu-2)/2)
r.ln <- r.ln - n * u * log(w)
r.ln <- r.ln + (nu - d) * (alpha - phi) + (phi/alpha) - 1
exp(r.ln)
}
## fn.filterAlpha
"fn.filterAlpha" <- function(u, beta){
n <- length(u)
filter.alpha <- .C('fnFilterAlphaC',
n = as.integer(n),
u = as.double(u),
beta = as.double(beta),
vstar = vector('double', n),
lstar = vector('double', n),
PACKAGE = 'bayesGARCH')
as.matrix(cbind(filter.alpha$lstar,
c(0, filter.alpha$vstar[1:(n-1)])))
}
## fn.filterW
"fn.filterW" <- function(u, alpha, beta){
n <- length(u)
.C('fnFilterWC',
n = as.integer(n),
u = as.double(u),
alpha = as.double(alpha),
beta = as.double(beta),
w = vector('double', n),
PACKAGE = 'bayesGARCH')$w
}
## fn.W
"fn.W" <- function(u, theta){
n <- length(u)
.C('fnQDiffC',
n = as.integer(n),
u = as.double(u),
theta = as.double(theta),
w = vector('double', n),
w0 = as.double(0),
PACKAGE = 'bayesGARCH')$w
}
## fn.Dd.D
## compute D and Dd
"fn.Dd.D" <- function(y, X, h, b0, iv.b0){
r <- list()
Xh <- X/h
D.inv <- t(Xh) %*% X + iv.b0
D <- solve(D.inv)
d <- t(Xh) %*% y + iv.b0 %*% b0
Dd <- D %*% d
list(D = D, Dd = Dd)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.