Nothing
#' HEAVY model estimation
#'
#' @description This function calculates the High frEquency bAsed VolatilitY (HEAVY) model proposed in Shephard and Sheppard (2010).
#'
#' @param data an \code{xts} object where the first column is a vector of returns
#' and the second column is a vector of realized stock market variation
#' @param startingValues a vector of alternative starting values: first three arguments for variance equation and last three arguments for measurement equation.
#'
#' @return The function outputs an object of class \code{HEAVYmodel}, a list containing
#' \itemize{
#' \item coefficients = estimated coefficients.
#' \item se = robust standard errors based on inverted Hessian matrix.
#' \item residuals = the residuals in the return equation.
#' \item llh = the two-component log-likelihood values.
#' \item varCondVariances = conditional variances in the variance equation.
#' \item RMCondVariances = conditional variances in the RM equation.
#' \item data = the input data.
#' }
#' The class HEAVYmodel has the following methods: plot.HEAVYmodel, predict.HEAVYmodel,
#' print.HEAVYmodel, and summary.HEAVYmodel.
#'
#' @details
#' Let \eqn{r_{t}} and \eqn{RM_{t}} be series of demeaned returns and realized measures of
#' daily stock price variation. The HEAVY model is a two-component model.
#' We assume \eqn{r_{t} = h_{t}^{1/2} Z_{t}} where \eqn{Z_t} is an i.i.d. zero-mean
#' and unit-variance innovation term. The dynamics of the HEAVY model are given by
#'
#' \deqn{
#' h_{t} = \omega + \alpha RM_{t-1} + \beta h_{t-1}
#' }
#' and
#' \deqn{
#' \mu_{t} = \omega_{R} + \alpha_{R} RM_{t-1} + \beta_{R} \mu_{t-1}.
#' }
#'
#' The two equations are estimated separately as mentioned in Shephard and Sheppard (2010).
#' We report robust standard errors based on the matrix-product of inverted Hessians and
#' the outer product of gradients.
#'
#' Note that we always demean the returns in the data input as we don't include a constant in the mean equation.
#'
#' @references Shephard, N. and Sheppard, K. (2010). Realising the future: Forecasting with high frequency based volatility (HEAVY) models. Journal of Applied Econometrics 25, 197--231.
#' @importFrom numDeriv jacobian hessian
#' @importFrom stats nlminb
#' @importFrom Rsolnp solnp
#' @author Onno Kleen and Emil Sjorup.
#'
#' @examplesIf !grepl("debian", sessionInfo()["platform"], fixed = TRUE)
#'
#' # Calculate returns in percentages
#' logReturns <- 100 * makeReturns(SPYRM$CLOSE)[-1]
#'
#' # Combine both returns and realized measures into one xts
#' # Due to return calculation, the first observation is missing
#' dataSPY <- xts::xts(cbind(logReturns, SPYRM$BPV5[-1] * 10000), order.by = SPYRM$DT[-1])
#'
#' # Fit the HEAVY model
#' fittedHEAVY <- HEAVYmodel(dataSPY)
#'
#' # Examine the estimated coefficients and robust standard errors
#' fittedHEAVY
#'
#' # Calculate iterative multi-step-ahead forecasts
#' predict(fittedHEAVY, stepsAhead = 12)
#' @seealso
#'
#' \code{\link{predict.HEAVYmodel}}
#'
#' @export
HEAVYmodel <- function(data, startingValues = NULL) {
ret <- as.numeric(data[,1])
rm <- as.numeric(data[,2])
ret <- ret - mean(ret)
if (is.null(startingValues)) {
startingValuesVar <- c(mean(ret^2) * (1 - 0.3 * mean(rm) / mean(ret^2) - 0.5), 0.3, 0.5)
startingValuesRM <- c(mean(rm) * (1 - 0.6 - 0.3), 0.6, 0.3)
} else {
startingValuesVar <-startingValues[1:3]
startingValuesRM <-startingValues[4:6]
}
rsolVar <- solnp(pars = startingValuesVar, fun = function(x) -sum(heavyLLH(x, ret = ret, rm = rm)),
ineqfun = function(x) c(x, x[2] + x[3]), ineqLB = c(0,0,0,0), ineqUB = c(Inf, Inf, 1, Inf), control = list(trace = 0))
rsolRM <- solnp(pars = startingValuesRM, fun = function(x) -sum(heavyLLH(x, rm = rm, RMEq = TRUE)),
ineqfun = function(x) c(x, x[2] + x[3]), ineqLB = c(0,0,0,0), ineqUB = c(Inf, 1, 1, 1), control = list(trace = 0))
varCondVariances <- calcRecVarEq(rsolVar$pars, rm)
RMCondVariances <- calcRecVarEq(rsolRM$pars, rm)
modelEstimates <- c(rsolVar$pars, rsolRM$pars)
names(modelEstimates) <- c("omega", "alpha", "beta",
"omegaR", "alphaR", "betaR")
invHessianVarEq <- try({
solve(-rsolVar$hessian[(nrow(rsolVar$hessian)-2):nrow(rsolVar$hessian),(nrow(rsolVar$hessian)-2):nrow(rsolVar$hessian)])
}, silent = TRUE)
invHessianRMEq <- try({
solve(-rsolRM$hessian[(nrow(rsolRM$hessian)-2):nrow(rsolRM$hessian),(nrow(rsolRM$hessian)-2):nrow(rsolRM$hessian)])
}, silent = TRUE)
if (class(invHessianVarEq)[1] == "try-error") {
warning("Inverting the Hessian matrix failed. No robust standard errors calculated for variance equation.")
robStdErrVarEq <- rep(NA, times = 3)
} else {
robStdErrVarEq <- sqrt(diag(invHessianVarEq %*%
crossprod(jacobian(function(theta) heavyLLH(theta, rm = rm, ret = ret), rsolVar$pars)) %*%
invHessianVarEq))
}
if (class(invHessianRMEq)[1] == "try-error") {
warning("Inverting the Hessian matrix failed. No robust standard errors calculated for RM equation.")
robStdErrRMEq <- rep(NA, times = 3)
} else {
robStdErrRMEq <-
sqrt(diag(invHessianRMEq %*%
crossprod(jacobian(function(theta) heavyLLH(theta, rm = rm, RMEq = TRUE), rsolRM$pars)) %*%
invHessianRMEq))
}
model <- list(
"coefficients" = modelEstimates,
"se" = c(robStdErrVarEq, robStdErrRMEq),
"residuals" = ret / sqrt(varCondVariances),
"llh" = c(llhVar = -rsolVar$value, llhRM = -rsolRM$value),
"varCondVariances" = varCondVariances,
"RMCondVariances" = RMCondVariances,
"data" = data
)
class(model) <- "HEAVYmodel"
model
}
#' Plotting method for HEAVYmodel objects
#' @param x an object of class \code{HEAVYmodel}.
#' @param ... extra arguments, see details.
#' @details
#' The plotting method has the following optional parameter:
#' \itemize{
#' \item{\code{legend.loc}}{ A string denoting the location of the legend passed on to \code{addLegend} of the \pkg{xts} package}
#' \item{\code{type}}{ A string denoting the type of lot to be made. If \code{type} is \code{"condVar"} the fitted values of the conditional variance of the returns
#' is shown. If \code{type} is different from \code{"condVar"}, the fitted values of the realized measure is shown. Default is \code{"condVar"}}
#' }
#'
#' @importFrom grDevices dev.interactive
#' @importFrom graphics plot panel.smooth points par
#' @importFrom stats residuals
#' @importFrom xts addLegend
#' @export
plot.HEAVYmodel <- function(x, ...){
options <- list(...)
#### List of standard options
opt <- list(legend.loc = "topleft", type = "condVar")
#### Override standard options where user passed new options
opt[names(options)] <- options
if(opt$type == "condVar"){
observed <- x$data[,1]^2
fitted <- x$varCondVariances
# dates <- x$data
# dates <- as.POSIXct(dates)
# type <- x$type
legend.loc <- opt$legend.loc
title <- "Observed squared returns and fitted variances based on HEAVY Model"
plot(cbind(fitted,observed), col = c('blue', 'red'), main = title, lty = c(1,2), yaxis.right = FALSE)
addLegend(legend.loc = legend.loc, on = 1,
legend.names = c("Fitted condidional variances", "Observed squared returns"),
lty = c(1, 2), lwd = c(2, 2),
col = c("blue", "red"))
} else {
observed <- x$data[,2]
fitted <- x$RMCondVariances
# dates <- x$data
# dates <- as.POSIXct(dates)
# type <- x$type
legend.loc <- opt$legend.loc
title <- "Observed realized measure and fitted realized measure based on HEAVY Model"
plot(cbind(fitted,observed), col = c('blue', 'red'), main = title, lty = c(1,2), yaxis.right = FALSE)
addLegend(legend.loc = legend.loc, on = 1,
legend.names = c("Fitted realized measure", "Observed realized measure"),
lty = c(1, 2), lwd = c(2, 2),
col = c("blue", "red"))
}
}
#' Iterative multi-step-ahead forecasting for HEAVY models
#'
#' Calculates forecasts for \eqn{h_{T+k}}, where \eqn{T} denotes the end of the estimation
#' period for fitting the HEAVYmodel and \eqn{k = 1, \dots, \code{stepsAhead}}.
#'
#' @param object an object of class HEAVYmodel.
#' @param stepsAhead the number of days iterative forecasts are calculated for (default 10).
#' @param ... further arguments passed to or from other methods.
#' @export
predict.HEAVYmodel <- function(object, stepsAhead = 10, ...) {
parVarEq <- object$coefficients[1:3]
parRMEq <- object$coefficients[4:6]
Bjs <- lapply(c(0:(stepsAhead)), function(s) Bj(s, parVarEq, parRMEq))
fcts <- t(vapply(c(1:stepsAhead), function(x) forecastsHEAVY(x, object, parVarEq, parRMEq, Bjs), numeric(2)))
rownames(fcts) <- paste0("T + ", c(1:stepsAhead))
colnames(fcts) <- c('CondVar', 'CondRM')
fcts
}
#' @export
print.HEAVYmodel <- function(x, digits = max(3, getOption("digits") - 3), ...) {
cat("")
summ <- summary(x)
print(summ$coefficients)
cat("\n")
cat(paste0("The log-likelihoods are: \n",
"Variance equation: ", sprintf("%.3f", summ$llhVar), "\n",
"RM equation:", sprintf("%.3f", summ$llhRM)), "\n")
}
#' @export
summary.HEAVYmodel <- function(object, ...) {
ans <- list()
ans$coefficients <-
cbind(Estimate = object$coefficients,
`Std. Error` = object$se,
`t value` = object$coefficients / object$se,
`Pr(>|t|)` = 2 * pnorm(abs(object$coefficients / object$se),
lower.tail = FALSE))
ans$llhVar <- object$llh[1]
ans$llhRM <- object$llh[2]
class(ans) <- "summary.heavy"
ans
}
#' @export
print.summary.heavy <- function(x, ...){
cat("")
print(x$coefficients)
cat("\n")
cat(paste0("The log-likelihoods are: \n",
"Variance equation: ", sprintf("%.3f", x$llhVar), "\n",
"RM equation:", sprintf("%.3f", x$llhRM)), "\n")
}
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.